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

权利要求书

1.一种片麻岩土石山区山坡尺度水文过程模拟方法，其特征在于，所述方法的计算过程如下：

山坡计算单元划分：采用等流时线法将山坡划分为若干基本计算单元；

计算单元垂直剖面划分：根据山坡植被、土壤和岩石特性，在基本计算单元内分为4层：植被冠层截留层、地表储留层、土壤层、基岩层；植被冠层截留层又细分为：高植被截留层和矮植被储留层、草地层和裸地；土壤层进一步分为均质土壤层、土石二元混合介质层；考虑到片麻岩基岩层内广泛分布着构造节理，且片麻岩本身具有一定持水性和透水性，进一步将基岩层分为基质流区和优先流区；

计算单元内状态变量包括：植被冠层截留量、洼地储留量、枯枝落叶储留量、土壤含水量；主要参数包括：植被最大截留深、洼地最大储留深、枯枝落叶干重、土壤导水系数、土壤水分特征曲线、土壤含水量、各层土壤厚度、土石二元混合介质碎石质量比系数、基岩层厚度、坡面糙率；

计算单元水文过程计算：包括水文气象数据展布、降雨期间产汇流过程计算、非降雨期间产汇流过程计算；

水文气象数据展布：包括水文气象过程空间尺度展布和降雨时间降尺度展布：

采用泰森多边形法和反距离加权平均法进行流域内气象数据的空间展布，包括降雨、气温、风速、空气湿度、净辐射，计算公式如下：

式中：D表示待插值点估计值；D_pi表示第pi个参证站点数据；pm表示参证站点个数；λ_pi表示第pi个参证站点数据权重；d_pi表示第pi个参证站点同待插值点的距离；pn表示权重指数；

由于日降雨过程的非稳定性，对日降雨数据进一步进行降尺度展布，具体公式如下：

S＝a·P+b (5)

式中：I为时段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为时段t_k内最大降水平均雨强；Rr表示植被冠层流出水量；LAI表示叶面积指数；

洼地储留计算，计算公式如下：

式中：I’为时段t_k内净雨强；H_u2为地表储留；H_umax2为土壤表层最大储留深；R_u2为土壤表面径流；f_in为入渗率；

枯枝落叶储留计算，计算公式如下：

H_umax＝z_maxG_k (11)

式中：G_k为枯枝落叶干重；z_max为枯枝落叶最大持水系数；

土壤及基岩层水分运动过程计算：

片麻岩分布区土壤和基岩均具有导水性和持水性，所有土壤层水分运动过程均采用理查兹公式计算：

式中：h为土壤水吸力；C为容水度；SS为源汇项；z为坐标轴；K(h)为导水系数；t为时间；考虑到土壤和基岩层均广泛分布着裂隙优先流，土壤水分运动过程中将计算单元分为基质区和优先流区；其中计算单元内基质区所占面积比例为W_m，优先流区所占面积比例为W_f；

根据碎石质量比系数R_v大小，将土壤层分为以下土层：

1)当R_v＝0时，土壤层为均质土壤；

2)当1>R_v>0时，土壤层为土石二元混合介质层；

3)当R_v＝1时，土壤层为基岩层；

此外，假设碎石内和土壤内水势相等，将基质区土壤水分运动过程进一步修正为：

其中，K_m(h)为基质区内平均导水系数，由基质区内碎石和土壤的导水系数加权平均得到，引入碎石形状系数ε和质量比系数R_v，得到公式：

K_m(h)＝(1-R_v)K_ss(h)+εR_vK_sr(h) (14)

其中：

式中：w_m为基质区面积比例；h为土壤水吸力；C_ms和C_mr分别为土壤和碎石的容水度；SS为源汇项；Γ为不同区间水量交换量；下标m表示基质区；下标i表示土壤层；K_ss(h)为土壤非饱和导水系数；K_sr(h)为碎石非饱和导水系数；K_ss为土壤饱和导水系数；K_sr为碎石饱和导水系数；α、vn和vm为参数，vm＝1-1/vn；下标1和2分别表示土壤和碎石；z为坐标轴；

土石二元混合介质达到饱和时，h＝0，基于公式(14)，计算单元内饱和导水系数为：

K_m＝(1-R_v)·K_ss+εR_v·K_sr (17)

土壤水分过程计算公式表示如下：

式中：下标f表示优先流区；w_f为优先流区所占面积比例；

上下边界条件均采用通量边界条件：

式中：q为各层土壤上下界面处水分通量；

基于Van Genuchten模型分别计算两者的土壤水分特征曲线：

式中：S_e为饱和度系数；θ为土石二元混合介质土壤含水量；θ_s为土石二元混合介质土壤饱和含水量；θ_r为土石二元混合介质残余含水量；α、vn和vm为参数，vm＝1-1/vn；h为土壤水吸力；下标1和2分别表示土壤和碎石；

壤中流计算：

R_sub,i＝q_m,i+q_f,i (25)

q_m,i＝w_m(q_m,i-1-q_m,i) (26)

q_f,i＝w_f(q_f,i-1-q_f,i) (27)

式中：Q为壤中过流断面流量；K_s为土壤饱和导水系数；W为计算单元宽度；Φ为壤中流水位高度，即土水势；x和z为坐标轴；e为土壤内孔隙度；R_sub为垂向入流量；q为各层土壤上下界面处水分通量；下标m和f分别表示基质区和优先流区；下标i表示土壤层；

地表汇流过程计算，计算公式如下：

连续方法：

运动方程：

S_f＝S₀ (29)

曼宁公式：

式中：Q₀表示地表过流断面流量；A表示过流断面面积；R_surf表示地表产流量；S_f表示摩擦坡降；S₀表示计算单元平均地面坡降或河道坡降；R_wr表示过流断面水力半径；kn表示曼宁糙率系数；

非降雨期土壤蒸散发过程计算，计算公式如下：

E＝E₁+E₂+E₃ (31)

式中：E表示计算单元总蒸散发量；E₁表示植被截留蒸发量；E₂表示植被蒸腾量；E₃表示裸地土壤蒸发量；

土壤潜在蒸发量由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 (36)

δ＝(Wr/Wr_max)^2/3 (39)

Wr_max＝0.2·Veg·LAI (40)

式中：Veg为裸地-植被域中植被的面积占计算单元的面积比例；δ为湿润叶面占植被叶面的面积比例；E_p为土壤潜在蒸发量，即最大蒸发量；Wr为植被截留水量；Wr_max为植被最大截留水量；I为时段t_k内最大降水平均雨强；Rr为植被冠层流出水量，即超出最大植被截留水量的部分；LAI表示叶面积指数；

植被蒸腾量E₂采用Penman-Monteith公式计算，土壤各层蒸散量采用雷志栋的根系吸水模型进行计算，公式如下：

E₂＝Veg·(1-δ)·E_PM (41)

式中：E_PM为采用Penman-Monteith公式计算的潜在蒸腾量；G为传入水中的热通量；r_c为植被群落阻抗；

根系吸水模型：

T_r＝E₂ (44)

式中：lr表示根系层厚度；T_r为实际植被蒸腾量；SS为源汇项；

植被群落阻抗计算：采用Dickinson公式计算植被群落阻抗：

σ₁^-1＝1-0.0016(25-T_a)² (46)

σ₂^-1＝1-VPD/VPD_c (47)

式中：r_c为植被群落阻抗；r_{s min}为最小气孔阻抗；LAI为叶面积指数；σ₁为温度影响函数；σ₂为大气水蒸气压饱和差影响函数；σ₃为光合作用有效放射的影响函数；σ₄为土壤含水量的影响函数；T_a为气温；VPD为饱和水蒸气压同实测值之间的差；VPD_c为气孔闭合时的VPD值；r_smax为最大气孔阻抗；PAR为光和作用有效放射；PAR_c为PAR的临界值；θ₃为根系层土壤含水量；θ_w为植被凋萎时的土壤含水量；θ_s为土壤饱和含水量；

裸地土壤蒸发量E₃由下式计算：

式中：β为土壤湿润函数；θ₄为表层土壤含水量；θ_m为土壤分子吸力对应的土壤含水量；θ_h为表层土壤田间持水量；

非降雨期土壤水分运动过程采用修正的Richards公式计算，计算公式同非降雨期；当非降雨期地表有地下水流出时，采用动力波方程进行汇流计算，计算公式同降雨期；非降雨期壤中流采用改进后的动力波方程，计算公式同降雨期。