第二章 土壤水分运动基本方程2

第二章 土壤水分运动基本方程
如前所述，达西定律是由达西（Darcy ，Henry 1856）通过饱和砂柱渗透试验得出，后由Richards （1931）将其扩伸至非饱和水流中，并规定导水率为土壤负压h 的函数，即
()H h k q ∇= （2-2-1）
式中：H ∇——为水势梯度；
k （h ）——为导水率，是土壤负压h 的函数； q ——为水流通量或流速。

Richards 方程垂向一维方程为
)
1)(( )
(±∂∂-=∂∂-=z
h
k z
H k q z θθ
注意：H=h ±z ，垂直坐标向上为“+”；向下时为“–”。

由于k （h ）受滞后影响较大，上式仅适用于单纯的吸湿或脱湿过程。

若将导水率作为容积含水率函数，即以k （θ）代替人k （h ），则可避免滞后作用的影响。

一般说来达西定律对饱和与非饱和水流均可适用，即水流通量与势能梯度成正比。

但在饱和土壤中，压力为正值，其总水头包括了由该点在地下水面以下深度来确定的静水压力（正值）和相对于基准面高度来确定的位置水头，总水头为压力水头和位置水头之和，水由总水头高处向低处流动。

在非饱和土壤中，基质势为负值，土水势在不考虑溶质势、温度势及气压势时，只包括重力势和基质势。

因此，总水头常以负压水头和位置水头之和来表示。

一维Richards 方程的几种形式：
根据()
()θθ
θD h
k =∂∂（K=C ×D ）得： x h k q x ∂∂-=)(θ x D q x ∂∂-=θ
θ)
( y h k q y ∂∂-=)
(θ y
D q y ∂∂-=θθ)( )1)(
(±∂∂-=z h k q z θ )]()([θθθk z
D q z ±∂∂-=
第一节 直角坐标系中土壤水分运动基本方程
一、基本方程的推导
土壤水分运动一般遵循达西定律，且符合质量守恒的连续性原理。

土壤水分运动基本方程可通过达西定律和连续方程进行推导。

如图2-2-1所示，从土壤中取出微分单元体abcdefgh ，其体积为z y x ∆∆∆，由于该立方体很小，在各个面上的每一点流速可以看成是相等的，设其流速为z y x v v v 、、，在t ～t+Δt 时段内，流入立方体的质量为（3个面流入）：
t y x v t z x v t z y v m z y x ∆∆∆+∆∆∆+∆∆∆=ρρρ入 （2-2-2）
流出立方体的质量为（3个面流出）：
t z y x x v v m x x ∆∆∆⎪⎭
⎫
⎝⎛∆∂∂+=ρ出
t y x z z v v t z x y y v v z z
y y ∆∆∆⎪⎭⎫ ⎝⎛∆∂∂++∆∆∆⎪⎪⎭
⎫ ⎝⎛∆∂∂++ρρ （2-2-3） 式中：ρ––––水的密度；
z y x ∆∆∆,,––––分别表示微分体x 、y 、z 方向长度；
x x v x ∆∂∂，y y v y ∆∂∂，z z
v
z ∆∂∂––––分别表示水流经微分体后，其流速在x 、y 、z 方向的变化值。

由式(2一2－2）、式（2－2－3）之差可求得流入和流出立方体的质量差：
出入m m m -=∆ ⎪⎪⎭⎫
⎝⎛∂∂+∂∂+∂∂-=z v y v x
v z y x ρt z y x ∆∆∆∆⨯ （2—2—4） 设θ为立方体内土壤含水率，则在Δt 时间内立方体内质量变化又可写为
t z y x t
m ∆∆∆∆∂∂=∆θ
ρ
（2—2—5） 根据质量平衡原理（流入量－流出量＝储存量变化量），式（3－2－4）、式（3—2—5）
应相等，即
⎪⎪⎭⎫ ⎝⎛∂∂+∂∂+∂∂-=∂∂z v y v x
v t z y x θ
（2-2-6） 根据达西定律得：
()
x H k v x ∂∂-=θ，()y H k v y ∂∂-=θ，()z
H
k v z ∂∂-=θ （2-2-7）
式中k （θ）––––土壤水力传导度，为含水率的函数；
H ––––总土水势，为基质势与重力势之和（H ＝h ＋z ）。

因此，式（2－2—6）可以写作以下形式：
()()()z
z H k y y H k x x H k t ∂⎥⎦⎤⎢⎣⎡∂∂∂+∂⎥⎦⎤⎢⎣⎡∂∂∂+∂⎥⎦⎤⎢⎣⎡∂∂∂=∂∂θθθθ （2-2-8）
上式可以简写为
()[]H k t
∇∇=∂∂θθ
（2-2-9） 式（2－2－8）或式（2－2－9）为土壤水分运动基本方程。

在饱和土壤中，含水量和基质势均为常量。

水力传导度也为常量，常称渗透系数，则方程（2－2－8）可写为
022222
2=∂∂+∂∂+∂∂z
H
y H x H （2－2－10） 或写作
02=∇H （2－2－10‘）
22
22222
z
y x ∂∂+∂∂+∂∂=∇ （2－2－11）
式中：▽2––––
拉普拉斯算子。

式（2－2－10）或式（2－2－10‘
）为饱和土壤水流的拉普拉斯方程。

二、基本方程的不同形式
为运用基本方程分析各种实际问题的方便，可将基本方程改写为多种表达形式。

为简便起见，以下均以一维垂向土壤水分运动为例，给出基本方程的不同表达形式。

（一）以含水率θ为变量的基本方程
由式（2－2－8）可得一维垂向土壤水分运动的基本方程为
()⎥⎦
⎤⎢⎣⎡∂∂∂∂=∂∂z H k z t θθ （2－2－12） 式中：H ––––总土水势；
z ––––为水流方向坐标，取z 向上为正。

因为H=h 十z ，所以上式可写作
()()z
k z h k z t ∂∂+⎥⎦⎤⎢⎣⎡∂∂∂∂=∂∂θθθ （2－2－13） 式(2－2－13）为以θ为变量的基本方程，将
z
h z h ∂∂∂∂=
∂∂θ
θ代入式（2－2－13）得： ()()z k z h k z t ∂∂+⎥⎦
⎤⎢⎣⎡∂∂∂∂∂∂=∂∂θθθθθ 令()
()θθ
θD h
k =∂∂，则式（2—2—13）可以写成（一维垂向土壤水分运动方程）： ()()z
k z D z t ∂∂+⎥⎦⎤⎢⎣⎡∂∂∂∂=∂∂θθθθ （2－2－14） 在水平运动的情况下，重力项等于0，所以()
x
D v x ∂∂-=θ
θ，其形式与Fick 扩散定律相同。

式（2－2－14）具有扩散方程的形式，故将D （θ）称为扩散度。

()⎥⎦
⎤⎢⎣⎡∂∂∂∂=∂∂x D x t θθθ （2－2－14‘
） Fick 定律：
自由水中溶质的分子扩散通量符合Fick 定律：
x
c D
J ∂∂-= 式中：J 为溶质的扩散通量； D 为溶质的扩散系数；
x
c
∂∂为溶质的浓度梯度。

（二）以基质势h 为变量的基本方程 由于
()t
h h c t h h t ∂∂=∂∂∂∂=∂∂θθ ，则式（2－2－14）可以写成： ()
()()z
h k z h h k z t h h c ∂∂+⎥⎦⎤⎢⎣⎡∂∂∂∂=∂∂ （2－2－15）
式中：c （h ）––––比水容量（也称容水度），c （h ）=h
∂∂θ
，表示单位基质势变化时含水率变化。

（三）以参数v 为因变量的基本方程
采用Kirchhoff 变换，令()()()⎰
⎰⎰
-∞
=
=
c
c
c
h h
h h
h d k V
d k d k v τττττ
τ1
则
()h k V
h v 1
=∂∂ ()⎰∞
=
c
h d k V ττ
由式（2-2-15）得：()()z
h k z z h h k t h h ∂∂+∂⎥
⎦⎤⎢⎣⎡
∂∂∂=∂∂∂∂θ
()
()z
v v h h h k z z v v h h k t v v h h ∂∂∂∂∂∂+∂⎥⎦⎤⎢⎣⎡
∂∂∂∂∂=∂∂∂∂∂∂θ
()()()()()z
v h k V h h k z z v h k V h k t v h k V h ∂∂∂∂-∂⎥⎦⎤
⎢⎣
⎡∂∂∂=∂∂∂∂θ
()
()z
v
v X z v t v v Y ∂∂+∂∂=∂∂22 （2－2－16） 式中h c ––––土壤的进气值，即土壤含水率开始小于饱和含水率时的负压值。

另外，()()()()()h k h c h D h h k v Y ==∂∂=
11θ；()()()h
h k h k v X ∂∂=1
在非饱和区：
()01
<=
⎰h
h c
d k V
v ττ
在饱和区：
()01
>=⎰h
h c
d k V
v ττ
且因为 ()0=∂∂=
h h c θ，()0=∂∂h
h k 所以 ()0=v Y ；()0=v X
则方程式（2－2－16）为：022=∂∂z
v
（四）以位置坐标z 为变量的土壤水运动方程
以z 为变量，则z 为θ、t 的函数，z （θ，t ）为未知函数。

已知θ=θ（z ，t ），当0≠∂∂z
θ
处，可以解出z= z （θ，t ），即[14]
()()0,,≡-t t z z z θ
对z ，t 分别求导数：
01=∂∂∂∂-
z z θθ，0=∂∂-∂∂∂∂-t
z
t z θθ
于是 θθ∂∂=∂∂z z 1及θ
θ∂∂∂∂-=∂∂z t z
t 将以上式子代入方程（2－2－14）得：
()()z k z D z t ∂∂+⎥⎦
⎤⎢⎣⎡∂∂∂∂=∂∂θθθθ ()z k z z D z t z
∂∂∂∂+∂∂⎪⎭⎫ ⎝⎛∂∂∂∂=∂∂∂∂-θ
θθθθθθ
()z
k z z D z t z ∂∂∂∂+∂∂⎪⎭⎫ ⎝⎛∂∂∂∂=∂∂∂∂-
θ
θθθθθθ ()θ
θθθθ∂∂+∂∂⎪⎭⎫ ⎝⎛∂∂∂∂=∂∂-
k
z z D t z （2－2－17）
（五）以参数u 为因变量的土壤水运动方程 定义()()()⎰⎰⎰=
=
θ
θθθ
θ
θθθθθθθi
s
i
i
d D U
d D d D u 1
式中：i θ––––初始含水率；
()⎰=
θ
θθθi
d D U ；
s θ—饱和含水率。

由式（2-2-14）得：
()()z
k z k z t ∂∂+⎥⎦⎤⎢⎣⎡∂∂∂∂=∂∂θθθθ
()t u
u k z u u D z t u u ∂∂∂∂∂∂+⎥⎦
⎤⎢⎣⎡∂∂∂∂∂∂=∂∂∂∂θθθθθ
将
()θθD U
u 1
=∂∂代入上式得：
()()()()z u
D U k z u D U D z t u D U ∂∂∂∂+
⎥⎦
⎤⎢⎣⎡∂∂∂∂=∂∂θθθθθ 所以 ()z u
k z
u D t u ∂∂∂∂+∂∂=∂∂θθ22 （2-2-18） 以上各式中式（2－2－14）、式（2－2－15）是二种经常采用的形式，形式的选定取决
于要解决问题的边界条件和初始条件。

以含水率θ为因变量的基本方程常用于求解均质土层或全剖面为非饱和流动问题，这种方程形式对于层状土壤或求解饱和—非饱和流问题不适用；以负压水头h 为因变量的基本方程是应用较多的一种形式，可适用于饱和—非饱和水流求解及层状土壤的水分运动分析计算，但由于非饱和土壤水的导水率k （h ）及容水度c(h)，受滞后影响较大，计算中参数选取不当会造成较大误差；以v ，u 为因变量基本方程实际上分别相当于以负压水头h 和含水率θ为因变量的基本方程，在某些情况下由于经代换后方程较为简单，易于求解；以坐标为因变量的基本方程根据定解条件需要求解较简单的土壤水分运动问题。

以上为直角坐标系中土壤水分运动的基本方程，求解某些土壤水分运动问题时，采用柱坐标系可能更方便。

第二节 柱坐标系中土壤水分运动基本方程
在推导柱坐标系中的基本方程时，方法同直角坐标系，同样可用达西定律与连续方程相结合的方法导出。

若以z 轴为轴的柱坐标系，根据达西定律，在此坐标系中可表示为：
()()⎪⎪⎪
⎩
⎪
⎪
⎪⎨⎧∂∂-=∂∂-=∂∂-=z H k q H k r q r H k q z r
θϕθθϕ1)( 式中：r 、φ、z ––––分别表示柱半径，角坐标和垂直坐标（图2—2—2）
q r 、q φ、q z ––––相应于r 、φ、z 三个方向的通量； H ——总水势。

下面利用质量守恒来推导连续方程。

Δt 时段内，在r 方向的流入量为r q r △φ△z △t ，流出量为
⎪⎭
⎫
⎝⎛∆∂∂+r r q q r r ()r r ∆+t z ∆∆∆ϕ，则流入与流出量之差（忽略高阶无穷小量）为
t z r r r q q r r ∆∆∆∆⎪⎭
⎫
⎝⎛∆∂∂+
-ϕ （2-2-19） 同理，在φ方向流入流出量之差为 t z r q ∆∆∆∆∂∂-
ϕϕ
ϕ （2-2-20）
在z 方向土壤水分流入流出量之差为 t z r q z
∆∆∆∆∂∂-
ϕϕ
（2-2-21） 上述三个方向流入和流出单元体的水量差总计为
t z r z q q r q
r q z
r r ∆∆∆∆⎪⎪⎭⎫
⎝
⎛∂∂+∂∂+∂∂+-ϕϕϕ （2-2-22） 单元体体积应为 z r r r ∆∆∆⎪⎭
⎫
⎝⎛
∆+ϕ2，略去高阶无穷小量后为r △φ△z △t ，在△t 时间内单元体内水分增量为
t z r t
∆∆∆∆∂∂ϕθ
（2-2-23） 根据质量守恒原理武（3－2－22）应与式（3－2－23）相等，即
⎪⎪⎭
⎫ ⎝⎛∂∂+∂∂+∂∂+-=∂∂z q q r r q q r t z r r ϕθ
ϕ11 （2-2-24） 式（2－2－24）为柱坐标系中土壤水分运动的连续方程。

将式（2－2－18）代入上式，
即得柱坐标系中土壤水分运动基本方程：
()()()()⎥⎦
⎤
⎢⎣⎡∂∂∂∂+⎥⎦⎤⎢⎣⎡∂∂∂∂+⎥⎦⎤⎢⎣⎡∂∂∂∂+∂∂=∂∂z H k z H k r r H k r r H r k t θϕθϕθθθ21 （2-2-25） 以总水势H=h+z ，水容度c=dh
d θ
，以及导水率k （θ），扩散度D （θ）等代入，基本方程可表示为
()()()()z
k z D z D r r rD r r t ∂∂+⎥⎦⎤⎢⎣⎡∂∂∂∂+⎥⎦⎤⎢⎣⎡∂∂∂∂+⎥⎦⎤⎢⎣⎡∂∂∂∂=∂∂θθθϕθθϕθθθ211 （2-2-26） 对于平面轴时称问题，上式可改写为
()⎥⎦
⎤⎢⎣⎡∂∂∂∂=∂∂r rD r r t θθθ1 （2-2-26‘
） 同理可推得以x （或y ）轴为轴的柱坐标系的基本方程：
()()()()()⎥⎦⎤
⎢⎣⎡∂∂∂∂+∂∂-∂∂+⎥⎦
⎤⎢⎣⎡∂∂∂∂+⎥⎦⎤⎢⎣⎡∂∂∂∂=∂∂x D x k r r k D r r rD r r t θθϕθϕθϕϕθθϕθθθsin cos 112
（2-2-27） 关于球坐标系中的基本方程应用较少，推导方法同上，这里不再论述。

第三节 土壤水分运动基本方程的定解条件
土壤水运动基本方程的定解条件包括初始条件和边界条件，为了简单起见，将以直角坐标系中基本方程常用形式为例进行论述。

（一）初始条件
相应于前式（２—２—1４）、式（２—２—15）的初始条件分别以下式表达：
()()z z i θθ=0, （2-2-28）
()()z h z h i =0, （2-2-29）
脚标“i ”表示初始已知量。

初始条件：t=0时剖面上θ、h 的分布已知。

（二）边界条件
边界条件一般有一类边界、二类边界、三类边界三种。

1．一类边界条件（变量已知的边界Γ1）
对干式（2－2－14）、式（2－2－15）的一类边界的表达式为
()()t t z 00,θθ= （2-2-30）
()()t h t z h 00,= （2-2-31）
脚标“0”均表示一类边界上的值；z 0为一类边界的坐标。

在一维垂向土壤水分运动中，一类边界的情况发生在压力入渗（地表形成水层）时，地表含水率达到饱和含水率，或当强烈蒸发时，表土达到风干土含水率的情况。

2．二类边界条件（边界Γ2上水流通量已知的情况） 相应于式（2－2－14）、式（2－2－15）表达式为
()
()()t k z D εθθ
θ=+∂∂-Γ2| （2-2-32） ()()()t z z h h k ε=∂-∂-Γ2| （2-2-33） 式（2－2—32）及式（2－2－33）中，均假设垂直坐标z 向下为正。

在一维垂向土壤水分运动中，这种情况常发生在降雨、灌水入渗或蒸发强度已知的边界。

在降雨或灌水入渗时，ε（t ）为正值，在蒸发时ε（t ）为负值。

在不透水边界和无蒸发入渗的边界，ε（t ）=0，则式（2－2－32）、式（2－2－33）分别为
()
()θθ
θk z D =∂∂ （2-2-32，） ()()h k z h
h k =∂∂ （2-2-33，）
3．三类运界条件［相当于水流通量随边界Γ3上的变量（θ或h ）值而变化的情况］
三类边界的一般形式为 321
a f a z
f
a =+∂∂ （2-2-34） 式中，f 为变量。

在土壤蒸发强度为表土含水率或表土负压的函数的情况下，式（2－2－14）、式（2－2－15）的三类边界条件表达式为
()
()b a k z D +=-∂∂θθθ
θ （2-2-35） ()()()b h a h k z
h
h k +=-∂∂ϕ （2-2-36）
方程式（2－2－35）右端的a θ+b 表示三类边界上水流通量为表土含水率θ的线性函数。

方程式（2－2－36）右端的a φ(h)+b 表示三类边界上水流通量为表土负压的函数。

例如、土壤的下部有弱透水层阻隔，边界受顶托补给，其补给强度决定于下部弱透水层的导水率k 2，弱透水层上、下的压力h 1、h 2，其厚度为δ，方程的三类边界条件可写成：
()
()()δ
δ--=-∂∂122h h k h k z h h k （2-2-36，
） 上述的二种边界条件是经常遇到的情况。

在野外实际情况下，有时还存在地下水位为已
知的运动边界，此时可将地下水面处h=0作为边界条件。

如在任一时间，地下水埋深为d （t ），则
()()0, ,==t d h t d z （2-2-37）
d （t ）表示 t 时刻的地下水面所在位置，如地下水位是等速下降的，则
()0d vt t d += （2-2-38）
式中：v ––––地下水位下降的速度。

如地下水位是由于水井抽水引起下降的，则
()⎪⎪⎭
⎫ ⎝⎛+=at r W T Q
d t d 4420π （2-2-39） 式中：d 0––––下水初始埋深；
Q ––––井的抽水流量；
T ––––地下水含水层的导水系数； a ––––地下水含水层的导压系数；
⎪⎪⎭
⎫ ⎝⎛at r W 42––––井函数（见泰斯公式）；
r ––––计算点离抽水井
若地下水位变化规律未知，不能作动边界处理。

第四节土壤水分运动参数
土壤水分运动中的主要参数有土壤水力传导度k（又称导水率），比水容量c（也称容水度）及扩散度D（也称扩散率）等。

这些参数的变化决定了土壤水分运动状态，所以了解和掌握这些参数的特性及其变化规律是十分重要的。

一、土壤水力传导度k
土壤水力传导度是反映土壤水分在压力水头差作用下流动的性能。

一般在饱和土壤中导水率称为渗透系数。

土壤水力传导度为在单位水头差作用下，单位断面面积上流过的水流通量。

它是土壤含水率或土壤负压的函数。

饱和土壤孔隙中都充满水，导水率达到最大值，且为常量。

在非饱和土壤中，因土壤孔隙中部分充气，导水孔隙相应减少，导水率低于饱和土壤水情况，而且导水率是负压或含水率的函数，随着含水率降低而减小。

由于在吸力作用下，土壤水首先从大孔隙中排除，随着吸力增加，水流仅能在小孔隙中流动。

所以，土壤从饱和到非饱和将引起导水率的急剧降低。

当吸力由零增至1×105 Pa时，导水率可能降低好几个数量级，有时降低到饱和导水率的1／100000。

对于不同结构土壤，饱和与非饱和土壤水导水性能的相对关系是不同的。

饱和土壤导水性能最好的是粗粒砂性土壤，导水最差的土壤是细质粘土，但非饱和土壤在较大负压情况下则情况可能相反。

具有大孔隙粗质土壤，在吸力作用下孔隙中水分很快排除，导水率迅速下降；而粘质细颗粒土壤，在较高吸力下，许多小孔隙仍充满水，仍具有一定的导水性能，因此，导水率下降较缓慢。

所以，细颗粒粘质土壤在同一吸力条件下可能较大孔隙粗质土壤具有较高导水率[15]。

所以，导水率与含水率（或负压）关系是较复杂的，目前还不能用理论分析方法推导它们之间关系式，需通过试验测定。

图2－2－3为非饱和土壤水在负压水头作用下流动的模型。

在水平土柱两端有多孔板，分别由平水箱保持一定水位，使其负压为h1和h2，在负压梯度△h/△x的作用下，立柱中土壤水从1端向2端运移。

土壤水通量q可由1端补给量或2端溢出量测得，两者相等时，水流处于稳定状态。

非饱和土壤水力传导度可由达西定律求得。

由于水平土柱沿程负压（或含水率）是变化的，求得的导水率k 也应是变化的，若距离较小，可以平均负压（或含水率）确定平均土壤水力传导度。

在不同的平均负压（吸力）值下，通量与负压梯度成正比，如图2－2－4所示[15]。

两者呈直线关系，但其斜率（即水力传导度）随平均负压而变。

此外，土壤水力传导度还与土质有关，如图2－2－5所示，砂性土壤饱和导水率高于粘性土壤，随着土壤吸力增加，砂性土壤导水率降低速率较粘性土壤快，所以吸力增大时，粘质土壤导水率反大于砂质土壤。

非饱和土壤水力传导度k 与土壤负压h 或含水率θ的关系通常由试验资料拟合成经验公式，一般有以下几种形式。

（1）土壤水力传导度与负压（吸力）h 的关系式：
()||h c s e k h k -= （2-2-40）
()n h a h k -=|| (2-2-41)
()b
h a h k n +=
|| (2-2-42)
式中：k s ––––饱和土壤水导水率，或称渗透系数；
a ，
b ––––经验常数；
c ，n ––––经验指数。

（2）土壤水力传导度与土壤含水率θ的关系式：
()()θθθ--=s c s e k k (2-2-43) ()m s k k θθ= (2-2-44)
()n
r s r
s k k ⎪⎪⎭
⎫ ⎝⎛--=θ
θθθθ (2-2-45)
()2
12
111⎪⎭
⎪
⎬⎫⎪⎩⎪
⎨⎧⎥⎥⎥⎦
⎤⎢⎢⎢⎣
⎡
⎪⎪⎭⎫
⎝⎛----⨯⎪⎪⎭
⎫ ⎝⎛--=m
m
r s r r s r
s k k θθθθθ
θθθθ (2-2-46) 式中：θ r ––––某一特征含水率，通常采用最大分子持水率；
θs ––––饱和含水率；
c ，m ，n ––––均为经验指数，在式（2－2－46）中，10,1
1<<-
=m n
m 其他符号意义同前。

k 值的量纲单位为［LT -1］。

由于土壤负压与含水率的关系曲线––––土壤水分特征曲线有滞后现象，所以，土壤水力传导度随负压的变化同样也存在滞后现象，即在同一负压下，干燥过程中的土壤水力传导度高于湿润过程中的土壤水力传导度。

但土壤水力传导度与含水率关系受滞后作用的影响较小。

二、土壤水分扩散度D
土壤水扩散度为单位含水率梯度下，通过单位面积的土壤水流量，其值为土壤含水率的函数，即()()
θ
θθ∂∂=h
k D 。

扩散度与土壤含水率的关系如图2－2－6所示，这种关系有时可用以下经验公式表示：
()θθb ae D = (2-2-47)
式中：a ，b ––––经验常数，D 的量纲为［L 2T -1］。

以上公式仅能用于含水率较高的阶段。

在土壤含水率很低时，由于土壤水汽扩散速度增大，使扩散度随土壤含水率降低而增大。

在土壤含水率很高的情况下，土壤接近于饱和，扩散率趋向于无限大。

三、容水度c
容水度表示压力水头减小一个单位时，自单位体积土壤中所能释放出来的水体体积，量
纲为［L
－１
］。

容水度可以用下式表示：
()dh
d h c θ=
它是负压的函数，为水分特征曲线上任一特定含水率θ值时的斜率（导数），并随土壤水分特征曲线而变化，所以它取决于土壤含水率和土壤质地等。

第五节 考虑水汽热耦合关系的土壤水分运动基本方程
长期以来，人们都采用等温模型研究土壤水分运动。

在自然条件下，日夜温差很大，地表以下不同深度处温度的差异和变化影响土壤水分的转化和运移，用等温模型来研究土壤水分运动常带来一定误差。

一些学者根据能量平衡和热传导理论，提出了考虑水、汽、热耦合关系的土壤水运动基本方程。

根据Philip 、De Vries 和Milly 的理论［１６~１８］
，多孔介质中的水分质量通量可以表示为
V L m q q q
+= (2-2-48)
式中：m q
––––水分总质量通量（kg/s ）；
L q
––––水流质量通量（kg/s ）； V q
––––水汽质量通量(kg/s)。

在入渗速率不很大的情况下，土壤中水和汽之间存在局部热动力学平衡时，L q 及V q
可表示为
［１７］
()1+∇-=h k R q L L ρ
(2-2-49) h D T D q HV V TV V V ∇+∇-=ρρ
(2-2-50)
式中：ρＬ––––土壤水密度（kg ／cm 3）；
h ––––土壤水压力水头（cmH 2O ，1cmH 2O=9．8×103Pa ）； T ––––绝对温度（C ）； D TV ––––热蒸汽扩散系数； D HV ––––蒸汽传导率；
Ｒ––––水蒸汽气体常数，R ＝4．615×１０６erg ／（g ℃）（1erg=10－７
J ）。

在水汽和多孔介质中水体达到局部平衡时，两者之间的自由能相等，则
()()⎪
⎪⎭
⎫
⎝⎛
+=273exp 0T R hg
T V ρρ （2-2-51） 式中：ρv ––––绝对湿度（g ／cm 3）；
ρ0（T ）––––水汽饱和状态下的绝对湿度（g ／cm 3），根据Camilo 等[19]研究表明：
()()[]273ex p 100+-=T R R T ρ (2-2-52)
Constantz 等对水力传导度与温度的关系进行了研究，提出了水力传导度可表示为
()()()T g h kk T h k L r μρ=, (2-2-53)
若忽略温度变化对流体密度影响，则k （h ，T ）可表示为
()()
()
()
T T T h k T h k μμ00,,= (2-2-54) 以上各式中：k ––––多孔介质内渗透率（cm 2）；
k r (h ）––––相对非饱和渗透率； g ––––重力加速度（cm ／s 2）；
k(h, T 0）––––在参考温度T 。

时的水力传导率（cm/d ） μ––––流体动力粘滞系数。

式（2－2－54）中μ（T 0）／μ（T ）可表示为
()()2
020000211.00384.01000211.00384.01T T T T T T ++++=
μμ (2-2-55)
根据Philip 和De Vries 等研究，热蒸汽扩散系数D TV 和蒸汽传导率D HV 分别可表示为
［１
６，１７］
()()
C s cm
T R h g T H fD D v a L TV ⋅⎪
⎪⎭
⎫ ⎝⎛+-∂∂Ω=-/2732
001
ρρθρ (2-2-56)
()
()s cm T R g
fD D v
a L HV /2731
+Ω=-θρ (2-2-57)
)
⎩⎨
⎧>++<='k
L k V k
L V
f θθθθθθθθθ1 (2-2-58)
()
()s cm T D a /2731229.02
75
.1-= (2-2-59)
()3
2L θθ-=Ω (2-2-60)
式中：f ––––在Philip 和 De Vries 公式中对水汽扩散引入的修正因子，f=f ′ξ，经验系数ξ=1．3～3．2；
θk ––––液体水流流动盯以忽略时的含水率； Da ––––空气中分于扩散系数（cm 2/s ）； Ω––––由气体所充填孔隙的弯曲率；
θL ––––孔隙中蒸汽的体积与土壤体积比； H ––––为相对湿度，H=exp （hg ／RT ）。

根据质量守恒定律得；
()()V L V V L L q q t
+-∇=∂+∂θρθρ (2-2-61)
将式（2－2－49）、式（2－2－50）、式（2－2－51）代人上式，得到水分运动方程：
()[][]z
k T D h D k t T S t h S TV HV T h
∂∂+∇∇+∇+∇=∂∂+∂∂ (2-2-62)
其中 T
V L V r L L
V
h h
h S ∂∂+
∂∂⎪⎪⎭
⎫ ⎝
⎛-
=ρρθθρρ1
H
T
T S V L V r L
L V T ∂∂+
∂∂⎪
⎪⎭⎫ ⎝
⎛-=ρρθθρρ1
在单位体积多孔介质中的热量为
()⎰-+-=L
Wd L T T C Q L V L L θθρθρ0
00 (2-2-63)
式中：C Ｌ=
∑=5
1
i i
i C θ
为水汽介质的热容量（cal ／cm 3，lcal ＝4．1868J ）；
C i ，θi ––––土壤中水、汽、石英和其他矿物质及有机质的体积热容量及占土体百分数；
L 0––––在参考温度T 下的蒸发潜热，cal ／g ，在T=20℃时，L0=585cal ／g 。

蒸发潜热L 可以下式表示：
()()00T T C C L L l v --+=
式中：C v ––––常压下水蒸汽比热；
C i ––––水的比热。

W 为微分吸湿热（cal ／g ），由热动力学原理，可表示为[17—１９]
：
()g cal T j T h g j W /1⎪
⎭⎫ ⎝
⎛
∂∂--=- (2-2-64)
式中：j ––––热功当量（erg ／cal ）。

多孔介质内任一点处的热通量为
()()m L V TV L h q T T C Lq T D L q 01-++∇-=ρλ
(2-2-65)
∑∑===51
5
1
/i i i i i i i k k θλθλ (2-2-66)
式中：λ––––土壤的热传导系数[cal ／（cm ·s ·℃）］，它是通过土壤中各种组分的热传导系数加权平均而求得的。

k i ––––固体颗粒温度梯度与水体温度梯度之比。

Chung [20]提出了一个热传导系数的经验式
()[]C s cm cal b b b L L ︒⋅⋅⎪⎭
⎫ ⎝⎛
++=69
.488/21
321θθλ (2-2-66,)
将式（2－2－63）及式（2－2－65）代入热量连续方程得：
()()()[]m L h HV L h r
q T T C LD T t
T C t T C 0-∇-∇∇+∇∇=∂∂+∂∂ρλ （2-2-67） 其中 h
L
h V
r T H T
H
C C ∂∂+∂∂+=θρ2
T
L
T V h h H h
H C ∂∂+∂∂=θρ2
1
()v v v L C T T H θθ001+-=
()()W L C C T T H L v v v L L ρρρρ----=02
C 为湿土壤的热容量［cal/（cm 3·℃）］。

式（2－2－62）和式（2－2－67）及其相应定解条件组成了水、汽、热耦合求解模型。

对一维问题，方程可简化为［２１］
：
()z k
z T D z z h D k z t T S t h S TV HV T h
∂∂+⎥⎦
⎤⎢⎣⎡∂∂∂∂+⎥⎦⎤⎢⎣⎡∂∂+∂∂=∂∂+∂∂ （2-2-68） ()[]m L HV v h T
q T T C z
z h LD z z T z t h C t T C 0-∂∂-⎥⎦⎤⎢⎣⎡∂∂∂∂+⎥⎦⎤⎢⎣⎡∂∂∂∂=∂∂+∂∂ρλ （2-2-69） 由于方程中系数也为待求函数的函数，这些非线性方程需通过迭代求解［７］。

在求解过
程中，常用的压力和热量（能量）单位有以下几种。

（1）压力单位常采用帕（Pa ）、巴（bar ）、标准大气压（atm ）、毫米汞柱（mmHg ）、毫个水柱（mmH 2O ）等。

（2）热量（能量）常用单位有焦尔（J ）、卡［路里］（cal ），尔格（erg ）。

以上单位的相互转换关系见下表所示。

第六节 土壤水分通量法
一、土壤水分通量法基本原理
土壤水分通量法是直接利用达西定律和质量守恒原理分析计算土壤水通量及潜水的入渗量或蒸发量的一种方法，这种方法在有定位点负压h （基膜势）和含水率观测资料地区应用简便。

根据质量守恒原理，一维垂向土壤水流连续方程（见2-2-6）可写作：
z
q
t ∂∂-=∂∂θ （2-2-70） 上式由z 1至z 2积分得：
dz t
z q z q z z ⎰
∂∂=-2
1
)()(21θ
（2-2-71） 式中：q （z 1）、q （z 2）––––分别表示高度为z 1和z 2处土壤水运动通量。

在t 1至t 2时段内（Δt ），上式可写作：
()()⎰⎰-=-2
1
21
1221,,)()(z z z z dz t z dz t z z Q z Q θθ （2-2-72）
式中：Q （z 1）、Q （z 2）––––分别表示高度为z 1和z 2处通过单位断面面积的水量（t 1至t 2时段内）。

式（2－2－72）表明，当已知时段前后两个瞬时土壤剖面上含水率分布时，仅需已知一断面上土壤水通量即可求得任一断面的通量或水量。

因此，称该方法为土壤水通量法。

由于这种方法是根据时段前后两个瞬时含水率剖面确定水流通量和水量的，在某些情况下，称之为瞬时剖面法。

通量法可分为零通量面法和已知通量法两种。

二、零通量面法
由前述可知，当不考虑溶质势、气压势和温度势时，土壤水的总水势为基膜势Ψm 及重力势Ψg 之和，常用负压水头h 和位置水头z 之和表示。

在测定土壤剖面上基质势和重力势后，可计算出土壤剖面上总水势分布曲线，如图2－2－7所示。

由达西定律，土壤水通量为：q ＝()
z H h k ∂∂-，当z
H
∂∂＝0时，q ＝0，即为零通量面。

图2－２—７中A 、B 两断面均为零通量面，但A 、B 两断面的水流状况是不同的。

对A 断面而
言，水流向上、下断面运移，也可称为发散型零通量面。

而B 断面的上、下断面水流分别向断面B 汇集，故称聚合型零通量面。

图2－2－7中有两个零通量面，这表明降雨（或灌水）入渗与蒸发是交替发生的。

一般土壤较长时间处于单一的蒸发或入渗状态，剖面上可能不存在零通量面。

当剖面上存在零通量面时，可根据式（2－2－72）求任一断面z 处土壤水通量，即：
()()⎰⎰-=z
z z z dz t z dz t z z Q 0
21,,)(θθ （2-2-72’）
也可自d 至z 0间土壤贮水量的变化，如图2－2－8中面积1234求自地表蒸发的水量，
()()⎰⎰-=d
z d z dz t z dz t z d Q 0
21,,)(θθ。

同样，也可根据面积154计算入渗补给地下水的水量。

三、已知通量法
如上所述，当长期处于蒸发或入渗状态时，土壤剖面上并不一定存在零通量面。

在这种
情况下，若能已知某一断面上土壤水通量，则可利用已知通量断面，推求其他断面通量，这
种方法称为已知通量法。

常用的已知通量法有表面通量法和定位通量法［８］。

表面通量法是已知地表入渗量或蒸发量，以地表为已知通量面，推求任一断面通量的方法。

若土壤表面在t 1至t 2时段内入渗量（或蒸发量）为Q m ，则根据式（2－2－72），任一断面z 处单位面积上流过水量为
()()⎰⎰-+=d
z
d z
m dz t z dz t z Q z Q 12,,)(θθ （2-2-73）
式中：d ––––地表距地下水面距离（以地下水面为基准面）。

地表的入渗量可以实际测定，蒸发量可近似地采用公式计算（如可用Penman 公式求蒸散发量）。

定位通量法是在作物根层以下某一特定位置（如地下水面以上一定位置处）上下z 1和z 2安装负压计，测定这两点负压。

如土壤水力传导度k （h ）已预先测定，则可计算这两点间平均断面z 1—2的通量为
⎪⎭
⎫
⎝⎛+∆--=-1)()(1221z h h h k z q （2-2-74）
式中：h 1、h 2––––分别为断面z 1和z 2处负压值。

2
,2
112h h h z z z +=
-=∆ 在已知断面z 1—2通量情况下，可求得t 1至t 2时段内流经断面z 1—2的单位面积土壤水量为Q （z 1—2）。

同样可由Q （z 1—2）求得任一断面流量Q （z ）：
()()⎰
⎰
---+=-2
12
11221,,)()(z z
z z
dz t z dz t z z Q z Q θθ （2-2-75）
使用上述已知通量法需有定位点含水率和负压资料，且需要预先测定土壤水力传导度k(h)值。

土壤的水力传导度常采用表面通量法，根据时段前后两个瞬时的土壤剖面含水率的变化进行测定。

在试验室条件下可用垂直土柱测定土壤水力传导度，在土柱的下边界利用供水装置瞬时供水，并保持一恒定水位。

在野外条件下可在田面灌水后，用塑料薄膜或其他不透水材料覆盖，在地表形成零通量面，每隔一定时间测定土壤剖面上各点的含水率和土壤负压，即可采用位于地表的零通量面法根据时段前后剖面含水率计算各断面的通量q 。

在已求得各断面通量q 后，即可根据任意两相邻断面z 1、z 2处的土壤负压值h 1、h 2自下式计算平均负压
2
2
1h h h +=
时的k(h)值： ()()
11
221±∆-=
-z
h h z q h k （2-2-76）
式中：△z= z 2－z 1。

当z 向下为正时，以上等式右端项分子中取＋1，在z 向上为正时取－1。