[发明专利]一种石灰岩土石山区山坡尺度水文过程模拟方法

摘要

本发明涉及一种石灰岩土石山区山坡尺度水文过程模拟方法，包括：山坡计算单元划分；水文气象数据展布；植被截留计算；基岩凹凸面储留计算；入渗及水分再分布过程计算；壤中流计算；地表汇流过程计算；蒸散发过程计算。本发明考虑了石灰岩土石山区广泛存在土石二元混合介质、大孔隙优先流和基岩凹凸面，一方面引入碎石形状系数ε和碎石质量比系数R_v量化了土石二元混合介质对山坡水循环过程的影响；另一方面基于山坡计算单元划分，将计算单元分为大孔隙优先流区和基质区，进行两流区水循环过程计算；同时，基于改进的洼地储留模型，计算了基岩凹凸面对水文过程影响。本发明使石灰岩土质的山坡水文过程模拟更加精确，更接近实际。

主权项

1.一种石灰岩土石山区山坡尺度水文过程模拟方法，其特征在于，所述方法的计算过程如下：山坡计算单元划分：采用等流时线法将山坡划分若干基本计算单元；计算单元内垂直剖面划分：根据山坡植被、土壤和岩石特性，在基本计算单元内分为4层：植被冠层截留层、地表储留层、土壤层、基岩层；植被冠层截留层又可细分为：高植被截留层和矮植被储留层、草地层和裸地；土壤层进一步分为均质土壤层、土石二元混合介质层、风化碎石层；基岩层进一步分为凹凸面储留层区和裂隙优先流区；计算单元内状态变量包括：植被冠层截留量、洼地储留量、基岩凹凸面储留量、枯枝落叶储留量、土壤含水量；主要参数包括：植被最大截留深、洼地最大储留深、基岩凹凸面最大储留深、枯枝落叶干重、土壤导水系数、土壤水分特征曲线、土壤含水量、土壤厚度、土石二元介质碎石质量比系数、坡面糙率；计算单元水文过程计算包括：水文气象数据展布、降雨期间产汇流过程计算、非降雨期间产汇流过程；水文气象数据展布包括水文气象过程空间降尺度展布和降雨时间降尺度展布；采用泰森多边形法和反距离加权平均法进行流域内气象数据的空间展布，包括降雨、气温、风速、空气湿度、净辐射，计算公式如下：式中：D表示待插值点估计值；D_pi表示第pi个参证站点数据；pm表示参证站点个数；λ_pi表示第pi个参证站点数据权重；d_pi表示第pi个参证站点同待插值点的距离；pn表示权重指数；由于日降雨过程的非稳定性，对日降雨数据进一步进行降尺度展布，具体公式如下：S＝a·P+b   (5)式中：为时段t_k内最大降水平均雨强；S表示暴雨参数；t为时间，其中，t_k‑1<t≤t_k；t_k为时段，其中，k∈0,1,...,N，N为时段数；μ表示暴雨衰减系数；P表示日降雨量；T表示日降雨总历时；a，b表示参数；降雨期间产汇流过程计算：降雨期间，土壤蒸散发量较小，将其忽略；计算单元内水文过程主要由降雨→植被冠层截留→入渗产流→汇流过程构成；植被冠层截留计算，计算公式如下：Wr_max＝0.2·Veg·LAI(8)式中：Veg表示裸地‑植被域中植被的面积占计算单元的面积比例；Wr表示植被冠层截留量；Wr_max表示植被冠层最大截留量；I为雨强；Rr表示植被冠层流出量；LAI表示叶面积指数；洼地储留计算，计算公式如下：式中：I’为净雨强，即：经植被冠层截留后的雨强；H_u2为土壤表层储留深度；H_umax2为土壤表层最大储留深度；R_u2为土壤表面径流；f_in为入渗率；枯枝落叶储留计算，计算公式如下：H_umax＝z_maxG_d(11)式中：G_d为枯枝落叶干重；z_max为枯枝落叶最大持水系数；基岩凹凸面储留计算：忽略降雨期间蒸散发条件下，基岩凹凸面储留采用下式计算：式中：H_u1为基岩凹凸面储留深度；H_umax1为基岩凹凸面最大储留深度；K_scr为风化壳导水系数；R_u1为基岩表面径流；L为计算单元长度；q为各层土壤上下界面处水分通量；降雨期间，基岩面上方下渗水分首先填充基岩凹凸面，然后产生基岩面壤中流；非降雨期，基岩凹凸面储留量主要消耗于植被蒸腾过程；土壤及基岩层水分动态过程计算：土壤层水分运动过程采用理查兹(Richards equation)公式计算：式中：h为土壤水吸力；C(h)为容水度；SS为源汇项，即根系吸水；z为坐标轴；K(h)为导水系数；t为时间；考虑到土壤和基岩层均广泛分布着裂隙优先流，土壤水分运动过程中将计算单元分为基质区和优先流区；其中计算单元内基质区所占面积比例为w_m，计算单元内大孔隙优先流区所占面积比例为w_f；对于石灰岩地区，其基质区土壤内含有大量石灰岩碎石，属于土石二元混合介质，但其中的碎石不具有持水性和导水性，影响了土壤水分动态过程；在对其基质区土壤水分动态过程模拟时，引入碎石质量比系数R_v描述石灰岩碎石对土壤含水量的影响；根据R_v大小，可以将土壤层分为以下土层：1)当R_v＝0时，土壤层为均质土壤；2)当1>R_v>0时，土壤层为土石二元混合介质层；3)当R_v＝1时，土壤层为基岩层；则基质区土壤水分运动过程进一步修正为：其中，w_m,i为第i层土壤的基质区所占面积比例；C_ms，i为第i层土壤的基质区内容水度；h_m，i为第i层土壤的基质区内土壤水吸力；SS_f，i为第i层土壤的大孔隙优先流区的源汇项；Γ_m，i为第i层土壤的基质区内不同区间水量交换量；K_m，i(h)为第i层土壤的基质区内导水系数，引入碎石形状系数ε和碎石质量比系数R_v，得到公式(16)；K_m，i(h)＝(1‑εR_v)K_ss(h)   (16)其中：式中；h_i为第i层土壤的土壤水吸力；K_ss(h)为土壤非饱和导水系数；K_ss为土壤饱和导水系数；α₁、vn和vm为参数，vm＝1‑1/vn；z为坐标轴；土石二元混合介质达到饱和时，h＝0，基于公式(16)，计算单元内饱和导水系数为：K_m＝(1－εR_v)·K_ss   (18)石灰岩地区计算单元内优先流区同样含有大量石灰岩碎石，但孔隙结构与基质区不同，导致水势、导水系数等有所不同；其土壤水分过程计算公式表示如下：式中：下标f表示大孔隙优先流区；w_f，i为第i层土壤的大孔隙优先流区所占面积比例；h_f，i为第i层土壤的大孔隙优先流区内的土壤水吸力；Γ_f，i为第i层土壤的大孔隙优先流区的不同区间水量交换量；其他参数同前；上下边界条件均采用通量边界条件：式中：q_m为基质区各层土壤上下界面处水分通量；q_f为大孔隙优先流区各层土壤上下界面处水分通量；w_f，i为第i层土壤的大孔隙优先流区所占面积比例；h_m，i为第i层土壤的基质区内土壤水吸力；下标m和f分别表示基质区和大孔隙优先流区；降雨期间，上边界通量为净雨强；非降雨期间为表层土壤蒸发量通量；下边界为石灰岩基岩面，由于石灰岩透水性差，故认为其是不透水层，但石灰岩基岩内分布的裂隙具有一定的导水性和持水性；因此，下边界大部分地区水分通量为0，只有小部分地区具有水分通量；壤中流计算：R_sub，i为第i层土壤中的垂向入流量；Φ_i为第i层土壤中的壤中流水位高度；Q_i为第i层土壤的壤中过流断面流量；W为计算单元宽度；e_i为第i层土壤的孔隙度；x为坐标轴；R_sub,i＝q_m,i+q_f,i(24)q_m，i为第i层土壤的基质区各层土壤上下界面处水分通量；q_f，i为第i层土壤的大孔隙优先流区各层土壤上下界面处水分通量；q_m,i＝w_m(q_m,i‑1‑q_m,i)   (25)q_f,i＝w_f(q_f,i‑1‑q_f,i)    (26)式中：Q_i为第i层土壤中的壤中过流断面流量；K_s为饱和导水系数；W为计算单元宽度；Φ为壤中流水位高度，即土水势；x和z为坐标轴；e为土壤内孔隙度；R_sub为垂向入流量；q为各层土壤上下界面处水分通量；下标m和f分别表示基质区和大孔隙优先流区；下标i表示土壤层；q_m，i‑1为第i‑1层土壤的基质区各层土壤上下界面处水分通量；q_f，i‑1为第i‑1层土壤的大孔隙优先流区各层土壤上下界面处水分通量；模型计算中，基质区和大孔隙优先流区土壤水分特征曲线采用相同曲线；土壤水分特征曲线采用考虑碎石Van Genuchten模型计算：式中：S_e为饱和度系数；R_v为碎石质量比系数；θ’为土石二元混合介质土壤含水量；θ_ss为土石二元混合介质土壤饱和含水量；θ_r为土石二元混合介质土壤残余含水量；θ_n为无碎石土壤含水量；θ_ns为无碎石土壤饱和含水量；θ_rs为无碎石土壤残余含水量；α、vn和vm为参数，vm＝1‑1/vn；h为土壤水吸力；碎石质量比系数地表汇流过程计算，计算公式如下：连续方法：运动方程：S_f＝S₀    (30)曼宁公式：式中：Q₀表示地表过流断面流量；A表示过流断面面积；R_surf表示地表产流量；S_f表示摩擦坡降；S₀表示计算单元平均地面坡降或河道坡降；R_wr表示过流断面水力半径；kn表示曼宁糙率系数；非降雨期间产汇流过程计算：非降雨期间，土壤水主要消耗于水分进行再分布过程，包括土壤蒸散发和深层渗漏过程；土壤蒸散发过程计算，计算公式如下：E＝E₁+E₂+E₃   (32)式中：E表示计算单元总蒸散发；E₁表示植被冠层截留蒸发量；E₂表示植被蒸腾量；E₃表示裸地土壤蒸发量；土壤潜在蒸发量E_p由Penman公式计算：式中：RN为净辐射量；G为传入水中的热通量；Δ为饱和水汽压对温度的导数；ρ_a为空气密度；C_p为空气的定压比热；δ_e为实际水蒸气压与饱和水蒸气压的差值；r_a为蒸发表面空气动力学阻抗；λ为水的气化潜热；γ为空气湿度常数；PR为大气压；空气动力学阻抗计算：其计算公式如下：式中：r_a为蒸发表面空气动力学阻抗；hz为气象站观测点离地面的高度；hd为置换高度；z_om表为水蒸气紊流扩散对应的地表粗度；z_ox为地表粗度；κ为von Karman常数；U为风速；h_c为植被高度；植被冠层截留蒸发量E₁使用Noilhan‑Planton模型计算：E₁＝Veg·δ·E_p   (37)δ＝(Wr/Wr_max)^2/3   (40)Wr_max＝0.2·Veg·LAI(41)式中：Veg为裸地‑植被域中植被的面积占计算单元的面积比例；δ为湿润叶面占植被叶面的面积比例；E_p为土壤潜在蒸发量，即最大蒸发量；Wr为植被冠层截留量；Wr_max为植被冠层最大截留量；I为雨强；Rr为植被冠层流出水量，即超出植被冠层最大截留量的部分；LAI表示叶面积指数；植被蒸腾量E₂采用Penman‑Monteith公式计算，土壤各层蒸散量采用雷志栋的根系吸水模型进行计算，具体公式如下：E₂＝Veg·(1‑δ)·E_PM   (42)式中：E_PM为采用Penman‑Monteith公式计算的潜在蒸腾量；G为传入水中的热通量；r_c为植被群落阻抗；根系吸水模型：T_r＝E₂(45)式中：lr表示根系层厚度；T_r为实际植被蒸腾量；SS为源汇项；植被群落阻抗计算：采用Dickinson提出公式计算植被群落阻抗：σ₁^‑1＝1‑0.0016(25‑T_a)²   (47)σ₂^‑1＝1‑VPD/VPD_c   (48)式中：r_c为植被群落阻抗；r_smin为最小气孔阻抗；LAI为叶面积指数；σ₁为温度影响函数；σ₂为大气水蒸气压饱和差影响函数；σ₃为光合作用有效放射的影响函数；σ₄为土壤含水量的影响函数；T_a为气温；VPD为饱和水蒸气压同实测值之间的差；VPD_c为气孔闭合时的VPD值；r_smax为最大气孔阻抗；PAR为光和作用有效放射；PAR_c为PAR的临界值；θ为表层土壤含水量；θ_w为植被凋萎时的土壤含水量；θ_s为饱和土壤含水量；裸地土壤蒸发量E₃由下式计算：式中：β为土壤湿润函数；θ为表层土壤含水量；θ_m为土壤分子吸力对应的土壤含水量；θ_h为表层土壤田间持水量；非降雨期土壤水分动态过程采用修正的Richards公式计算，计算公式同非降雨期；当非降雨期地表有地下水流出时，采用运动波方程进行汇流计算，计算公式同降雨期；非降雨期壤中流采用改进后的运动波方程，计算公式同降雨期。