第五章 降雨和灌水入渗条件下土壤水分运动2.doc

第五章降雨和灌水入渗条件下土壤水分运动第一节 水向土中入渗过程一、概述降雨和灌水入渗是田间水循环的重要环节，与潜水蒸发一样，是水资源评价和农田水分状况调控的重要依据。水渗入土壤的强度主要取决于降雨或灌水的方式和强度以及土壤渗水性能。如果土壤渗水性能较强，大于外界供水强度，则入渗强度主要决定于外界供水强度，在入渗过程中土壤表面含水率随入渗而逐渐提高，直至达到某一稳定值。如果降雨或灌水强度较大，超过了土壤渗水能力，入渗强度就决定于土壤的入渗性能，这样就会形成径流或地表积水。这两种情况可能发生在入渗过程的不同阶段，如在稳定灌溉强度（例如喷灌）下，开始时灌溉强度小于土壤入渗能力，入渗率等于灌溉强度；但经过一定时间后，土壤入渗能力减少，灌水强度大于土壤入渗能力，于是产生余水，如图251所示的降雨或灌水条件下的入渗过程。开始时入渗速率较高，以后逐渐减小。土壤的入渗能力随时间而变化，与土壤原始湿度和土壤水的吸力有关，同时也与土壤剖面上土质条件、结构等因素有关。一般来说，开始入渗阶段，土壤入渗能力较高，尤其是在入渗初期，土壤比较干燥的情况，然后随土壤水的入渗速率逐渐减小，最后接近于一常量，而达到稳定入渗阶段。在较干旱的条件下，土壤表层的水势梯度较陡。所以，入渗速率较大，但随着入渗水渗入土中，土壤中基模吸力下降。湿润层的下移使基模吸力梯度减小。在垂直入渗情况下，如供水强度较大，使土壤剖面上达到饱和，当入渗强度等于土壤饱和水力传导度时，将达到稳定入渗阶段。如供水强度较小，小于饱和土壤水力传导度时，达到稳定入渗阶段的入渗强度将等于该湿度条件下的非饱和土壤水力传导度。入渗过程中，土壤剖面上水分分布与土表入渗条件有关。根据Coleman和Bodman的研究，当均质土壤地表有积水入渗时，典型含水率分布剖面可分为四个区，即表层有一薄层为饱和带，以下是含水率变化较大的过渡带，其下是含水率分布较均匀的传导层，以下是湿润程度随深度减小的湿润层，该层湿度梯度越向下越陡，直到湿润锋。随着入渗时间延续，传导层会不断向深层发展，湿润层和湿润锋也会下移，含水率分布曲线逐渐变平缓。二、影响入渗过程的条件降雨或灌水条件下的入渗过程和初始土壤剖面上水分分布与地下水位条件有关，一般入渗问题的定解条件有以下几种情况。（一）初始条件入渗过程的初始条件一般为初始剖面含水率或负压分布已知的条件，即 （2-5-1）（二）边界条件1地表边界条件（1）通过降雨或灌水使地表湿润，但不形成积水，表土达到某一接近饱和的含水率，即（一类边界） （2-5-2）（2）降雨和喷灌强度已知，且不超过土壤入渗强度，地表不形成积水，即（二类边界） 或 （2-5-3）式中：R（t）降雨或灌水入渗强度。（3）当降雨或灌水强度大于土壤入渗强度，地表形成积水，成为压力入渗。即（一类边界） （2-5-4）式中：H（t）地表积水深度。当地表积水而没有产生径流时，地表水深为H（t）；若产生地表径流，积水深度H（t）可根据来水强度R（t）、土壤入渗强度i（t）及地表径流量Q（t）求得。2下边界条件（1）地下水埋深较小，以地下水位作边界。当地下水位变化很小或基本保持不变时，则地下水面处土壤含水率为饱和含水率（地下水面离地面距离为d），故 （2-5-5）当地下水面随时间而变化时，即地下水埋深d为时间t函数d（t），则地下水面处负压为零，即 （2-5-6）（2）地下水埋深较大的情况，在计算范围内，下边界土壤剖面含水率保持初始含水率，即 (2-5-7)在上述条件下，如初始含水率上下一致，得则 （2-5-8）式中：k(i)离地表距离d处断面通量。（3）不透水边界。下边界为流量等于零的边界，即 （2-5-9）上述表明，研究入渗时边界条件是较为复杂的，所以，计算方法也较为复杂。第二节 土壤水入渗线性化方程的近似解在垂直入渗情况下，一维土壤水分运动的基本方程可写作： （2-5-10）如降雨或灌水前的初始含水率（在土壤剖面上含水率均匀分布）为i，则初始条件为 （2-5-11）在地表有一薄水层时，表层含水率等于饱和含水率S；在地下水埋深较大时，计算时段内入渗水不会到达下边界。为此，下边界土壤含水率不变，等于初始含水率，则边界条件可以写作以下形式： （2-5-12）由于式（2-5-10）为非线性方程（因为扩散度D（）及水力传导度k（）均为待求含水率的函数），求解比较困难，为了简化计算，近似地以平均扩散度代替D（），并以代替，则式中（2-5-10）可简化为 （2-5-13）式中（2-5-13）为常系数线性方程，可以用拉普拉斯变换求解。对式（2513）采用拉普拉斯变换后可得象函数方程： （2-5-14）式（2-5-14）的通解为 （2-2-15）式（2-5-12）经拉氏变换后，得： （2-2-16） （2-2-17）根据边界条件式（2516）、式（2一517）确定常数：代入式（2-5-15），得象函数的表达式为 （2-5-18）进行逆变换后，得含水率的表达式为 （2-5-19）补余误差函数可自表122查得。式（2519）中可用下式计算： （2-5-20）若已知D与的关系式，代入式（2520）积分，即可求得。采用式（2-5-19）求得的土壤剖面上含水率分布如示意图2-5-2所示。由于地表的入渗强度，为了推求入渗强度，首先根据的象函的表达式求： （2-5-51）地表处，z=0，则 （2-5-21）在入渗初期，相当于时，则式（2-5-21）可近似写成： （2-5-21）经逆变换得： （2-5-22）入渗初期：当D（）取平均值D时 （2-5-23）入渗时间较久，即当之，相当干时ot (min)i (cm/min)Ks图2-5-3 入渗率随时间变化图，代入（2-5-21）式，则，所以则入渗时间久时，入渗强度（iks）为 （2-2-24）自式（2523）、式（2524）得入渗速度在时间上的变化过程如图（2-5-3）所示。第三节 GreenAmpt模型的入渗解GreenAmpt模型50是1911年提出的一种简化的入渗模型，它是建立在毛管理论基础上的一种入渗模型。假定土壤是由一束直径不相同的毛管组成，水在土壤入渗过程中，湿润锋面几乎是水平锋面，且在锋面上各点的吸力水头均为Sm。锋面后面的土壤含水率为均一的，如图（2-5-4）所示。所以k（）也为常数，这种模型又称活塞模型。根据达西定律： (2-5-25)式中：H地面以上水层厚度；Sm锋面处土壤负压；z一锋面推进距离。式（2-5-25）为单位时间，单位面积流入土体的水量。根据水量平衡原理，应等于土体内增加的水量，即 （2-5-26）式（2-5-26）积分：所以 （2-5-27）式（2527）为zt关系式，原则上可以求得任何时刻t时入渗锋面所达到的位置，当然也就不难求得该时刻的累计入渗量： （2-5-28）H0时，式（2-5-27）可写作： (2-2-27)或由式（2-5-26），当t很小时，该式的H+Sm+z项中z略去，所以。此时积分得 （2-2-27）t时入渗总量 （2-5-29）式（2-5-27）表明，入渗初期，入渗深度z与成正比。将I对t求导，得： （2-5-30）当t大时，式（2 -5-26）中Z H+Sm，因此，则由式（2525）可知： （2-5-31）即入渗强度近似等于土壤饱和渗透系数。第四节 水平入渗条件下的Philip解法水平入渗条件下的Philip解51是一种半解析法，即前半部用解析法，利用博茨曼（Boltzmann）变换，将偏微分方程转换为常微分方程；后半部采用迭代计算，求解常微分方程。由于求解过程中未作过分简化，求得结果较为严密。一、水平入渗的常微分方程推求水平入渗的基本方程 （2-5-32）将式（2-5-32）中基本方程改写为以坐标x（，t）为变量的方程，根据第二章中方程（2-2-17），在水平入渗时应为 （2-5-33）采用Boltzmann变换，引入变量从（），且令 （2-5-34）则 （2-5-35）故 （2-5-36） （2-5-37）将式（25-36）、式（2537）代入式（2533）得：经整理后得微分方程 （2-5-38）由边界条件已知：为了求得的关系式，将式（2538）常微分方程自i至积分得： （2-5-39）式（2-5-39）为关系式，若已知D（）关系式代入上式即可求得关系式，因，即可由求得（x，t），从而求得剖面上任何时间，任何距离的含水率分布。若能实测（），则代入上式可求得D（）关系曲线。二、迭代计算法在D（）已知时，式（2539）为关系式，它的关系曲线如图255所示，Philip的送代法是将0i等份成n份，步长为=（0-i）/n，若等分点按顺序其含水率为1，2，n-1，n。相应各等分点相应值为1，2，n，两个相邻等分点中点的含水率为1/2，1+1/2，（n-2）+1/2，其相应扩散率值为D1/2，D1+1/2，D（n-2）+1/2，任一点含水率r，可用下式计算 （2-5-40）Dr+1/2取平均值，定义 （2-5-41）写出1/2，1+1/2，（n-2）+1/2各点的式（3-5-39）： （2-5-42）令 （2-5-43）当r=0时， （2-5-44）式（2543）表示关系曲线所包含的面积。式（2-5-44）表示=r+1/2时，r+1/2n之间曲线所包含面积的近似值。由图可知面积为 （2-5-45）式中：r+1为的曲线与r，r+1/2之间所包含面积减去矩形面积后的一小块面积。由式（25一43）及式（2545）得：当取得足够小时，则r+1可以忽略，所以 （2-5-45）则式（2一5一42）可以改写为 （2-5-46）由图形可知：式中：r矩形面积与曲线与r1/2、r之间包含的面积之差。上式可写为 （2-5-47）当足够小时，同样，（r-r）可忽略。 （2-5-48）根据式（2547），我们可以得到： （2-5-49）根据式（2546）和式（2549），若已知J1/2及D（）则可交替使用此二式，求得1，2r，n-1。因此，为了求得关系值必须首先求得J1/2及D关系曲线。由上述可知J1/2代表了曲线所包含面积的近似值，由式（2549）可知：若将面积看作一系列的长为（一n），高为d的水平方向的矩形，则 （2-5-50）式中：和D（）的函数。对于初始值J1/2，我们可将D（）考虑为在整个含水率变化范围内的平均值，由式（2550）可知：若用加权平均值来代替上式中D值，以D表示： （2-5-51）式中：的加权因子。其中等号右端系数2是为使等式衡等，D=D（若将D代替D（）即可算得该系数）。将式（2-5-51）写成有限差的形式： （2-5-52）因D为常数，直接利用入渗问题线性化解，以代人得： （2-5-53）将式（2-5-53）代入式（2-5-50）得： （2-5-54）式（2-5-54）为的初始近似值第一次迭代值，因此可写作。将=/2（D）1/2代人上式，且令积分上限为0=以n代替（0n）则式（2-5-54）可写作： （2-5-55）为了求解式（2555）积分式，利用erfc=1erf，则已知 得 (2-5-56)当时,，所以式（2556）为 (2-5-57)将式（2-5-57）代入式（2-5-55）得： (2-5-58)有了J1/2就可以由式（2一556）和式（2549）重复交替计算求得关系曲线，但由于J1/2是估算值，存在偏差，求得的曲线必然也有偏差，为了使迭代计算尽快收敛，采用两种不同方法计算J n1/2，不断地修正J1/2，以使两种方法求得的J n1/2值足够的接近，以便得到一个比较精确的曲线。第一种方法采用 （2559）即式（2549）第二种方法采用 （2559）即式（2549）JS式中，由于变化于n，n-1之间，若取得足够小，则D可以看作为一常量Dn1/2： (2-5-61)式（2-5-61）表明，在n和n1范围内将D视作常量，当取得足够小时，不会引起明显误差。上述问题可为求解半无限长水平土柱的入渗问题，即 (2-5-62)根据土壤水入渗的线性化方程的近似解法，式（2562）定解问题的解为 (2-5-63)式中A（y）定义为 (2-5-64)其中 (2-5-65)A(y)的函数表见表（2-5-1）。对y较大时，A（y）可以渐近级数表示： （2-5-66）n-1由式（2559）求得，对某一个曲线n-1及Dn-1/2为常量，所以，A（y）也为常数，可采用式（2564）计算，根据y求A（y）。首先，已知 （2-5-67）式中的为虚拟变量，可以改写为 （2-5-68）由式（2-5-68）求得erfc（y），代入式（2-5-64）即可求得A（y）、A（y）y已列于表（2-5-1）中，所以根据某Dn1/2，n1求得y，查表251求得。将式（2563）代人式（2560）得： （2-5-69）当足够小时，且假定J1/2是正确的，则 JFn-1/2与JSn-1/2应该是相等的，但由于计算中采用了近似值两者误差为；，若在所要求精度范围内即可，否则进行第二、第三、，更多次的选代。由与（J1/2）1；可以得到一个更接近于实际的初始值（J1/2）2。重复上述过程进行第二个循环的计算。（J1/2）2一般可写作： 迭代P次后： （2-5-71） 式中：N式（2570）两边相等（或误差在所要求的精确度范围内）时送代次数。为了加速收敛，J1/2还有另一改进公式： （2-5-71）由于在其他著作中采用了与上述不同符号，这里介绍一下用其他符号表示的Philip的迭代公式。Philip定义I r+1/2为 （2-5-72）所以，其他以I为参数的各式中均比J为参数的多除一个。即 （2-5-73）迭代计算步骤：（1）确定0至n的分段数n，步长=（0-n）。（2）以计算平均扩散度。（3）根据公式J1/2=2n（D/）1/2计计算第一次送代初值（J1/2）1，并由公式（2546）计算1，由（J1/2）1和1即可计算（J1+1/2），再代入公式（2-5-49）求2，依次类推，最后求得。（4）由以上计算同时可求得n-1和Dn-1/2，代入公式（2-5-60）计算。（5）计算误差1=，若1满足要求，以上计算结束，否则要进行第二、三、次的送代，直至满足要求为止。（6）进行多次送代时，重复以上（2）（5）步骤，直至迭代误差很小，且1，2，、不再随J1/2改变时为止。若计算结果不理想，可能是选用不合适，则可重新选用值，重复以上步骤进行迭代计算。（7）计算出关系曲线后，因，故，对于给定的时间t，由可计算x值，即得x关系曲线。也即求得了在剖面上的分布。Philip水平入渗解的理论可用于水平土柱法测定土壤水扩散度。土壤水分在较长的（理论上说是半无限边界）均质土柱中发生水平运动时，一维土壤水平运动方程如式（2532）所示，经Boltzmann变换可得如前述方程（2539）： （2-5-39）式中：（）x，t 的函数。当水平土柱首端供水，根据观测水分运移资料，即可绘出曲线，（），为由初始含水率i至x之间曲线所包含的面积，如图2-5-5所示。为在曲线上对应含水率x点的斜率。因此，通过水平土柱试验即可求得土壤水扩散度D值。第五节 垂直入渗条件下的Philip解法一维垂直入渗基本方程写成z（，t）为函数时为 （2-5-74）对一类边界定解条件为 （2-5-75）Philip级数解的形式52 （2-5-76）相应边界条件为 （2-5-77） （2-5-78）若求得式（2-5-77）中系数i（）值，则可求得在一定时的位置z。以下我们以待定系数法求解i（）。式（2574）可改写为 （2-5-79）将z的解式（2-5-76）分别代入以上各项： （2-5-80）将式（2-5-80）代入式（2-5-79）归并后得： （2-5-81）式中前四项系数1，2，3，4分别为 （2-5-82）由于式（258l）左端等于零，所以1=2=3=4=01=0，即全式除以得即 式（2583）表明，1（）即为（）为水平入渗时解。同理，以2=0求得： （2-5-84）3=0得： （2-5-85）4=0得： （2-5-86）根据以上各式，若已知1（）则可逐步求得外2（），3（），即可求得方程（2574）的解。由于收敛较快，求得前四项就足够精确了。上述1（），2（），3（），4（）都较复杂，1（）即为水平入渗时的，为了求解这些值并不容易，为此Philip 采用了递推公式求解。以下以求得2为例说明求解方法。为了求解系数，我们可以写出与式（2-5-39）类似的近似系数一般方程，则 （2-5-87）式中：，均为的已知函数。式（2587）左边可以写作式（2588），如图256所示。 （2588）式（2一588）中右端第二个积分为 由图 256可知，该积分近似等于，其中h为四边形BCDE的平均高度。取得足够小时，DE可以看作为一直线，在梯形DFAC中h可以下式表示： （2-5-89） （2-5-90）式（2-5-87）中，若=r+1/2，则方程为 （2-5-91）r+1/2为r+1r范围内平均值，且以来近似（df/d）r+1/2（其负号表示值增加的方向与f值增加方向相反），代入式（2591）得： （2-5-92）将式（2589）、式（2591）代入式（2588），经整理后得 令 （2-5-94）将式（2一594）代入式（2593）经整理后得： （2-5-95）式（2595）为Philip的第一个递推公式。第二个递推公式，将r改为rl，则式（2594）为 （2-5-96）将式（2594t）4）与式（2596）两式相减，得： （2-5-97）根据式（2590）得： 将式（2598）代入式（2597）经整理得：当足够小时，fr，fr+1/2，fr-1/2，fr-1非常接近，则上式可写作：或 （2-5-99）式（2599）即为PhiliP的第二个速推公式。根据该公式，令 （2-5-100）由式（2594），若r=nl，则 （2-5-101）当r=n时，fn并不趋于无穷大，则由式（25101）和式（2599）可知，In=0。以上两递推公式如果已知初始和边界条件即可交替使用，求得f曲线。如果边界条件为=0时，f=0，且f0=0，则可由式（2595）求得I1/2，由式（2599）又可求得I，当 r= 1时，代入式（2594）求得f2。重复以上过程即可求得f曲线。当f=，等不同项的系数时，都可求得对应的与关系曲线。在某一t时，可以求得不同的z值，则可求得对应于不同t时的z关系曲线，即求得任何时间剖面上含水率分布（z，t）。入渗量计算如下：渗入土壤的水量应等于流出水量与剖面上含水量变化之和，在半无限情况下，底部含水率将保持初始含水率，即为n，导水率也为常量k，因此由下部流出水量 （2-5-102）如将土柱剖分为n等份，步长为z，取横断面为1个单位时，经t时土柱顶部含水率为0，对半无限土壤剖面含水量增量近似表示为精确值可用积分表示：上式表示=n以上z曲线所包含的面积，这面积也可表示为 （2-5-103）式（25102）与式（25103）相加即为t时间内渗入剖面总水量，即 （2-5-104）将Philip解的表达式代入式（25104），并加以整理得， （2-5-105）式中： A、B、C、D等均为常数。式（25105）中右边第一项相当于水平入渗，所以水平入渗量为 （2-5-106）入渗率的计算如下：水进入土壤的速率即为z=0处的水流通量，即为I以式（25105）代入： （2-5-107）按式（25107）求得的入渗率当t较小时误差较小，但时间长后，计算值与实验资料偏离较大，Philip对其进行了修正，但修正式并不适用于t小时情况。Miller和Klute用Hagan等人的结果，用下式计算与实验值接近。 （2-5-108）式中：A，if均为土壤特性常数。A值大小决定于土壤初始含水率，一般称为吸水率。if为稳定入渗速度，相当于土壤饱和时的水力传导度（即渗透系数）ks。入渗初期，入渗速度很大，ks远较1/2At-1/2为小，可忽略不计。随着时间增大，入渗速度迅速减小，在入渗时间很久，t，式（25108）中右端第一项趋近于0，i=if