【课程】04 土壤水动力学

《土壤水动力学》1987.3 第三章1~2 节及第四章

第三章 土壤水分的入渗

3.1 土壤水分入渗过程概述

入渗是指水分进入土壤的过程，这是自然界水循环中的一个重要环节。

3.1.1 入渗过程的一般分析

初始含水率${\theta }_t$较低的土壤，在积水条件下入渗一定时间后，土壤剖面中含水率分布可分为四个区：饱和区、含水率有明显降落的过渡区、含水率变化不大的传导区和含水率迅速减少至初始值的湿润区。

干土积水后，土壤含水率分布随时间的变化：

（1）在水施加于土壤表面后的很短时间内，表土的含水率$\theta (0,t)$将很快地由初始值${\theta }_i$增大到某一最大值${\theta }_0$、由于完全饱和在自然条件下一般是不可能的，故${\theta }_0$值较饱和含水率${\theta }_s$为小。

（2）随着入渗的进行，湿润锋不断前移，含水率的分布曲线由比较陡直逐渐变为相对缓平。

（3）在地表$z=0$处，含水率梯度的绝对值逐渐由大变小，当$t$足够大时，即接近地表处含水率不变。

累积入渗量：
$$I(t)={\int }_0^L[\theta (z,t)-\theta (z,0)]dz$$
入渗率为单位时间内通过地表单位面积入渗到土壤中的水量：
$$i(t)=q(0,t)=[-D(\theta )\frac{\partial \theta }{\partial z}+K(\theta ){]}_{z=0}$$

3.1.2 地表入渗条件

为了求出入渗过程中土壤含水率分布，以及入渗率等随时间变化的定量结果，可以在一定的初始含水率分布条件下，根据入渗边界条件，求解水流运动的基本方程。对于入渗边界，可有下列三种简化模型：

（1）地表含水率已知

灌溉模型，第一类边界条件：
$$\theta ={\theta }_{0}t>0z=0$$
（2）地表通量已知

降水模型，第二类边界条件：
$$-D(\theta )\frac{\partial \theta }{\partial z}+K(\theta )=R(t)t>0z=0$$
（3）地表积水

积水模型，土壤中水分运动为饱和-非饱和流动。

3.2 土壤水分运动线性化方程的入渗解

用解析方法求解非饱和土壤水分运动方程，主要困难之一在于该微分方程的非线性，即土壤水分运动参数——扩散率、导水率及比水容量本身是含水率或基质势的函数。最简单的线性化方法是用某种意义上的平均值来代替$D(\theta )$和$K(\theta )$，使土壤水分运动方程近似为线性偏微分方程。这种简化显然与实际相差较大，仅适用于含水率$\theta$变化不大的情况，但只有这种简化才能得出解析解。

3.2.1 水平入渗

对水平半无限长均质土柱来说，第一类边界条件，且扩散率为常数$\overline{D}$的定解问题
$$\{\begin{array}{lr}\frac{\partial \theta }{\partial t}=\overline{D}\frac{{\partial }^2\theta }{\partial x^2}& \\ \theta (0,x)={\theta }_i& \\ \theta (t,0)={\theta }_0& \\ \theta (t,\infty )={\theta }_i& \end{array}$$
上述定解问题可用拉普拉斯变换求解，即水平入渗时线性化方程的入渗解为
$$\theta (x,t)=({\theta }_0-{\theta }_i)erfc(\frac{x}{2\sqrt{\overline{D}t}})+{\theta }_i$$

3.2.2 垂直入渗

定解问题
$$\{\begin{array}{lr}\frac{\partial \theta }{\partial t}=\overline{D}\frac{{\partial }^2\theta }{\partial z^2}-\overline{K}\frac{\partial \theta }{\partial z}& \\ \theta (0,z)={\theta }_i& \\ \theta (t,0)={\theta }_0& \\ \theta (t,\infty )={\theta }_i& \end{array}$$
垂直入渗条件下。线性方程入渗解的解析表达式
$$\theta (z,t)=\frac{{\theta }_0-{\theta }_i}{2}[erfc(\frac{z-\overline{K}t}{2\sqrt{\overline{D}t}})+e^{K_z/D}erfc(\frac{z+\overline{K}t}{2\sqrt{\overline{D}t}})]$$

第四章 蒸发条件下土壤水分的运动

4.1 概述

陆地的水分蒸发可发生在土壤表面，也可发生在植物的叶面。前者称为土面蒸发，后者称为植物蒸腾。陆地上的水分蒸发总量约占降水总量的70%左右，因此它是陆地水分循环中的主要纽成部分。
在潜水位较高的地区，随着水分的腾发盐分将向表土积累，若降水量较少，则土壤势必受到盐渍化的威胁。蒸发作用的强弱常以蒸发强度表示，即单位时间内单位面积地面上所蒸发的水量$(mm/d)$。

土面蒸发的形成及蒸发强度的大小主要取决于两方面的因素。一是受辐射、气温、湿度和风速等气象因素的影响。即大气蒸发能力。二是受土壤中含水率的大小和分布的影响，即土壤供水能力。根据大气蒸发能力和土壤供水能力所起的作用、土面蒸发所呈现的特点及规律，将士面蒸发过程区分为两个或三个阶段。表土蒸发强度保持稳定的阶段、表土蒸发强度随含水率变化的阶段、水汽扩散阶段。

蒸发第一阶段地表处蒸发边界条件可表示为
$$E=E_0\theta \le {\theta }_k$$
蒸发第二阶段，线性关系
$$\frac{E}{E_0}=a\theta +b{\theta }_e<\theta <{\theta }_k$$

4.2 定水位条件下均质土壤的稳定蒸发

4.2.1 稳定蒸发条件下土壤含水率及吸力分布

定解问题
$$\{\begin{array}{lr}-D(\theta )\frac{d\theta }{dz}-K(\theta )=E& \\ \theta (t,0)={\theta }_s& \end{array}$$
积分可得
$$z={\int }_{\theta }^{{\theta }_s}\frac{D(\theta )}{K(\theta )+E}d\theta$$

$$z=\int \frac{1}{a{s}^m+\beta }ds$$

当土壤导水率的参数$a_1$、$a_2$和$m$已知时，由以上各式可求出在给定稳定蒸发强度$E$时，潜水位以上的土壤水吸力$s$ 的分布。

4.2.2 定水位下潜水极限蒸发强度

地表处土壤水吸力$S_H$和稳定蒸发强度$E$的关系
$$H=\frac{1}{\sqrt{\frac{E}{a_1}(\frac{E}{a_1}{a}_2+1)}}arctan\sqrt{\frac{E{S}_H^2}{a_1+a_{2}E}}$$
使用下列两式之一进行积分计算，可求出潜水极限蒸发强度$E_{max}$和地下水位埋深$H$的关系
$$H={\int }_{\theta }^{{\theta }_s}\frac{D(\theta )}{K(\theta )+E_{max}}d\theta$$

4.2.3 潜水稳定蒸发及其经验公式

潜水稳定蒸发强度$E$，受土壤供水能力和大气蒸发能力两者的制约，两者之间正是通过地表含水率${\theta }_H$；联系在一起的。对于一种确定的土壤，潜水稳定蒸发强度$E$完全可由水面蒸发强度$E_0$和地下水位埋深$H$来决定。

稳定蒸发条件下的潜水蒸发强度$E$、水面蒸发强度$E_0$和地下水位埋深$H$之间存在一定关系，经验公式应满足以下条件：
$$\{\begin{array}{lr}E=0E_0=0& \\ E<E_{0}andE<E_{max}{E}_0>0& \\ E=E_{max}and\frac{dE}{d{E}_0}=0E_0\to \infty & \end{array}$$
为了避免使用极限埋深$H_{max}$，下列指数型公式在实际应用时也有较好的结果，公式形式如下
$$E=E_0{e}^{-aH}$$
式中，参数$a$需由实测资料定出。

4.3 定水位条件下层状土壤的稳定蒸发

地下水位一定，非均质层状土壤的稳定蒸发问题的分析方法基本上同均质土壤稳定蒸发问题。并有以下假定：

（1）大气蒸发能力足够大，蒸发主要受土壤输水能力控制。潜水蒸发强度、各层的主壤水分运动通量和地表蒸发强度相等。

（2）层间界面处土壤水吸力$s$是连续的。

（3）每层土壤本身是均质的。

（4）各层土壤的导水率$K(s)$均以同一经验公式表示，但参数$a_1$、$a_2$和$m$各层可不相同。

凭借关系曲线，可推求层状土壤稳定蒸发强度或土壤吸力分布。

已知稳定蒸发强度，求地下水位以上的土壤水吸力分布；已知地表处吸力（或含水率），求稳定蒸发强度。采用试算法，即先假定一个稳定蒸发强度$E$，按前述方法由下而上求得表土的吸力值$S_H$，和已知值比较，两者不相等时,则假定另一$E$重算，直到表土吸力的计算结果和已知值相等为止。对于非均质层状分布的土壤，可由上述方法求出稳定蒸发强度$E$和地表吸力$S_H$的关系。进而可推算出表土吸力$S_H\to \infty$时的潜水极限蒸发强度$E_{max}$。对不同的地下水位埋深$H$进行计算，可得出$E_{max}H$关系。

4.4 蒸发条件下土壤水分的非稳定运行

稳定蒸发只在地下水位埋深较浅、有侧向补给维持水位不变的情况下才有可能。但通常发生的是因蒸发使土壤不断变干的情况，即蒸发时土壤水分运动是非稳定的。

4.4.1 表土蒸发强度保持稳定阶段的土壤水分运动求解

均质土壤，初始含水率${\theta }_i$分布均匀，蒸发时，吸力梯度在数值上远大于1，定解问题为
$$\{\begin{array}{lr}\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial z}[D\frac{\partial \theta }{\partial z}]& \\ \theta (0,z)={\theta }_i& \\ D\frac{\partial \theta }{\partial z}=Ez=0t>0& \\ \frac{\partial \theta }{\partial z}=0z=It\ge 0& \end{array}$$
1.Covey的求解分析方法

引入无因次变量，将定解问题无因次化，相应定解问题转化为
$$\{\begin{array}{lr}\frac{\partial D^{\ast }}{\partial \tilde{t}}=D^{\ast }\frac{{\partial }^2{D}^{\ast }}{\partial {\tilde{z}}^2}& \\ D^{\ast }=1\tilde{t}=00\le \tilde{z}<\infty & \\ \frac{\partial D^{\ast }}{\partial \tilde{z}}=1\tilde{z}=0\tilde{t}>0& \\ \frac{\partial D^{\ast }}{\partial \tilde{z}}=0\tilde{z}\to \infty \tilde{t}\ge 0& \end{array}$$
解析解为
$$D^{\ast }=1-2\sqrt{D^{\ast }\tilde{t}}ierfc(\frac{\tilde{z}}{2\sqrt{D^{\ast }\tilde{t}}})$$
2.Gardner与Hillel关于土柱均匀蒸发的近似解
θ c ( z ) = θ k + 1 β l n { 1 + β E L 2 D k [ 1 − ( 1 − z / L ) 2 ] } \theta_c(z)=\theta_k + \frac{1}{\beta}ln\{1+\frac{\beta EL}{2D_k}[1-(1-z/L)^2]\} θc​(z)=θk​+β1​ln{1+2Dk​βEL​[1−(1−z/L)2]}

4.4.2 表土蒸发强度随地表含水率递减阶段土壤水分运动的近似解

在蒸发的第二阶段，地表蒸发强度$E$受上壤供水限制，将随地表含水率的降低而减小，并可近似为直线关系。对于半无限长均质土柱，初始含水率${\theta }_i$假定为均勾分布，其定解问题可表述为
$$\{\begin{array}{lr}\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial z}[D\frac{\partial \theta }{\partial z}]-\frac{\partial K(\theta )}{\partial z}& \\ \theta (0,z)={\theta }_i& \\ D(\theta )\frac{\partial \theta }{\partial z}=E=a\theta +bz=0t>0& \\ \theta ={\theta }_{i}z=\infty t>0& \end{array}$$
解得
$$\theta ={\theta }_i-({\theta }_i+\frac{b}{a})\eta (z,t)$$
表土蒸发强度$E$的计算式可写为
$$E=a\theta (0,t)+b=E_i[1-\eta (0,t)]$$

4.4.3 表土瞬时变干时土壤水分运动的求解

研究范围均属半无限均质土柱、初始含水率均匀分布的理想状况。

（1）表土瞬时变干的Gardner解
$$\{\begin{array}{lr}\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial z}[D(\theta )\frac{\partial \theta }{\partial z}]& \\ \theta ={\theta }_{i}t=0z\ge 0& \\ \theta ={\theta }_{e}z=0t>0& \\ \theta ={\theta }_{i}z=\infty t>0& \end{array}$$
经Boltzmann变换后，最后结果为
$${\theta }^{\ast }=\frac{{\int }_0^{\zeta }\frac{1}{D}{e}^{-{\int }_0^{\zeta }\frac{2\zeta D_c}{D}}d\zeta }{{\int }_0^{\infty }\frac{1}{D}{e}^{-{\int }_0^{\zeta }\frac{2\zeta D_e}{D}}d\zeta }$$
选用迭代法求解。

表土蒸发强度计算式为
$$E=({\theta }_i-{\theta }_s)\sqrt{\frac{\overline{D}}{\pi t}}$$
（2）表土瞬时变干的线性化解

表土瞬时变干时，土壤水分运动的线性化解即是将土壤水分运动参数$D(\theta )$及$K(\theta )$近似为常数。以下分别讨论忽略重力项和考虑重力项两种情况的线性化解。

1）忽略重力项的线性化解

定解问题可表述为
$$\{\begin{array}{lr}\frac{\partial u}{\partial t}=D(u)\frac{{\partial }^{2}u}{\partial z^2}& \\ u=1t=0z\ge 0& \\ u=0z=0t>0& \\ u=1z=\infty t>0& \end{array}$$
解得
$${\theta }^{\ast }=\frac{\theta -{\theta }_c}{{\theta }_i-{\theta }_c}=\frac{1}{{\beta }^{\prime }}ln[(e^{{\beta }^{\prime }}-1)erf(\frac{1}{2\sqrt{\overline{D}t}})+1]$$
2）考虑重力作用的线性化解

将一维垂直方向土壤水分运动的基本方程改换成以$v$为因变量的基本方程及相应的定解条件为
$$\{\begin{array}{lr}\frac{1}{D(v)}\frac{\partial v}{\partial t}=\frac{{\partial }^{2}v}{\partial z^2}-M(v)\frac{\partial v}{\partial z}& \\ v=1t=0z\ge 0& \\ v=0z=0t>0& \\ v=1z=\infty t>0& \end{array}$$
应用拉氏变换解得
$$v(z,t)=\frac{1}{2}[erfc(-\frac{z}{2\sqrt{\overline{D}t}}-\frac{m}{2}\sqrt{\overline{D}t})-e^{-mz}erfc(\frac{z}{2\sqrt{\overline{D}t}}-\frac{m}{2}\sqrt{\overline{D}t})]$$
已知$\overline{D}$和$m$后，$v(z,t)$的关系由上式解出。