由于地质、水文、风化及生物过程等的复合作用，边坡土体多表现为明显的层状结构。不少学者对边坡降雨入渗过程进行过大量的研究工作，并建立了各自的求解方法。姚海林等^[1]研究降雨入渗对非饱和膨胀土边坡稳定性的影响，从降雨强度、饱和渗透系数、降雨持时等方面对膨胀土边坡稳定性进行参数研究。钱纪芸等^[2]通过对降雨条件下的边坡进行离心模型试验发现，随着降雨量的增大，边坡的位移逐渐发展，发生明显变形的区域也逐渐变大，且变形区主要集中在边坡表面。陈善雄等^[3]根据极限平衡分析方法建立了考虑雨水入渗下的非饱和边坡稳定性的方法。孙冬梅等^[4]通过建立水–气二相流数学模型对土质边坡降雨入渗过程进行了分析。蒋忠明等^[5]研究了降雨入渗条件下3维边坡渗流场的变化过程。但上述研究均未考虑坡面产流的影响。张家发^[6]、张书函^[7]、Prudic^[8]等通过有限元法修正降雨入渗边界，实现了降雨条件下坡面径流和入渗耦合数值模拟，可应用于一些简单边坡，但对复杂边界条件，如连续坡面等情况适用性不足。张培文等^[9-10]采用有限元法，通过考虑降雨过程中坡面边界条件的变化，假定入渗率迭代求解入渗流量，将降雨–坡面径流–入渗视为一个系统，建立了相应的耦合方程并求解，但该方法需要分别计算入渗量和入渗速率，迭代计算量大，难以保证较高的计算效率和计算精度。童富果等^[11]推导了有限元方法的降雨入渗与坡面径流直接耦合计算模型与求解方法，解决了容水度不连续问题，避免了迭代计算，提高了计算速度与精度，但却忽略降雨入渗率参数值，评价方法不够精确。田东方等^[12]通过有限元法，将滑坡渗流计算域缩小为滑体，以避免因滑体与滑床（滑带）渗透性差异巨大引起的数值计算困难，但忽略了滑体与滑床之间的流量交换，误差偏大。传统层状土边坡降雨入渗研究^[13-15]，通过改进GA入渗模型，提出适合多层非饱和土边坡降雨入渗的计算方法，但忽略了坡面径流带来的影响，结果偏小。

上粗下细型层状土可增加土壤入渗、抑制水分蒸发、降低水分深层渗漏，从而提高土壤的持水能力。不同结构土层之间的渗透系数存在较大差异，降雨条件下边坡的渗流过程也较为复杂，导致其稳定性的变化规律变得更加复杂。因此，上粗下细型层状土的降雨入渗分析对土壤水分运动研究具有重要意义。作者在前人研究的基础上，基于水量平衡原理，建立坡面径流和上粗下细型层状土地下渗流耦合模型，并采用Python编制相应的程序对某边坡降雨过程中径流和渗流变化进行计算模拟。

1 耦合模型控制方程

降雨入渗过程是时空变化的动态过程，因为降雨强度、降雨量及边坡的物质构成与土层的不同，在强降雨作用下会产生坡体饱和–非饱和入渗及坡面地表产流现象。坡体饱和非饱和入渗遵循坡体土壤水分运动控制方程，坡面地表产流遵循坡面径流控制方程，耦合方程由坡体土壤水分运动控制方程与坡面径流的控制方程组成。

雨水入渗及重分布表现为土体剖面体积含水率随时间的变化及湿润锋前移的过程，如图1所示。图1中：表层土体厚度为L₁，表层土初始含水率为θ_1i，饱和含水率θ_1s；次层土假设向下无限延伸，初始含水率为θ_2i，饱和含水率θ_2s。在0～t₁降雨入渗初期，t₁时刻湿润锋入渗到表层土体Z₁处，湿润锋Z₁以上部分土体达到饱和状态。在t₁～t₂时间段，持续降雨入渗，t₂时刻降雨停止，湿润锋跨过分界面到达下层土体Z₂处，湿润锋Z₂以上部分土体达到饱和状态。当降雨停止后，原入渗雨水在湿润锋下方非饱和区域土体的基质吸力和自身重力作用下继续向下运动，土体中的雨水处于重分布阶段，即t₂～t₃时间段；t₃时刻湿润锋到达Z₃处，由于入渗雨水重分布，湿润锋上方土体和下方土体均处于非饱和状态，表层土的含水率为θ₁，次层土的含水率为θ₂。

1.1 入渗阶段控制方程 1.1.1 入渗阶段控制方程

降雨条件下双层土入渗控制方程可用Moore入渗模型^[16]描述：

$ f = \left\{\! \begin{array}{l} {K_{1{\rm{s}}}}\bigg(1 + \dfrac{{\Delta {\theta _{1{\rm{s}}}} \cdot \left( {{S\!_{1{\rm{f}}}} + h} \right)}}{F} \bigg) , {\textit{Z}} \le {L_1};\\ {K_{2{\rm{s}}}} \cdot \dfrac{{H + F - {F_1}}}{{E + F - {F_1}}} ,{\textit{Z}} > {L_1} \end{array} \right. $

(1)

式中： $f$ 为入渗速率； ${K_{1{\rm{s}}}}$ 、 ${K_{2{\rm{s}}}}$ 分别为表层土和次层土的饱和渗透系数；Z为概化湿润锋深度； ${L_1} $ 为表层土厚度；h为径流深度； $H = \Delta {\theta _{2{\rm{s}}}}({L_1} + {S\!_{2{\rm{f}}}})$ ； $E = {L_1}\Delta {\theta _{2{\rm{s}}}}({K_2}/{K_1}) $ ； $F $ 为累计入渗量； ${F_1} = {L_1} \cdot \Delta {\theta _{1{\rm{s}}}} $ ； $\Delta {\theta _{{\rm{1s}}}}{\rm{ = }}{\theta _{1{\rm{s}}}} - {\theta _{1{\rm{i}}}} $ ， $\Delta {\theta _{2{\rm{s}}}}{\rm{ = }}{\theta _{2{\rm{s}}}} - {\theta _{2{\rm{i}}}}$ ， ${\theta _{1{\rm{i}}}}$ 、 ${\theta _{1{\rm{s}}}}$ 分别为表层土的初始含水率和饱和含水率， ${\theta _{2{\rm{i}}}}$ 、 ${\theta _{2{\rm{s}}}}$ 分别为次层土的初始含水率和饱和含水率； ${S\!_{1{\rm{f}}}}$ 、 ${S\!_{2{\rm{f}}}}$ 分别为表层土和次层土在湿润锋处的平均基质吸力。

坡面产流前，坡体入渗速率等于降雨强度，即f=I，记坡面产流时刻为t_p，该时刻累计入渗量为F_p，则有F_p=I·t_p。根据式（1）可以求出坡面产流时刻入渗量为：

$ {F_{\rm{p}}} = \left\{\! \begin{array}{l} \dfrac{{\Delta {\theta _{1{\rm{s}}}} \cdot {S\!_{1{\rm{f}}}}}}{{I/{K_1} - 1}} , {\textit{Z}} \le {L_1};\\ \dfrac{{H - E \cdot I/{K_2}}}{{I/{K_2} - 1}} {\rm{ + }}{F_1} ,{\textit{Z}} > {L_1} \end{array} \right. $

(2)

由t_p=F_P/I_p可求得积水时刻t_p。

Green–Amp假设湿润锋至入渗面间的土体完全饱和，根据累计入渗量和入渗率的关系，即 $f{\rm{ = d}}F/{\rm{d}}t$ ，联立式（1）可得出累计入渗量F与时间t的隐式方程：

$ \left\{\! \begin{array}{l} {K_1}(t - {t_{\rm{p}}}) = F - \Delta {\theta _{1{\rm{s}}}}({S\!_{1{\rm{f}}}} + h) \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\ln \bigg[1 + \dfrac{F}{{\Delta {\theta _{1{\rm{s}}}}({S\!_{1{\rm{f}}}} + h)}}\bigg]{\kern 1pt} {\kern 1pt} ,{\rm{ }}{\textit{Z}} \le {L_1};\\ {K_2}\left( {t - {t_{\rm{p}}}} \right) = F - (H - E) \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\ln \bigg[1 + \dfrac{F}{{H - {F_1}}}\bigg]{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ,{\kern 1pt} {\kern 1pt} {\textit{Z}} > {L_1} \end{array} \right. $

(3)

再根据累计入渗量 $F $ 等于土层含水量的变化量 $(({\theta _{\rm{s}}} - {\theta _{\rm{i}}})Z)$ ，可求得湿润锋深度。

1.1.2 重分布阶段控制方程

降雨停止后，坡体内水分在自身重力和湿润锋下方土体基质吸力的双重作用下继续下渗，记停雨时刻为t_s，停雨时刻的入渗量为F_s，次层细质土由饱和状态变为非饱和状态时所对应的时刻为t_c。重分布过程分为两个阶段。

第1阶段，表层粗质土处于非饱和状态而次层细质土处于饱和状态，即 ${t_{\rm{s}}} < t \le {t_{\rm{c}}}$ 。根据式（4）可确定表层土含水率分布：

$\frac{{{\rm{d}}{\theta _1}}}{{{\rm{d}}t}} = - \frac{{{K_{2{\rm{s}}}}}}{{{L_1}}} \cdot \bigg[1 + ({\theta _{2{\rm{s}}}} - {\theta _{2{\rm{i}}}}) \cdot \frac{{{S\!_{2{\rm{f}}}} - {\psi _1}({L_1})}}{{{F_{\rm{s}}} - ({\theta _1} - {\theta _{1{\rm{i}}}}) \cdot {L_1}}}\bigg]$

(4)

式中， ${\psi _1}({L_1})$ 为两层土接触面处的表层土基质吸力， ${\theta _1}$ 为表层土的含水率，其他符号同式（1）。

第2阶段，表层土与次层土均处于非饱和状态，即 $t > {t_{\rm{c}}}$ 。根据式（5）～（7）可确定土体的含水率分布，具体参见文献[17]。

当 $t = {t_{\rm{c}}}$ 时，表层土的含水率为：

$\theta _1' = {\theta _{1{\rm{r}}}} + ({\theta _{1{\rm{s}}}} - {\theta _{1{\rm{r}}}}) \cdot {\bigg(\frac{{{\psi _{1{\rm{b}}}}}}{{{\psi _{2{\rm{b}}}}}}\bigg)^{{\lambda _1}}}$

(5)

式中， $\theta _1'$ 为t_c时刻表层粗质土表面的体积含水率， ${\psi _{1{\rm{b}}}}$ 、 ${\psi _{2{\rm{b}}}}$ 分别为上、下层土进气值， ${\theta _{1{\rm{r}}}}$ 为表层土残余饱和度， $\lambda $ 为Brooks（1964）模型拟合公式^[18]中用于反映土体孔径分布特征的经验系数。

当 ${\textit{Z}} \le {L_1}$ 时，湿润锋以上的含水量分布：

$\begin{aligned}[b] \frac{{{\rm{d}}{\theta _1}}}{{{\rm{d}}t}} =& - \frac{{({\theta _1} - {\theta _{1{\rm{i}}}})}}{{{F_{\rm{s}}}}} \cdot\\ & \bigg[K({\theta _{1{\rm{i}}}}) + K({\theta _1}) + \frac{{({\theta _1} - {\theta _{1{\rm{i}}}}) \cdot {K_{1{\rm{s}}}} \cdot G({\theta _{1{\rm{i}}}},{\theta _1})}}{{{F_{\rm{s}}}}}\bigg] \end{aligned} $

(6)

当 ${\textit{Z}} > {L_1}$ 时，湿润锋以上的含水量分布：

$ \begin{aligned}[b] \frac{{{\rm{d}}{\theta _2}}}{{{\rm{d}}t}} =& - \frac{1}{{\dfrac{{{\theta _1} - {\theta _{1{\rm{r}}}}}}{{{\theta _2} - {\theta _{2{\rm{r}}}}}} \cdot \dfrac{{{\lambda _1}}}{{{\lambda _2}}} \cdot {L_1} + \dfrac{{{F_{\rm{s}}} - ({\theta _1} - {\theta _{1{\rm{i}}}}) \cdot {L_1}}}{{({\theta _2} - {\theta _{2{\rm{i}}}})}}}} \cdot \\ &\left[ {K({\theta _{2{\rm{i}}}}) + K({\theta _2}) + \frac{{({\theta _2} - {\theta _{2{\rm{i}}}}) \cdot {K_{2{\rm{s}}}} \cdot G({\theta _{2{\rm{i}}}},{\theta _2})}}{{{F_{\rm{s}}} - ({\theta _1} - {\theta _{1{\rm{i}}}}) \cdot {L_1}}}} \right] \end{aligned} $

(7)

式（6）～（7）中： $G$ 为毛细驱动力； ${\theta _{2{\rm{r}}}} $ 为次层土残余含水量； ${K({\theta )}}$ 为土层含水为 $ \theta$ 时，对应的渗透系数；其他符号同前。

1.2 坡面径流控制方程

Saint–Venant坡面流1维微分方程^[19]可以很好地描述坡面径流过程。因此，采取1维运动波理论来描述坡面径流过程：

$ \left\{\! \begin{array}{l} \dfrac{{\partial h}}{{\partial t}} + \dfrac{{\partial q}}{{\partial x}} = I(t)\cos \;a - f(x,t),\\ q = \dfrac{1}{n}{h^{5/3}}\sin \;{a^{1/2}} \end{array} \right. $

(8)

式中，n为Manning粗糙系数^[20]，a、I、 $f$ 分别为坡角、降雨强度和入渗率，q为单宽流量，v、h为坡长x处的流速和水深。

坡面径流简图如图2所示。

2 控制方程离散和模型耦合过程

假设坡面为均匀降雨，降雨强度在坡面上处相同，仅随时间发生变化。坡面径流和坡面渗流通过坡面的流深和入渗相互作用，方程（1）、（4）、（6）、（7）和（8）共同构成耦合模型的控制方程。坡面是否产流决定了坡面渗流过程的上边界条件，而坡面是否积水是由边坡入渗率决定的。两者通过入渗率进行耦合，并采用交替迭代的方法进行坡面径流和渗流的耦合求解。Pressimann差分法^[21]作为计算流体力学一种经典算法，具有结构简单、稳定性较好、精度较高等优点，应用广泛。Pressimann的四点隐式差分格式为：

$\left\{\! \begin{array}{l} f{\bigg|_M} = \dfrac{{\theta (f_{i + 1}^{j + 1} + f_i^{j + 1}) + (1 - \theta )(f_{i + 1}^j + f_i^j)}}{2}, \\ \dfrac{{\partial f}}{{\partial x}}\bigg|_M \bigg. = \theta \dfrac{{f_{i + 1}^{j + 1} - f_i^{j + 1}}}{{\Delta x}} + (1 - \theta )\dfrac{{f_{i + 1}^j - f_i^j}}{{\Delta x}}, \\ \dfrac{{\partial f}}{{\partial t}}\bigg|_M \bigg. = \dfrac{{f_{i + 1}^{j + 1} + f_i^{j + 1} - f_{i + 1}^j - f_i^j}}{{2\Delta t}} \end{array} \right.$

(9)

式中： $\theta $ 为加权系数， $0 \le \theta \le 1$ ； $f$ 代表需要求解的量，即流量 $q$ 或水位 $h$ 或流速 $v$ ；上标 $j$ 为时间节点，下标 $i$ 为空间节点。采用式（9）形式的Pressimann差分格式对径流方程（8）进行离散，可得：

$({h^{5/3}})_{i + 1}^{j + 1} - ({h^{5/3}})_{i}^{j + 1} + {D_1}(h_{i + 1}^{j + 1} + h_i^{j + 1}) = {D_2}$

(10)

式中：

$ \left\{\! \begin{array}{l} {D_1} = \dfrac{n}{{\sqrt {{S\!_0}} }}\dfrac{{\Delta x}}{{2\Delta t\theta }}{\text{，}}\\ {D_2} = \dfrac{{n\Delta x}}{{2\theta \sqrt {{S\!_0}} }}\bigg[(Q)_i^j + (Q)_{i + 1}^j\bigg] - \dfrac{{1 - \theta }}{\theta }\bigg[{(h_{i + 1}^j)^{5/3}}-\\ \;\;\;\;\;\;\;\;\; {(h_i^j)^{5/3}}\bigg] + \dfrac{n}{{\sqrt {{S\!_0}} }}\dfrac{{\Delta x}}{{2\Delta t\theta }}(h_{i + 1}^j + h_i^j) \end{array} \right. $

(11)

$Q = I(t)\cos \;a - f(x,t)$ ，为坡面径流方程源项， $S\!_0$ 为坡面坡度。

初始条件（ $j = 0$ ）和边界条件（ $i = 0$ ）分别为：

$\left\{\! \begin{array}{l} h_i^0 = 0, \\ h_0^j = 0 \\ \end{array} \right.$

(12)

将边界条件和初始条件代入方程（10）中得：

${(h_1^1)^{5/3}} + {D_1}h_1^1 = {D_2}|_0^0 + {(h_1^1)^{5/3}} - {D_1}h_0^1$

(13)

当 $j = 0$ 时，入渗速率和降雨均为0，则源项 $N_0^0 = 0$ ， $N_1^0 = 0$ 。采用牛顿迭代法可求出上述一元非线性方程的近似实根 $h_1^1$ ，将边界条件 $h_0^1$ 和 $h_1^1$ 求均值代入渗流方程（3）可求得相邻两个空间节点间的土体的入渗量 $F$ ；然后，将入渗量 $F $ 与 $(h_0^1 + h_1^1)/2$ 代入方程（1），进而求得两个空间节点间的土体的入渗速率 $f_{}^{}$ 。根据GA入渗模型均匀入渗假设 $f_0^1 = f_1^1 = f$ ，可求得源项 $N_0^1 = I_0^1 - f_0^1$ ， $N_1^1 = I_1^1 - f_1^1$ 。同理，根据 $h_1^1$ 、 $N_0^1$ 、 $N_1^1$ 可求得 $h_1^2$ ，如此反复可求得 $i = 1$ 的空间节点的各时刻的径流深度，同理可一次求得 $i=2,3,4\cdots$ 各节点的各个时刻径流深度及入渗过程。当 $t > {t_{\rm{s}}}$ 时，将入渗量 $F $ 代入重分布模块，计算次层细质土由饱和状态变为非饱和状态时所对应的时刻 ${t_{\rm{c}}}$ ，判断重分布状态，通过4阶Runge–Kutta法求解重分布微分方程，计算土体的含水率，从而计算土体的湿润锋深度。耦合流程图如图3所示。

3 模型验证和算例数值分析

为验证提出的耦合模型的正确性，根据文献[22]中的多层土降雨入渗试验，当降雨时间在500 min以内时，土层分布为砂土–黏土–壤土的降雨入渗试验中土体湿润锋仍然处于第2层土中。取多层土中表面两层土入渗过程进行数值模拟，计算域与试验装置保持一致，如图4所示。图4中：左右两侧为不透水边界；上层土层为砂土，厚度为15 cm，下层为黏土，土层参数如表1所示。试验坡长为29 cm，坡度为0.01°，粗糙系数n为0.025；降雨时间为400 min，降雨强度随时间发生变化取自文献[22]中土层分布为砂土–黏土–壤土的降雨入渗试验，如图5所示。

为确保差分格式计算结果收敛，取计算空间步长Δx=1.5 cm，时间步长Δt=1 min。当t=4 min，降雨强度变为1.3 cm/min时，坡面产生径流；t=22 min时，坡面径流强度达到峰值，坡脚处的径流速率如图6所示。

由图6可知，本文耦合方法计算的坡面径流与文献模拟^[22]计算结果基本吻合。

为进行对比研究，分别对双层边坡耦合模型与不考虑坡面径流、简化坡面边界条件的Moore双层入渗模型（以下称为简化模型）进行数值模拟，典型剖面及土层概况如图7所示。

假设该模型向下无限延伸，上层土层为砂质黏性壤土，厚度为0.4 m，下层为粉质黏土，为典型的粗–细型边坡。坡长L为50 m，粗糙率n为0.015，坡度a=6°，降雨强度I=2.5 cm/h，降雨历时为12 h且为均匀降雨。土体的力学参数和渗透系数根据Rawls等^[23]的统计数据确定，具体见表2。

双层边坡模型坡面分为10个单元、11个节点。降雨历时12 h，重分布阶段历时120 h，时间步长为0.1 h。

图8为简化模型与耦合模型计算所得的入渗速率随时间的变化趋势。由图8可以看出，耦合模型计算结果与不考虑坡面径流影响的简化模型模拟结果基本一致。根据式（2）可以计算出坡面在降雨初期（t=0.45 h）开始产流，坡面产流之前入渗速率等于降雨速率。产流后，由于坡面径流的出现，坡面土体进入暂态饱和状态，土体的含水率增加，基质吸力不断减小，降雨速率随着降雨历时的增加而逐渐减小，并趋于稳定。由于坡面存在径流深度，耦合模型比简化模型先进入次层土体，当湿润锋进入次层土时，入渗速率快速减小，最终次层土体达到完全饱和，入渗速率减小到稳定值，即次层土体的饱和入渗速率。

图9、10分别为体积含水率在土层剖面的分布及简化模型与耦合模型计算所得的湿润锋深度随时间的变化趋势。

由图9、10可见：当t=6.5 h左右时，由于当湿润锋到达次表面时，上层土的流量与次层土保持相同，次层土的基质吸力比上层土大，所以湿润锋短时间内在次层土的基质吸力作用下会快速向下推进一段距离。同时，上层土的基质势远大于次层土的基质势；为了使上下层土壤间的基质势相同，交界面的入渗通量必然会大量减少，次层土体快速达到饱和状态，从而导致了入渗速率大幅度降低，最终减小至次层土饱和入渗速率。因此可知，对于上粗下细型层状边坡，入渗速率主要由次层细质土控制。当t=12 h，降雨停止后，土体中已入渗的雨水在自身重力和下方土体基质吸力作用下继续向下运动，在12～100 h时段向下推移深度小于0～12 h时段，而100～132 h时段，湿润锋只往下渗透了5 cm左右，重分布基本趋于稳定。

沿坡面不同位置处径流随时间变化如图11所示。由图11可知，坡面产流后，径流深度迅速上升，约在6 h左右基本达到稳定。坡面较高位置的点，稳定时水深较浅，并率先达到稳定水位。沿坡面向下，水深逐渐增加，且达到稳定水位的时间相对滞后。当t=6.5 h，湿润锋到达下层土体时，由于入渗速率降低，降雨强度不变，径流深度再一次上升后很快达到稳定。降雨停止后，坡面径流在很短的时间内衰减为0。

沿坡面方向不同位置处湿润锋深度随时间变化如图12所示。由图12可见，沿坡面向下，同一时刻、不同径流深度坡面处的湿润锋深度相差不大，湿润锋深度随径流深度增大而增大。

综上所述，双层结构边坡入渗过程与均质土边坡入渗过程不同，当湿润锋到达土层交界面时，入渗速率和湿润锋深度会产生向下突变，导致径流深度发生变化，从而影响土体的入渗速率。层状土边坡实际入渗过程不仅与时间有关，还与坡面径流深度有关，湿润锋深度随径流深度增大而增大。

4 结　论

通过算例验证和数值对比分析表明，坡面产流时，层状土边坡降雨入渗过程采用简化入渗边界进行计算，结果偏小。采用本文耦合模型与计算方法能够较好地反映实际层状土边坡产流和降雨入渗过程，可为类似的层状土边坡稳定、水土流失和山地洪水分析提供计算依据。