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

摘要

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

主权项

1.一种片麻岩土石山区山坡尺度水文过程模拟方法，其特征在于，所述方法的计算过程如下：山坡计算单元划分：采用等流时线法将山坡划分为若干基本计算单元；计算单元垂直剖面划分：根据山坡植被、土壤和岩石特性，在基本计算单元内分为4层：植被冠层截留层、地表储留层、土壤层、基岩层；植被冠层截留层又细分为：高植被截留层和矮植被储留层、草地层和裸地；土壤层进一步分为均质土壤层、土石二元混合介质层；考虑到片麻岩基岩层内广泛分布着构造节理，且片麻岩本身具有一定持水性和透水性，进一步将基岩层分为基质流区和优先流区；计算单元内状态变量包括：植被冠层截留量、洼地储留量、枯枝落叶储留量、土壤含水量；主要参数包括：植被最大截留深、洼地最大储留深、枯枝落叶干重、土壤导水系数、土壤水分特征曲线、土壤含水量、各层土壤厚度、土石二元混合介质碎石质量比系数、基岩层厚度、坡面糙率；计算单元水文过程计算：包括水文气象数据展布、降雨期间产汇流过程计算、非降雨期间产汇流过程计算；水文气象数据展布：包括水文气象过程空间尺度展布和降雨时间降尺度展布：采用泰森多边形法和反距离加权平均法进行流域内气象数据的空间展布，包括降雨、气温、风速、空气湿度、净辐射，计算公式如下：式中：D表示待插值点估计值；Dpi表示第pi个参证站点数据；pm表示参证站点个数；λpi表示第pi个参证站点数据权重；dpi表示第pi个参证站点同待插值点的距离；pn表示权重指数；由于日降雨过程的非稳定性，对日降雨数据进一步进行降尺度展布，具体公式如下：S＝a·P+b    (5)式中：I为时段tk内最大降水平均雨强；S表示暴雨参数；t为时间，其中，tk‑1Wrmax＝0.2·Veg·LAI    (8)式中：Veg表示裸地‑植被域中植被的面积占计算单元的面积比例；Wr表示植被截留水量；Wrmax表示植被最大截留水量；I为时段tk内最大降水平均雨强；Rr表示植被冠层流出水量；LAI表示叶面积指数；洼地储留计算，计算公式如下：式中：I’为时段tk内净雨强；Hu2为地表储留；Humax2为土壤表层最大储留深；Ru2为土壤表面径流；fin为入渗率；枯枝落叶储留计算，计算公式如下：Humax＝zmaxGk    (11)式中：Gk为枯枝落叶干重；zmax为枯枝落叶最大持水系数；土壤及基岩层水分运动过程计算：片麻岩分布区土壤和基岩均具有导水性和持水性，所有土壤层水分运动过程均采用理查兹公式计算：式中：h为土壤水吸力；C为容水度；SS为源汇项；z为坐标轴；K(h)为导水系数；t为时间；考虑到土壤和基岩层均广泛分布着裂隙优先流，土壤水分运动过程中将计算单元分为基质区和优先流区；其中计算单元内基质区所占面积比例为Wm，优先流区所占面积比例为Wf；根据碎石质量比系数Rv大小，将土壤层分为以下土层：1)当Rv＝0时，土壤层为均质土壤；2)当1>Rv>0时，土壤层为土石二元混合介质层；3)当Rv＝1时，土壤层为基岩层；此外，假设碎石内和土壤内水势相等，将基质区土壤水分运动过程进一步修正为：其中，Km(h)为基质区内平均导水系数，由基质区内碎石和土壤的导水系数加权平均得到，引入碎石形状系数ε和质量比系数Rv，得到公式：Km(h)＝(1‑Rv)Kss(h)+εRvKsr(h)    (14)其中：式中：wm为基质区面积比例；h为土壤水吸力；Cms和Cmr分别为土壤和碎石的容水度；SS为源汇项；Γ为不同区间水量交换量；下标m表示基质区；下标i表示土壤层；Kss(h)为土壤非饱和导水系数；Ksr(h)为碎石非饱和导水系数；Kss为土壤饱和导水系数；Ksr为碎石饱和导水系数；α、vn和vm为参数，vm＝1‑1/vn；下标1和2分别表示土壤和碎石；z为坐标轴；土石二元混合介质达到饱和时，h＝0，基于公式(14)，计算单元内饱和导水系数为：Km＝(1‑Rv)·Kss+εRv·Ksr    (17)土壤水分过程计算公式表示如下：式中：下标f表示优先流区；wf为优先流区所占面积比例；上下边界条件均采用通量边界条件：式中：q为各层土壤上下界面处水分通量；基于Van Genuchten模型分别计算两者的土壤水分特征曲线：式中：Se为饱和度系数；θ为土石二元混合介质土壤含水量；θs为土石二元混合介质土壤饱和含水量；θr为土石二元混合介质残余含水量；α、vn和vm为参数，vm＝1‑1/vn；h为土壤水吸力；下标1和2分别表示土壤和碎石；壤中流计算：Rsub,i＝qm,i+qf,i        (25)qm,i＝wm(qm,i‑1‑qm,i)    (26)qf,i＝wf(qf,i‑1‑qf,i)      (27)式中：Q为壤中过流断面流量；Ks为土壤饱和导水系数；W为计算单元宽度；Φ为壤中流水位高度，即土水势；x和z为坐标轴；e为土壤内孔隙度；Rsub为垂向入流量；q为各层土壤上下界面处水分通量；下标m和f分别表示基质区和优先流区；下标i表示土壤层；地表汇流过程计算，计算公式如下：连续方法：运动方程：Sf＝S0       (29)曼宁公式：式中：Q0表示地表过流断面流量；A表示过流断面面积；Rsurf表示地表产流量；Sf表示摩擦坡降；S0表示计算单元平均地面坡降或河道坡降；Rwr表示过流断面水力半径；kn表示曼宁糙率系数；非降雨期土壤蒸散发过程计算，计算公式如下：E＝E1+E2+E3      (31)式中：E表示计算单元总蒸散发量；E1表示植被截留蒸发量；E2表示植被蒸腾量；E3表示裸地土壤蒸发量；土壤潜在蒸发量由Penman公式计算：式中：RN为净辐射量；G为传入水中的热通量；Δ为饱和水汽压对温度的导数；ρa为空气密度；Cp为空气的定压比热；δe为实际水蒸气压与饱和水蒸气压的差值；ra为蒸发表面空气动力学阻抗；λ为水的气化潜热；γ为空气湿度常数；PR为大气压；空气动力学阻抗计算：其计算公式如下：式中：ra为蒸发表面空气动力学阻抗；hz为气象站观测点离地面的高度；hd为置换高度；zom表为水蒸气紊流扩散对应的地表粗度；zox为地表粗度；κ为von Karman常数；U为风速；hc为植被高度；植被截留蒸发量E1使用Noilhan‑Planton模型计算：E1＝Veg·δ·Ep     (36)δ＝(Wr/Wrmax)2/3       (39)Wrmax＝0.2·Veg·LAI       (40)式中：Veg为裸地‑植被域中植被的面积占计算单元的面积比例；δ为湿润叶面占植被叶面的面积比例；Ep为土壤潜在蒸发量，即最大蒸发量；Wr为植被截留水量；Wrmax为植被最大截留水量；I为时段tk内最大降水平均雨强；Rr为植被冠层流出水量，即超出最大植被截留水量的部分；LAI表示叶面积指数；植被蒸腾量E2采用Penman‑Monteith公式计算，土壤各层蒸散量采用雷志栋的根系吸水模型进行计算，公式如下：E2＝Veg·(1‑δ)·EPM       (41)式中：EPM为采用Penman‑Monteith公式计算的潜在蒸腾量；G为传入水中的热通量；rc为植被群落阻抗；根系吸水模型：Tr＝E2     (44)式中：lr表示根系层厚度；Tr为实际植被蒸腾量；SS为源汇项；植被群落阻抗计算：采用Dickinson公式计算植被群落阻抗：σ1‑1＝1‑0.0016(25‑Ta)2      (46)σ2‑1＝1‑VPD/VPDc       (47)式中：rc为植被群落阻抗；rs min为最小气孔阻抗；LAI为叶面积指数；σ1为温度影响函数；σ2为大气水蒸气压饱和差影响函数；σ3为光合作用有效放射的影响函数；σ4为土壤含水量的影响函数；Ta为气温；VPD为饱和水蒸气压同实测值之间的差；VPDc为气孔闭合时的VPD值；rsmax为最大气孔阻抗；PAR为光和作用有效放射；PARc为PAR的临界值；θ3为根系层土壤含水量；θw为植被凋萎时的土壤含水量；θs为土壤饱和含水量；裸地土壤蒸发量E3由下式计算：式中：β为土壤湿润函数；θ4为表层土壤含水量；θm为土壤分子吸力对应的土壤含水量；θh为表层土壤田间持水量；非降雨期土壤水分运动过程采用修正的Richards公式计算，计算公式同非降雨期；当非降雨期地表有地下水流出时，采用动力波方程进行汇流计算，计算公式同降雨期；非降雨期壤中流采用改进后的动力波方程，计算公式同降雨期。