地表径流-土壤水流模型的耦合是水文循环模型的重点之一，由于耦合模型的复杂性，部分工程实践和案例忽略了地表积水造成的影响，Fu等^[1]考虑在地表径流产生时直接将边界条件设置成饱和入渗边界，Xiong等^[2]考虑将土体的饱和渗透系数作为判断边界条件的标准，袁俊平等^[3]考虑裂隙条件下的边坡稳定，将降雨边界简化成零水头边界条件.

以上降雨入渗研究，未考虑地表径流作用. 为了反映这一作用，目前已经有不少学者在理论和数值上进行地表径流与土壤水流的模拟，且大致分为2种假设，一种是假设当土体饱和时，地表产生积水，另一种是假设当土体的入渗强度小于降雨强度时，地表将产生积水. 张培文等^[4]考虑了地表积水的现象，但简化了模型，未考虑坡面水流的动量守恒；Zhao^[5]采用第1种假设计算坡面流；Wang等^[6]用有限元法将Navier-Stokes方程同Richards方程耦合，假设饱和积水的条件，考虑复杂地面的积水情况；Zhang等^[7]同样采用有限元法耦合Navier-Stokes方程和Richards方程，也设置了饱和积水的条件，计算边坡的降雨入渗；刘育田等^[8-9]对地表径流方程采用有限差分法求解，土壤水流方程采用有限元法求解，假设饱和积水条件，对比考虑与不考虑地表径流条件下的坡体入渗. Tan等^[10]根据第2种假设，在理论上计算地表径流的质量守恒式，并采用Green-Ampt模型模拟土壤水流，用有限差分法考虑坡面流对无限长边坡造成的影响，该模型并未考虑地表水的动量守恒；Jahnson等^[11]使用Richards方程推导一维、半无限土层的理论解，并依据土壤水流场和地表径流场之间的质量交换将地表径流与土壤水流模型耦合，用有限差分法求解该模型，此模型同样忽略了地表水在坡面的流动；Zhu等^[12]采用假设2，将地表径流方程简化成扩散波方程，并通过有限元法计算边坡降雨入渗的问题；Pato等^[13]采用有限体积法模拟地表径流条件下饱和土体的入渗；Pato等^[14-15]采用有限体积法模拟两者耦合作用，土壤水流模型采用Green-Ampt模型.

分析已有研究成果可知，许多学者在模拟雨水入渗过程中采用Green-Ampt模型，数值方法采用有限元法. G-A模型要求入渗土体中水的初始体积分数是均匀分布的，且在入渗过程中将湿润锋简化为一条线，将入渗土体分为湿润锋上部饱和区域和下部初始体积分数区域2个部分，即湿润锋由初始体积分数变为饱和体积分数的土层厚度被忽略. Green-Apmt入渗模型是对复杂的水入渗过程的简化. Richards方程是在达西定律和质量守恒的基础上，应用数学物理方法推导出的非饱和土壤水分流动方程，其初始条件和边界条件的选择更加灵活，且各单元的水头值，以及各单元表面的通量大小都可以计算得出，只是求解较困难. 在数值方法的选取上，本研究采用有限体积法求解耦合方程，相比于有限单元法，有限体积法要求因变量的积分守恒对任意一组控制体积都得到满足，对整个计算区域，自然也得到满足，而有限元法则是满足了整体的质量守恒，局部守恒无法保证，且有限体积法在计算方面也更为简单.

本研究采用有限体积法耦合求解地表径流和土壤水流方程，用Navier-Stokes方程模拟地表径流，用Richards方程模拟土壤水流，通过对比数值耦合解与SEEP/W计算的解，验证耦合模型的有效性，通过计算土坡降雨入渗的案例，对比分析坡顶和坡底的入渗深度情况.

地表径流模型采用不可压缩流体的Navier-Stokes^[15]方程，其中有2个前提假设：1）相对于水平面的尺寸而言，水头高度是微小的，因此土体表面水流的流速在水的深度范围内是定值，且以水流在深度方向的平均流速作为该定值；2）水流所产生的压力是静水压力，因此方程在积分时采用的竖向的压力场为 $ p\left( {t,x,y} \right) = {{gh{{\left( {t,x,y} \right)}^2}}}/{2} $，其中 $ g $为重力加速度常量，h为水头高度. 一维的Navier-Stokes方程为

(1) $ \left. {\begin{array}{*{20}{l}} {{\partial _t}h + {\partial _x}\left( {hu} \right) = R - I}, \\ {{\partial _t}\left( {hu} \right) + {\partial _x}( {h{u^2} + {{g{h^2}}}/{2}} ) = gh\left( {{S_{0{{x}}}} - {S_{{\text{f}}x}}} \right)} . \end{array}} \right\} $

式中： $ u $为 $ x $方向的流速； $ R $、 $ I $分别为降雨强度和入渗强度； $ {S_{0{{x}}}} $、 $ {S_{{\text{f}}x}} $分别为地表坡度和摩擦项. 水头分布 $ h $以及地表高程 $ {Z_{\text{b}}} $分布注解如图1所示.

一般将式（1）化简成向量形式的欧拉方程组：

(2) $ {\partial _t}{\boldsymbol{W}} + {\partial _x}F\left( {\boldsymbol{W}} \right) = {\bf{0}} ,\; {\boldsymbol{W}} = \left[ h ,\; {hu} \right]^{\rm{T}} . $

式中： $F\left( {\boldsymbol{W}} \right) = \left[ {hu},\; {h{u^2} + {{g{h^2}}}/{2}} \right]^{\rm{T}}$.

进一步化简得到

(3) $ {\partial _t}{\boldsymbol{W}} + A\left( {\boldsymbol{W}} \right){\partial _x}{\boldsymbol{W}} = {\bf{0}} . $

式中： $A\left( {\boldsymbol{W}} \right) = F'\left( {\boldsymbol{W}} \right) = \left[ {\begin{array}{*{20}{c}} 0&1 \\ { - {u^2} + {{g{h^2}}}/{2}}&{2u} \end{array}} \right]$.

可以根据 $ A\left( {\boldsymbol{W}} \right) $是否可对角化将式（3）划分为严格的双曲型方程和非双曲型方程，当 $ h > 0 $时， $ A\left( {\boldsymbol{W}} \right) $是可对角化的，当 $ h = 0 $时， $ A\left( {\boldsymbol{W}} \right) $不可对角化. 在 $ h > 0 $的条件下，矩阵 $ A\left( {\boldsymbol{W}} \right) $的2个特征值分别为 $ u - \sqrt {gh} $和 $ u + \sqrt {gh} $. 由于单波方程的特性，可以通过这2个特征值将流动状态进行划分. 当 $ \left| u \right| < \sqrt {gh} $时，2个特征值符号相反，信息在上游和下游双向传递，该流动状态被称为亚音速状态，此时须上下游各给一个边界条件；当 $ \left| u \right| > \sqrt {gh} $时，信息只向下游传递，该流动状态被称为超音速状态，此时只须在上游给定2个边界条件. 因此，可以用 ${F_{\rm{r}}} ={{\left| u \right|}}/{{\sqrt {gh} }}$来判断流动状态. 地表径流方程的时间步和空间步的离散方式如图2所示.

式（1）在求解时被划分为了2部分，第1部分为去除源项后的方程即式（2），第2部分为源项即式（1）等号右边项，因此在求解方程组时须进行2步计算.

第1步须用有限体积法化简式（2）：

(4) $ W_i^{k + 1} = W_i^k - \frac{{\Delta t}}{{\Delta x}}\left( {F_{i + {1}/{2}}^k - F_{i - {1}/{2}}^k} \right) . $

式中： $W_i^k $为W在第i个单元的第k时间步的变量值； $F_{i + {1}/{2}}^k$和 $F_{i - {1}/{2}}^k$表示单元体 $ i $左右两侧的数值通量，上标 $ k $表示第k时间步； $ \Delta t $为时间步长； $ \Delta x $为空间步长. 目前有较多方法计算该通量的值，例如HLL、Rusanov、VFRoe-ncv等^[16-17]，本研究采用HLL法，通过该方法可以得到一阶精度的解，再采用线性重构技术^[16]（例如MUSCL、ENO），可以使方程的空间项具有二阶精度. 该方法的核心是捕捉由间断引起的激波，使计算出来的结果不发生振荡. 可以通过化简将式（4）变成与时间步相关的形式：

(5) $ W_i^{k + 1} = W_i^k - \Delta t{\mathit{\Phi}} \left( {{W_i^k}} \right) . $

可以使用heun法（也被称为改进的欧拉法）保证时间项具有二阶精度，该方法须对时间项进行3步计算以修正式（5）. 修正计算公式如下：

(6) $\left. \begin{array}{*{20}{l}} W_i^ * = W_i^k - \Delta t{\mathit{\Phi}} \left( {{W_i^k}} \right) ,\;\; W_i^{ * * } = W_i^ * - \Delta t{\mathit{\Phi}} \left( {{W^ * }} \right), \\ W_i^{k + 1} = {{(W_i^k + W_i^{ * * })}}/{2} . \end{array} \right\}$

式中： $ W_i^{ * * } $、 $ W_i^ * $为预测值. 通过该方法可以使时间项具有二阶精度. 待求时间项如果需要迭代求解则称为隐式格式，如果可以通过方程组直接求出则称为显式格式. 对显式格式的方程而言，在将微分方程离散成差分方程时须考虑到方程解的一致性，因此只有在方程满足相容性、收敛性、稳定性的情况下才能得出收敛于微分方程的解，在该条件下CFL^[17]数要满足特定的范围才能得出正确的解. CFL数的计算方法如下：

(7) $ c\frac{{\Delta x}}{{\Delta t}} \geqslant \mathop {\sup }\limits_{(t,x)} \;\left\{ {\left| {u\left( {t,x} \right)} \right| + \sqrt {gh\left( {t,x} \right)} } \right\} . $

式中： $ {c} $为可以控制求解精度的参数.

第2步对源项进行处理，源项包括 $ R $、 $ I $、 ${S_{0{{x}}}}$、 $ {S_{{\text{f}}x}} $共4个部分，针对摩擦项有Manning-Strickler摩擦准则和Darcy-Weisbach准则^[18]. 本研究采用Darcy-Weibach准则，对应的表达式为 ${S_{{\text{f}}x}} = \dfrac{{f {{u}}^2}}{{8gh}}$. 其中 $ f $是一个无量纲系数； ${S_{0{{x}}}} = {\partial _x}{Z_{\rm{b}}}$. 为了保证数值格式的和谐性^[19]，采用静水重构技术^[20]，目的是保证水头的非负性，以及保证在计算静止状态时，计算的水头能保持稳定性. 对于地表摩阻项，一般显式处理会导致模型不稳定，尤其是当水头较小时该问题更为突出，而采用全隐格式则会产生较大的计算负担，因此采用半隐式^[21]的方法. 降雨强度R和入渗强度I则不会涉及到数值计算的稳定性，因此基于不增加计算难度的考虑，将它们作为显式处理即可.

变饱和孔隙介质流动模型来用古典的非线性Richards方程^[22]描述，该方程由质量守恒和达西定律推导出来，并且可以被表示成3种不同的形式，如下所示.

1）以 $ h $为基本变量的方程：

(8) $ \left[ {C\left( h \right) + {S_{\rm{s}}}{S_{\rm{w}}}} \right]\frac{{\partial h}}{{\partial t}} = \nabla \cdot \left[ {{\boldsymbol{K}}\left( h \right) \nabla \left( {h + z} \right)} \right] + {Q_{\rm{s}}} . $

式中： $ C\left( h \right) $为容水度，表达式为 ${{\partial \varphi }}/{{\partial h}}$， $\varphi$为水的体积分数； ${Q_{\rm{s}}}$为体积源汇项； ${S_{\rm{s}}}$为单位储水量； ${S_{\rm{w}}}$为饱和系数，S_w= $\varphi / \eta $， $ \eta $为孔隙度； $ z $为高程； $ {\boldsymbol{K}}\left( h \right) $为水力传导系数.

2）以 $\varphi$为基本变量的方程：

(9) $ \frac{{\partial \varphi }}{{\partial t}} + {S_{\rm{s}}}{S_{\rm{w}}}\frac{{\partial h}}{{\partial t}} = \nabla \cdot \left[ {{\boldsymbol{D}}\left( \varphi \right) \nabla \varphi } \right] + \nabla \cdot \left[ {{\boldsymbol{K}}\left( h \right) \nabla z} \right] + {Q_{\rm{s}}} . $

式中： $ {\boldsymbol{D}}\left( \varphi \right) $为水力扩散率.

3）混合形式的方程：

(10) $ \frac{{\partial \varphi}}{{\partial t}} + {S_{\rm{s}}}{S_{\rm{w}}}\frac{{\partial h}}{{\partial t}} = \nabla \cdot \left[ {{\boldsymbol{K}}\left( h \right) \nabla \left( {h + z} \right)} \right] + {Q_{\rm{s}}} . $

在确保每一时间步中 $ h $和 $ C\left( h \right) $足够小的情况下，第1）种方程能够有较好的质量守恒性，但其在高度非线性条件下会产生较大的质量误差^[23]，这是由于 $ \dfrac{{\partial \varphi }}{{\partial t}} $只在连续空间下等于 $ C\left( h \right)\dfrac{{\partial h}}{{\partial t}} $，而在离散条件下不是严格相等的，且 $\varphi \left( h \right) $本身也是非线性的函数，又加大了误差. 第2）种方程在连续和非连续空间条件下都有完全的质量守恒性，然而在接近饱和状态时水力扩散系数 $ {\boldsymbol{D}}\left( \varphi \right) $将趋近于无穷，饱和态的水头与体积分数的关系并不存在^[24]，因此在接近饱和时该方程将不再适用. 第3）种方程结合了前2种方程，避免了前2种方程出现的问题，具有更好的质量守恒性，但是在处理底部边界条件为自由排水的灌溉实践工程中出现了较大的质量守恒误差^[25]. 由于3个方程各有优劣，模拟土壤水流的程序使用Hao等^[25]提出的切换算法，它将式（8）、（10）结合起来，收敛性和质量守恒性更好，且应用范围更广.

以式（8）为例，通过有限体积法，将控制方程离散成如下形式:

(11) $\begin{split} \frac{\partial }{{\partial t}}\;&\int_{{V_{\text{p}}}} {\left[ {C\left( h \right) + {S_{\rm{s}}}{S_{\rm{w}}}} \right]} h{\rm{d}}V = \int_{{V_{\text{p}}}} {\nabla \cdot \left[ {{\boldsymbol{K}}\left( h \right) \nabla h} \right]} {\rm{d}}V + \\ &\int_{{V_{\text{p}}}} {\nabla \cdot \left[ {{\boldsymbol{K}}\left( h \right) \nabla z} \right]} {\rm{d}}V + \int_{{V_{\text{p}}}} {{Q_{\rm{s}}}} {\rm{d}}V . \end{split} $

式中： $ {V_{\text{p}}} $为控制单元的体积.

通过假定在 $ {V_{\text{p}}} $控制体内的变量为线性变化，来使方程计算时在空间离散上保持二阶精度. 用高斯迭代法将体积分方程转换成面积分方程，对含有 $ z $的项，将其合并到源项，对含有 $ h $的项，则将其转换成面积分形式：

(12) $ \int_{{V_{\text{p}}}} {\nabla \cdot \left[ {{\boldsymbol{K}}\left( h \right) \nabla h} \right]{\rm{d}}V = \sum\nolimits_{{f}} {\boldsymbol{S}} } \cdot {\left[ {{\boldsymbol{K}}\left( h \right) \nabla h} \right]_{{f}}} . $

式中： $ {\boldsymbol{S}} $为面积法向量，f表示相邻2个单元的公共面f.

对f面的定义如图3所示，该面是相邻2个单元的公共面，P、N分别为相邻单元的中心点，由P点指向N点构成向量 $ {\boldsymbol{d}} $. 对于正交的网格，向量 $ {\boldsymbol{d}} $和向量 ${\boldsymbol{S}}$是平行的，而对于非正交的向量， $ {\boldsymbol{d}} $和 $ {\boldsymbol{S}}$之间将会有一定夹角，如图4所示. 在计算式（12）右边项时须估算该面的渗透系数 $ \nabla h $，对于正交网格而言，可以直接列出如下等式：

(13) $ {\boldsymbol{S}} \cdot {\left[ {{\boldsymbol{K}}\left( h \right) \nabla h} \right]_{{f}}} = \left[ {{\boldsymbol{K}}{{\left( h \right)}_{{f}}} {\boldsymbol{S}}} \right] \cdot \frac{{\boldsymbol{S}}}{{\left| {\boldsymbol{S}} \right|}}\frac{{{h_{\text{P}}} - {h_{\text{N}}}}}{{\left| {\boldsymbol{d}} \right|}} . $

式中：表面的 $ {\boldsymbol{K}}{\left( h \right)_{{f}}} $是由相邻单元体中心点处 $ {\boldsymbol{K}}\left( h \right) $的值插值得到的； $ \left| {\boldsymbol{d}} \right| $、 $ \left| {\boldsymbol{S}} \right| $分别为向量 $ {\boldsymbol{d}} $、 $ {\boldsymbol{S}} $的模. 当网格不是正交时，为了确保空间离散的二阶精度，Jasak^[26]提出将 ${\boldsymbol{S}} \cdot {\left[ {{\boldsymbol{K}}\left( h \right) \nabla h} \right]_{{f}}}$项分成2部分（见图4），分别为 ${\boldsymbol{\varDelta }} \cdot {\left[ {{\boldsymbol{K}}\left( h \right) \nabla h} \right]_{{f}}}$、 ${\boldsymbol{k}} \cdot {\left[ {{\boldsymbol{K}}\left( h \right) \nabla h} \right]_{{f}}}$，第1项为正交网格影响的部分，第2项为非正交网格影响的部分. 因此有 ${\boldsymbol{S}} = {\boldsymbol{\varDelta }} + {\boldsymbol{k}}$， ${\boldsymbol{\varDelta }}$与 $ {\boldsymbol{d}} $的方向相同.

与时间相关的积分控制方程如下：

(14) $ \frac{\partial }{{\partial t}}\int_V {\phi h{\rm{d}}V} . $

式中： $ \phi $为 $ \dfrac{{\partial h}}{{\partial t}} $通用的系数， $ V $为控制体积.

式（14）的2种隐式离散方式表达式如下.

1）一阶精度的欧拉法：

(15) $ {\partial _t}\int_V {\phi h{\rm{d}}V = \frac{{{{\left( {\phi {h_{\text{P}}}V} \right)}^{k + 1}} - {{\left( {\phi {h_{\text{P}}}V} \right)}^k}}}{{\Delta t}}} . $

2）二阶精度的向后欧拉法：

(16) $ {\partial _t}\int_V {\phi h{\rm{d}}V = \frac{{3{{\left( {\phi {h_{\text{P}}}V} \right)}^{k + 1}} - 4{{\left( {\phi {h_{\text{P}}}V} \right)}^k} + {{\left( {\phi {h_{\text{P}}}V} \right)}^{k - 1}}}}{{8\Delta t}}} . $

对于瞬态问题，由于须处理空间项的倒数以及确定时间项精度，对式（14）进行时间上的积分：

(17) $ \int_t^{t + \Delta t} {\left[ {\int_V {\forall h{\rm{d}}V} } \right]} {\rm{d}}t = \int_t^{t + \Delta t} {{{\forall}^ * }h} {\rm{d}}t . $

式中： $ \forall h $表示含 $ h $的通用项， $ {\forall ^ * } $表示 $ h $对应的空间离散算子，对该式的时间积分有3种表示方法：

1）隐式欧拉法：

(18) $ \int_t^{t + \Delta t} {{\forall ^ * }h} {\rm{d}}t = {\forall ^ * }{h^{n + 1}}\Delta t . $

2）显式法：

(19) $ \int_t^{t + \Delta t} {{\forall ^ * }h} {\rm{d}}t = {\forall ^ * }{h^n}\Delta t . $

3）Crank Nicolson法:

(20) $ \int_t^{t + \Delta t} {{\forall ^ * }h} {\rm{d}}t = {\forall ^ * }\left( {\frac{{{h^n} + {h^{n + 1}}}}{2}} \right)\Delta t . $

方法1）和2）在时间上都是一阶精度，方法1）保证了 $ h $的有界性且是无条件稳定的，方法2）只有在库朗数大于1的情况下才能稳定，方法3）在时间上有二阶精度且无条件稳定，但是不能保证 $ h $的有界性. 该软件将各自的积分式子优点结合起来，通过时间步的动态控制实现了计算的准确性.

由于材料本构关系的非线性，对控制方程的直接求解是不可能的，使用古典的Celia等^[27]提出的迭代法，并结合式（10）可以得到：

(21) $ \begin{split} &\frac{{{\varphi ^{k + 1,m + 1}} - {\varphi ^k}}}{{\Delta t}} + {S_{\rm{s}}}S_{\rm{w}}^k \frac{{{h^{k + 1,m + 1}} - {h^k}}}{{\Delta t}} =\\ &\qquad \nabla \cdot \left[ {{\boldsymbol{K}}\left( {{h^k}} \right) \nabla \left( {{h^{k + 1,m + 1}} + z} \right)} \right] + Q_{\rm{s}}^{k + 1} . \end{split} $

式中：上标 $ m $表示迭代步.

${\varphi ^{k + 1,m + 1}}$可以表示为以 $ h $为变量的公式：

(22) $ {\varphi ^{k + 1,m + 1}} \approx {\varphi ^{k + 1,m}} + {\left( {\frac{{{\rm{d}}\varphi }}{{{\rm{d}}h}}} \right)^{k + 1,m}} \left( {{h^{k + 1,m + 1}} - {h^{k + 1,m}}} \right) . $

将式（22）代入到式（21），可以得到

(23) $ \begin{split} &\left( {{{\left( {\frac{{{\rm{d}}\varphi }}{{{\rm{d}}h}}} \right)}^{k + 1,m}} + {S_{\rm{s}}}S_{\rm{w}}^k} \right)\frac{{{h^{k + 1,m + 1}} - {h^k}}}{{\Delta t}} =\\ &\qquad \nabla \cdot \left[ {{\boldsymbol{K}}\left( {{h^k}} \right) \nabla \left( {{h^{k + 1,m + 1}} + z} \right)} \right] + S . \end{split} $

式中： $S = Q_{\rm{s}}^{k + 1} + \dfrac{{{\varphi ^k} - {\varphi ^{k + 1,m}}}}{{\Delta t}} + {S_{\rm{s}}}S_{\rm{w}}^k\dfrac{{{h^k} - {h^{k + 1,m}}}}{{\Delta t}}$.

通过将Navier-Stokes方程中的入渗项用Richards方程中计算的边界流通量表示以实现2个方程的耦合. 具体方法如图5所示. 图中， $ I \times {\rm{d}}t $为该时刻入渗量， $ T $为计算结束时间.

偏微分方程有3类边界条件，分别为Dirichlet边界、Neumann边界、Cauchy边界. 计算土壤水流的程序中有4种类型的边界，分别为常水头边界、常总水头边界、特定通量边界、内部位置固定值边界. 它们都可以是时间的函数. 其他边界的模拟都可以用这4种边界组合来实现. 常水头边界和常总水头边界可以分别表示为 $ h = {h_{\text{D}}}\left( {x,y,z,t} \right) $、 $ h + z = {h_{\text{D}}}\left( {x,y,z,t} \right) $. 其中 $ {h_{\text{D}}}\left( {x,y,z,t} \right) $为水头分布函数， $ x $、 $ y $、 $ z $为边界表面的中心点坐标. 2种边界都是 $ z $的函数，因此可以化简成 $ h = {h_{\text{G}}}\left( {x,y,z,t} \right) $，在有限体积法的计算框架中，边界面的水力梯度 ${\left( {\nabla h} \right)_{{f}}}$可以通过边界面上中心点处的水头同边界相邻单元的中心点的水头差与两者之间距离的比值表示，即 ${{\boldsymbol{S}}} \cdot {\left( {\nabla h} \right)_{{f}}} = \left| {{{\boldsymbol{S}}}} \right|\dfrac{{{h_{\text{G}}} - {h_{\text{p}}}}}{{\left| {\boldsymbol{d}} \right|}}$，其中 $ {h_{\text{G}}} $、 $ {h_{\text{p}}} $分别为边界面中心水头高度和边界面对应的单元中心点水头高度. 通量边界条件可以表示为

(24) $ - \frac{{{{\boldsymbol{S}}}}}{{\left| {{{\boldsymbol{S}}}} \right|}} \cdot \left[ {{{\boldsymbol{K}}_{\rm{s}}} {k_{\rm{r}}}\left( {\nabla h + \nabla z} \right)} \right] = {q_{\rm{c}}}\left( {x,y,z,t} \right) . $

式中： ${{\boldsymbol{S}} }$为面积向量， ${{\boldsymbol{K}}_{\rm{s}}}$为饱和渗透系数张量， ${q_{\rm{c}}}\left( {x,y,z,t} \right)$为边界通量， ${k_{\rm{r}}}$为相对渗透系数. 在降雨条件下，土体表面的边界条件并非是一成不变的，它将随着土体的入渗能力和降雨两者的相互作用而发生改变. 例如，当初始时刻土体中水的体积分数较低时，表面的入渗强度会较大，若入渗强度大于降雨强度，边界条件变为通量边界条件，而随着时间的推移，入渗强度会随着土壤中水的体积分数的增大而减弱，因此将会出现入渗强度小于降雨强度的现象，此时将产生地表径流，且水流将通过重力作用向地势低的地方流动，若在单元上的降雨量与通量变化达到动态平衡，则从该时刻开始，单元的水头高度将不再变化，此时瞬态边界将变为水头边界条件，当降雨减小或结束时，土壤入渗强度又将大于降雨强度，但由于入渗的时间效应，若地表水未全部入渗至土体内，地表边界条件将仍为水头边界，若地表水全部入渗，则水头边界条件将变成通量边界条件，当降雨最终为零时，边界条件将变成零通量边界条件. 因此，可以通过判断当前时间步单元边界的水头高度是否超过上一时刻土体入渗能力来决定边界条件的选取. 若有 $h_{{f}}^{k + 1} - I_{\text{s}}^k \times {\rm{d}}t > 0$（ $h_{{f}}^{k + 1}$表示 $ k + 1 $时刻的边界水头高度， $h_{{f}}^k$表示 $ k $时刻的边界水头高度，两者均可以由地表径流方程计算得出， $I_{\rm{s}}^k$表示 $ k $时刻土边界的入渗强度），表明土体入渗强度小于单元边界的水头高度，土体在该时间步下无法完全入渗表面积水，当前时间步在迭代开始时水头高度为恒定值 $h_{{f}}^{k + 1}$，式中 $h_{{f}}^{k + 1}$的值由程序计算得出， $ I_{\text{s}}^k $由程序在该时刻以水头高度 $h_{{f}}^k$为边界条件计算的入渗强度表示， $ {I}_{\mathrm{s}} $的具体计算公式如下：

(25) $ {I_{\text{s}}} = \frac{{{\boldsymbol{S}} \cdot {{\left[ {{\boldsymbol{K}}\left( h \right) \nabla h} \right]}_{{f}}}}}{{\left| {\boldsymbol{S}} \right|}} = \left[ {{\boldsymbol{K}}{{\left( h \right)}_{{f}}} {\boldsymbol{S}}} \right] \cdot \frac{{\boldsymbol{S}}}{{{{\left| {\boldsymbol{S}} \right|}^2}}}\frac{{{h_{\text{G}}} - {h_{\text{P}}}}}{{\left| {\boldsymbol{d}} \right|}} . $

若有 $h_{{f}}^{k + 1} - I_{\text{s}}^k \times {\rm{d}}t < 0$，表示该时间步的边界单元的水头高度小于土体的入渗能力，地表积水将全部入渗，此时边界条件将成为通量边界条件，边界入渗强度等于降雨强度.

采用一维土柱入渗试验来验证耦合模型的算法，土柱的高度为1 m，在土柱顶部的固定水头为0 m，土柱底部的固定水头为−8 m，水的饱和体积分数为0.363，残余水的体积分数为0.186，饱和渗透系数为1×10⁻⁶ m/s，土柱的初始水头均匀分布，设置为−7.848 m，VG模型参数 $ \alpha $=1 m⁻¹，n=1.53，入渗时间为13 h.

分别采用SEEP/W和本研究耦合模型计算t=6.5、13.0 h时刻土柱的水头分布情况，如图6所示，由结果可以看出，两者差别较小.

本例中土柱的高度也为1 m，土柱顶部的通量为2.78×10⁻⁷ m/s，VG模型参数值同2.1节，土柱的初始水头为静水条件，如图7所示，底部边界为不透水边界，入渗时间为12 h. 图中，x表示水平位置，y表示土柱高度，土条内部数字代表土柱水头高度. 采用SEEP/W和耦合模型计算的t=6、12 h时刻土柱的水头分布情况如图8所示. 可以看出，两者相差较小，基本一致，说明在通量边界条件下，耦合模型计算的结果是有效的.

计算结果表明，土壤水流和地表径流耦合模型在固定水头边界和通量边界条件下计算的结果是有效的.

非饱和土渗透模型采用Van Genuchten^[28]提出的模型，该模型下体积分数 $ \varphi $与水头 $ h $的关系如下：

(26) $ \varphi \left(h\right)=\left\{\begin{array}{l}{\varphi }_{{\rm{r}}}+\dfrac{{\varphi }_{{\rm{s}}}-{\varphi }_{{\rm{r}}}}{{\left[1+{\left(\alpha \left|h\right|\right)}^{n}\right]}^{m}},\quad \; h < 0;\hfill \\ {\varphi }_{{\rm{s}}},\qquad \qquad \qquad \qquad 其他. \end{array} \right.$

相对渗透系数 $ {k_{\rm{r}}} $与体积分数 $ \varphi $的关系如下：

(27) $ {k_{\rm{r}}}\left( \varphi \right) = {\left( {\frac{{\varphi - {\varphi _{\rm{r}}}}}{{{\varphi _{\rm{s}}} - {\varphi _{\rm{r}}}}}} \right)^{0.5}}{\left\{ {1 - {{\left[ {1 - {{\left( {\frac{{\varphi - {\varphi _{\rm{r}}}}}{{{\varphi _{\rm{s}}} - {\varphi _{\rm{r}}}}}} \right)}^{{1}/{m}}}} \right]}^m}} \right\}^2} . $

式中： $ \alpha $、 $ n $为VanGenuchten模型的2个参数， $m = 1 - {1}/{n}$， ${\varphi _{\text{s}}}$为饱和体积分数， ${\varphi _{\text{r}}}$为残余体积分数.

模型所采用的裂隙土为火山碎屑土^[29]，其VanGenuchten模型参数如下：α=1.8 m⁻¹，n=1.34， $ \varphi_{\rm{r}} $=0.1， $\varphi_{\rm{s}} $=0.7，k_s=9.8×10⁻⁶ m/s.

模型将模拟3.5 m高、坡度为26°的边坡，模型几何尺寸如图9所示.

模型假定坡体是均质的且入渗是各向同性的，初始时刻地表径流场的初始水头为零，时间步长和空间步长分别为0.05 h和0.5 m，对土坡的地表径流场的左边界采用不透水边界条件，对土坡的地表径流场的右边界采用自由排水边界. 土壤水流场中水的初始体积分数为0.4，侧面边界和底面边界则采用与初始水头值相同的水头边界，其值为−3.9 m，每0.05 h记录一次运算结果. 初始降雨强度为0.2 m/h，在半小时后降至0.1 m/h，持续2.5 h后降雨停止，直至5.0 h结束计算，降雨强度随时间的变化如图10所示.

如图9所示的A-A剖面处和B-B剖面处的水头随时间的变化如图11所示. 可以看出，在初始时刻，土的入渗能力较强，坡面积水较小，坡顶和坡底的水头较接近，坡面水的流动较小；随着地表入渗强度的降低，坡体表面积水深度开始增加，由于坡顶水的流动，坡底水出现从亚临界流到超临界流再到亚临界流的变化，坡底水出现先增大再减小又增大又减小的波动特征，坡顶水头变化较为单一；到降雨结束时，由于入渗和重力作用，坡顶和坡底积水开始迅速减少，坡顶水头和坡底水头大约在t=4.65 h时刻减小为零. 如图12所示为降雨结束后坡体表面水头的分布情况. 图中，H为总水头. 可以看出，在降雨结束时坡底部分位置水头并没有很快降低，这主要是由坡面水的流动对坡底水头的补充所导致的，且在降雨结束后地表径流大约持续了1.65 h.

为了表示坡面径流对土壤入渗水流的影响，如图13所示为坡顶处（A-A剖面）和坡脚处（B-B剖面）在时间t=3 h时的入渗雨水分布情况. 图中，l为入渗深度. 可以看出，坡脚处入渗深度大于坡顶的入渗深度，说明水的流动导致的坡面水头分布的不均对土体入渗的影响是存在的.

（1） 以边界水头作为自变量将地表径流和土壤水流方程耦合后，可以模拟地表水在坡面的流动和确定坡面水头的动态变化，以及相应的土壤中水分流动.

（2） 根据计算结果，在降雨停止后，由于地表水的径流，坡底水头仍会在一段时间内维持在一定高度.

（3） 对比2种工况下坡底的水头分布可知，地表径流引起地表水流动，造成地表水头在坡面上分布的不均匀性，对坡底和坡顶入渗的影响是不同，坡底的入渗量相对较高，说明考虑地表径流对边坡入渗的影响是有必要的.

须指出的是，本研究只考虑了均质各向同性坡体，对于坡体渗透各向异性、入渗土体分层情况还须进一步研究.