如图1，将原点取在坡面处，轴以垂直坡面向下为正。根据达西定律，地表处的入渗率为：

(1)

式中：i为入渗率(m/s)；K_s为饱和渗透系数(m/s)；z^*_f为湿润锋深度(m)；H为积水深度(m)；γ为边坡坡度(°)；s_f为湿润锋处基质吸力(m)。

根据水量平衡原理，有：

I=(θ_s-θ_i)z^*_f(2)

根据入渗率与累积入渗量之间的关系，有：

(3)

将式(1)代入式(2)后化简，得：

(4)

对式(4)积分，并代入初始条件z^*_f |t=0化简，得：

(5)

将式(2)代入式(5)化简，得：

(6)

当边坡表面积水很浅或者雨水入渗的时间很长时，积水深度H相对于湿润锋深度z^*_f可以忽略，则式(1)可以写为：

(7)

联立式(2)和式(7)化简，得：

(8)

在K_s、H、s_f、γ、θ_s和θ_i已知的条件下，由式(5)求得z^*_f(t)，然后再由式(6)求得I(t)，最后再由式(7)求得i(t)。但由式(5)一般得不到z^*_f(t)的显式表达式，需要采用数值方法进行计算。

积水入渗阶段(t>t_p)，由于式(6)中的I是表示积水入渗以后的累积入渗量，而积水入渗不是从t=0开始的，故需要对式(6)进行修正：

(9)

式中：t_s为由t=0开始积水到I=I_p时所需的时间(s)；t_p为地表实际开始积水的时间(s)；I_p为地表实际开始积水入渗以后的累积入渗量。

式(9)中的I_p、t_s、t_p可由式(10)计算得到：

(10)

非积水入渗阶段，入渗率可表示为：

i=Rcosγ t≤t_p(11)

式中：R为降雨强度(m/s)。

根据非饱和土渗流达西定律有：

(12)

根据坐标转换关系有：

z=x^*sinγ+z^*cosγ (13)

联立式(12)和式(13)，得：

(14)

根据Green-Ampt入渗模型的假定，湿润锋上部土体的含水率分布均匀，则式(14)中的第二项为0。联立式(11)和式(14)，可得到：

R=K(h_c)(15)

根据Van Genuchten^[12]非饱和渗透系数拟合模型有：

(16)

联立式(15)和式(16)，得：

(17)

式中：α、m、n均为拟合参数

利用式(17)求得h_c，然后代入Van Genuchten土水特征曲线拟合方程^[13]计算湿润锋以上土体的体积含水率θ_w：

(18)

式中：θ_w为非饱和土含水率；θ_r为非饱和土残余含水率。

根据水量均衡原理，可得：

I=Rt=(θ_w-θ_r)z^*f(19)

故降雨持续时间t的入渗深度为：

(20)

通过式(20)便能求得z^*_f(t)，而i(t)和I(t)可直接得到，分别等于R和Rt。

改进的降雨入渗模型如图2所示，引入的土体初始含水率反比例函数如下：

(21)

式中：A为土体初始含水率分布参数；B为间接参数；d为地下水埋深(m)；θ₀为地表处初始含水率。

图2 改进降雨入渗模型计算图
Fig.2 Calculation diagram of improved rainfall infiltration model

降雨强度取R=0.45 m/d，采用式(17)计算负压水头h_c=0.000 75 m，然后再用式(18)计算湿润锋以上土体的体积含水率θ_w=0.399 983，最后将A和θ_w代入式(32)便能得到基于改进降雨入渗模型的降雨入渗深度-时间关系曲线，如图8红线所示。同时将θ_i=0.25和θ_w=0.399 983代入式(20)便能得到基于Mein-Larson入渗模型的降雨入渗深度-时间关系曲线，如图8蓝线所示。

在数值模型中施加降雨边界条件，如图6所示，设置降雨时长为2.5 d，时间步长为0.02 d，然后在SEEP/W模块进行瞬态渗流分析。由于篇幅有限，仅列出降雨入渗达到1 m和3 m的瞬态分析结果，见图7，图中蓝色区域为体积含水率为0.4的土体，蓝色虚线为降雨入渗浸润线和地下水位线。提取每个时刻的降雨入渗深度绘制降雨入渗深度-时间关系曲线，如图8黑线所示。

将改进模型理论解、Mein-Larson入渗模型理论解和数值解绘于图8，由图可知，降雨入渗在小于2.5 m时，改进模型理论解和数值解存在一定的误差，降雨入渗在2.5~5 m时，改进模型理论解和数值解基本一致。而且降雨入渗深度一定时，离地表越近，入渗时间的误差越大，最大误差约为7.9 h，产生误差的原因是理论解在计算渗透率时忽略了负压水头随深度变化的一项，而负压水头与深度呈递增关系，斜率逐渐减小到最后趋于零，所以理论解与数值解的误差也是随着深度的增加逐渐减少最后趋于零。改进模型理论解较Mein-Larson入渗模型理论解更接近数值解，因此考虑土体的初始含水率分布能够让降雨入渗模型更加完善和准确。

图6 降雨边界布置图
Fig.6 Rainfall boundary layout

图7 瞬态渗流数值模拟结果(0.45 m/d)
Fig.7 Numerical simulation results of transient seepage(0.45 m/d)

图8 降雨入渗深度理论和数值解对比(0.45 m/d)
Fig.8 Comparison of theoretical and numerical solutions of rainfall infiltration depths(0.45 m/d)

降雨强度取0.60 m/d，湿润锋处的平均基质吸力取s_f=0.06 m^[14]，采用式(28)、式(30)、式(27)和式(31)分别求得I_p=0.097 376 m、t_p=0.187 406 s、z^*_f=0.346 42 m和t_s=0.145 250 s，然后采用式(23)计算t≤t_p时的降雨入渗深度，采用式(29)计算时的降雨入渗深度，将计算得到的各时刻降雨入渗深度绘成降雨入渗深度-时间关系曲线，如图 10红线所示。同时将θ_i=0.25、θ_n=0.4和s_f=0.06 m代入式(10)分别求得I_p=0.075 m、t_p=0.125 d、z^*_f=0.300 m和 t_s=0.096 25 s，然后采用式(2)计算t≤t_p时的降雨入渗深度，采用式(9)计算t≤t_p时的降雨入渗深度，将计算得到的各时刻降雨入渗深度绘成降雨入渗深度-时间关系曲线，如图 10蓝线所示。

在数值模型中施加降雨边界条件，如图6所示，设置降雨时长为2.5 d，时间步长为1 h，然后在SEEP/W模块进行瞬态渗流分析。由于篇幅有限仅列出降雨入渗达到2 m和4 m的瞬态分析结果，见图9，图中蓝色区域为体积含水率为0.4的土体，蓝色虚线为降雨入渗浸润线和地下水位线。提取每个时刻的降雨入渗深度绘制降雨入渗深度-时间关系曲线，如图 10黑线所示。

将改进模型理论解、Mein-Larson入渗模型理论解和数值解绘于图 10，由图可知，降雨入渗在小于2.5 m时，理论解和数值解存在一定的误差，降雨入渗在2.5~5 m时，理论解和数值解基本一致。而且降雨入渗深度一定时，离地表越近，入渗时间的误差越大，最大误差约为8.2 h，产生误差的原因是理论解在求解过程中假设湿润锋处的基质吸力为定值 。非饱和土的基质吸力 是随深度的增加而不断减少最后趋于零(地下水位处)，而 0.06 m是土体接近饱和状态的基质吸力，也就是说基质吸力 0.06 m时，越接近地面基质吸力误差就越大，于是便导致了图 10中的误差。改进模型理论解较Mein-Larson入渗模型理论解更接近数值解，因此考虑土体的初始含水率分布能够让降雨入渗模型更加完善和准确。

图9 瞬态渗流数值模拟结果(0.60 m/d)
Fig.9 Numerical simulation results of transient seepage(0.60 m/d)

图 10 降雨入渗深度理论和数值解对比(0.60 m/d)
Fig.10 Comparison of theoretical and numerical solutions of rainfall infiltration depths(0.60 m/d)