一种石灰岩土石山区山坡尺度水文过程模拟方法

技术领域

本发明涉及一种石灰岩土石山区山坡尺度水文过程模拟方法，是一种水文模拟方法，是一种针对石灰岩山区的水循环过程的模拟方法。

背景技术

研究石灰岩区山坡水文过程，可为石灰岩区水资源评价、管理和水土流失治理提供基础。石灰岩土石山区山坡水文过程异于平原区。石灰岩地区山坡土层较薄，一般可以分为均质土壤层、土石二元混合介质层、风化碎石层和基岩层。石灰岩导水性和持水性均较差，当土壤含有石灰岩碎石时，其碎石可以被认为非透水区。由于石灰岩具有可蚀性，土壤下部下覆巨厚基岩经常年累月水蚀，广泛发育着地下暗河或优先流通道。同时，基岩表面凹凸不平，这些凹凸面是储存水分的地下水库。石灰岩土石山区水文特征异常复杂尤其是在北方干旱半干旱地区。当前，山坡水文模型均未考虑石灰岩地区土石水文特征对山坡水文过程影响，导致模型出现模拟失真、精度不高等问题。

发明内容

为了克服现有技术的问题，本发明提出了一种石灰岩土石山区山坡尺度水文过程模拟方法。本发明考虑了石灰岩山区广泛存在土石二元混合介质、大孔隙优先流和基岩凹凸面，一方面引入碎石形状系数ε和碎石质量比系数R_v量化了土石二元混合介质对山坡水循环过程的影响；另一方面基于山坡计算单元划分，将计算单元分为大孔隙优先流区和基质区，进行两流区水循环过程计算；同时，基于改进的洼地储留模型，计算了基岩凹凸面对水文过程影响。

本发明的目的是这样实现的：一种石灰岩地区山坡尺度水文过程模拟方法，所述方法的计算过程如下：

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

计算单元内垂直剖面划分：根据山坡植被、土壤和岩石特性，在基本计算单元内分为4层：植被冠层截留层、地表储留层、土壤层、基岩层；植被截留层又可细分为：高植被截留层和矮植被储留层、草地层和裸地；土壤层进一步分为均质土壤层、土石二元混合介质层、风化碎石层；基岩层进一步分为凹凸面储留层区和裂隙优先流区；

计算单元内状态变量包括：植被冠层截留量、洼地储留量、枯枝落叶储留量、土壤含水量；主要参数包括：植被最大截留深、洼地最大储留深、基岩凹凸面最大储留深、枯枝落叶干重、土壤导水系数、土壤水分特征曲线、土壤含水量、土壤厚度、土石二元介质碎石质量比系数、坡面糙率等。计算单元水文过程计算：包括水文气象数据展布、降雨期间产汇流过程计算、非降雨期间产汇流过程计算等。

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

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

（1）

（2）

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

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

（3）

（4）

（5）

式中：I为时段t_k内最大降水平均雨强；S表示暴雨参数；t为时间>t_k-1<t≤t_k)；t_k为时段()区间的时间，N为时段数；表示暴雨衰减系数；P表示日降雨量；T表示日降雨总历时；a，b表示参数。

降雨期间产汇流过程计算：降雨期间，土壤蒸散发量较小，可以忽略。计算单元内水文过程主要由降雨→植被截留→入渗产流→汇流过程构成。

植被截留计算：

植被截留计算，计算公式如下：

（6）

（7）

（8）

式中：Veg表示植被的面积；Wr表示植被截留水量；Wr_max表示最大植被截留水量；I表示雨强；Rr表示植被冠层流出水量；LAI表示叶面积指数。

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

（9）

（10）

式中：为净雨强（经植被截留后的雨强）； H_u2为地表储留；H_umax2为土壤表层最大储留深；R_u2为土壤表面径流；f_in为入渗率。地表储留量分别为洼地储留量和枯枝落叶储留量构成，枯枝落叶储留量消耗于后期蒸散发过程，而洼地储留消耗于土壤入渗。

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

H_umax=z_maxG>

式中：G为枯枝落叶干重；z为枯枝落叶最大持水系数。

基岩凹凸面储留计算：

忽略降雨期间蒸散发条件下，基岩凹凸面储留采用下式计算：

（12）

（13）

式中：H_u1为基岩凹凸面储留；H_umax1为基岩凹凸面最大储留深；K_scr为风化壳导水系数；R_u1为基岩表面径流；L为计算单元长度；q为上层土壤下边界水分通量。降雨期间，基岩面上方下渗水分首先填充基岩凹凸面，然后产生基岩面壤中流；非降雨期，基岩凹凸面储留量主要消耗于植被蒸腾过程。

土壤及基岩层水分动态过程计算：土壤层水分运动过程采用理查兹（Richards equation）公式计算：

（14）

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

对于石灰岩地区，其基质区土壤内含有大量石灰岩碎石，属于土石二元混合介质，但其中的碎石不具有持水性和导水性，影响了土壤水分运动过程。在对其基质区土壤水分动态过程模拟时，引入碎石质量比系数R_v描述石灰岩碎石对土壤含水量的影响。根据R_v大小，可以将土壤层分为以下土层：1）当Rv=0时，土壤层为均质土壤；2）当1>R_v>0时，土壤层为土石二元混合介质层；3）当R_v=1时，土壤层为基岩层。则基质区土壤水分运动过程进一步修正为：

（15）

其中，为基质区内导水系数，引入碎石形状系数ε和碎石质量比系数R_v，得到公式（16）。

（16）

其中：

（17）

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

土石二元混合介质达到饱和时，h=0，基于公式（16），计算单元内饱和导水系数为：

K_m=>εR_v)·K_ss>

石灰岩地区计算单元内优先流区同样含有大量石灰岩碎石，但孔隙结构与基质区不同，导致水势、导水系数等有所不同。其土壤水分过程计算公式表示如下：

（19）

式中：下标f表示优先流区；w_f为优先流区面积比例；其他参数同前。

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

（20）

（21）

式中：q为水分通量；下标m和f分别表示基质区和大孔隙优先流区；下标i表示土壤层。降雨期间，上边界通量为净雨强；非降雨期间为表层土壤蒸发量通量。下边界为石灰岩基岩面，由于石灰岩透水性差，可以认为是不透水层，但石灰岩基岩内分布的裂隙具有一定的导水性和持水性。因此，下边界大部分地区水分通量为0，只有小部分地区具有水分通量。

壤中流计算：

（22）

（23）

（24）

（25）

（26）

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

模型计算中，基质区和大孔隙优先流区土壤水分特征曲线采用相同曲线。考虑到土石二元混合介质中含有大量碎石，改变了土壤水分特征参数。模型计算中，土壤水分特征曲线采用考虑碎石Van Genuchten模型计算：

（27）

（28）

式中：S_e为饱和度系数；R_v为碎石碎石质量比系数；θ为土石二元混合介质土壤含水量；θ_s为土石二元混合介质土壤饱和含水量；θ_r为土石二元混合介质土壤残余含水量；θ_n为无碎石土壤含水量；θ_ns为无碎石土壤饱和含水量；θ_nr为无碎石土壤残余含水量；α、vn和vm为参数，vm=1-1/vn；h为土壤水吸力。

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

连续方法：

（29）

运动方程：

S_f=S_0>>

曼宁公式：

（31）

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

非降雨期间产汇流过程计算：非降雨期间，土壤水主要消耗于水分再分布过程，包括土壤蒸散发和深层渗漏过程等。

土壤蒸散发过程计算：

计算单元内的土壤蒸散发量是植被截留蒸发、土壤蒸发和植被蒸腾的之和，计算公式如下：

（32）

式中：表示计算单元总蒸散发（mm）；下标1表示植被截留蒸发；下标2表示植被蒸腾；下标3表示裸土蒸发。

土壤潜在蒸发能力由Penman公式计算（最大蒸发强度）：

（33）

（34）

式中：为净辐射量；为传入水中的热通量；为饱和水汽压对温度的导数；为空气密度；为空气的定压比热；为实际水蒸气压与饱和水蒸气压的差值；为蒸发表面空气动力学阻抗；为水的气化潜热；为空气湿度常数；为大气压。

空气动力学阻抗计算：其计算公式如下：

（35）

（36）

式中：为空气动力学阻抗（s/m）；为气象站观测点离地面的高度（m）；为置换高度（m）；表为水蒸气紊流扩散对应的地表粗度（m）；为地表粗度；为von Karman 常数；为风速；为植被高度。

植被截留蒸发量（）使用Noilhan-Planton模型计算：

（37）

（38）

（39）

（40）

（41）

式中：为裸地-植被域中植被的面积占计算单元的面积比例；为湿润叶面占植被叶面的面积比例；E_p为潜在蒸发量，即最大蒸发量；Wr为植被截留量；Wr_max为植被最大截留水量；I为雨强；Rr>LAI表示叶面积指数。

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

（42）

（43）

式中：T_r为实际植被蒸腾量（mm）；为采用Penman-Monteith公式计算的潜在蒸腾量；为热通量；为植被阻抗（s/m）。

根系吸水模型：

（44）

T_r=E₂>

式中：lr表示根系层厚度；T_r为植被蒸腾量；SS为根系吸水。

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

（46）

（47）

（48）

（49）

（50）

式中：r_c为植被群落阻抗（s/m）；r_smin为最小气孔阻抗（s/m）；LAI为叶面积指数；为温度影响函数；为大气水蒸气压饱和差影响函数；为光合作用有效放射的影响函数；为土壤含水量的影响函数；为气温（℃）；为饱和水蒸气压同实测值之间的差（kPa）；为气孔闭合时的值（大约等于4kPa）；r_smax为最大气孔阻抗（5000s/m）；为光和作用有效放射（W/m²）；为的临界值（高植被：30W/m²；低植被：100W/m²）；为根系层土壤含水量（如不加说明，本文所有土壤含水量均指体积含水量）；为植被凋萎时的土壤含水量；为饱和土壤含水量。

裸地土壤蒸发量（）由下式计算：

（51）

（52）

式中：为土壤湿润函数；为表层土壤含水量；为土壤分子吸力（约1000-10000个大气压）对应的土壤含水量；为表层土壤田间持水量；其他参数意义同前。

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

本发明产生的有益效果是：计算单元采用等流时线法划分，计算单元内根据山坡剖面特征，在垂向上依次化分为植被层、土壤层和基岩层，每一层内根据其特性进一步分为亚层。计算单元内水分通量包括降雨、植被截留、洼地储留、入渗及水分动态过程、壤中流、地表流和蒸散发过程等方面。降雨期产汇流过程模拟中，主要对降雨数据在时间和空间上进行离散；植被截留是叶面积指数的函数；考虑到山区土壤入渗存在优先流，土壤降雨入渗过程计算中，将山区划分为大孔隙优先流区和基质流区，基质流和优先流土壤入渗及水分再分布均采用理查兹方程（Richards）计算；土壤壤中流采用改进后的动力波方程（Saint-Venant equation）计算；地表汇流采用动力波方程（Saint-Venant equation）计算。非降雨期间，土壤水分消耗于土壤蒸散发和水分再分布，土壤蒸散发包括植被截留蒸发、植被蒸腾和裸土蒸发，计算中分别计算；土壤水分动态过程、壤中流过程和地表汇流过程均同于降雨期间。本发明使石灰岩分布区土石山坡水文过程模拟更加精确，更接近实际。

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1是本发明的实施例一所述方法的流程图。

具体实施方式

实施例一：

本实施例是一种石灰岩土石山区山坡尺度水文过程模拟方法，流程见图1所示。本实施例所述方法的计算过程如下：

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

计算单元内垂直剖面划分：根据山坡植被、土壤和岩石特性，在基本计算单元内（等高带）分为4层，植被冠层截留层、地表储留层（枯枝落叶储留层和洼地储留层）、土壤层、基岩层。土壤层进一步分为均质土壤层、土石二元混合介质层、风化碎石层和基岩层。计算单元内状态变量包括：植被冠层截留量、洼地储留量、枯枝落叶储留量、土壤含水量等。主要参数包括：植被最大截留深、洼地最大储留深、基岩凹凸面最大储留深、枯枝落叶干重、土壤导水系数、土壤水分特征曲线、土壤含水量、各层土壤厚度、土石二元介质碎石质量比系数、坡面糙率等。植被截留层又可细分为：高植被截留层（分为针叶林、针阔混交林、阔叶林和常绿阔叶林）和矮植被储留层（灌木林）、草地层和裸地。

计算单元内状态变量包括：植被冠层截留量、洼地储留量、基岩凹凸面储留量、枯枝落叶储留量、土壤含水量；主要参数包括：植被最大截留深、洼地最大储留深、基岩凹凸面最大储留深、枯枝落叶干重、土壤导水系数、土壤水分特征曲线、土壤含水量、土壤厚度、土石二元介质碎石质量比系数、坡面糙率等。计算单元水文过程计算：包括水文气象数据展布、降雨期间产汇流过程计算、非降雨期间产汇流过程计算等。

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

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

（1）

（2）

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

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

（3）

（4）

（5）

式中：I为时段t_k内最大降水平均雨强；S表示暴雨参数；t为时间>t_k-1<t≤t_k)；t_k为时段()区间的时间，N为时段数；表示暴雨衰减系数；P表示日降雨量；T表示日降雨总历时；a，b表示参数。

降雨期间产汇流过程计算：降雨期间，土壤蒸散发量较小，可以忽略。计算单元内水文过程主要由降雨→植被截留→入渗产流→汇流过程构成。

植被截留计算：

植被截留计算，计算公式如下：

（6）

（7）

（8）

式中：Veg表示植被的面积；Wr表示植被截留水量；Wr_max表示最大植被截留水量；I表示雨强；Rr表示植被冠层流出水量；LAI表示叶面积指数。

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

（9）

（10）

式中：为净雨强（经植被截留后的雨强）；H_u2为地表储留；H_umax2为土壤表层最大储留深；R_u2为土壤表面径流；f_in为入渗率。地表储留量分别为洼地储留量和枯枝落叶储留量构成，枯枝落叶储留量消耗于后期蒸散发过程，而洼地储留消耗于土壤入渗。

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

H_umax=z_maxG>

式中：G为枯枝落叶干重；z为枯枝落叶最大持水系数。

基岩凹凸面储留计算：

忽略降雨期间蒸散发条件下，基岩凹凸面储留采用下式计算：

（12）

（13）

式中：H_u1为基岩凹凸面储留；H_umax1为基岩凹凸面最大储留深；K_scr为风化壳导水系数；R_u1为基岩表面径流；L为计算单元长度；q为上层土壤下边界水分通量。降雨期间，基岩面上方下渗水分首先填充基岩凹凸面，然后产生基岩面壤中流；非降雨期，基岩凹凸面储留量主要消耗于植被蒸腾过程。

土壤及基岩层水分动态过程计算：土壤层水分运动过程采用理查兹（Richards equation）公式计算：

（14）

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

对于石灰岩地区，其基质区土壤内含有大量石灰岩碎石，属于土石二元混合介质，但其中的碎石不具有持水性和导水性，影响了土壤水分动态过程。在对其基质区土壤水分动态过程模拟时，引入碎石质量比系数R_v描述石灰岩碎石对土壤含水量的影响。根据R_v大小，可以将土壤层分为以下土层：1）当R_v=0时，土壤层为均质土壤；2）当1>R_v>0时，土壤层为土石二元混合介质层；3）当R_v=1时，土壤层为基岩层。则基质区土壤水分运动过程进一步修正为：

（15）

其中，为基质区内导水系数，引入碎石形状系数ε和碎石质量比系数R_v，得到公式（16）。

（16）

其中：

（17）

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

土石二元混合介质达到饱和时，h=0，基于公式（16），计算单元内饱和导水系数为：

K_m=>εR_v)·K_ss（18）

石灰岩地区计算单元内优先流区同样含有大量石灰岩碎石，但孔隙结构与基质区不同，导致水势、导水系数等有所不同。其土壤水分过程计算公式表示如下：

（19）

式中：下标f表示优先流区；w_f为优先流区面积比例；其他参数同前。

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

（20）

（21）

式中：q为水分通量；下标m和f分别表示基质区和大孔隙优先流区；下标i表示土壤层。降雨期间，上边界通量为净雨强；非降雨期间为表层土壤蒸发量通量。下边界为石灰岩基岩面，由于石灰岩透水性差，可以认为是不透水层，但石灰岩基岩内分布的裂隙具有一定的导水性和持水性。因此，下边界大部分地区水分通量为0，只有小部分地区具有水分通量。

壤中流计算：

（22）

（23）

（24）

（25）

（26）

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

模型计算中，基质区和大孔隙优先流区土壤水分特征曲线采用相同曲线。考虑到土石二元混合介质中含有大量碎石，改变了土壤水分特征参数。模型计算中，土壤水分特征曲线采用考虑碎石Van Genuchten模型计算：

（27）

（28）

式中：S_e为饱和度系数；R_v为碎石碎石质量比系数；θ为土石二元混合介质土壤含水量；θ_s为土石二元混合介质土壤饱和含水量；θ_r为土石二元混合介质土壤残余含水量；θ_n为无碎石土壤含水量；θ_ns为无碎石土壤饱和含水量；θ_nr为无碎石土壤残余含水量；α、vn和vm为参数，vm=1-1/vn；h为土壤水吸力。

碎石质量比系数地表汇流过程计算：计算公式如下：

连续方法：

（29）

运动方程：

S_f=S₀>

曼宁公式：

（31）

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

非降雨期间产汇流过程计算：非降雨期间，土壤水主要消耗于水分进行再分布过程，包括土壤蒸散发和深层渗漏过程等。

土壤蒸散发过程计算：

计算单元内的土壤蒸散发量是植被截留蒸发、土壤蒸发和植被蒸腾的之和，计算公式如下：

（32）

式中：表示计算单元总蒸散发（mm）；下标1表示植被截留蒸发；下标2表示植被蒸腾；下标3表示裸土蒸发。

土壤潜在蒸发能力由Penman公式计算（最大蒸发强度）：

（33）

（34）

式中：为净辐射量；为传入水中的热通量；为饱和水汽压对温度的导数；为空气密度；为空气的定压比热；为实际水蒸气压与饱和水蒸气压的差值；为蒸发表面空气动力学阻抗；为水的气化潜热；为空气湿度常数；为大气压。

空气动力学阻抗计算：其计算公式如下：

（35）

（36）

式中：为空气动力学阻抗（s/m）；为气象站观测点离地面的高度（m）；为置换高度（m）；表为水蒸气紊流扩散对应的地表粗度（m）；为地表粗度；为von Karman 常数；为风速；为植被高度。

植被截留蒸发量（）使用Noilhan-Planton模型计算：

（37）

（38）

（39）

（40）

（41）

式中：为裸地-植被域中植被的面积占计算单元的面积比例；为湿润叶面占植被叶面的面积比例；E_p为潜在蒸发量，即最大蒸发量；Wr为植被截留量；Wr_max为植被最大截留水量；I为雨强；Rr为植被冠层流出水量，即超出最大植被截留水量的部分（mm）；LAI表示叶面积指数。

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

（42）

（43）

式中：T_r为实际植被蒸腾量（mm）；>

根系吸水模型：

（44）

T_r=E₂>

式中：lr表示根系层厚度；T_r为植被蒸腾量；SS为根系吸水。

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

（46）

（47）

（48）

（49）

（50）

式中：r_c为植被群落阻抗（s/m）；r_smin为最小气孔阻抗（s/m）；LAI为叶面积指数；为温度影响函数；为大气水蒸气压饱和差影响函数；为光合作用有效放射的影响函数；为土壤含水量的影响函数；为气温（℃）；为饱和水蒸气压同实测值之间的差（kPa）；为气孔闭合时的值（大约等于4kPa）；r_smax为最大气孔阻抗（5000s/m）；为光和作用有效放射（W/m²）；为的临界值（高植被：30W/m²；低植被：100W/m²）；为根系层土壤含水量（如不加说明，本文所有土壤含水量均指体积含水量）；为植被凋萎时的土壤含水量；为饱和土壤含水量。

裸地土壤蒸发量（）由下式计算：

（51）

（52）

式中：为土壤湿润函数；为表层土壤含水量；为土壤分子吸力（约1000-10000个大气压）对应的土壤含水量；为表层土壤田间持水量；其他参数意义同前。

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

最后应说明的是，以上仅用以说明本发明的技术方案而非限制，尽管参照较佳布置方案对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案（比如数据的采集、各种公式的运用、计算的先后顺序等）进行修改或者等同替换，而不脱离本发明技术方案的精神和范围。