第三节 蒸发条件下土壤水分运动

一、土壤水蒸发及其特点

土壤水的蒸发，发生在土壤的表层，其强度一般取决于两个因素：一是外界蒸发能力，即气象条件所限定的最大可能蒸发强度；二是土壤自下部土层向上的输水能力，其数值随含水率的降低而减小。据此，常将土壤蒸发分为三个阶段，如图1-8所示。蒸发开始时，由于土壤水较丰富，导水率大，在大气蒸发力作用下，表层土壤源源不断地从土体内部得到水分补给。这时土壤蒸发主要受大气蒸发力（常以水面蒸发量表示）的控制（Ⅰ）；随着土壤含水量的降低，地表与下面湿润土层的土壤水吸力梯度逐渐加大，以致蒸发速率愈来愈小，这时土壤蒸发进入第二个阶段（Ⅱ），其蒸发速率主要由土壤导水率控制；当表土含水率很低，土壤输水能力很弱，不能补充表土蒸发损失的水分，土壤表明形成干土层。土壤水分蒸发发生在干土层底部，然后以水汽扩散的方式穿过干土层而进入大气。在此阶段内，蒸发面不是在地表，而是在土壤内部。其蒸发强度的大小受干土层内水汽扩散能力控制，主要取决于干土层的厚度。该过程可持续数周或数月，是蒸发过程的最后阶段（Ⅲ）。

二、蒸发条件下土壤水分运动的定解问题

地表土壤蒸发属于垂直一维土壤水分运动，可用式（1-22）描述，并与初始条件和边界条件一起构成了定解问题。常见的初始条件和边界条件如下。

1.初始条件

有两种情况：

（1）土壤含水量沿剖面均匀分布，θ（z，0）=θ_i。

（2）土壤含水量沿剖面非均匀分布，θ（z，0）=θ（z）。

2.上边界条件

有四种情况：

（1）表土蒸发速率E₀。已知，且为常数。

（2）表土蒸发速率E已知，且考虑日变化，如按正弦周期变化。

式中：E_max为白天中午发生的最大蒸发率；t是以日出为起始的时间，s。

（3）表土蒸发率E随表土含水量的改变发生变化。

式中：θ_a、θ_k分别为风干含水率和临界含水率，cm³/cm³；a、b为待定系数。

（4）表土含水量一定。如在强烈蒸发作用下，表土很快降低到风干含水率θ_a，即

3.下边界条件

（1）对于半无限蒸发土柱，表土的蒸发不能影响到无穷深处的含水量，其下边界含水量不变，为初始含水量θ_i。

（2）若为有限土柱，如，底部为不透水层，即下边界处水分通量为零。

再如，底部为浅层地下水，地下水面处土壤基质势为0：

以上定解问题需用解析法或数值法求解，由于土壤水分运动的复杂性，实用中常采用数值法，但解析法有助于土壤水分运动机理的认识。以下对几种常见的解析法做介绍。

三、潜水稳定蒸发

在潜水埋深较浅，且相对稳定的条件下，潜水蒸发可看作是稳定蒸发，即地表处的蒸发强度与任一断面处的土壤水分通量相等，此通量即为潜水蒸发强度，如图1-9所示。其中z坐标原点在潜水面处，向上为正。

潜水蒸发定解问题可写为

对式（1-46a）积分，并利用条件式（1-46b），可得出

式中：θ_s为土壤饱和含水率，当非饱和土壤扩散率D（θ）和导水率K（θ）已知时，式（1-47）给出了稳定蒸发强度为E时土壤含水率分布z～θ。

若用土壤水吸力s表示未知函数，则有

对式（1-48a）积分，并利用条件式（1-48b），可得出：

式中：当非饱和土壤导水率K（s）已知时，式（1-49）给出了稳定蒸发强度为E时土壤含水率分布z～s。

令式（1-47）中的z等于潜水埋深，即z=H，可得到

由式（1-50）可计算确定不同潜水埋深条件下的蒸发强度E和θ_H，进而可给出以H为参数的E～θ_H关系。

如某重壤土，其饱和含水率θ_s=0.48cm³/cm³，非饱和土壤扩散率D（θ）=92.16（θ/θ_s）^2.391（cm²/h），导水率K（θ）=0.346（θ/θ_s）^9.133（cm/h）。可采用数值积分的方法由式（1-50）确定E～θ_H关系，见图1-10。其计算过程如下。

（1）针对某一潜水埋深（H=50cm），可假定一表土含水率，如θ_H=0.08，见表1-1，再假定一潜水蒸发速率E。

（2）采用数值方法对式（1-50）右端做积分计算。将区间［θ_H，θ_s］等分为n个小区间，相应的区间长度为，Δθ=（θ_s-θ_H）/n。

（3）采用梯形求积方法。令

式（1-47）可写为：

由式（1 52）可求得z，若z=H，则假定的E即为所求值；若z≠H，则需改变E，再做积分计算，直到z=H。数值积分是一种近似计算，计算误差是不可避免的，只要z和H的绝对值|z-H|小于给定误差ε（如ε=0.0001）即可。该计算过程实际上是一个优化计算问题，其决策变量为潜水蒸发强度E，目标函数为min|z-H|。可采用一般的优化方法对该问题求解，如可以采用Excel软件中规划求解功能计算。n的取值越大，计算结果越精确。对本例，取n=100，其计算结果见表11。如此进行可求得潜水埋深H=50cm时的潜水蒸发强度与表土含水率的关系曲线E～θ_H，见图1 10。

由表1-1可以看出，当表土含水率较小时，如小于0.22时，潜水蒸发量较大，可以达到16mm/d以上，但这仅仅是该种土壤在这一潜水埋深条件下可能的潜水蒸发速率。潜水蒸发速率还受大气蒸发力的控制，当大气蒸发力小于可能的潜水蒸发速率时，表土含水率会相应地增大，这就是为什么潜水埋深较小时，地表面常常潮湿的缘故。

按照同样方法，可求得其他潜水埋深条件下，潜水蒸发强度与表土含水率的关系曲线E～θ_H，见图1-10。

曲线与纵轴的交点，表示表土含水率θ_H=0时的潜水蒸发速率，称为该埋深下的潜水极限蒸发强度E_max，与潜水埋深的关系可用式（1-53）表示：

式中：A、m为随土壤而异的参数。

在分析潜水埋深H、潜水蒸发强度E和地表含水率θ_H关系的基础上，雷志栋（1983年）提出了如下计算潜水蒸发速率的经验公式

式中：E为土壤蒸发强度，mm/d；E₀为水面蒸发强度，mm/d；η为经验常数，与土质和地下水位埋深H有关。由式（1-53）可知，式（1-54）隐含了潜水埋深H，表示了表土蒸发量E、潜水埋深H和大气蒸发力E₀之间的关系。

对于潜水蒸发强度与潜水埋深的关系，国内外广泛采用的经验公式还有式（1-55）：

式中：H_max为极限蒸发深度，即当地下水位埋深H≥H_max时，潜水蒸发强度E=0；n为与土壤质地和植被情况有关的经验常数，一般为1～3。轻质地土壤或植物根系吸水深度大、蒸腾旺盛时，则n值较小；土质黏重或有黏土夹层时，n值较大。

极限蒸发深度在理论上是不存在的，但当地下水位埋深较大（如大于3～4m）时，潜水蒸发已十分微弱，可近似为零。为了避免使用极限埋深H_max，下列指数型公式也常会用到。其公式形式为：

式中：参数α需由实测资料定出。如根据山西省汾河灌溉试验站资料分析得出类似模型E/E₀=ae^-bH的参数a和b，见表1-2。

四、蒸发条件下土壤水分非稳定运动问题的解析法

稳定蒸发只在一些特定的情况下才有可能。通常发生的是因蒸发使土壤不断变干的情况，即蒸发时土壤水分运动是非稳定的。对土壤水分非稳定运动问题的求解包括解析法和数值法。采用解析法求解过程中，通常假定大气蒸发能力不变，并以水面蒸发量E₀表示；不考虑地下水位的情况，或者说地下水位埋藏足够深，对土壤水分运动没有影响。土壤蒸发分为三个阶段，不同阶段土壤水分运动定解问题的数学模型基本一致，只是边界条件略有差异，其求解方法也相应的有所不同。现以表土瞬时变干情况为例，简要介绍土壤水分运动问题求解的解析法。

大气蒸发能力无限大时，可以认为地表迅速变干，含水率由初始值θ_i瞬时地降低为某一风干含水率θ_c。此时可近似地认为蒸发不经过第一和第二阶段，而直接按第三阶段处理，地表为第一类边界。对此有多种解法，下面仅介绍Gardner的解法。研究范围属半无限均质土柱、初始含水率均匀分布的理想状况。

当忽略重力作用时，表土瞬时变干的土壤水分运动定解问题为：

对此问题，Gardner（1959年）提出了一种半解析求解方法。使用Boltzmann变换，即令：

式中：D_c为风干含水率θ_c相应的扩散率值，为一常量。此时，含水率θ可表示为ξ的函数

θ=f（ξ）=f［ξ（z，t）］

由式（1-58）可知，。因此，定解问题式（1-57）经变换后为：

为方便起见，设

将式（1-60）用于定解问题式（1-59），则转化为：

由方程式（1-61a）可得：

上式是关于的一阶齐次微分方程，其通解为

其中：A为任意常数。对上式积分，并应用式（1-73c）所示的条件，得

由式（1-61b）的条件，可得

故，最后结果为：

式（1-62）只在理论上得出了θ^*～ξ的关系，但由于被积函数中包含D，而D又是θ亦即θ^*的函数，因而还不能直接积分。为此，拟选用迭代解法。首先取D为θ的指数函数

迭代计算步骤如下：

（1）取θ^*=0，即θ=θ_c作为迭代初值。相应有D=D_c，将其代入式（1-62），可求出第1次迭代结果

已知误差函数erfc（u）=du，且erfc（∞）=1，故：

此即第1次迭代结果。由误差函数表或误差函数的级数表达式，可查得或算出θ^*～ξ关系。

（2）已知θ^*～ξ的第1次迭代值，由式（1-62）可得出相应的D～ξ关系，代入式（1-62），利用数值积分，则可求得θ^*～ξ的第2次迭代值。

（3）重复上一步骤，直到前后二次θ^*～ξ迭代值之差小于所规定的误差要求为止。

求得θ^*～ξ的关系后，由式（1-60）可转化为θ～ξ的关系。利用式（1-58），对任一给定的时间t，可求得土壤剖面的含水率分布θ～z。

表土的蒸发强度E是随表土的干燥逐渐减小的。表土蒸发强度E可表示为：

在z=0的地表处，θ=θ_c，D=D_c，ξ=0，所以

由于θ^*～ξ的关系已经求得，为已知，由式（1-65）可求得随时间变化的蒸发强度。

为了近似地分析蒸发强度E与时间t的关系，可将扩散率D（θ）简化为常数。此时，由式（1-64）可知：

因表土蒸发强度，从而得出其计算式为：

相应地累积蒸发量W_E为

式（166）和式（167）表明蒸发强度E和成反比，累积蒸发量W_E和成正比。对于，根据Crank（1956年）的研究，在蒸发脱水过程中，扩散率D（θ）的加权平均值可表示为

五、蒸发条件下土壤水分非稳定运动问题的数值法求解

由上述求解过程可见，解析法均需要做一些简化，如忽略重力项、将D（θ）和简化为常数等，由此会导致求解结果与实际情况不一致，有时其偏离程度达到无法接受的程度，因此实践中更多地倾向于采用数值法求解土壤水分运动问题。现以田间均质一维非饱和土壤水分运动为例，对数值法做简要介绍。其定解问题可表示如下：

式中：θ_a为均匀分布的初始含水率；L为下边界深度；E为地表蒸发强度；其余符号意义同前。

为了对问题求解，常将全部区域z=0～L离散化为n个单元，共有n结点，编号为i=0，1，2，…，n，其中i=0和i=n为边界结点（图1-11），距离长为ΔZ，时间步长为Δt，时间步长常采用变步长，开始时取小值，以后逐步加大。

对于任一内结点，按隐式差分格式写出式（1-81）的差分方程：

令，，上式经整理后可写为：

式中

当i=1时，即上边界条件，由此可写出方程：

其中

式（1-72）为方程组式（1-71）中的第一个方程。对于本定解问题式（1-69），上式中的。

当i=n-1时，即下边界条件，由此可写出方程：

式（1-70）、式（1-72）和式（1-73）形成了三对角方程组式（1-75），可按照追赶法求解。

由上可知，方程组式（1-75）中的系数包含了待求变量，必须对方程组进行线性化，才能对方程组求解。线性化的方法有多种，如显式线性化、预报校正法、迭代法，包括参数取值方法等，可参见《土壤水动力学》（雷志栋等，1988年）。

对本定解问题，雷志栋等（1988年）以某一轻壤土的灌溉和蒸发过程为例进行了模拟，其初始土壤含水量为0.15cm³/cm³，采用喷灌，喷灌强度为20mm/h，喷灌90min后，以5mm/d的强度蒸发，当表土达风干含水率（0.02）后，表土含水率维持不变。其土壤剖面含水率动态模拟结果见图1-12。