本发明涉及一种片麻岩土石山区山坡尺度水文过程模拟方法，是一种水文模拟方法，是一种针对片麻岩山区的水循环过程的模拟方法。
背景技术：
分析片麻岩土石山区山坡水文过程内在关系，实现山坡水文过程模拟，是认识片麻岩土石山区山坡水资源评价、管理、水土保持等领域的基础。片麻岩土石山区山坡水文过程异于平原区，也异于石灰岩地区、花岗岩地区等坚硬基岩分布区域山坡水文过程。片麻岩土石山区山坡地区土层较薄、土壤存在分层且土壤内广泛分布着土石二元混合介质；土壤下部下覆巨厚片麻岩基岩，片麻岩基岩具有很强导水性和持水性（相较于花岗岩、石灰岩等坚硬岩石），导致片麻岩土石山区山坡在连续降雨或强降雨条件下，入渗、产流、蒸散发和土壤水分运动等过程异于平原地区，异于石灰岩等坚硬岩石分布区。当前，山坡水文模型均未考虑片麻岩土石山区水文特征对山坡水文过程影响，导致模型出现模拟失真、精度不高等问题。
技术实现要素：
为了克服现有技术的问题，本发明提出了一种片麻岩土石山区山坡尺度水文过程模拟方法。考虑到土壤内存在大量优先流孔隙，将土壤层分为基质流区和优先流区，两区土壤入渗及水分再分布均采用改进的理查兹方程（Richards equation）计算，土壤壤中流采用改进后的动力波方程（Saint-Venant equation）计算，地表汇流采用动力波方程（Saint-Venant equation）计算。鉴于土壤内含有大量片麻岩分化碎石，基质流区计算中，引入碎石质量比系数Rv来描述碎石对土壤含水量的影响；引入碎石质量比系数Rv和形状系数ε来描述碎石对土壤导水系数的影响。
本发明的目的是这样实现的：一种片麻岩土石山区山坡尺度水文过程模拟方法，所述方法的计算过程如下：
山坡计算单元划分：采用等流时线法将山坡划分为若干基本计算单元；
计算单元垂直剖面划分：根据山坡植被、土壤和岩石特性，在基本计算单元内分为4层：植被冠层截留层、地表储留层、土壤层和基岩层；植被截留层又可细分为：高植被截留层和矮植被储留层、草地层和裸地；土壤层进一步分为均质土壤层、土石二元混合介质层；考虑到片麻岩基岩层内广泛分布着构造节理，且片麻岩本身具有一定持水性和透水性，进一步将基岩层分为基质流区和优先流区。
计算单元内状态变量包括：植被冠层截留量、洼地储留量、枯枝落叶储留量、土壤含水量；主要参数包括：植被最大截留深、洼地最大储留深、枯枝落叶干重、土壤导水系数、土壤水分特征曲线、土壤含水量、各层土壤厚度、土石二元混合介质碎石质量比系数、基岩层厚度、坡面糙率；
计算单元水文过程计算，包括水文气象数据展布、降雨期间产汇流过程计算、非降雨期间产汇流过程计算等。
水文气象数据展布，包括水文气象过程空间尺度展布和降雨时间降尺度展布：
采用泰森多边形法和反距离加权平均法进行流域内气象数据的空间展布，包括降雨、气温、风速、空气湿度、净辐射，计算公式如下：
（1）
（2）
式中：D表示待插值点估计值；Dpi表示第pi个参证站点数据；pm表示参证站点个数；λpi表示第pi个参证站点数据权重；dpi表示第pi个参证站点同待插值点的距离；pn表示权重指数；
由于日降雨过程的非稳定性，对日降雨数据进一步进行降尺度展布，具体公式如下：
（3）
（4）
（5）
式中：I为时段tk内最大降水平均雨强；S表示暴雨参数；t为时间 (tk-1<t≤tk)；tk为时段()区间的时间，N为时段数；表示暴雨衰减系数；P表示日降雨量；T表示日降雨总历时；a，b表示参数。
降雨期间产汇流过程计算：降雨期间，土壤蒸散发量较小，可以忽略；计算单元内水文过程主要由降雨→植被截留→入渗产流→汇流过程构成；
植被截留计算，计算公式如下：
（6）
（7）
（8）
式中：Veg表示植被的面积；Wr表示植被截留水量；表示最大植被截留水量；I为时段tk内雨强；Rr表示植被冠层流出水量；LAI表示叶面积指数。
洼地储留计算，计算公式如下：
（9）
（10）
式中：为时段tk内净雨强（经植被截留后的雨强）； Hu2为地表储留；Humax2为土壤表层最大储留深；Ru2为土壤表面径流；fin为入渗率。地表储留量分别为洼地储留量和枯枝落叶储留量构成，枯枝落叶储留量消耗于后期蒸散发过程，而洼地储留消耗于土壤入渗。
枯枝落叶储留计算，计算公式如下：
Humax=zmaxG （11）
式中：G为枯枝落叶干重；z为枯枝落叶最大持水系数；
土壤及基岩层水分运动过程计算：片麻岩分布区土壤和基岩均具有导水性和持水性，所有土壤层水分运动过程均采用理查兹（Richards equation）公式计算：
（12）
式中：h为土壤水吸力；C为容水度；SS为源汇项（即根系吸水）；z为坐标轴；K(h)为导水系数；t为时间。考虑到土壤和基岩层均广泛分布着裂隙优先流，土壤水分运动过程中将计算单元分为基质区和优先流区。其中计算单元内基质区所占面积比例为，优先流区所占面积比例为。
对于片麻岩土石山区，其基质区土壤内含有大量片麻岩碎石，属于土石二元混合介质，碎石具有一定持水性和导水性，影响了土壤水分运动过程。在对其基质区土壤水分运动过程模拟时，引入碎石质量比系数Rv描述片麻岩碎石对土壤含水量的影响。根据Rv大小，可以将土壤层分为以下土层：1）当Rv=0时，土壤层为均质土壤；2）当1>Rv>0时，土壤层为土石二元混合介质层；3）当Rv=1时，土壤层为基岩层。此外，假设碎石内和土壤内水势相等，可将基质区土壤水分运动过程进一步修正为：
（13）
式中：为基质区内平均导水系数，由基质区内碎石和土壤的导水系数加权平均得到，引入碎石形状系数ε和体积系数Rv，得到公式（14）。
（14）
其中：
（15）
（16）
式中：为基质区面积比例；h为土壤水吸力；和分别为土壤和碎石的容水度；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）
片麻岩土石山区计算单元内优先流区同样含有大量片麻岩碎石，但孔隙结构与基质区不同，导致水势、导水系数等有所不同。其土壤水分过程计算公式表示如下：
（18）
式中：下标f表示优先流区；wf为优先流区面积比例；其他参数同前。
上下边界条件均采用通量边界条件：
（19）
（20）
式中：q为上下边界通量；其它参数同前。
模型计算中，基质区和优先流区土壤水分特征曲线采用相同曲线。考虑到土石二元混合介质中含有大量片麻岩碎石，且碎石具有一定持水性和导水性，进而改变了土壤水分特征参数。考虑到土壤和碎石的土壤水分特征参数变化的差异，基于Van Genuchten模型分别计算两者的土壤水分特征曲线。
（21）
（22）
式中：Se为饱和度系数；θ为土石二元混合介质土壤含水量；θs为土石二元混合介质土壤饱和含水量；θr为土石二元混合介质残余含水量；α、vn和vm为参数，vm=1-1/vn；h为土壤水吸力；下标1和2分别表示土壤和碎石。
壤中流计算：
（23）
（24）
（25）
（26）
（27）
式中：Q为壤中过流断面流量；Ks为土壤饱和导水系数；W为计算单元宽度；为壤中流水位高度，即土水势；x和z为坐标轴；e为土壤内孔隙度；为垂向入流量（即壤中流产流量）；为各层土壤上下界面处水分通量；下标m和f分别表示基质区和优先流；下标i表示土壤层。
地表汇流过程计算：计算公式如下：
连续方法：
（28）
运动方程：
Sf=S0 （29）
曼宁公式：
（30）
式中：Q0表示地表过流断面流量；A表示过流断面面积；Rsurf表示地表产流量；Sf表示摩擦坡降；S0表示计算单元平均地面坡降或河道坡降；Rwr表示过流断面水力半径；kn表示曼宁糙率系数；
非降雨期土壤蒸散发过程计算步骤如下：
计算单元内的土壤蒸散发量是植被截留蒸发、土壤蒸发和植被蒸腾的之和，计算公式如下：
（31）
式中：表示计算单元总蒸散发（mm）；下标1表示植被截留蒸发；下标2表示植被蒸腾；下标3表示裸土蒸发。
土壤潜在蒸发能力由Penman公式计算（最大蒸发强度）：
（32）
（33）
式中：为净辐射量；为传入水中的热通量；为饱和水汽压对温度的导数；为空气密度；为空气的定压比热；为实际水蒸气压与饱和水蒸气压的差值；为蒸发表面空气动力学阻抗；为水的气化潜热；为空气湿度常数；为大气压。
空气动力学阻抗计算：其计算公式如下：
（34）
（35）
式中：为空气动力学阻抗（s/m）；为气象站观测点离地面的高度（m）；为置换高度（m）；表为水蒸气紊流扩散对应的地表粗度（m）；为地表粗度；为von Karman 常数；为风速；为植被高度。
植被截留蒸发量（）使用Noilhan-Planton模型计算：
（36）
（37）
（38）
（39）
（40）
式中：为裸地-植被域中植被的面积占计算单元的面积比例；为湿润叶面占植被叶面的面积比例；Ep为潜在蒸发量，即最大蒸发量；Wr为植被截留量；Wrmax为植被最大截留水量；I为雨强；Rr为植被冠层流出水量，即超出最大植被截留水量的部分（mm）；LAI表示叶面积指数。
植被蒸腾量（）采用Penman-Monteith公式计算。土壤各层蒸散量采用雷志栋的根系吸水模型进行计算。具体公式如下：
（41）
（42）
式中：Tr为实际植被蒸腾量（mm）；为采用Penman-Monteith公式计算的潜在蒸腾量；为热通量；为植被阻抗（s/m）。
根系吸水模型：
（43）
Tr=E2 （44）
式中：lr表示根系层厚度；Tr为植被蒸腾量；SS为源汇项（即根系吸水）。
植被群落阻抗计算：采用Dickinson等提出公式计算植被群落阻抗：
（45）
（46）
（47）
（48）
（49）
式中：rc为植被群落阻抗（s/m）；rs min为最小气孔阻抗（s/m）；LAI为叶面积指数；为温度影响函数；为大气水蒸气压饱和差影响函数；为光合作用有效放射的影响函数；为土壤含水量的影响函数；为气温（℃）；为饱和水蒸气压同实测值之间的差（kPa）；为气孔闭合时的值（大约等于4kPa）；rs max为最大气孔阻抗（5000 s/m）；为光和作用有效放射（W/m2）；为的临界值（高植被：30W/m2；低植被：100W/m2）；为根系层土壤含水量（如不加说明，本文所有土壤含水量均指体积含水量）；为植被凋萎时的土壤含水量；为土壤饱和含水量。
裸地土壤蒸发量（）由下式计算：
（50）
（51）
式中：为土壤湿润函数；为表层土壤含水量；为土壤分子吸力（约1000-10000个大气压）对应的土壤含水量；为表层土壤田间持水量；其他参数意义同前。
非降雨期土壤水分运动过程采用修正的Richards公式计算，计算公式同非降雨期；当非降雨期地表有地下水流出时，采用动力波方程进行汇流计算，计算公式同降雨期；非降雨期壤中流采用改进后的动力波方程，计算公式同降雨期。
本发明产生的有益效果是：本发明的计算单元采用等流时线法划分，计算单元内根据山坡剖面特征，在垂向上依次划分为植被层、土壤层和基岩层，每一层内根据其特性进一步分为亚层。计算单元内水分通量包括降雨、植被截留、洼地储留、入渗及水分再分布过程、壤中流、地表流和蒸散发过程等方面。降雨过程模拟中，主要对降雨数据在时间和空间上进行离散；植被截留是叶面积指数的函数；考虑到山区土壤入渗存在优先流，土壤降雨入渗过程计算中，将山区划分为优先流区和基质流区。非降雨期间，土壤水分消耗于土壤蒸散发和水分再分布，土壤蒸散发包括植被截留蒸发、植被蒸腾和裸土蒸发，计算中分别计算；土壤水分运动过程、壤中流过程和地表汇流过程均同于降雨期间。本发明使片麻岩土质的山坡水文过程模拟更加精确，更接近实际。
附图说明
下面结合附图和实施例对本发明作进一步说明。
图1是本发明的实施例一所述方法的流程图。
具体实施方式
实施例一：
本实施例是一种片麻岩土石山区山坡尺度水文过程模拟方法，流程如图1所示。本实施例所述方法的计算过程如下：
山坡计算单元划分：采用等流时线法将山坡划分为若干基本计算单元。
计算单元内垂直剖面划分：根据山坡植被、土壤和岩石特性，在基本计算单元内（等高带）分为4层：植被冠层截留层、地表储留层（枯枝落叶储留层和洼地储留层）、土壤层、基岩层。土壤层进一步分为均质土壤层、土石二元混合介质层、风化碎石层和基岩层。计算单元内状态变量包括：植被冠层截留量、洼地储留量、枯枝落叶储留量、土壤含水量等。主要参数包括：植被最大截留深、洼地最大储留深、枯枝落叶干重、土壤导水系数、土壤水分特征曲线、土壤含水量、各层土壤厚度、土石二元介质碎石质量比系数、基岩层厚度、坡面糙率等。植被截留层又可细分为：高植被截留层（分为针叶林、针阔混交林、阔叶林和常绿阔叶林）和矮植被储留层（灌木林）、草地层和裸地。
计算单元内状态变量包括：植被冠层截留量、洼地储留量、枯枝落叶储留量、土壤含水量；主要参数包括：植被最大截留深、洼地最大储留深、枯枝落叶干重、土壤导水系数、土壤水分特征曲线、土壤含水量、各层土壤厚度、土石二元混合介质碎石质量比系数、基岩层厚度、坡面糙率。
计算单元水文过程计算：包括水文气象数据展布、降雨期间产汇流过程计算、非降雨期间产汇流过程计算等。
水文气象数据展布：包括水文气象过程空间尺度展布和降雨时间降尺度展布：
采用泰森多边形法和反距离加权平均法进行流域内气象数据的空间展布，包括降雨、气温、风速、空气湿度、净辐射等，计算公式如下：
（1）
（2）
式中：D表示待插值点估计值；Dpi表示第pi个参证站点数据；pm表示参证站点个数；λpi表示第pi个参证站点数据权重；dpi表示第pi个参证站点同待插值点的距离；pn表示权重指数。
由于日降雨过程的非稳定性，对日降雨数据进一步进行降尺度展布，具体公式如下：
（3）
（4）
（5）
式中：I为时段tk内最大降水平均雨强；S表示暴雨参数；t为时间 (tk-1<t≤tk)；tk为时段()区间的时间，N为时段数；表示暴雨衰减系数；P表示日降雨量；T表示日降雨总历时；a，b表示参数。
降雨期间产汇流过程计算：降雨期间，土壤蒸散发量较小，可以忽略。计算单元内水文过程主要由降雨→植被截留→入渗产流→汇流过程构成。
植被截留计算，计算公式如下：
（6）
（7）
（8）
式中：Veg表示植被的面积；Wr表示植被截留水量；Wrmax表示最大植被截留水量；I为时段tk内雨强；Rr表示植被冠层流出水量；LAI表示叶面积指数。
洼地储留计算，计算公式如下：
（9）
（10）
式中：为时段tk内净雨强（经植被截留后的雨强）；Hu2为地表储留；Humax2为土壤表层最大储留深；Ru2为土壤表面径流；fin为入渗率。地表储留量分别为洼地储留量和枯枝落叶储留量构成，枯枝落叶储留量消耗于后期蒸散发过程，而洼地储留消耗于土壤入渗。
枯枝落叶储留计算，计算公式如下：
Humax=zmaxG （11）
式中：G为枯枝落叶干重；z为枯枝落叶最大持水系数。
土壤及基岩层水分运动过程计算：
片麻岩分布区土壤和基岩均具有导水性和持水性，所有土壤层水分运动过程均采用理查兹（Richards）公式计算，公式如下：
（12）
式中：h为土壤水吸力；C为容水度；SS为源汇项（即根系吸水）；z为坐标轴；K(h)为导水系数；t为时间。考虑到土壤和基岩层均广泛分布着裂隙优先流，土壤水分运动过程中将计算单元分为基质区和优先流区。其中计算单元内基质区所占面积比例为，优先流区所占面积比例为。
对于片麻岩土石山区，其基质区土壤内含有大量片麻岩碎石，属于土石二元混合介质，碎石具有一定持水性和导水性，影响了土壤水分运动过程。在对其基质区土壤水分运动过程模拟时，引入碎石质量比系数Rv描述片麻岩碎石对土壤含水量的影响。根据Rv大小，可以将土壤层分为以下土层：1）当Rv=0时，土壤层为均质土壤；2）当1>Rv>0时，土壤层为土石二元混合介质层；3）当Rv=1时，土壤层为基岩层。此外，假设碎石内和土壤内水势相等，可将基质区土壤水分运动过程进一步修正为：
（13）
式中：为基质区内平均导水系数，由基质区内碎石和土壤的导水系数加权平均得到，引入碎石形状系数ε和体积系数Rv，得到公式（14）。
（14）
其中：
（15）
（16）
式中：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）
片麻岩土石山区计算单元内优先流区同样含有大量片麻岩碎石，但孔隙结构与基质区不同，导致水势、导水系数等有所不同。其土壤水分过程计算公式表示如下：
（18）
式中：下标f表示优先流区；wf为优先流区面积比例；其他参数同前。
上下边界条件均采用通量边界条件：
（19）
（20）
式中：q为上下边界通量；其它参数同前。
模型计算中，基质区和优先流区土壤水分特征曲线采用相同曲线。考虑到土石二元混合介质中含有大量片麻岩碎石，且碎石具有一定持水性和导水性，进而改变了土壤水分特征参数。考虑到土壤和碎石的土壤水分特征参数变化的差异，基于Van Genuchten模型分别计算两者的土壤水分特征曲线。
（21）
（22）
式中：Se为饱和度系数；θ为土石二元混合介质土壤含水量；θs为土石二元混合介质土壤饱和含水量；θr为土石二元混合介质残余含水量；α、vn和vm为参数，vm=1-1/vn；h为土壤水吸力；下标1和2分别表示土壤和碎石。
壤中流计算：
（23）
（24）
（25）
（26）
（27）
式中：Q为壤中过流断面流量；Ks为饱和导水系数；W为计算单元宽度；Φ为壤中流水位高度，即土水势；x和z为坐标轴；e为土壤内孔隙度；Rsub为垂向入流量（即壤中流产流量）；q为各层土壤上下界面处水分通量；下标m和f分别表示基质区和优先流；下标i表示土壤层。
地表汇流过程计算：计算公式如下：
连续方法：
（28）
运动方程：
Sf=S0 （29）
曼宁公式：
（30）
式中：Q0表示地表过流断面流量；A表示过流断面面积；Rsurf表示地表产流量；Sf表示摩擦坡降；S0表示计算单元平均地面坡降或河道坡降；Rwr表示过流断面水力半径；kn表示曼宁糙率系数；
非降雨期土壤蒸散发过程计算。计算单元内的土壤蒸散发量是植被截留蒸发、土壤蒸发和植被蒸腾之和，计算公式如下：
（31）
式中：表示计算单元总蒸散发（mm）；下标1表示植被截留蒸发；下标2表示植被蒸腾；下标3表示裸土蒸发。
土壤潜在蒸发能力由Penman公式计算（最大蒸发强度）：
（32）
（33）
式中：为净辐射量；为传入水中的热通量；为饱和水汽压对温度的导数；为空气密度；为空气的定压比热；为实际水蒸气压与饱和水蒸气压的差值；为蒸发表面空气动力学阻抗；为水的气化潜热；为空气湿度常数；为大气压。
空气动力学阻抗计算：其计算公式如下：
（34）
（35）
式中：为空气动力学阻抗（s/m）；为气象站观测点离地面的高度（m）；为置换高度（m）；表为水蒸气紊流扩散对应的地表粗度（m）；为地表粗度；为von Karman 常数；为风速；为植被高度。
植被截留蒸发量（）使用Noilhan-Planton模型计算：
（36）
（37）
（38）
（39）
（40）
式中：为裸地-植被域中植被的面积占计算单元的面积比例；为湿润叶面占植被叶面的面积比例；Ep为潜在蒸发量，即最大蒸发量；Wr为植被截留量；Wrmax为植被最大截留水量；I为雨强；Rr 为植被冠层流出水量，即超出最大植被截留水量的部分（mm）；LAI表示叶面积指数。
植被蒸腾量（）采用Penman-Monteith公式计算。土壤各层蒸散量采用雷志栋的根系吸水模型进行计算。具体公式如下：
（41）
（42）
式中：Tr为实际植被蒸腾量（mm）；为采用Penman-Monteith公式计算的潜在蒸腾量；为热通量；为植被阻抗（s/m）。
根系吸水模型：
（43）
Tr=E2 （44）
式中：lr表示根系层厚度；Tr为植被蒸腾量；SS为源汇项（即根系吸水）。
植被群落阻抗计算：采用Dickinson等提出公式计算植被群落阻抗：
（45）
（46）
（47）
（48）
（49）
式中：rc为植被群落阻抗（s/m）；rs min为最小气孔阻抗（s/m）；LAI为叶面积指数；为温度影响函数；为大气水蒸气压饱和差影响函数；为光合作用有效放射的影响函数；为土壤含水量的影响函数；为气温（℃）；为饱和水蒸气压同实测值之间的差（kPa）；为气孔闭合时的值（大约等于4kPa）；rs max为最大气孔阻抗（5000s/m）；为光和作用有效放射（W/m2）；为的临界值（高植被：30W/m2；低植被：100W/m2）；为根系层土壤含水量（如不加说明，本文所有土壤含水量均指体积含水量）；为植被凋萎时的土壤含水量；为土壤饱和含水量。
裸地土壤蒸发量（）由下式计算：
（50）
（51）
式中：为土壤湿润函数；为表层土壤含水量；为土壤分子吸力（约1000-10000个大气压）对应的土壤含水量；为表层土壤田间持水量；其他参数意义同前。
非降雨期土壤水分运动过程采用修正的Richards公式计算，计算公式同非降雨期；当非降雨期地表有地下水流出时，采用动力波方程进行汇流计算，计算公式同降雨期；非降雨期壤中流采用改进后的动力波方程，计算公式同降雨期。
最后应说明的是，以上仅用以说明本发明的技术方案而非限制，尽管参照较佳布置方案对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案（比如数据的采集、各种公式的运用、步骤的先后顺序等）进行修改或者等同替换，而不脱离本发明技术方案的精神和范围。