20世纪中叶，Mein和Larson两位学者针对传统G-A模型未考虑边坡倾角和弱降雨—自由入渗的不足进行推导，提出了Mein-Larson入渗模型（以下简称M-L入渗模型）。它主要建立在以下几个假定的前提下：①边坡为均质土体；②雨水入渗过程中，土体分为湿润层和初始层两层，两层的含水率均匀分布，初始层含水率为初始体积含水率；③水分的传导方向为垂直坡面向下即湿润锋面垂直坡面向下平行推进；④湿润锋面处基质吸力不变，为常量。

Mein-Larson入渗模型如图1所示，先分析弱降雨—自由入渗阶段即降雨强度小于土体入渗率阶段，此时湿润锋以上的含水率为湿润含水率θ_w，湿润锋以下的含水率为土体的初始含水率θ_i，坡面倾角为β，坡面上z^*方向的入渗率为i，以坡面作为参考面，以沿坡面垂直向下作为z^*轴正方向，θ轴为土体的含水率轴，根据达西定律可得式：

式（1）中：h_w为压力水头；K ( )h_w 为与压力水头h_w相关的土体渗透系数。

图1 Mein-Larson入渗模型（自由入渗阶段） Fig.1 Mein-Larson infiltration model（free infiltration stage）

根据坐标的转换关系，可得式：

z=x^∗sinβ + z^∗cosβ（2）

联立式和式，整理可得：

根据M-L降雨入渗模型的基本假定：湿润区的体积含水率均匀分布，故上式第2项可以当成0。此时，式化为式。

i=K(h_w ) cosβ（4）

自由入渗阶段，水分全部渗入土体，坡面不产生径流，沿垂直坡面向下方向的入渗率等于降雨强度在该方向上的分量，即：

i=Rcosβ（5）

在非饱和土的入渗分析中，表征土体含水率及渗透系数分别与基质吸力关系的土水特征曲线（SWCC）和水力传导方程（HFC）是两个十分重要的水力特性。VG模型对于上述两个水力特性的描述分别如式（6）与式（7）所示：

式（6）中：θ：体积含水率（%）；θ_r：残余体积含水率（%）；θ_s：饱和体积含水率（%）；α、m：均为土水特征曲线拟合参数；n：孔隙分布系数；h：土体负压水头（m）。

式（7）中：K：压力水头为h时，土体的渗透系数（m/s）；K_s为饱和渗透系数（m/s）。

将式（4）、式（5）及式（7）联立，可求得：

通过式（8）可以求得压力水头的大小，将其代入式（6），可以求出边坡土体的体积含水率表达式如（9）：

根据降雨量守恒原理，可得：

R·t cos β=( θ _w - θ _i ) _f cos β（10）

将式（10）进行整理，可得降雨入渗深度随时间的变化规律：

在有压入渗阶段即土体入渗能力控制阶段，降雨入渗模型如图2所示，其中湿润锋以上的部分土体体积含水率为饱和含水率θ_s，坡面上形成了径流h₀，根据达西定律，此时的土体入渗率如式（12）所示：

根据降雨量守恒原理，可得：

I=(θ_s - θ_i) Z_f（13）

入渗雨量对时间进行求导可得土体的入渗率，可表示为：

联立上述式（12）和式（14），并求积分，代入初始条件 |Z_{f t=0} =0，可以求得：

图2 Mein-Larson入渗模型（有压入渗阶段） Fig.2 Mein-Larson infiltration model （with pressure infiltration stage）

将式（13）代入式（15）可得入渗雨量随时间的变化关系：

在入渗时间很长的情况下，坡面径流的深度h_o可以略去不计，此时土体入渗率可以表示为：

将式（13）代入上式（17）可得：

式（18）中I为降雨入渗量，S_f为基质吸力水头，按照下式（19）进行计算：

降雨时刻是从非积水入渗开始的，而上述的推导是以积水时刻作为起始时刻的，中间还有t_p时间段，且降雨入渗量I实际上是累计入渗量，故必须对上式进行相应的修正：

上述式中，t_p代表积水时刻。由下式给出：

传统M-L入渗模型假设土体初始含水率均匀分布，这显然与自然土体中的含水率分布相差甚远，难以准确描述土体中水分的分布情况，且为了简化研究，将湿润锋的发展规律假设为直线均匀推进，与实际情况不符。本文针对M-L入渗模型存在的上述不足，基于分层假定引入一种新的指数型分布函数表征土体的含水率对其进行改进。1.2.1 土体含水率分布函数修正

假定土体初始含水率分布函数为指数型分布，湿润层内土体的含水率为湿润含水率θ_w，初始层内的土体含水率为θ_i，改进模型如图3所示，两层区域的土体含水率分布函数如式（22）和式（23）：

式（23）中，Z_f为湿润锋的入渗深度（m）；Z^*为土体深度（m）；θ_s、θ_w分别为土体饱和含水率与湿润层土体含水率；d为地下水埋深（m）；θ₀为土体含水率拟合参数，物理意义为表层土体含水率；a为土体含水率拟合参数。

图3 改进M-L降雨入渗模型Fig.3 Improved M-L rainfall infiltration model

1.2.2 土体累计入渗量修正

上述的假定给出了两层土体含水率的分布情况，根据式式（22）和式式（23）可以算出土体累计入渗量I随着湿润锋入渗深度Z_f变化的关系，累计入渗量I由式（24）式给出：

将式（22）与式（23）代入上式（24），进行整理化简，可得累计入渗总量I的表达式如式（25）：

本文提出的改进降雨入渗模型分为弱降雨和强降雨两个情况分别对应自由入渗和有压入渗，基于上述的假定来推导湿润锋入渗深度Z_f与时间t的关系，先分析弱降雨—自由入渗阶段，改进降雨入渗模型如上图3所示，由于此阶段降雨强度恒小于土体入渗率，此时坡体法向的入渗率即为降雨强度在此方向的分量，表达式见式（26）。

i_(t)=Rcosβ（26）

式（26）中入渗率i_{( t )}又可表示为累计入渗总量I对时间t的导数，即：

将式（25）代入式（27），整理化简可得：

对式（28）进行积分运算，并带入初始条件|

Z_{f t=0} =0，积分可得：

式（29）即为改进之后，自由入渗阶段湿润锋入渗深度随时间增长的变化关系式。

再来推导强降雨—有压入渗阶段，湿润锋随时间的变化关系式。此阶段降雨强度恒大于土体入渗率，由土体入渗能力控制水分入渗，由达西定律可以得知：

考虑边坡入渗时间很长即形成稳定的坡面径流时的情况或坡面径流高度较小时，此时可以将h_o忽略不计。

同理，将累计入渗总量I对时间t求导，即为土体入渗率i_{( )}t ，联立式（27）及式（30），化简整理可得：

式（31）中C为积分常数，代入初始条件 |Zf _t=0 =0便可求得。式（31）即为改进之后，有压入渗阶段湿润锋入渗深度随时间增长的变化关系式。

当降雨强度略大于土体的饱和入渗率时，降雨前期会发生自由入渗，经过积水时刻后，发生有压入渗，此时的水分入渗经过自由入渗到有压入渗的转化。现进行此入渗过程中湿润锋入渗深度随时间的变化关系的推导。

两个阶段的转化时刻即为坡面开始积水时刻，设为t_p，累计入渗雨量设为I_p，湿润锋入渗深度设为Z_p，此时，土体入渗率在数值上就等于降雨强度在垂直坡面方向的分量。

联立式（2 6）及式（3 0），可求解得到Z_p如下所示：

坡面开始积水时刻的累计入渗雨量I_p表达式如式（33）：

将式（22）及式（23）代入式（33），化简整理可得：

同时，积水时刻t_p可以表达为式（35）：

将积水时刻的累计入渗雨量I_p代入式（35），可得：

当t<t（p 即自由入渗阶段）时，湿润锋入渗深度随时间的变化关系表达式同式（29）。

当t≥t（p 即有压入渗阶段）时，湿润锋入渗深度随时间的变化关系表达式需要对式（31）进行修正，式（31）是以 t=0 作为入渗的开始时刻点，而实际上的有压入渗阶段发生在积水点之后也即 t=t _p之后的阶段，并考虑将初始条件 |Z _f t=t _p =Z _p，则修正之后的表达式（37）如下式所示，其中累计入渗总量见式（25）：

式（29）及式（37）即为降雨强度略大于土体饱和入渗率时，全降雨阶段的湿润锋入渗深度随时间的变化关系表达式。

本文选取华南地区花岗岩残积土边坡来进行改进M-L模型的试验验证。该试验边坡坡高10 m，坡长20 m，坡角为30°，地下水埋深7 m，其基质吸力水头取为8.14×10^{− 2} m，土体的具体参数列于下表1所示。

表1 花岗岩残积土计算参数Table 1 Calculation parameters of granite residual soil

表中：θ_r为土体残余含水率；θ_s为土体饱和含水率；K_s为土体饱和渗透系数( m•d^-1)；α、n为土体特征曲线拟合参数。

根据该边坡在不同深度处的初始含水率数据，利用式（24）进行拟合，该拟合结果见图4。

图4 土体初始含水率拟合结果Fig.4 Fitting results of soil initial moisture content

由上图可知，R²=0.974拟合优度接近于1，拟合效果具有可信度，依此便可以得到土体初始含水率随深度的变化规律曲线，具体表达式见下式（38）。

本文改进M-L模型主要有弱降雨—自由入渗阶段和强降雨—有压入渗阶段，依据该地区降雨数据，针对本试验边坡模型选取两个降雨强度，分别为R=8.0和R=50.8，另外在降雨强度略大于土体饱和入渗率时，存在一个弱降雨向强降雨转化的阶段，故增加一个降雨强度为R=26.0，进行验证。水分在入渗过程中，湿润锋上部的土体体积含水率通常不饱和，故此湿润锋上部土体体积含水率取为θ_w =0.9,θ_s =0.378，将a=10.686，d=7m， β=30^o，θ₀ =0.23，R=0.192代入式（29），便可以求得自由入渗阶段降雨强度为8.0的湿润锋入渗深度随时间的变化曲线，绘出该曲线如图5所示。将R=0.624， S_f =0.0814 m， K_s =0.594与上述参数分别代入式（34）、式（36），便可求得I_p =0.2643 m， t_p =0.489 d，将上述数据分别代入式（29）及 式（37）便可以求得降雨强度R=26.0时，在整个降雨入渗过程中，湿润锋入渗深度随时间的变化曲线。绘出该湿润锋入渗深度-时间变化曲线如图6所示。将R=1.22与上述参数代入式（37），便可求得降雨强度R=50.8时，在整个降雨入渗过程中，湿润锋入渗深度随时间的变化曲线。绘出该湿润锋入渗深度-时间变化曲线如图7所示。

文献^[6]中的改进降雨入渗模型是基于M-L模型，考虑倾角和初始含水率呈反比例函数分布推导得出的，该反比例型初始含水率分布见式（39），式（40）为自由入渗和有压入渗两个阶段的表达式，本文将其作为对比模型1。文献^[7]中的改进降雨入渗模型则是基于G-A模型，考虑倾角的影响推导得出的，式（41）为该模型的表达式，本文将其作为对比模型2。式（40）及式（41）中各物理量意义同本文。

式（39）中： A、 B为拟合参数； d为地下水埋深（m）；θ₀为坡面表层土初始含水率；θ_s为饱和含水率。针对本试验边坡，各参数分别取值为A=0.8306，B=3.0709，d=7 m、θ₀ =0.232，θ_s =0.42， h₀ =0 。式（4 1）中K_w取为土体的饱和渗透系数。将上述两个模型得出的解以及本文M-L改进模型的解与数值解进行对比分析，分析结果列于表2。

表2 不同降雨条件下湿润锋入渗深度随时间的变化关系Table 2 Variation of wetting front infiltration depth with time under different rainfall conditions

工况2中：本文改进模型积水时刻为0.536 d，此时湿润锋入渗深度2.183 m；对比模型1积水时刻为0.722 d，湿润锋入渗深度3.523 m；对比模型2积水时刻为0.569 d，入渗深度2.182 m。

由上表可以看出，本文改进降雨入渗M-L模型的计算结果与模型1、模型2及数值解均具有较好的一致性，按照本文的模型来考虑初始含水率非线性的指数型分布是可靠的。分析上表，即可看出在工况1和工况3两种情况下，本文提出的模型计算结果与模型-1的计算结果相比，湿润锋入渗深度较小；在工况-2情况下，本文的改进模型计算结果较模型-1的计算结果相比较大，这主要的原因在于：模型-1中的初始含水率分布是假定为均匀分布与实际土体中含水率分布的不均匀所导致的。这得出的结论与文献^[6]得到的结论也是一致的。由下图5与图6对比分析可以看出模型-2的计算结果与本文模型变化规律是大致相同的，这主要是因为两种模型都是基于指数型初始含水率的分布形式推导得出的。

图5 工况1湿润锋-时间变化曲线Fig.5 Wetting front-time curve of working condition 1

图6 工况2湿润锋-时间变化曲线Fig.6 Wetting front-time curve of working condition 2

图7 工况3湿润锋-时间变化曲线Fig.7 Wetting front-time curve of working condition 3

对比上述三图可以得出：在弱降雨-自由入渗阶段，湿润锋的入渗规律随时间基本上呈线性变化的，这大致上有两方面的原因：其一是该阶段的入渗主要由降雨强度控制；其二是降雨强度过小，累积入渗量过小，不足以推进湿润锋向深处发展。在自由入渗—有压入渗转化阶段，本文改进模型的积水时刻较对比模型1早了4.5 h，在工程上来讲是更偏于安全的，完全可以接受。而对比模型2与其他三个模型的解在同一时刻，湿润锋入渗深度差值较大，吻合度较低，主要是因为对比模型2将初始含水率假定为均匀分布，过于依赖初始含水率的取值，敏感性太大，这充分说明初始含水率非线性分布的必要性。积水时刻之后，本阶段大致以湿润锋入渗深度为4m时为界，变化曲线逐渐呈现出指数变化的形式，其原因大致有两个，其一是越靠近地下水位，其土体初始含水率便越大，入渗所需的雨量便越小，则在降雨强度不变的条件下，湿润锋的入渗速率便越快。其二是坡面上方形成径流，导致湿润锋下渗的动力增加，下渗的速率加快。同理，在强降雨-有压入渗阶段，大致以3.5 m为界，变化曲线呈现指数变化规律，也可由此解释。