降雨和灌水入渗条件下土壤水分运动2.docx

第五章降雨和灌水入渗条件下土壤水分运
动
第一节水向土中入渗过程
一、概述
降雨和灌水入渗是田间水循环的重要环节，与潜水蒸发一样，是水资源评价和农田水分
状况调控的重要依据。

水渗入土壤的强度主要取决于降雨或灌水的方式和强度以及土壤渗水性能。

如果土壤渗水性能较强，大于外界供水强度，则入渗强度主要决定于外界供水强度，在入渗过程中土壤表面含水率随入渗而逐渐提高，直至达到某一稳定值。

如果降雨或灌水强度较大，超过了土壤渗水能力，入渗强度就决定于土壤的入渗性能，这样就会形成径流或地表积水。

这两种情况可能发生在入渗过程的不同阶段，如在稳定灌溉强度（例如喷灌）下，开始时灌溉强度小于土壤入渗能力，入渗率等于灌溉强度；但经过一定时间后，土壤入渗能力减少，灌水强度大于土壤入渗能力，于是产生余水，如图2－5－ 1所示的降雨或灌水条件下的入渗过程。

开始时入渗速率较高，以后逐渐减小。

土壤的入渗能力随时间而变化，与土壤原始湿度和土壤
水的吸力有关，同时也与土壤剖面上土质条件、结构等因素有关。

一般来说，开始入渗阶段，土壤入渗能力较高，尤其是在入渗初期，土壤
比较干燥的情况，然后随土壤水的入渗速率逐
渐减小，最后接近于一常量，而达到稳定入渗
阶段。

在较干旱的条件下，土壤表层的水势梯度
较陡。

所以，入渗速率较大，但随着入渗水渗
入土中，土壤中基模吸力下降。

湿润层的下移
使基模吸力梯度减小。

在垂直入渗情况下，如
供水强度较大，使土壤剖面上达到饱和，当入
渗强度等于土壤饱和水力传导度时，将达到稳
定入渗阶段。

如供水强度较小，小于饱和土壤
水力传导度时，达到稳定入渗阶段的入渗强度将等于该湿度条件下的非饱和土壤水力传导
度。

入渗过程中，土壤剖面上水分分布与土表入渗条件有关。

根据 Coleman和Bodman 的研
究，
当均质土壤地表有积水入渗时，典型含水率分布剖面可分为四个区，即表层有一薄层为饱和带，以下是含水率变化较大的过渡带，其下是含水率分布较均匀的传导层，以下是湿润程度随深度减小的湿润层，该层湿度梯度越向下越陡，直到湿润锋。

随着入渗时间延续，传导层
会不断向深层发展，湿润层和湿润锋也会下移，含水率分布曲线逐渐变平缓。

二、影响入渗过程的条件
降雨或灌水条件下的入渗过程和初始土壤剖面上水分分布与地下水位条件有关， 渗问题的定解条件 有以下几种情况。

（一）初始条件
入渗过程的初始条件一般为
初始剖面含水率或负压分布已知
的条件，即
一般 入
( z,0) i ( z)
(t 0, z 0) h(z,0)
h i (z) (t
0, z （ 2-5-1）
0)
（二）边界条件 1．地表边界条件
（ 1）通过 降雨或灌水使地表湿润 ，但 不形成积水 ，表土达到某一 接近饱和的含水率 ，即（一类边界）
(0, t )
t 0, z 0
（ 2-5-2）
（ 2）降雨和喷灌强度已知 ，且不超过土壤入渗强度， 地表不形成积水 ，即（二类边界）
或
k ( )(
h
1)( ) t
0, z
（ 2-5-3 ）
R t
z
式中： R （ t ）——降雨或灌水入渗强度。

（ 3）当降雨或灌水强度 大于土壤入渗强度 ，地表形成积水 ，成为 压力入渗 。

即（一类
边界）
h(0,t )
H (t ) t 0, z 0
（ 2-5-4）
式中： H （ t ）—— 地表积水深度 。

当地表积水而没有产生径流时，地表水深为H （ t ）；若产 生地表径流，积水深度 H （ t ）可根据来水强度 R （ t ）、土壤入渗强度 i （ t ）及地表径流量 Q
（t ）求得。

2．下边界条件
（ 1）地下水埋深较小，以 地下水位作边界 。

当地下水位变化很小或基本保持不变时， 则地下水面处 土壤含水率为饱和含水率 （地下水面离地面距离为 d ），故
(d , t)
s ,
z d , t
（2-5-5 ）
h(d , t) 0, z
d , t
当地下水面随时间而变化
时，即地下水埋深 d 为时间 t 函数 d （ t ），则地下水面处负压
为零，即
h(d (t), t) 0, z d(t), t 0
（ 2-5-6）
（ 2）地下水埋深较大 的情况， 在计算范围内， 下边界土壤剖面含水率保持初始含水率
，
即
(d ,t ) i ( d) z d , t 0
( 2- 5- 7)
在上述条件下，如 初始含水率上下一致 ， i (z)
i ，得
i ( z)
0 则
z
qD ( )
i
k ( i ) k( i ) z d , t 0
（ 2-5-8）
z
式中： k(θ i )––––离地表距离
d 处断面通量。

（ 3）不透水边界 。

下边界为流量等于零的边界 ，即
q
k(h)(
h
1) 0, h
1, z d , t 0
（2-5-9 ）
z
z
上述表明，研究入渗时边界条件是较为复杂的，所以，计算方法也较为复杂。

第二节 土壤水入渗线性化方程的近似解
在垂直入渗情况下， 一维土壤水分运动的基本方程
可写作：
D
k
（ 2-5-10）
z
z
z
z
如降雨或灌水前的初始含水率（在土壤
剖面上含水率均匀分布 ）为 θ i ，则 初始条件 为
( z,0)
i
（ 2-5-11 ）
在地表 有一薄水层时，表层含水率等于
饱和含水率 θS ；在地下水埋深较大时，计算时
段内入渗水不会到达下边界。

为此， 下边界 土壤含水率不变 ，等于初始含水率 ，则边界条件
可以写作以下形式：
(0, t)
( ,t )
s i
z 0, t 0
（ 2-5-12 ）
z
, t 0
由于式（ 2-5-10）为 非线性方程 （因为扩散度
D （θ ）及水力传导度 k （ θ ）均为待求
含水率 θ 的函数），求解比较困难，为了简化计算，近似地以
平均扩散度 D 代替 D （ θ ），
并以 N
k ( s )
k ( 0 )
代替
k ( )
，则式中（ 2-5-10 ）可简化为
s
2
t
D
N
z
（ 2-5-13）
z 2
式中（ 2-5-13）为 常系数线性方程 ，可以用 拉普拉斯变换求解 。

对式（ 2－5－ 13）采用拉普拉斯变换后可得象函数方程：
2
d N d P
式（ 2-5-14 ）的通解为
i
（ 2-5-14 ）
N 1 N 2 N 1 N 2 2D
D
P z 2D
D
P z
4 D
4 D
z, P C 1 e
C 2 e
（ 2-2-15）
式（ 2-5-12 ）经拉氏变换后，得：
(0, P)
s
（2-2-16 ）
P
( , P)
i
（ 2-2-17 ）
P
根据边界条件式（ 2－ 5－ 16）、式（ 2 一 5－ 17）确定常数：
代入式（ 2-5-15），得象函数的表达式为
N 1 N 2
z
2 D
D
P
z, P 4 D
（2-5-18 ）
i s
P i e
P
进行逆变换后，得 含水率 的表达式为
Nz
z, t
i
s
i
erfc z Nt
e D erfc z
Nt （ 2-5-19 ）
2
2 D t
2 D t
补余误差函数可自表
1－ 2－2 查得。

式（ 2－ 5－ 19）
中 D 可用下式计算：
5/ 3
2
D
D 3 d
（ 2-5-20 ）
5 i
3
s
若已知 D 与 θ的关系式，代入式（ 2－ 5－ 20）积分，
即可求得 D 。

采用式（ 2-5-19 ）求得的 土壤剖面上含水率分布如示
意图 2-5-2 所示 。

由于地表的入渗强度
i
D
k
，为了推
z
求入渗强度，首先根据
的象函 的表达式求
z ：
N
1 N 2
N 1 N
2
P z
s
i
P
2D
D
4 D
（ 2-5-51）
z P
2D
D
4D
e
地表处， z=0，则
s i
N
1 N
2 P
（ 2-5-21 ’）
z P
2D
D 4D
在 入 渗 初 期 ， t 0.2
D
， 相 当 于 P 20
N 2 时 ， P N 2 P ，
N 2
4D 4D
N P P
，则式（ 2-5-21 ）可近似写成：
2D
D
D
s
i
P si
1 z
P D
D
P
经逆变换得：
1
s
i
t
2
z
D
入渗初期： ［当 D （ θ ）取平均值 D 时］
1 i
D
k s
s
i
D
t 2
k s
z
入渗时间较久，即当
P
1 N 2
之，相当干 t 80
D
20 4D N 2 N 2
N 2 1 N 2 P
N 代 入
P
，
D 4D
2D
4D
4D
（2-5-21 ’）式，则
0 ，所以
z
z
则入渗时间久时 ，入渗强度（ i → k s ）为
iD
k
s
k
s
（ 2-2-24）
z
自式（ 2－ 5－ 23）、式（ 2－5－ 24）得入渗速度在时间上的变化过程如图（ 2-5-3）所示。

（ 2-5-21 ’’）
（ 2-5-22 ）
（ 2-5-23）
时
i (cm/min)
Ks
o
t (min)
图 2-5-3 入渗率随时间变化图
第三节 Green － Ampt 模型的入渗解
Green － Ampt 模型 ［50］
是 1911 年提出的一种 简化的入渗模型 ，它是建立在 毛管理论基
础上的一种入渗模型。

假定土壤是由一束直径不相同的毛管组成， 水在土壤入渗过程中， 湿
润锋面几乎是水平锋面 ，且在锋面上各点的吸力水头均为
S m 。

锋面后面的土壤含水率为均
一的 ，如图（ 2-5-4 ）所示。

所以 k （ θ ）也为常数，这种模
型又称 活塞模型 。

根据达西定律：
qk J
H
S m z
k
(2-5-25)
z
式中： H ——地面以上水层厚度；
S m —— 锋面处土壤负压； z 一锋面推进距离。

式（ 2-5-25 ）为 单位时间，单位面积流入土体的水量。

根据水量平衡原理，应等于
土体内增加的水量 q
dz ，即
dt
dz k
s
H S m z （ 2-5-26 ）
dt
z
s
式（ 2-5-26 ）积分：
所以
k
s
z H S m ln
H S m z （2-5-27 ）
t
H
S m
s
式（ 2－ 5－27）为 z ～ t 关系式，原则上可以求得任何时刻
t 时入渗锋面所达到的位置，
当然也就不难求得该时刻的
累计入渗量 ：
W i
s
0 z
（ 2-5-28 ）
H → 0 时，式（ 2-5-27 ）可写作：
t
s
z S m ln
S m z
(2-2-27 ’)
k s
S m
或由式（ 2-5-26）
dz
k
s
H
S m z
H+S m +z 项中 z 略去，所
dt
z
，当 t 很小时，该式的
s
以 H
S m z 。

此时
积分得
t 时入渗总量
2k
s
H S m t
（ 2-2-27 ’’）
z
s
I is0
z
2k s s 0 H S m t
（ 2-5-29 ）
式（ 2-5-27 ’’）表明，入渗初期，入渗深度 z 与
t 成正比。

将 I 对 t 求导，得：
dI
k s H
S m
1
t 2
（ 2-5-30 ）
i
s0
dt 2
s
当 t 大时 ，式（ 2 -5-26 ）中 Z >>H+S m ，因此
H
S m z 1 ，则由式（ 2－ 5－25）可
z
知：
i k s
（ 2-5-31 ）
即入渗强度近似等于土壤饱和渗透系数。

第四节水平入渗条件下的Philip解法
水平入渗条件下的Philip 解[51]是一种半解析法，即前半部用解析法，利用博茨曼（B oltzmann ）变换，将偏微分方程转换为常微分方程；后半部采用迭代计算，求解常微
分方程。

由于求解过程中未作过分简化，求得结果较为严密。

一、水平入渗的常微分方程推求
水平入渗的基本方程
t x D
x
i t0, x0,（ 2-5-32 ）
0t0, x0,
lim i t0, x,
x
将式（ 2-5-32 ）中基本方程改写为以坐标x（θ，t）为变量的方程，根据第二章中方程（2-2-17 ），在水平入渗时应为
x x
）
D（ 2-5-33
t
采用 Boltzmann 变换，引入变量从λ （θ），且令
x( , t )t 则x ,t 1
2（ 2-5-34 ）
1
t 2（ 2-5-35 ）
1
故x 1
t2
x
t 2（ 2-5-36 ）
1
t 2（ 2-5-37 ）
将式（ 2－5-36）、式（ 2－ 5－ 37）代入式（ 2－ 5－33）得：
经整理后得微分方程
1 d D d
（ 2-5-38）2d d
由边界条件已知：
为了求得λ～θ的关系式，将式（2－ 5－ 38）常微分方程自θi至θ 积分得：
i d2D
d
（ 2-5-39 ）
d
式（ 2-5-39 ）为λ～θ关系式，若已知 D（θ ）关系式代入上式即可求得λ～θ 关系式，
1
因xt 2，即可由λ ～θ求得θ（x，t），从而求得剖面上任何，任何距离的含水率
分布。

若能λ（θ ），代入上式可求得D（θ）关系曲。

二、迭代计算法
在 D（θ）已知，式（ 2－ 5－ 39）λ ～θ关系式，它的关系曲如2－ 5－ 5 所示，
Philip 的送代法是将θ0→θ i 等份成 n 份，步θ =（θ0-θ i）/n，若等分点按序其含
水率θ 1，θ 2，⋯，θ n-1，θ n。

相各等分点相λλ 1，λ 2，⋯λ n，两个相等分
点中点的含水率θ 1/2，θ1+1/2，⋯，θ （n-2）+1/2，其相散率D
1/2
，D
1+1/2
，⋯D（
n-2
）
+1/2 ，任一点含水率θr，可用下式算
r0r
（ 2-5-40）
D r+1/2取平均，定
D r r
1Dd d（2-5-41）
r
2r 1r 1
写出θ1/2，θ1+1/2，⋯θ（n-2）+1/2各点的式（ 3-5-39 ）：
1
点
2
因为
所以
所以
1
点
1
2
所以
任一点1
d
n
1
2d
1
2D 1
2
11
n
2d
21
r r 1
2D 1
201
2D 1
21
1
d
2
n
2D1
1
21
2
2D 111d
2
1n
2
2D1
n 2令
r
2
1
点
2
n 1
J1
r
2
r 11d
r 1
22
n
n 22D11
n 2d
2
r 2（ 2-5-42）
2
n
r
d r
2
（ 2-5-43 ）
n
当r=0 时，J
1
2
d（ 2-5-44 ）n
式（ 2－ 5－43）表示λ～θ关系曲线所包含的面积。

式（ 2-5-44 ）表示θ =θr+1/2时，θr+1/2～θn之间λ～θ 曲线所包含面积的近似值。

由图可知面积为
1r
r d d 2
n n r r 1
（ 2-5-45）2
式中：δ r+1—––为λ ～θ 的曲线与λ r，λ r+1/2之间所包含面积减去r矩形面积后的一小
2
块面积。

由式（ 2－5 一 43）及式（ 2－ 5－ 45）得：
当θ 取得足够小时，则δ r+1可以忽略，所以
J 11d
r（ 2-5-45 ’）
2
r n
2
则式（ 2 一 5 一 42）可以改写为
1
2D
12
J 1
2
2
2D
1
1
2
J 1
1
2
（ 2-5-46 ）
r 1
n 1
2D
2D
1
r
2 J
1
2
n 2
1
2 J 1
2
由 形可知：
式中：
r ––––矩形面
r
与 λ～ θ 曲 与 λ r －1/2、 λ r 之 包含的面 之差。

2
上式可写
J
1 J
1
r
rr
r
2 r
1
2
当
θ 足 小 ，同 ，（
r -δ r ）可忽略。

J r
1 J r
1 1
r
2
2
根据式（ 2－ 5－ 47），我 可以得到：
J 1
d
2
n
J
1 J 1
1
1
2 2
J
1 J 1
2
2
2
2
2
J 1 J
1
r
r
r 1
2 2
J
1
J
1
n 1
n 1
n 2
2
2
（ 2-5-47）
（ 2-5-48）
（ 2-5-49 ）
根据式（ 2－5－ 46）和式（ 2－ 5－ 49），若已知 J 1/2及 D （ θ） 可交替使用此二式，求
得λ 1，λ 2⋯ λ r ⋯，λ n-1。

因此， 了求得 λ～ θ 关系 必 首先求得
J 1/2及 D ～ θ 关系曲 。

由上述可知 J 1/2 代表了 λ～ θ 曲 所包含面 的近似 ，由式（ 2－ 5－ 49）可知：
若将面 看作一系列的 （
θ 一 θ n ），高 d λ的水平方向的矩形，
J 1
n d
（ 2-5-50 ）
2
式中： λ ––––θ和 D （ θ）的函数。

对于初始值 J 1/2 ，我们可将 D （ θ ）考虑为在整个含水率变化范围内的平均值，由式（
2
－5－ 50）可知：
若用加权平均值来代替上式中
D 值，以 D ’表示：
D
0 n D d
d
（ 2-5-51 ）
2
n
n
n
式中：
n —— D
的加权因子。

n
其中等号右端系数
2 是为使等式衡等， D ’=D ’（若将 D ’代替 D （ θ ）即可算得该系
数）。

将式（ 2-5-51）写成 有限差 的形式：
2
n
D
r
n
D
r
（ 2-5-52 ）
2
n
r 0
1
因 D ’为常数，直接 利用入渗问题线性化解
，以
xt 2 代人得：
n 0 n
erfc
1
（2-5-53 ）
2D 2
将式（ 2-5-53）代入式（ 2-5-50）得：
J 1
0 0
n erfc
1 d
（ 2-5-54 ）
2
2D
2
式（ 2-5-54 ）为 J 1 的初始近似值第一次迭代值，因此可写作
(J 1 )1 。

2
2
将 α =λ /2（ D ’） 1/2 代人上式，且令积分上限为 β 0=β →∞以 n · Δθ 代替（ θ 0－ θ n ）
则式（ 2-5-54 ）可写作：
( J 1 )1
n
2 D
1
erfc d
2 0
（ 2-5-55 ）
2
为了求解式（ 2－ 5－55）积分式，利用
erfc α =1－ erf α ，则
1
1
erf d
erf
2
已知
2 e
2
1 1
erfc d
erf
2
得
2
e
2
(2-5-56)
1 2
当
时 ,
1 erf
2
e
，所以式（ 2－5－ 56）为
1
erfc
d
2 ,
(2-5-57)
将式（ 2-5-57）代入式（ 2-5-55）得：
D
1
( J 1 )1 2n
2
(2-5-58)
2
有了 J 1/2 就可以由式（ 2 一 5－ 56）和式（ 2－ 5－ 49）重复交替计算求得 λ ～θ 关系曲
线，但由于 J 1/2 是估算值，存在偏差，求得的 λ ～ θ曲线必然也有偏差，为了使迭代计算尽 快收敛，采用两种不同方法计算 J n － 1/2，不断地修正 J 1/2，以使两种方法求得的 J n － 1/2 值足够
的接近，以便得到一个
比较精确的 λ ～ θ曲线 。

第一种方法采用
J F 1
J
1 1 n 1
（ 2－5－ 59）［即式（ 2－ 5－ 49）］
n
2
n
2
第二种方法采用
S
J
1
2
n 1
（ 2－ 5－ 59）［即式（ 2－ 5－ 49）
d n 1
n
2
J S 式中，由于 θ变化于 θn ，θ n-1 之间，若 θ 取得足够小，则 D 可以看作为一常量
D n
－
1/2：
1
n 1 d (2-5-61)
D
n
式（ 2-5-61 ）表明，在 θ n 和 θ n － 1 范围内将 D 视作常量，当 θ 取得足够小时，不会引
起明显误差。

上述问题可为求解半无限长水平土柱的入渗问题，即
2
t
D
1 x 2
n 2
i
t 0, x
(2-5-62)
1
n 1
i
t
0, x n 1t
2
t
0, x
根据土壤水入渗的线性化方程的近似解法，式（
2－ 5－ 62）定解问题的解为
2D
1
n 1
n
2 A y
d
n 1
(2-5-63)
n
n 1
式中 A （ y ）定义为
2ye y 2
2 y 2
A y 1
(2-5-64)
2
erf y
1
其中
y
n 1 /[ 2( D 1 )
2
]
n
2
A(y) 的函数表见表（ 2-5-1）。

对 y 较大时， A （ y ）可以渐近级数表示：
A y
1
2 10 74
706
2 y
2
2 y
2 2
2 y
2 2
2y
2 2
(2-5-65)
（ 2-5-66 ）
λn-1由式（ 2－ 5－ 59）求得，某一个λ～θ曲λn-1及 D n-1/2常量，所以， A （ y）也常数，可采用式（ 2－ 5－64）算，根据 y 求 A （ y）。

首先，已知
erfc y 1erf y12y2
（ 2-5-67）1
e d
2
式中的β 虚量，可以改写
2
e2 1d（ 2-5-68）
y
2
由式（ 2-5-68 ）求得 erfc（ y），代入式（ 2-5-64 ）即可求得 A （ y）、 A （ y）～ y已列于表（ 2-5-1）中，所以根据某 D n－1/2，λn－1求得 y，表 2－ 5－ 1求得。

将式（ 2－5－ 63）代人式（ 2－5－ 60）得：
12D1
J S1
n
2 A y
n 1
（2-5-69 ）
n
2n1
2
当θ 足小，且假定J1/2是正确的， J F n-1/2与 J S n-1/2是相等的，但由于算
中采用了近似两者差 1 ；，若在所要求精度范内即可，否行第二、第三、⋯⋯，
更多次的代。

由1与（ J1/2）1；可以得到一个更接近于的初始（J1/2）2。

重复上述程行第二个循的算。

（J1/2）2一般可写作：
迭代 P 次后：
J 1J 1P 1（ 2-5-71 ）
2
2P2P 1
式中： N––––式（ 2－ 5－ 70）两相等（或差在所要求的精确度范内）送代次数。

了加速收，J1/2有另一改公式：
P2,3, , N（ 2-5-71 ）由于在其他著作中采用了与上述不同符号，里介一下用其他符号表示的Philip 的迭代公式。

Philip 定 I r+1/2
1
r
1
（ 2-5-72 ）
I 1d2r
2
n
r
所以，其他以 I 参数的各式中均比J 参数的多除一个θ。

即I 1J 1 /（ 2-5-73 ）
r r
22
迭代算步：
（ 1）确定θ0至θn的分段数 n，步θ =（θ0 -θn）。

2n
D r
（ 2）以D2rn算平均散度。

0n r
’ 1/2
算第一次送代初 （
J 1/2）1，并由公式（ 2
（ 3）根据公式 J 1/2 =2n ·Δ θ（ D /π ）
－5－ 46） 算 λ 1，由（ J 1/2） 1 和λ 1 即可 算（ J 1+1/2），再代入公式（
2-5-49 ）求 λ 2，依次
推，最后求得 ( J F
1 )1 。

n 2
（ 4）由以上 算同 可求得
λ n-1和 D n-1/2 ，代入公式（ 2-5-60 ） 算 ( J S
1 )1 。

n 2
（ 5） 算 差
1=
( J
F 1
)
1 －
( J S 1 )1 ，若 1 足要求，以上 算 束，否 要 行第
n 2
n 2
二、三、⋯次的送代，直至 足要求 止。

（ 6） 行多次送代 ，重复以上（
2）～（ 5）步 ，直至迭代 差
很小，且 λ 1， λ
2，⋯、不再随 J 1/2 改 止。

若 算 果不理想，可能是 用
θ不合适， 可重新 用
Δθ ，重复以上步 行迭代 算。

1
1
（ 7） 算出 λ ～ θ关系曲 后，因
xt 2 ，故 x t 2 ， 于 定的 t ，由 λ 可
算 x ，即得 θ ～x 关系曲 。

也即求得了 θ 在剖面上的分布 。

Philip 水平入渗解的理 可用于 水平土柱法 定土壤水 散度 。

土壤水分在 的（理
上 是半无限 界） 均 土柱中 生水平运 ， 一 土壤水平运 方程如式 （ 2－ 5－ 32）
所示， Boltzmann
可得 [如前述方程（ 2－ 5－ 39）] ：
x
d
D x
1 i
（ 2-5-39 ）
2
d
d
x
式中： λ （θ ） —––x ， t 的函数。

1
当水平土柱首端供水， 根据 水分运移 料， 即可 出 λ ～θ 曲 ，（
( ) xt 2 ），
x )d 由初始含水率 至 θ x
2-5-5 所示。

(
θ i 之 λ ～ θ 曲 所包含的面 ，如
i
( d
)
x
在 λ ～ θ曲 上 含水率 θx 点的斜率。

因此，通 水平土柱 即可求得土壤
d
水 散度 D 。

第五节 垂直入渗条件下的
Philip 解法
一 垂直入渗
基本方程写成 z （θ ， t ） 函数
z
D
dk t
z
（ 2-5-74 ）
d
一 界 定解条件
Philip 级数解的形式
i 0 i [52]
z 0, t
0,
z 0, t 0
（ 2-5-75 ）
z
, t
1
3
i
z ,t
1 t
2
2 t
1
3
t 2
4
t 2
i
t 2
（2-5-76 ）
i
1
相应边界条件 为
1
020
（ 2-5-77 ） i
i
（ 2-5-78）
若求得式（ 2-5-77 ）中系数 η i （ θ ）值，则可求得在一定 θ 时的位置 z 。

以下我们以待定
系数法求解 η i （ θ）。

式（ 2－ 5－74）可改写为
2
z z 2
2
z dk z
dD z
D 2
d
t
d
将 z 的解式（ 2-5-76 ）分别代入以上各项：
z
1
3
i
1
t
2
2 t
3
t
2
i
t 2
i 1
2
z
1
3
i
1
t 2
2
t 3
t 2
i
t 2
2
i 1
2
3
z
1
2 t
2
1
2
t 2 2
1
3
2
2
t 2
5
2
1 4
2
3
t 2
z
1
3
1
1
t
2 2
3 t 2
2 4
t t
2
2
1 i
2
i
t
2
i 1 2
将式（ 2-5-80）代入式（ 2-5-79）归并后得：
1
3
i
1
t 2
2
t
3
t 2 i
t 2
i 1
式中前四项系数 η 1，η 2 ，η 3， η 4分别为
（ 2-5-79 ）
（ 2-5-80 ）
（ 2-5-81 ）
1D11
1
2D
1 2
D k22 22112
112
D 2
D2k 332
3312
21
2 212113（ 2-5-82 ）
1
22
D3
2
D k 2 12
4432
223
14123
22
13222
141231
D 4
由于式（ 2－ 5－ 8l）左端等于零，所以η 1=η 2=η3=η 4⋯=0η1=0，即
全式除以1 2 得
即
d D1
1 ( )
d d ()2
d
式（ 2－ 5－83）表明，η1（θ）即λ（θ ）水平入渗解。

同理，以η2=0求得：
d
2
d 2
k k i（ 2-5-84 ）
i dD d
i d 1
η3=0得：
22 3
dD
d d 3d d
2
3
d 1d
D
d
i
d 1
η4=0得：
3
1
（ 2-5-85）
22
2
d d 2d d
4dD d D d
i d 1 d 1
2
1
d
2
d
3
2
d
d
2
1
（ 2-5-86 ）根据以上各式，若已知η 1（θ ）可逐步求得外η 2（θ ），η 3（θ ），⋯，即可求得
方程（ 2－ 5－ 74）的解。

由于收快，求得前四就足精确了。

上述η 1（θ），η2（θ ），η 3（θ ），η 4（θ ）都复，η 1（θ ）即水平入渗
的λ，了求解些并不容易，此 Philip 采用了推公式求解。

以下以求得η2例明求解方法。

为 了 求 解 系 数 ， 我 们 可 以 写 出 与 式
（2-5-39 ）类似的近似系数一般方程，则
df fd
（ 2-5-87 ）
n
d
式中： α ，β — 均为 θ 的已知函数。

式（ 2 － 5－ 87）左边可以写作式（ 2 － 5
－88），如图 2－ 5－ 6 所示。

r 1
r r 1
fd （ 2 － 5
2
fd
fd
2
n n
r
－ 88）
式（ 2 一 5－ 88）中右端第二个积分为 由图 2 － 5 － 6可知，该积分近似等于
h
，其中 h 为四边形 BCDE 的平均高度。

2
θ取得足够小时， DE 可以看作为一直线，在梯形
DFAC 中 h 可以下式表示：
h r
f r
1 f r
1 f r 1
3 f r （ 2-5-89 ）
f r 1
4
4
1
1
f r 1
3 f r
r 2 fd
1 r
n
4
2
r
（ 2-5-90 ）
1
1
3 f
r
f r 2
4
式（ 2-5-87 ）中，若 θ =θ r+1/2，则方程为
r 1 df
n 2 fd
d
r
1
（ 2-5-91 ）
1 r
2
α r+1/2 为 θ r+1 <θ<θ r 范围内平均值，且以
f r 1
f r
来近似（ df/d θ ）r+1/2 （其负号
表示 θ 值增加的方向与 f 值增加方向相反），代入式（
2－ 5－91）得：
r 1 fd
r 1
f
r 1
f r
r 1
（2-5-92 ）
n 2
2
2
将式（ 2－5－ 89）、式（ 2－5－ 91）代入式（ 2－ 5－ 88），经整理后得
I 1 1 r
1
f r
令
fd
（ 2-5-94）
2
n
2
r
将式（ 2 一 5－ 94）代入式（ 2－ 5－ 93）经整理后得：
1
r
2
I
1
r
2
f
r 1
f r
（ 2-5-95 ）
r 1 1 2
2
8
式（ 2－ 5－95）为 Philip 的第一个递推公式。

第二个递推公式，将
r 改为 r — l ，则式（ 2－ 5－ 94）为
1
r 1 1
I
r 1
n
2
fd
f r 1
（ 2-5-96）
2
2
将式（ 2－594t ） 4）与式（ 2－ 5－ 96）两式相减，得：
1
1 r
1
f r f r 1
r 2
fd
fd
（ 2-5-97 ）
r 1 r
2
1
2
根据式（ 2－ 5－ 90）得：
将式（ 2－5－ 98）代入式（ 2－5－ 97）经整理得：
当 θ 足够小时， f r ， f r+1/2 ， f r-1/2 ， f r-1 非常接近，则上式可写作：
或
I
1I
1f
r
（ 2-5-99 ）
r
r
2 2
式（ 2－ 5－99）即为 PhiliP 的第二个速推公式。

根据该公式，令
I n I 1 1
f n
（ 2-5-100 ）
2
2
n
由式（ 2－ 5－ 94），若 r=n — l ，则
1
n 1
1
（ 2-5-101）
I 1 fd
f n 1
n
n
当 r=n 时， f n 并不趋于无穷大，则由式（ 2－ 5－ 101）和式（ 2－ 5－ 99）可知， I n =0。

以上两递推公式如果 已知初始和边界条件 即可 交替使用 ，求得 f ～θ 曲线。

如果边界条
件为 θ =θ 0 时， f=0 ，且 f 0=0 ，则可由式（ 2－ 5－ 95）求得 I 1/2 ，由式（ 2－5－ 99）又可求得
I ，当 r= 1 时，代入式（ 2－ 5－ 94）求得 f 2。

重复以上过程即可求得 f ～θ 曲线。

当 f= λ ， η， ψ ，ω 等不同项的系数时，都可求得对应的与 θ 关系曲线。

在某一 t 时，可以求得不同
θ的 z 值，则可求得对应于不同
t 时的 θ～ z 关系曲线，即求得任何时间剖面上含水率分布
θ（ z ， t ）。

入渗量计算如下：
渗入土壤的水量应等于流出水量与剖面上含水量变化之和， 在半无限情况下， 底部含水率将保持初始含水率，即为 θ n ，导水率也为常量 k ，因此由下部流出水量
W 出 k n t
（2-5-102）
如将土柱剖分为
n 等份，步长为
z ，取横断面为 1 个单位时，经 t 时土柱顶部含水率
为θ 0，对半无限土壤剖面含水量增量近似表示为
精确值可用积分表示：
上式表示 θ=θ n 以上 θ ～ z 曲线所包含的面积，这面积也可表示为
zd
（2-5-103 ）
W
n
式（ 2－ 5－102）与式（ 2－ 5－103）相加即为 t 时间内渗入剖面总水量，即
I
0 k n t
zd
（2-5-104 ）
n
1 2 3 4
将 Philip 解的表达式 z
,t
1t
2
2t
2
3 t
2
4t 2
代入式 （ 2－ 5－ 104），并
加以整理得，
1
3
I
At 2
k n t Bt
Ct 2
Dt 2
（ 2-5-105）
式中： A
0 B 0 d
1d
2
n
n
A 、
B 、
C 、
D 等均为常数。

式（ 2－ 5－105）中右边第一项相当于水平入渗，所以水平入渗量为
1
I k At 2
（ 2-5-106）
入渗率的计算如下： 水进入土壤的速率即为
z=0 处的水流通量，即为
I 以式（ 2－ 5－ 105）代入：
i
1
At
2
1 3 1
2
k n B
Ct 2 2Dt （ 2-5-107）
2
按式（ 2－5－ 107）求得的入渗率当 t 较小时误差较小，但时间长后，计算值与实验资
料偏离较大， Philip 对其进行了修正，但修正式并不适用于 t 小时情况。

Miller 和 Klute 用
Hagan 等人的结果，用下式计算与实验值接近。

1
i
1
i f
（2-5-108 ）
At 2
2
式中： A ， i f —— 均为土壤特性常数。

A 值大小决定于土壤初始含水率，一般称为吸水率。

i f 为稳定入渗速度，相当于土壤饱
和时的水力传导度（即渗透系数） k s 。

入渗初期，入渗速度很大，
k s 远较 1/2At -1/2 为小，可
忽略不计。

随着时间增大，入渗速度迅速减小，在入渗时间很久，
t →∞，式（ 2－ 5－ 108）
中右端第一项趋近于
0， i=i f ，即为稳定入渗速度。

入渗速度在时间上变化如图
2-5-3 所示。

A 值可通过初期入渗总量 I h 确定：
1
I h
At 2
（ 2-5-109）
1
故
A
I h t 2
（ 2-5-109）
入渗率理论公式中的常数需通过试验确定。

在生产中常采用经验公式计算入渗率
i 和入
渗量 I 。

第六节土壤水入渗的经验公式
前节介绍的理论公式，既可以计算土壤的入渗速度也可以计算入渗条件下土壤含水率的分布，但这些公式都是在简化的条件下求得的，公式计算比较复杂，且常数也往往需要通过实验确定，应用比较困难。

在生产实践中为了计算入渗强度和入渗总量常采用经验公式，以下介绍两种应用较多的经验公式。

一、考斯恰阔夫（Костяков.H）经А验公式
Костяков入渗公式的形式
为
[53]
i i1t（ 2-5-110）
式中： i1––––第一单位时间的入渗速度，决定于土壤质地和初始土壤含水率；
α ––––经验指数，变化干 0．3～0．8之间，轻质土α值较小，重质土α值较大，初始含水率愈高，α值愈小，一般土壤α值常取 0.5。

在时间 t 内入渗总量I（以水层深度表示）为
I t t i1i1
t1（ 2-5-111）
idt
0 t a dt
1
当α =0.5 时
1
i i1t 2（2-5-110 ’）
1
I2i1t 2（ 2-5-111 ’）以上两式与入渗初期的线性化方程入渗解（2－ 5－23）、 Green－Ampt 入渗解式（ 2－ 5－29）和 Philip 入渗解式（ 2－ 5－ 108）均具有相同形式。

公式（2－ 5－ 110）在入渗时间较久时，入渗强度i 趋近于 0，不符合地下水埋深较大时入渗情况。

这一公式多用于入渗初期
或入渗时间较短的农田灌水时入渗计算。

二、 Horton 入渗公式[50]
Horton 的入渗公式形式为
i i f i 0 i f e t（2-5-112）
式中： i f––––时间较久时的稳定入渗率；
i 0––––初始入渗速度；
β ––––反映土壤特性的常数。

Horton 公式表明，入渗初期入渗速度i不是趋于无穷大，而是一有限入渗速度i 0，当入渗时间久时，趋近于一个稳定入渗强度i f。

这与实际的入渗过程比较符合，但计算中需要有三
个经验常数。

本章所介绍的理论公式都是在表层边界条件简单，均质土壤的情况，且初始土壤剖面上含水率分布均匀，地下水埋深较大的条件下求得的。

在入渗边界随时间而变化，非均质层状
土壤，初始土壤剖面上含水率非均匀分布，地下水埋深较小对入渗有显著影响的情况下，需采用第六章介绍的数值方法求解。

本节所介绍的入渗计算经验公式各有其适用条件，在生产实践中，必须根据具体情况选用。