降雨和灌水入渗条件下土壤水分运动2.docx
《降雨和灌水入渗条件下土壤水分运动2.docx》由会员分享，可在线阅读，更多相关《降雨和灌水入渗条件下土壤水分运动2.docx（31页珍藏版）》请在冰点文库上搜索。
降雨和灌水入渗条件下土壤水分运动2
第五章降雨和灌水入渗条件下土壤水分运动
第一节水向土中入渗过程
、概述
降雨和灌水入渗是田间水循环的重要环节，与潜水蒸发一样，是水资源评价和农田水分
状况调控的重要依据。
水渗入土壤的强度主要取决于降雨或灌水的方式和强度以及土壤渗水性能。
如果土壤渗水性能较强，大于外界供水强度，则入渗强度主要决定于外界供水强度，在入渗过程中土壤表面含水率随入渗而逐渐提高，直至达到某一稳定值。
如果降雨或灌水强度较大，超过了土壤渗水能力，入渗强度就决定于土壤的入渗性能，这样就会形成径流或地表积水。
这两种情况可能发生在入渗过程的不同阶段，如在稳定灌溉强度（例如喷灌）下，开始时灌溉强度小于土壤入渗能力，入渗率等于灌溉强度；但经过一定时间后，土壤入渗能力减少，灌水强度大于土壤入渗能力，于是产生余水，如图2-5-1所示的降雨或灌水条件下的入渗过程。
开
始时入渗速率较高，以后逐渐减小。
土壤的入渗能力随时间而变化，与土壤原始湿度和土壤
结构等因素有关。
一般来说，开始入渗阶段,
图2-5-1入滲率对时间的茏祭曲线
水的吸力有关，同时也与土壤剖面上土质条件、土壤入渗能力较高，尤其是在入渗初期，土壤比较干燥的情况，然后随土壤水的入渗速率逐渐减小，最后接近于一常量，而达到稳定入渗阶段。
在较干旱的条件下，土壤表层的水势梯度较陡。
所以，入渗速率较大，但随着入渗水渗入土中，土壤中基模吸力下降。
湿润层的下移使基模吸力梯度减小。
在垂直入渗情况下，如供水强度较大，使土壤剖面上达到饱和，当入渗强度等于土壤饱和水力传导度时，将达到稳定入渗阶段。
如供水强度较小，小于饱和土壤水力传导度时，达到稳定入渗阶段的入渗强度将等于该湿度条件下的非饱和土壤水力传导度。
入渗过程中，土壤剖面上水分分布与土表入渗条件有关。
根据Coleman和Bodman的研究,
当均质土壤地表有积水入渗时，典型含水率分布剖面可分为四个区，即表层有一薄层为饱和
带，以下是含水率变化较大的过渡带，其下是含水率分布较均匀的传导层，以下是湿润程度
随深度减小的湿润层，该层湿度梯度越向下越陡，直到湿润锋。
随着入渗时间延续，传导层
会不断向深层发展，湿润层和湿润锋也会下移，含水率分布曲线逐渐变平缓。
、影响入渗过程的条件
般入
降雨或灌水条件下的入渗过程和初始土壤剖面上水分分布与地下水位条件有关,渗问题的定解条件有以下几种情况。
（一）初始条件
的条件，即
入渗过程的初始条件一般为初始剖面含水率或负压分布已知
（2-5-1）
乜（z,0）=d（z）（t=0,Z〉0）
、h（z,0）=hi（z）（t=0,z>0）
（二）边界条件
1.地表边界条件
（1）通过降雨或灌水使地表湿润，但不形成积水，表土达到某一接近饱和的含水率，即（一类边界）
二（0,t）-v0t•0,z=0
边界）
式中：
生地表径流，积水深度H（t）可根据来水强度R（t）、土壤入渗强度i（t）及地表径流量Q（t）求得。
2.下边界条件
（1）地下水埋深较小，以地下水位作边界。
当地下水位变化很小或基本保持不变时，则地下水面处土壤含水率为饱和含水率（地下
水面离地面距离为d）,故
（2-5-5）
$（d,t）=乞，z=d,t=0
h（d,t）=0,z=d,t>0
当地下水面随时间而变化时，即地下水埋深d为时间t函数d（t）,则地下水面处负压为零，即
（2-5-6）
h（d（t）,t）=0,z=d（t）,t0
（2）地下水埋深较大的情况，在计算范围内，下边界土壤剖面含水率保持初始含水率，
即
在上述条件下，如初始含水率上下一致，耳（z）-耳,得「i（z）=0则
胡i
q=—D（旳【k（m）二k（R）z=d,tO（2-5-8）
cz
式中：
k（0i）--离地表距离d处断面通量。
（3）不透水边界。
下边界为流量等于零的边界，即
hh
q--k（h）
（1）=0,1,z=d,t，0（2-5-9）
dzdz
上述表明，研究入渗时边界条件是较为复杂的，所以，计算方法也较为复杂。
第二节土壤水入渗线性化方程的近似解
在垂直入渗情况下，一维土壤水分运动的基本方程可写作：
化行d（日严（2-5-10）
:
z:
z:
z:
z
如降雨或灌水前的初始含水率（在土壤剖面上含水率均匀分布）为0i，则初始条件为
二（乙0）-（2-5-11）
在地表有一薄水层时，表层含水率等于饱和含水率0s；在地下水埋深较大时，计算时
段内入渗水不会到达下边界。
为此，下边界土壤含水率不变，等于初始含水率，则边界条件
可以写作以下形式：
E（0,t）“sZ=O,t>0丿（2-5-12）
日（°°,t）=©ZT旳，tA0
由于式（2-5-10）为非线性方程（因为扩散度D（0）及水力传导度k（0）均为待求
平均扩散度D代替D（0）,
含水率0的函数），求解比较困难，为了简化计算，近似地以并以⑴孚严代替字，则式中（25）可简化为
2•
（2-5-13）
（2-5-14）
otj—C0少
D「-N-:
z:
z
式中（2-5-13）为常系数线性方程，可以用拉普拉斯变换求解对式（2-5-13）采用拉普拉斯变换后可得象函数方程：
2——
dNZP-
—.—U
dz2DdzD
式（2-5-14）的通解为
（2-2-15）
（2-2-16）
（2-2-17）
根据边界条件式（2—5—16）、式（2一5—17）确定常数:
代入式（2-5-15），得象函数的表达式为
N~
--iE
（2-5-18）
二Z,P's-
PP
进行逆变换后，得含水率的表达式为
eJ皿〕+e畔erfc
7十Nt
2
l2^Dt丿
l2』Dt/
1—2—2查得。
式（2—5-19）
（2-5-19）
补余误差函数可自表
中D可用下式计算:
5/3-
D」D…f®
入-71。
3
（2-5-20）
若已知D与0的关系式，代入式（2—5—20）积分,
即可求得D。
采用式（2-5-19）求得的土壤剖面上含水率分布如示意图2-5-2所示。
由于地表的入渗强度i二-Dkb〕，为了推
cz
圉2-^-2人邃条件下，剖血上含水率分布示总图
求入渗强度，首先根据二的象函二的表达式求—
（2-5-51）
-Z
2D
'2
NP
D（4D
在入渗初期，t:
：
0.2-Dy,
N2
N2
相当于P204D
（2-5-21'
P4：
4D
，则式（2-5-21）可近似写成:
宀基7PU1
；z
P■.DD、P
经逆变换得:
1
胡_（兀-3）严
；z、、二D
入渗初期：
［当D（0）取平均值D时］
（2-5-21'）'
（2-5-22）
i=-Dk忑…7
Sz
（2-5-23）
入渗时间较久，即当
P」也之,
204D
t时
N2
4D
N2
4D
cQcQ
（2-5-21'）式，贝U=0,所以=0
.1-..
:
zz
则入渗时间久时，入渗强度（iTks）为
i—Dk忑二k忑（2-2-24）
cz
自式（2—5-23）、式（2—5-24）得入渗速度在时间上的变化过程如图（2-5-3）所示。
第三节Green—Ampt模型的入渗解
Green—Ampt模型5^是1911年提出的一种简化的入渗模型，它是建立在毛管理论基础上的一种入渗模型。
假定土壤是由一束直径不相同的毛管组成，水在土壤入渗过程中，湿
润锋面几乎是水平锋面，且在锋面上各点的吸力水头均为
Sm。
锋面后面的土壤含水率为均
一的，如图（2-5-4）所示。
所以k（0）也为常数，这种模型又称活塞模型。
根据达西定律：
H+S+7
q=—krJ=m——（2-5-25）
z
式中：
H――地面以上水层厚度；
Sm――锋面处土壤负压；
z一锋面推进距离。
图2-5-4Green---AmpiA.
根据水量平衡原理，应等于土体内增加的水量q=dA-,
dt
式（2-5-26）积分:
所以
式（2-5-27）为z〜t关系式，原则上可以求得任何时刻当然也就不难求得该时刻的累计入渗量：
（2-5-28）
W\~-^0z
Ht0时，式（2-5-27）可写作:
t时入渗总量
Z=,2门厂〒0—厂Smt
（2-5-29）
I对t求导，得:
式（2-5-27）'表明，入渗初期，入渗深度z与•..t成正比。
将
（2-5-30）
当t大时，式（2-5-26）中
Z>>H+Sm,因此H—Sm一Z:
'1，则由式（2-5-25）可
知:
（2-5-31）
第四节水平入渗条件下的Philip解法
水平入渗条件下的Philip解【51］是一种半解析法，即前半部用解析法，利用博茨曼
（Boltzmann）变换，将偏微分方程转换为常微分方程；后半部采用迭代计算，求解常微分方程。
由于求解过程中未作过分简化，求得结果较为严密。
、水平入渗的常微分方程推求
水平入渗的基本方程
c9
ft
limy-补
t=0,x0,
t0,x=0,
t0,x—',
（2-5-32）
将式（2-5-32）中基本方程改写为以坐标x（0,t）为变量的方程，根据第二章中方程
（2-5-34）
（2-5-35）
（2-5-36）
（2-5-37）
（2-2-17），在水平入渗时应为
（2-5-33）
采用Boltzmann变换，引入变量从入（0）,且令
1
'v-X（V,t）t2
1
则XV,t-'Vt2
m丄V门
.:
t2
将式（2—5-36）、式（2—5—37）代入式（2-5-33）得:
经整理后得微分方程
（2-5-38）
由边界条件已知：
为了求得入〜0的关系式，将式（2—5—38）常微分方程自0i至0积分得：
八dV
'-二dv--2D*L（2-5-39）
因•二Xt3，即可由入〜B求得B（x,t），从而求得剖面上任何时间，任何距离的含水率分布。
若能实测入（B）,则代入上式可求得D（0）关系曲线。
、迭代计算法
在D（0）已知时，式（2-5—39）为入〜0关系式，它的关系曲线如图2—5—5所示,
00^0i等份成n份，步长为=（00-0i）/n,若等分点按顺序其含
0n-1,0n。
相应各等分点相应入值为入1,入2,…入n,两个相邻等分
1/2,D1+1/2，…D（n-2）
亠_入_r宀（2-5-40）
Dr+1/2取平均值，定义
式（2-5-43）表示入〜0关系曲线所包含的面积。
式（2-5-44）表示0=0r+1/2时，0r+1/2〜0n之间入〜0曲线所包含面积的近似值。
由图可知面积为
'r矩形面积后的一小
2
（2-5-45'
式中：
S「+1—-为入〜0的曲线与入r,入r+1/2之间所包含面积减去块面积。
由式（2—5一43）及式（2-5-45）得：
当△0取得足够小时，则Sr+1可以忽略，所以
—J"
2
则式（2一5一42）可以改写为
J1:
2
J1丄
2
J21
2
=J1-'1L.：
：
l
2
7丄
2
由图形可知:
A9
1
2J1
2
（2-5-46）
（2-5-49）
J1
r-
2
2"
根据式（2—5—46）和式（2—5—49），若已知J1/2及D（0）则可交替使用此二式，求得入1,入2…入r…，入n-1。
因此，为了求得入〜0关系值必须首先求得J1/2及D〜0关系曲线。
由上述可知J1/2代表了入〜0曲线所包含面积的近似值，由式（2—5—49）可知：
若将面积看作一系列的长为（0一0n）,高为d入的水平方向的矩形，则
od
J1=v一片d'（2-5-50）
2
式中：
入---和D（0）的函数。
对于初始值J1/2，我们可将D（0）考虑为在整个含水率变化范围内的平均值，由式（2
—5—50）可知：
若用加权平均值来代替上式中D值，以D'表示：
（2-5-51）
其中等号右端系数2是为使等式衡等，D'=D'（若将D'代替D（0）即可算得该系数）。
将式（2-5-51）写成有限差的形式：
1
因D'为常数，直接利用入渗问题线性化解，以，二xt^代人得:
1
（2-5-58）
（J）=2n.cD2
2/
有了Ji/2就可以由式（2一5—56）和式（2-5-49）重复交替计算求得入〜0关系曲线，但由于J/2是估算值，存在偏差，求得的入〜0曲线必然也有偏差，为了使迭代计算尽快收敛，采用两种不同方法计算Jn-1/2,不断地修正Ji/2，以使两种方法求得的Jn-1/2值足够
的接近，以便得到一个比较精确的入〜0曲线。
第一种方法采用
jFj^J1-■nJ^（2-5-59）[即式（2-5-49）]
勺n一21
第二种方法采用
.:
t
2y2
式中A（y）定义为
（2-5-64）
二2erfy
A（y）的函数表见表（2-5-1）。
对y较大时，A（y）可以渐近级数表示:
（2-5-68）
当厶0足够小时，且假定J1/2是正确的，则JFn-1/2与JSn-1/2应该是相等的，但由于计算
中采用了近似值两者误差为.訂；，若在所要求精度范围内即可，否则进行第二、第三、……,
式中：
N--式-（2—5-70）两边相等（或误差在所要求的精确度范围内）时送代次数。
为了加速收敛，Jl/2还有另一改进公式：
（2-5-71）
P=2,3,,N
由于在其他著作中采用了与上述不同符号,代公式。
这里介绍一下用其他符号表示的Philip的迭
Philip
定义Ir+1/2为
所以,
即
（2-5-72）
1
'r：
2
其他以I为参数的各式中均比J为参数的多除一个A0。
—r丄小
22
Dr计算平均扩散度。
（3）根据公式J”2=2n•△B（D'n）1/2计计算第一次送代初值（Ji/2）1,并由公式（2—5—46）计算入1,由（J1/2）i和入i即可计算（Ji+i/2），再代入公式（2-5-49）求入2,依次类推，最后求得（丿尸“人。
n
2
（4）由以上计算同时可求得入n-1和Dn-1/2，代入公式（2-5-60）计算（JS1）1。
n一
2
（5）计算误差△i=（JF1）1-（JS1）1，若△“满足要求，以上计算结束，否则要进行第
匕nP
二、三、…次的送代，直至满足要求为止。
（6）进行多次送代时，重复以上
（2）〜（5）步骤，直至迭代误差△很小，且入1,入2,…、不再随Ji/2改变时为止。
若计算结果不理想，可能是选用△B不合适，则可重新选用
△B值，重复以上步骤进行迭代计算。
11
（7）计算出入〜B关系曲线后，因-=xt2，故X='t2，对于给定的时间t，由入可计算x值，即得0〜x关系曲线。
也即求得了B在剖面上的分布。
Philip水平入渗解的理论可用于水平土柱法测定土壤水扩散度。
土壤水分在较长的（理
论上说是半无限边界）均质土柱中发生水平运动时，一维土壤水平运动方程如式（2-5-32）
所示，经Boltzmann变换可得[如前述方程（2-5—39）]:
屮ed日fde
——
1心丿
式中：
入（0）―--,t的函数。
1当水平土柱首端供水，根据观测水分运移资料，即可绘出入〜0曲线，（■C）二xt_2）,
Q
吕-（Rdr为由初始含水率0i至0x之间入〜0曲线所包含的面积，如图2-5-5所示。
d■'
（）d为在入〜0曲线上对应含水率0X点的斜率。
因此，通过水平土柱试验即可求得土壤
d，
水扩散度D值。
第五节垂直入渗条件下的Philip解法
一维垂直入渗基本方程写成z（0,t）为函数时为
D®）
&
」
对一类边界定解条件为
■---iz■?
0,t=0,
炉=日0Z=0,t>0（2-5-75）
0=3ZT°O,t：
>0
Philip级数解的形式［52］
13
=送・（日）2（2-5-76）
i二
ZRt「二t2•2二t1•3一厂t2■■-
（2-5-77）
相应边界条件为
1航=2卞一…一%=0
9时的位置z。
以下我们以待
（2-5-78）
若求得式（2-5-77）中系数ni（9）值，则可求得在一定定系数法求解ni（9）。
式（2-5-74）可改写为
将式（2-5-80）代入式（2-5-79）归并后得：
丄m：
：
丄
1二t22二t3二存八i二存=0（2-5-81）
im
式中前四项系数n1,n2,n3,n4分别为
D
d（旳
n3=0得:
（2-5-85）
n4=0得:
閉（日）丿d日
「d日徑凹丫,2屮3（日）_屮2（日）［23（0）八护1（日）丿］小2（日）d=（0）」
（2-5-86）
根据以上各式，若已知n1（0）则可逐步求得外n2（9）,n3（B）,…，即可求得方程（2—5—74）的解。
由于收敛较快，求得前四项就足够精确了。
上述n1（0）,n2（0）,n3（0）,n4（0）都较复杂，n1（0）即为水平入渗时的入，为了求解这些值并不容易，为此Philip采用了递推公式求解。
以下以求得n2为例
说明求解方法。
2-5-6关系福线
DFAC中h可以下式表示:
为了求解系数，我们可以写出与式
（2-5-39）类似的近似系数一般方程，则
I■-df
:
fdv-:
—一：
（2-5-87）
直d日
式中：
a,3—均为9的已知函数。
式（2—5-87）左边可以写作式（2—5
—88），如图2—5—6所示。
1；fd=「fdv‘2fd：
（2—5
r'~n雹r
—88）
式（2一5—88）中右端第二个积分为
由图2—5—6可知，该积分近似等于
h，其中h为四边形BCDE的平均高度。
△
2
9取得足够小时，DE可以看作为一直线，在梯形
（2-5-89）
1FX1f
hr二f「[fri—f「fr13fr
44
将式（2—5—89）、式（2—5—91）代入式（2—5—88）,经整理后得
（2-5-94）
将式（2一5—94）代入式（2—5—93）经整理后得:
式（2-5-95）为Philip的第一个递推公式。
第二个递推公式，将r改为r—l，则式（2-5-94）为
A0
fd—1f
由式（2—5—94），若r=n—I，则
（2-5-101）
当r=n时，fn并不趋于无穷大，则由式（2—5-101）和式（2-5—99）可知，一=0。
以上两递推公式如果已知初始和边界条件即可交替使用，求得f〜0曲线。
如果边界条件为0=00时，f=0,且f0=0，则可由式（2-5-95）求得11/2，由式（2—5—99）又可求得I，当r=1时，代入式（2-5-94）求得f2。
重复以上过程即可求得f〜0曲线。
当f=入，n,^,3等不同项的系数时，都可求得对应的与0关系曲线。
在某一t时，可以求得不同
0的z值，则可求得对应于不同t时的0〜z关系曲线，即求得任何时间剖面上含水率分布0（乙t）。
入渗量计算如下：
渗入土壤的水量应等于流出水量与剖面上含水量变化之和，在半无限情况下，底部含水
率将保持初始含水率，即为0n,导水率也为常量k，因此由下部流出水量
如将土柱剖分为n等份，步长为△z,取横断面为1个单位时，经t时土柱顶部含水率
为B0,对半无限土壤剖面含水量增量近似表示为
精确值可用积分表示：
上式表示0=0n以上B〜z曲线所包含的面积，这面积也可表示为
0
Wzd）（2-5-103）
式（2-5-102）与式（2-5—103）相加即为t时间内渗入剖面总水量，即
I二.zdrknt（2-5-104）
r
1234
将Philip解的表达式•2t八3t"4住•…代入式（2—5—104）,并
加以整理得，
13
2
I=At2kntBtCt2Dt（2-5-105）
也」6
式中：
A二q「drB「冷2rdr
A、B、C、D等均为常数。
式（2—5—105）中右边第一项相当于水平入渗，所以水平入渗量为
1
（2-5-106）
丨k=At2
入渗率的计算如下：
水进入土壤的速率即为z=0处的水流通量，即为
I以式（2—5—105）代入：
按式（2—5—107）求得的入渗率当t较小时误差较小，但时间长后，计算值与实验资
料偏离较大，Philip对其进行了修正，但修正式并不适用于t小时情况。
Miller和Klute用
Hagan等人的结果，用下式计算与实验值接近。
1丄
iAt2if（2-5-108）
2
式中：
A,if——均为土壤特性常数。
A值大小决定于土壤初始含水率，一般称为吸水率。
if为稳定入渗速度，相当于土壤饱
和时的水力传导度（即渗透系数）ks。
入渗初期，入渗速度很大，ks远较1/2At-1/2为小，可
忽略不计。
随着时间增大，入渗速度迅速减小，在入渗时间很久，上78,式（2—5—108）
中右端第一项趋近于0,i=if,即为稳定入渗速度。
入渗速度在时间上变化如图2-5-3所示。
A值可通过初期入渗总量Ih确定：
1
lh=At2（2-5-109）
故A=1ht2（2-5-109）
入渗率理论公式中的常数需通过试验确定。
在生产中常采用经验公式计算入渗率i和入
渗量I。
第六节土壤水入渗的经验公式
前节介绍的理论公式，既
升级会员 