1 模拟方法

1.1　基于Green-Ampt修正模型和王文焰等提出模型的改进模型

1.1.1　积水入渗情况

对于积水入渗情况，Green-Ampt模型采用的水分剖面如图1所示，其将实际入渗深度为Z_f的浸润体简化为长度为Z_s的饱和体，水力梯度i采用式（1）计算，根据达西定律得到入渗量增量，该值等于饱和体长度增量下土体水量的增量如式（2）所示，由式（1）和式（2）得到饱和体深度Z_s和时间的关系如式（3）所示。

式中： H为积水水头，Z_s为Green-Ampt模型的浸润锋深度（也即饱和体深度），S_m为Green-Ampt模型浸润锋Z_s处的吸力水头（取为正值），K_s为饱和渗透系数，θ_i为土体初始体积含水率，θ_s为土体饱和体积含水率。

图1 积水入渗下Green-Ampt模型、王文焰等提出的以及实际的水分剖面示意图

Fig.1 The water content profiles adopted by Green-Ampt model, proposed by Wang W Y et al and the actual water content profile during ponding infiltration

该模型浸润锋处吸力水头S_m并非土体初始含水率θ_i对应的吸力水头，Philip^[14]认为S_m是一个数学量，没有实际的物理意义，该值远小于浸润锋下土体初始含水率对应的吸力水头S_i。Bourwer^[15]最早给出了S_m的计算公式如式（4）所示，并认为S_m可近似取为增湿过程中进水值即增湿过程中气不连续时对应的吸力值，然而该值不容易测定，其建议取为饱和样（对应于抽真空饱和样）脱湿过程中进气值的50%，该确定方法比较简单可行^[16]。在浸润锋以上为饱和的假设下，Mein和Larson^[11]应用式（5）获得S_m，该值在初始含水率变化较大范围内相差很小，其对不同初始含水率和不同降雨强度下采用相同的值。Mein和Farrell^[17]经过推导认为应用式（4）获得S_m是合理可靠的。Neuman^[18]通过推导得到了浸润锋处吸力值表达式与式（4）相同。Brakensiek^[19]通过比较，表明上述各种方法可得相近的结果，对于采用上述等效吸力的模型称为Green-Ampt修正模型。已有大量研究表明， Green-Ampt修正模型在模拟入渗量和时间的关系上具有较高的准确度^[9,24]，说明了上述确定S_m的方法是有足够精度的。从公式（4）可见，S_m小于初始含水率对应的吸力水头S_i。

式中： S_i和K_i分别为初始含水率对应的吸力水头和渗透系数，S和 K分别为吸力水头和渗透系数。

Green-Ampt模型及修正模型仅能计算如图1所示等效饱和体深度，王文焰等^[10]根据大量黄土区实测水分剖面，假定入渗深度为Z_f的区域近似由上部长度为Z_f/2的饱和区（严格来讲是传导区，因为饱和区和过渡区的区间都很小^[20]，其中传导区和过渡区为接近饱和的区域，所以这里仍用饱和区称之）和下部长度为Z_f/2、水分剖面为椭圆的非饱和浸润区组成如图1所示，提出了能计算不同时刻水分剖面的模型，然而在计算浸润锋处等效吸力S_m上却没有采用上述确定方法，其采用下部Z_f/2非饱和浸润区平均初始含水率对应的吸力S_i施加于上部饱和区底部，得到上部Z_f/2饱和区的水力梯度如式（6）所示，根据达西定律和水量平衡原理，得到式（7），将式（6）代入式（7）整理可得入渗深度Z_f和时间t的关系如式（8）所示；如上文所述作用于长度Z_s饱和体的等效吸力S_m小于S_i，而将S_i作用于长度为Z_f/2的饱和区，可见过大估计了下部吸力的作用。该模型^[10]计算所得现场积水入渗试验中的水分剖面运移过程较为符合实测，可能在于原状黄土常存在节理和大孔隙，节理^[21]或大孔隙^[22,23]存在导渗的作用，加快了浸润锋迁移速率，对此问题应采用相应模型^[21,22,23]进行分析，而不应单独在模型中增加吸力作用。

彭振阳等^[24]借助王文焰等^[10]提出的分层假定，采用求解Richard方程的数值方法获得了饱和区和非饱和区的比例变化；在计算浸润锋锋面吸力上，其建立了非饱和浸润区等效渗透系数K¯和整个入渗区等效渗透系数K¯w的计算公式，如式（9a）和式（9b）所示，采用式（10）计算浸润锋Z_f处吸力S′_m。张杰等^[25]沿用式（10）并给出了更严格的Z_f-t的关系表达式，并在此基础上提出了考虑气阻效应下的改进模型。式（10）的理论基础在文献[24]中并没有明确给出，在不考虑气阻效应情况下，本文认为应该用式（11）计算S′_m，因为Green-Ampt修正模型在模拟入渗率方面具有很高的精度^[9,24]，式（11）能保证与Green-Ampt修正模型具有相同的入渗率。为了求解方便，本文并不用式（11）求解S′_m，而是应用该式的右端项来计算入渗量的增量也即Green-Ampt修正模型来求出等效饱和体长度Z_s，然后根据等效饱和体水量如图1所示阴影面积与饱和区及浸润区水量相等，来得到实际浸润锋深度Z_f和饱和体深度Z_s的关系如式（12a）所示。当应用王文焰等^[10]所采用饱和区和非饱和浸润区长度相等的假定，也即ε=0.5，此时Z_f和Z_s的关系如式（12b）所示。下文采用式（12b）来求解实际浸润锋位置并结合浸润区为椭圆形求解水分剖面，并称之为本文改进模型。

式中：ε为浸润区占入渗区的比例，其随Z_f的增加而线性增加，可表示为ε=aZ_f+b。

式中：Z_s为如图1所示的饱和体长度，按浸润区为椭圆过渡，经过换算Z_s等于(1-ε+0.25πε)Z_f。

1.1.2　降雨入渗情况

对于不同降雨强度下的水分入渗模拟，Mein和Larson应用Green-Ampt修正模型进行了降雨入渗分析 ^[11]，得到与求解Richard方程的数值计算方法相近的入渗率和时间关系，但其并没有引入水分剖面形状来预测实际水分剖面。本文应用上述改进模型，获得实际水分剖面。当降雨强度小于饱和渗透系数时，不会发生径流，入渗率f(t)一直等于降雨强度P。当降雨强度大于饱和渗透系数时，在地表发生径流前，入渗率等于降雨强度；发生径流时，降雨强度P等于入渗强度K_si，据此可得径流发生时等效饱和体长度Z_s，进而根据式（12b）得到此时的实际入渗深度Z_f；根据此时入渗量F_p等于等效饱和体范围内土体水量变化，可得F_p如式（13）所示，进一步可得径流发生的时刻t_p如式（14）所示。径流发生后，假定水及时排走没有积水水头，即发生积水水头为0的积水入渗，入渗量F和时间的关系可按式（15）计算，入渗深度Z_f由式（16）和式（12b）计算得到。

1.2　基于求解Richard方程的数值方法

Hydrus-1d软件应用迦辽金有限元法求解给定初始条件和边界条件下的一维Richard方程，获得非饱和土水分迁移过程。一维Richard方程如式（17）所示。对于本文积水和降雨入渗中初始含水率沿深度均相等情况，初始条件如式（18）所示，下边界条件为自由排水边界如式（19）所示。对于积水入渗情况，上边界条件如式（20a）所示；对于降雨入渗情况，径流发生前，上边界条件如式（20b）所示，当计算到上边界水头为0的时刻，即为径流时刻，从该时刻开始，将上边界条件定义为水头为0的边界，即应用式（20a）所示边界条件并将H设置为0；另外还分析了停雨后1小时的水分再分布，该过程将上边界设置为流量为0的边界即将式（20b）所示P设置为0。求解时需要h和θ的关系，以及K与h或θ的关系；然而有时不能采用间接法如VG模型获得的非饱和渗透系数参数而需要反演^[4,5]，本文采用改进的浸润锋前进法^[13]进行求取，具体参数确定过程见下节。

式中：z坐标取向上为正，地表处z=0；此处h为吸力水头，当孔隙水压力为负时取为负值，H代表积水水头，P为降雨强度。

2 两类模拟方法的验证与分析

2.1　积水入渗情况水分运移计算的验证与结果分析

2.1.1　试验情况及参数确定

为检验上述两类模拟方法在模拟重塑黄土积水情况下水分迁移的可靠性，以延安治沟造地工程建设中开挖边坡土料制取的两种干密度重塑黄土土柱积水入渗试验^[27]为对象进行模拟，两种干密度分别为1.35g/cm³和1.53g/cm³，土柱初始质量含水率均约为12.5%，初始体积含水率分别为0.174和0.194，积水水头为2cm，土柱高220cm，在不同高度处布置水分仪获得了积水过程中的含水率变化，土柱底部为透气板，且入渗速率相对松散土壤或砂土较慢，浸润锋下的气压阻渗作用可以不考虑。

为了能用上述两类方法模拟试验中水分运移过程，本文测试了饱和渗透系数和持水曲线，获得了各分析方法所需参数如表1所示，其中间接法1采用室内制取相同初始状态的土样进行增湿或脱湿测试的持水曲线如图2，应用VG模型如式（21a）拟合；对制成相同初始状态的试样进行浸水饱和，应用变水头法测得的饱和渗透系数K_s，应用VG-Mualem模型即式（21b）间接获得非饱和渗透系数函数，进而获得S_m。间接法2除饱和渗透系数应用改进的浸润锋前进法根据积水入渗土柱试验获得的饱和渗透系数外，其他参数与间接法1相同。直接法则采用改进的浸润锋前进法^[13]根据积水入渗土柱试验获得的非饱和渗透系数函数。

图2 持水曲线

Fig.2 Soil water retention curve

（21a）

式中： θ_w为体积含水率，ψ为基质吸力，单位为kPa；a、n、θ_r、θ_s为拟合参数; a单位为kPa；n为曲线斜率有关的参数，m=1-1/n；θ_s选取为土柱积水入渗后试样的饱和含水率；Θ=θw-θrθs-θr=1+(ψa)n-m。

其中根据改进的浸润锋前进法获得非饱和渗透参数的过程如下，根据土柱上各测点含水率开始变化的时间即浸润锋达到该处的时间，按幂函数关系拟合得到浸润锋迁移距离和时间的关系，对于两种干密度样，所得结果分别如式（22a）和（22b）所示，相关系数均在0.997以上。应用该关系很容易得到不同时刻浸润锋迁移速率vZf，以其中一个测点A为例，按式（23a）计算得到t₁~t₂时间段通过A断面的水流速度v，浸润锋前进法^[12]采用式（23b）计算相应水力梯度i，依据改进的浸润锋前进法^[13]相应采用式（23c）计算t₁~t₂时间段A断面的水力坡降i，根据达西定律便可得到该时间段平均吸力下的渗透系数，选取不同的时间段便可得到不同吸力下的渗透系数。图3给出了应用改进的浸润锋前进法所得不同吸力下的渗透系数结果，这里选取了水分变化过程比较合理的测点。从结果上来看虽然有一定的离散，但每个土柱所用的两个测点获得的不同吸力下的渗透系数都比较接近。应用VG模型即式（20b），对所得非饱和渗透系数和基质吸力的关系进行拟合，拟合得到的饱和渗透系数在抽真空饱和样和浸水饱和样所测渗透系数之间，可见虽然浸水饱和样和土柱入渗试验后的饱和度接近，在渗透系数上仍存在试样不对等性，文献^[28]则建议采用入渗试验后的土柱进行渗透试验获得饱和渗透系数。由于应用式（20b）拟合非饱和渗透系数时，所得a和n与原持水曲线拟合参数不同，因此需要在新的a和n下，重新对持水曲线拟合，得到残余含水率θ_r分别为0.114和0.138，虽然该参数有较大改变，但前后两次所得持水曲线在大于初始含水率时都很接近测试点，而入渗过程是增湿过程，所得持水曲线和非饱和渗透系数函数在含水率大于初始含水率时具有足够的精度。

式中：θzA,t2代表t₂时刻，深度z_A处测点A的体积含水率；ψ表示基质吸力(kPa)，由前面所测持水曲线拟合函数根据含水率反算得到；θ_i表示初始体积含水率；vZf表示t₁~t₂时间段浸润锋迁移速度。

（a） 干密度1.35g/cm³土柱

（b） 干密度1.53g/cm³土柱

图3 改进的浸润锋前进法所得非饱和渗透系数

Fig.3 Unsaturated hydraulic conductivity obtained by improved wetting front advancing method

2.1.2　模拟结果分析

图4给出了各方法所得浸润锋入渗深度和时间的关系以及与实测结果的对比。总体上，王文焰等提出的模型^[10]过快地预测了入渗过程，Green-Ampt修正模型较慢地预测了入渗过程，而采用相同参数情况下改进模型较Green-Ampt修正模型提高了浸润锋迁移速率，相比前两种方法均更加接近实测值。针对采用不同方法确定的参数，采用间接法1所得参数误差较大，而采用间接法2或直接法所得参数预测结果均较为接近实测值，说明饱和渗透系数的准确性在很大程度决定了预测的准确度。另外整体上，对干密度1.53 g/cm³土柱预测较为准确，这可能与1.35g/cm³土柱具有导渗的大孔隙有关。综合以上分析，下文仅给出改进模型以及Hydrus采用直接法所得参数的分析结果。

（a） 干密度1.35g/cm³土柱

（b） 干密度1.53g/cm³土柱

图4 积水入渗过程中浸润锋入渗深度与时间的关系

Fig.4 The relationship between the depth of wetting front and time during ponding infiltration

图5给出了离土柱顶面不同深度测点处体积含水率随时间的变化，其中由于深200cm处测点水分

（a） 干密度1.35g/cm³土柱

（b） 干密度1.53g/cm³土柱

图5 积水入渗过程中各测点含水率随时间变化预测与

Fig.5 Comparison of predicted and measured water content change at measuring points during ponding infiltration

注：实测值的对比

变化过程中，浸润锋到达了底部，因此这里便没给出该情况的模拟结果。从结果上看，改进模型和Hydrus软件都能模拟出浅部测点体积含水率变化快，深部测点体积含水率变化稍慢的特点。总体上，对于干密度1.53g/cm³土柱预测较好，而对干密度1.35g/cm³土柱预测稍差。

图6给出了浸润锋到达各测点时的水分分布预测值和实测结果的对比。从结果上来看，改进模型和Hydrus软件所得结果很接近，只是在浸润锋位置处稍有差异。有限的实测点结果接近于两种方法所得的水分分布线，说明两种方法的可靠性。另外根据Hydrus软件所得结果，分析水分剖面从饱和到非饱和的过渡点，得到饱和区（严格来讲为传导区）所占比例随着入渗深度的增加由0.3变化到0.6。

（a） 干密度1.35g/cm³土柱

（b） 干密度1.53g/cm³土柱

图6 积水入渗过程浸润锋到达各测点时的水分分布预测值和实测结果对比

Fig.6 Comparison of the predicted and measured water content profile when wetting front arrived at measurement points during ponding infiltration

2.2　降雨入渗情况的验证与分析

为了获得积水入渗和降雨入渗的差别以及改进模型相对于Green-Ampt修正模型在降雨入渗过程分析的适宜性，对比了Green-Ampt修正模型、本文改进模型和Hydrus降雨情况下浸润锋入渗规律和水分剖面并与积水情况进行了对比。

根据近年来陕北地区7月份经常遭受暴雨，24小时降雨量可达200mm以上，结合取土地区的降雨量资料^[29]，这里给出特大暴雨强度1.17cm/h，大雨强度0.208cm/h和中雨强度0.104cm/h下降雨1天和停雨1天的分析结果。对于Green-Ampt修正模型和本文改进模型，按式（14）计算得到干密度1.35g/cm³和干密度1.53g/cm³黄土地层仅特大暴雨情况下分别在3.4h和0.73h发生径流，其他两种雨强在1天内均未发生径流；Hydrus分析结果表明1天内也仅特大暴雨下发生了径流，径流时间分别为4.5h和0.8h，这与Green-Ampt修正模型和本文改进模型所得结果相近。

图7给出了积水和三种雨强条件下浸润锋入渗深度和时间的关系。对于Green-Ampt修正模型和本文改进模型，径流之前，由于地表处于非饱和状态不能应用改进模型预测，图中仅给出了径流发生后的入渗深度和时间的关系，从结果上可见本文改进模型比Green-Ampt修正模型更加接近Hydrus所得结果。对比三种雨强和积水入渗结果，特大暴雨情况比较接近于积水入渗情况，特别是干密度为1.53g/cm³地层；这也说明了积水入渗可以作为降雨入渗的极端条件予以研究。

（a） 干密度1.35g/cm³地层

（b） 干密度1.53g/cm³地层

图7 降雨24h过程中入渗深度和时间的关系及与积水入渗的对比

Fig.7 The relationship between the infiltration depth and time during 24h rainfall compared with that under ponding condition

图8给出了特大暴雨情况下，在降雨24小时和停雨24小时时间段，两种干密度地层水分入渗及水分再分布过程水分剖面图。Hydrus分析结果表明初始入渗时地表由非饱和状态逐渐过渡到饱和状态，Green-Ampt修正模型和本文改进模型不能预测径流发生前的水分分布变化；发生径流时及降雨24h时，本文改进模型所得水分分布与Hydrus所得水分分布接近，特别是干密度1.53g/cm³地层，而Green-Ampt修正模型预测水分剖面精度较差。相对于降雨24h，发生径流时Green-Ampt修正模型和本

（a） 干密度1.35g/cm³地层

（b） 干密度1.53 g/cm³地层

图8 降雨一天及停雨一天过程中的水分剖面变化

Fig.8 The change of water content profile during one day of rain and following one day without rain

文改进模型所得结果与Hydrus分析结果差异较大，

这主要是径流发生的时间不同导致两者入渗量不同，另外如前面分析表明入渗深度浅时，饱和区所占比例小，改进模型在此阶段采用了相对较大的饱和区比例。降雨结束后，改进模型不能预测水分重分布，Hydrus得到地表含水率减少并且更深处水分增加，这与已有实测结果^[6,7]相符；需要说明的是，该分析过程宜用增湿后脱湿的持水曲线函数，另外停雨后一般涉及到蒸发，需要测试或计算蒸发强度E，应用将式（20b）中P替换为-E的上边界条件进行计算。

3 结论

在Green-Ampt修正模型和王文焰等提出的模型^[10]基础上，提出了可较为合理预测入渗过程中水分剖面及迁移过程的改进模型，应用各模型和求解Richard方程的数值方法模拟了室内重塑黄土积水入渗试验，并对比了特定雨强下改进模型和数值方法所得的水分剖面。主要结论如下：

（1）改进模型合理考虑了浸润锋处吸力作用和水分剖面形状，相比Green-Ampt修正模型和王文焰等提出的模型，能更好地模拟室内非饱和重塑黄土积水入渗试验中水分迁移过程。

（2）根据土柱入渗试验各测点水分变化，采用改进的浸润锋前进法^[13]获得的饱和渗透系数或非饱和渗透系数函数时，能较好地预测水分迁移过程，而采用室内制样获得的饱和渗透系数和间接法获得的非饱和渗透系数函数时预测误差较大。

（3）根据Mein和Larson提出的降雨入渗理论，将改进模型推广到降雨情况下径流发生后的水分剖面预测，所得结果与求解Richard方程的数值方法分析结果接近，特别是在试样干密度比较大的情况，改进模型相对数值方法计算简便、易于应用。另外应用Hydrus软件可以分析得到径流发生前和停雨后的水分迁移变化过程，结果与已有试验规律相符，显示出较强的模拟能力。

参考文献