中国西部山区地形地质条件复杂，地质灾害频繁发生。这些地质灾害不仅本身具有极大的危害性，而且其产生的松散岩土碎屑物堆积体在为后续降雨滑坡、泥石流等次生灾害的发育提供了充足的物源条件^[1]。相较于原坡面土体，这类松散堆积体一般属于宽级配弱固结土体，具有结构松散，孔隙率大，表面植被覆盖率低的特点；在降雨条件下，这类堆积土中细小颗粒极易随雨水渗流发生迁移，诱发斜坡失稳^[2-3]（如图1所示）。由这类松散堆积体降雨失稳引发的次生灾害影响范围广、持续时间长，成为灾区恢复重建过程中的重大隐患。

针对松散堆积体中的细颗粒运移现象，国内外研究者开展了大量研究，通过野外调查^[3-4]、现场人工降雨试验^[5]和室内水槽降雨模型试验^[6-8]描述了降雨条件下堆积土边坡中细颗粒伴随雨水入渗由坡顶和坡面向坡脚和土体深处运移的趋势，定性指明了细颗粒运移架空粗颗粒骨架，改变土体水力特性，诱发坡体局部失稳的作用机制。相对来说，针对堆积土边坡降雨条件下细颗粒运移效应的定量化研究较少。张磊^[9]、Zhang^[10]、Cividini^[11]等提出的饱和土渗流潜蚀模型，分析了降雨入渗侵蚀细颗粒引发的土体孔隙率及渗透性增大效应对边坡稳定性的影响；Lei等^[12]进一步提出了可考虑细颗粒运移–沉积堵塞效应的渗流潜蚀模型，分析了细颗粒沉积堵塞引发的坡面浅层破坏机理。

堆积体边坡降雨失稳是一个包含雨水非饱和入渗、细颗粒运移、土体水力特性动态演化等多物理现象的降雨–渗流潜蚀耦合瞬态过程。总体来说，国内外研究者已在宏观上充分认识到了松散堆积土边坡中的细颗粒运移现象及其对坡体稳定性的潜在作用效应，但缺乏针对堆积土降雨–渗流潜蚀耦合过程的深入分析，针对细颗粒运移对边坡降雨稳定性影响的定量化研究依然比较薄弱。既有研究主要考虑了细颗粒运移引发土体渗透性演化对降雨入渗及边坡稳定性的影响，尚没有涉及潜蚀对土体持水性的影响效应研究。为此，作者对既往研究中相关试验现象进行概化分析，建立能考虑堆积土水土特征曲线、力学性质演化的降雨–渗流潜蚀模型及无限边坡稳定性分析方法，旨在定量分析细颗粒迁移引发的堆积土渗透性、持水性及强度演化对降雨非饱和入渗过程及边坡稳定性的影响规律。

1 物理模型构建 1.1 松散堆积土边坡细颗粒迁移现象概化

松散堆积土一般由粗细不均的各粒级土石颗粒混杂堆积而成，其中的粗大颗粒经常相互接触构成土体的骨架结构，而黏粒（<0.005 mm）、粉粒（0.005～0.075 mm）等细小土颗粒则充填于粗大颗粒构成的大孔隙间^[2]。这类土体级配宽，固结程度低，其细小颗粒极易在水流作用下发生迁移^[2]。本文所概化的模型如图2（a）所示。根据既往研究，图2中在降雨作用下，松散堆积土中细小颗粒在壤中流溶蚀、拖拽作用下脱离土体基质，并伴随雨水入渗在粗粒间的大孔隙中由坡顶和坡面向坡脚和土体深处运移^[9,13]。细颗粒在松散堆积土坡体中随雨水入渗向坡内的运移，使坡面浅表层土体的细颗粒含量不断减少，增大浅表层堆积土的孔隙率，使其土体颗粒逐渐粗化^[2,6]。此外，细颗粒在运移过程中，可能在其渗流路径发生沉积，堵塞土体内部孔隙，堵塞土体内部孔隙，导致局部土体渗透性降低^[6,14]。在这种细颗粒随雨水侵蚀–运移–沉积过程的反复作用下，堆积体入渗面以下的若干分米处往往会逐渐形成一层大致平行于坡面的相对不透水层，阻碍雨水的继续入渗，加速浅层土体基质吸力的消散，可能继续发展形成滞水并产生地表径流，诱发浅层失稳破坏^[5-7,13]。

细颗粒迁移的侵蚀–运移–沉积堵塞机制与单个细颗粒在粗颗粒骨架间孔隙内的迁移过程密切相关（图2（b））。目前，堆积土边坡领域尚没有针对孔隙尺度细颗粒迁移机制的研究报道。在这方面，水利工程领域土石坝潜蚀破坏^[15]及石油工程领域储层破坏^[16-17]的相关研究值得借鉴。参考Civan^[16],堆积土体中孔隙可抽象为一系列连通的孔室；各孔隙体通过孔喉（孔径狭小处）相互连通。如图2（b）所示，当流体流经土体时，在胶体力或水动力的作用下孔壁表层的细颗粒（状态1）逐渐被释放或脱离出来^[16]；细颗粒一旦被释放，将作为孔隙水的悬浮物（状态2）随入渗水流运移；在其运移过程中，细颗粒悬浮物可能在孔隙中孔喉、孔壁等其他位置重新沉积或截留（状态3）或运移出孔隙^[17]。细颗粒从某一孔壁处启动又在另一孔壁处发生沉积通常对孔隙介质的渗透性影响较小。但细颗粒在孔喉处的沉积、截留将严重损伤孔隙介质的渗透性^[18]：细颗粒在孔喉的截留及桥接有可能引发孔喉的关闭，从而阻止水流的通过；更甚的是，细颗粒一旦在孔喉发生截留，后续流经的悬浮液中的细颗粒将不断在该处积累，类似于“雪球效应”，细颗粒不断被过滤截留最终形成宏观坡体尺度的不透水薄层（图2（a））。

1.2 三相多组分孔隙介质体积分数

自然状态下的松散堆积土可视为一典型的非饱和多孔介质：堆积土基质构成孔隙介质的固相土骨架，其孔隙则由液体和气体共同填充。对于易遭受渗流潜蚀的堆积土，其土体基质可假定由粗颗粒及细颗粒两大类颗粒组成。在满足一定的水力条件及几何条件下，细颗粒可能从土体骨架脱离成为液体悬浮物，并随水流在粗颗粒间的孔隙中运移。为对这一系列过程进行数学描述，将非饱和堆积土概化为三相（固 ${\rm{S}}$ 、液 ${\rm{L}}$ 、气 ${\rm{G}}$ ）多组份（粗颗粒 ${\rm{c}}$ 、细颗粒 ${\rm{e}}$ 、水 ${\rm{w}}$ 、空气 ${\rm{g}}$ ）孔隙介质，并由下列集合表示：

${\rm{S}} = \left\{ {{\rm{c}}',{\rm{e}}'} \right\};{\rm{L}} = \left\{ {{\rm{w}},{\rm{e}}} \right\};{\rm{G}} = \left\{ {\rm{g}} \right\}$

(1)

其中，加上标的组分归属于固相，而没有上标的组分归属于流体（液、气）相。细颗粒的侵蚀及沉积过程，可视为固相细颗粒组分 ${\rm{e}}'$ 和液相细颗粒组分 ${\rm{e}}$ 相互交换质量的相变过程。

如图3所示，在初始时刻 ${t_0}$ ，初始构形中一初始体积为 ${\rm{d}}{\varOmega _0}$ 的三相多孔介质被用作非饱和堆积土的代表体积单元。在任意时刻 $t$ ，经过雨水入渗后该单元体积变为 ${\rm{d}}{\varOmega _{\rm{t}}}$ 。其土骨架变形可由固相构形坐标转换的雅可比行列式表示 $J={\rm{d}}{\varOmega }_{{\rm{t}}}/{\rm{d}}{\varOmega }_{0}$ 。根据混合物理论，代表体积单元的总体积为其所包含的各组分体积之和，也等于其所包含的各相体积之和（ ${\rm{d}}{\varOmega _{\rm{t}}} = \displaystyle\sum\limits_{j} {{\rm{d}}\varOmega _{\rm{t}}^{j}} = \displaystyle\sum\limits_{\alpha} {{\rm{d}}\varOmega _{\rm{t}}^{\alpha} }$ ， ${j}$ 、 ${\alpha}$ 分别为任意组分及相）^[19]。任意组分及相的拉格朗日形式体积分数可表示为 ${\phi }_{{j}}={\rm{d}}{\varOmega }_{{\rm{t}}}^{{j}}/ {\rm{d}}{\varOmega }_{0}$ 和 ${\phi _{\alpha} } = {{{\rm{d}}\varOmega _{\rm{t}}^{\alpha} } / {{\rm{d}}{\varOmega _0}}}$ ；土力学中的孔隙率 $\phi $ 在此可表示为所有孔隙空间的拉格朗日体积分数^[20]；而液相饱和度 ${S_{\rm{L}}}$ 为液体所占体积相对整个空隙体积的体积分数；液相细颗粒体积浓度 ${c_{\rm{e}}}$ 为液相细颗粒组分在液体中的体积分数。由此，单位多孔介质体积 ${\rm{d}}{\varOmega _0}$ 内液相体积分数 ${\phi _{\rm{L}}}$ 及细颗粒组分体积分数 ${\phi _{\rm{e}}}$ 可表示为：

${\phi _{\rm{L}}} = {S_{\rm{L}}}\phi ;{\phi _{\rm{e}}} = {c_{\rm{e}}}{S_{\rm{L}}}\phi $

(2)

在雨水入渗过程中，饱和度 ${S_{\rm{L}}}$ 和孔隙率 $\phi $ 均随时间变化。在任意时刻 $ t$ ，液相及气相所占据的体积分别从其初始时刻 $ t_0 $ 的初始体积 $ {S_{{\rm{L}}0}}{\phi _0}{\rm{d}}{\varOmega _0} $ 及 $\left( {1 - {S_{{\rm{L}}0}}} \right){\phi _0}{\rm{d}}{\varOmega _0} $ 变化为 ${S_{\rm{L}}}\phi {\rm{d}}{\varOmega _0} $ 及 $ \left( {1 - {S_{\rm{L}}}} \right)\phi {\rm{d}}{\varOmega _0} $ 。与其对应，土体固相颗粒在 $t_0 $ 及 $ t$ 时刻可分别表示为 $\left( {1 - \phi } \right){\rm{d}}{\varOmega _0} $ 及 $\left( {J - \phi } \right){\rm{d}}{\varOmega _0} $ ^[21]。假设在雨水入渗过程中，由细颗粒迁移直接引发的孔隙体积变化为 ${\phi _{\rm{er}}}{\rm{d}}{\varOmega _0}$ ，则由图3可得 $ t $ 时刻的拉格朗日孔隙率 $ \phi $ ：

$\phi = {\phi _0} + (J - 1) + {\phi _{\rm{er}}}$

(3)

式中： ${\phi _0}$ 为初始孔隙率； $J$ 为土体骨架变形转换的雅可比行列式， $(J - 1)$ 代表土体骨架变形对孔隙率的贡献； ${\phi _{\rm{er}}}$ 为固相细颗粒质量分数 ${m_{{\rm{e}}'}}$ 变化对土体孔隙率的贡献（ ${{\partial {\phi _{\rm{er}}}}/ {\partial t}} = - {{\partial {m_{{\rm{e}}'}}} / {\left( {{\rho _{{\rm{e}}'}}\partial t} \right)}}$ ，其中 ${\rho _{{\rm{e}}'}}$ 为固相细颗粒的内在密度），其演化由细颗粒潜蚀本构模型控制。式（3）指出，孔隙率的变化是由土骨架变形和细颗粒潜蚀共同作用造成的。只针对性地研究细颗粒的潜蚀效应，将假定雨水入渗过程中（失稳前）土体骨架变形对孔隙率的贡献可以忽略不计（ $J \approx 1$ ）。

1.3 质量守恒法则

基于连续介质力学，对于发生渗流潜蚀的土体，其骨架内占据同一空间区域的任一组分 $j$ 的质量变化 ${{{{\rm{d}}}^{\rm{S}} m_{{j}}} / {{\rm{d}}t}}$ （ ${{\rm{d}}}^{\rm{S}}/{\rm{d}}t$ 为相对于固相的物质导数的运算符^[20]）事实上是两种物理现象造成的：1）由细颗粒潜蚀（侵蚀及沉积）所引发的组分质量变化 ${{\partial {m_{j}}} / {\partial t}}$ ；2）液体流经土体骨架所在空间携带引入的该组分质量变化 $ \nabla \cdot {{{I}}^{j}} $ （ $ {{{I}}^{j}} $ 为组分 ${j}$ 的质量通量^[19]）。在小变形条件下，其连续性方程可表示为 ${{{{\rm{d}}}^{\rm{S}}{m_{j}}}}/ {{{\rm{d}}t}} + \nabla \cdot {{{J}}^{j}} = {{\partial {m_{j}}} / {\partial t}}$ ^[19]。为简化模型，固、液相中所有组分均假定具有内在不可压缩性（内在密度 ${\;\rho _{j}}$ 为常数），其连续性方程可进一步简化为 ${{{{\rm{d}}}^{\rm{S}}{\phi _{j}}} / {{\rm{d}}t}} + \nabla \cdot {{{J}}^{j}} = {{\partial {m_{j}}} / {\left( {{\rho _{j}}\partial t} \right)}}$ （ $ {{{J}}^{j}} $ 为组分 ${j}$ 的体积通量^[19]）。对于本文中的堆积土而言，只有液相悬浮细颗粒组分 $ {\rm{e}} $ 可与其对应的固相可侵蚀细颗粒组分 ${\rm{e}}' $ 相互交换质量（ ${\dot m_{{\rm{e}}'}} = - {\dot m_{\rm{e}}} $ ）。假定液相悬浮细颗粒 ${\rm{e}} $ 与固相可侵蚀细颗粒 $ {\rm{e}}' $ 具有相同的内在密度（ ${\;\rho _{\rm{e}}} = {\rho _{{\rm{e}}'}}$ ），在小变形条件下，液相 ${\rm{L}}$ 及液相细颗粒组分 ${\rm{e}}$ 的质量守恒方程可表示为：

$\frac{{\rm{d}}^{\rm{S}} {\phi _{\rm{L}}}}{{{\rm{d}}t}} + \nabla \cdot {{{J}}^{\rm{L}}} = - \frac{1}{{{\rho _{{\rm{e}}'}}}}\frac{{\partial {m_{{\rm{e}}'}}}}{{\partial t}}$

(4)

$\frac{{\rm{d}}^{\rm{S}} {\phi _{\rm{e}}}}{{{\rm{d}}t}} + \nabla \cdot {{{J}}^{\rm{e}}} = - \frac{1}{{{\rho _{{\rm{e}}'}}}}\frac{{\partial {m_{{\rm{e}}'}}}}{{\partial t}}$

(5)

式中， $ {{{J}}^{\rm{L}}} $ 为液相整体的体积通量 ${{{J}}^{\rm{L}}} = \displaystyle\sum\limits_{j \in L} {{{{J}}^{j}}}$ ， $ {{{J}}^{\rm{e}}} $ 为液相细颗粒组分的体积通量。本文忽略液相细颗粒的扩散运移效应，只考虑随水流平流运移；细颗粒相的体积通量可计算为 ${{{J}}^{\rm{e}}} = {c_{\rm{e}}}{{{J}}^{\rm{L}}} $ ^[19]。

结合式（2）中液相体积分数 ${\phi _{\rm{L}}} = {S_{\rm{L}}}\phi $ ，液相 ${\rm{L}}$ 的质量守恒方程式（4）可最终展开为：

$\frac{{\left( {1 - {S_{\rm{L}}}} \right)}}{{{\rho _{{\rm{e}}'}}}}\frac{\partial {m_{{\rm{e}}'}}}{{\partial t}} + \phi \frac{{\partial {S_{\rm{L}}}}}{{\partial {p_{\rm{L}}}}}\frac{{\partial {p_{\rm{L}}}}}{{\partial t}} + \nabla \cdot {{{J}}^{\rm{L}}} = 0$

(6)

式中， $p_{\rm{L}} $ 为液相压强，饱和度 $ S_{\rm{L}} $ 假定为一个只依赖于基质吸力 $ p_{\rm{c}} $ 的函数（忽略水土特征曲线的滞回效应），并假定气体压强与大气压保持平衡取为0，故 $ {p_{\rm{c}}} = - {p_{\rm{L}}} $ 。

将式（2）中液相细颗粒体积分数 $ {\phi _{\rm{e}}} = {c_{\rm{e}}}{S_{\rm{L}}}\phi $ 及液相细颗粒相的体积通量 $ {{{J}}^{\rm{e}}} = {c_{\rm{e}}}{{{J}}^{\rm{L}}} $ 代入液相细颗粒组分 $ {\rm{e}} $ 的质量守恒方程（5），并借助于连续性方程（6），可得出液相细颗粒组分的质量守恒方程：

$\frac{{1 - {c_{\rm{e}}}}}{{{\rho _{{\rm{e}}'}}}}\frac{\partial {m_{{\rm{e}}'}}}{{\partial t}} + {S_{\rm{L}}}\phi \frac{\partial {c_{\rm{e}}}}{{\partial t}} + {{{J}}^{\rm{L}}} \cdot \nabla {c_{\rm{e}}} = 0$

(7)

1.4 潜蚀及非饱和渗流法则

细颗粒迁移诱发的渗透系数演化是影响边坡降雨稳定性最主要的原因之一^[9]。由于细颗粒孔喉截留对渗透性的影响远大于细颗粒孔壁沉积^[16]，参考Liu等^[18]将孔壁视为提供固相可侵蚀细颗粒的主要源区，而孔喉则被认为是液相悬浮细颗粒发生截留重新变为固相沉积细颗粒的场所（图2（b））。鉴于此，将固相细颗粒质量分数 ${m_{{\rm{e}}'}}$ 视为由孔壁处可侵蚀细颗粒 ${m_{\rm{er}}}$ 及孔喉处沉积细颗粒 ${m_{\rm{de}}}$ 两部分组成，式（6）及（7）中源项可表示为 ${{\partial {m_{{\rm{e}}'}}} / {\partial t}} = {{\partial {m_{\rm{er}}}} / {\partial t}} + {{\partial {m_{\rm{de}}}} / {\partial t}}$ 。

1.4.1 细颗粒侵蚀法则

文中细颗粒的侵蚀假定完全由雨水入渗造成的水动力机制所引发。在孔隙尺度，土骨架中细颗粒的侵蚀与孔隙中流体所施加的剪切力、土体的抗侵蚀能力及孔壁表层可供侵蚀的细颗粒含量有关^[15]。在此，采用Cividini等^[11]所提出的侵蚀法则描述非饱和堆积土中细颗粒的侵蚀过程：

$\left\{ \begin{array}{l} \dfrac{\partial }{{\partial t}}\dfrac{{{m_{\rm{er}}}}}{{{\rho _{{\rm{e}}'}}}} = - {k_{\rm{e}}}\left( {{x_{{\rm{e}}'}} - {x_{{\rm{e}}'\infty }}} \right)\left\| {{{{v}}^{\rm{LS}}}} \right\|; \\ {x_{{\rm{e}}'\infty }} = {x_{{\rm{e}}'0}} - \left( {{x_{{\rm{e}}'0}} - x_{{\rm{e}}'\infty }^ * } \right)\dfrac{{{{{v}}^{\rm{LS}}}}}{{{{{v}}^ * }}},\;{{\rm{if}}}\;\;\;{0 \le {{{v}}^{\rm{LS}}} \le } {{{v}}^ * } ; \\ {x_{{\rm{e}}'\infty }} = x_{{\rm{e}}'\infty }^ * - {\alpha _{\rm{er}}}\lg \left( {\dfrac{{{{{v}}^{\rm{LS}}}}}{{{{{v}}^ * }}}} \right),\;{{\rm{if}}}\;\;\;{{{{v}}^{\rm{LS}}} > {{{v}}^ * }} \end{array}\right. $

(8)

式中： ${k_{\rm{e}}}$ 为一非负材料常数，反映土体对流水动力的敏感性； $\left\| {{{{v}}^{\rm{LS}}}} \right\|$ 为水流相对于土骨架的渗透速度 ${{{v}}^{\rm{LS}}}$ （ ${{{v}}^{\rm{LS}}} = {{{{{J}}^{\rm{L}}}} / \phi }$ ）的范数； ${x_{{\rm{e}}'}}$ 为当前土骨架中可供侵蚀的细颗粒的质量分数（ ${x_{{\rm{e}}'}} = {m_{{\rm{e}}'}}/{m_{{\rm{S}}0}}$ ）； ${x_{{\rm{e}}'}}_\infty $ 为当前流速下土骨架最终剩余细颗粒含量； ${{{v}}^ * }$ 为细颗粒发生线性侵蚀时孔隙水的参考临界流速； $x_{{\rm{e}}'\infty }^ * $ 为参考流速下土骨架最终剩余细颗粒含量； ${x_{{\rm{e}}'0}}$ 为侵蚀前土体初始细颗粒含量（ ${x_{{\rm{e}}'0}} = {m_{{\rm{e}}'0}}/{m_{{\rm{S}}0}}$ ）； ${\alpha _{\rm{er}}}$ 为控制 ${x_{{\rm{e}}'\infty }}$ 随流速 ${{{v}}^{\rm{LS}}}$ 增加而减少速率的侵蚀参数^[9]。

1.4.2 细颗粒沉积法则

被侵蚀后的细颗粒作为液相悬浮物会随水流一起运移，随后可能在孔壁或孔喉处发生沉积。对于孔壁沉积，Civan等^[16]指出当水流小于某一临界速度时，只有颗粒孔壁沉积现象，而大于某一流速时则颗粒沉积与携带现象会同时发生。细颗粒在孔壁处的沉积对土体水力特性的影响主要表现为轻微减小其过水通道，对整个土体水力特性影响较小；在此，不单独对孔壁处发生的沉积效应进行讨论，而是假定其效应可通过调节孔壁侵蚀过程予以考虑。采用Schaufler等^[22]提出的细颗粒沉积法则，其假定孔喉处被拦截的概率随着单位时间内流经的细颗粒数量成正比：

$\frac{\partial }{{\partial t}}\frac{{{m_{\rm{de}}}}}{{{\rho _{{\rm{e}}'}}}} = {k_{\rm{d}}}{c_{\rm{e}}}\left\| {{{{v}}^{\rm{LS}}}} \right\|$

(9)

式中： $ k_{\rm{d}} $ 为一非负材料常数，用以衡量孔喉对悬浮细颗粒的截留能力； $c_{\rm{e}} $ 为液相悬浮细颗粒的体积浓度； $ \left\| {{{{v}}^{\rm{LS}}}} \right\| $ 为水流相对于土骨架的渗透速度的范数。

1.4.3 水土特征曲线

式（6）中饱和度 $ S_{\rm{L}} $ 及其导数 ${{\partial {S_{\rm{L}}}} / {\partial {p_{\rm{L}}}}}$ 须根据水土特征曲线函数计算。对于发生渗流潜蚀的非饱和堆积土而言，其水土特征曲线不仅为基质吸力的函数，而且可能受细颗粒迁移引发的孔隙率、细颗粒含量演化影响。为确定细颗粒含量对堆积土水土特征曲线的影响，本研究对采自都江堰某泥石流源区的松散堆积土进行了筛分重塑，配制了细颗粒含量（粒径 <0.075 mm）分别为20%、15%和10%的3种重塑土，其级配曲线如 图4（a）所示。这3种土体被制成干密度为1.70、1.62和1.53 g/m³的试样，其孔隙率分别测定为0.34，0.38及0.41（以期反演土体在渗流潜蚀不同阶段的土体水土特征曲线演化情况）。随后，通过压力膜仪测定的这3组土样的土体水土特征关系如图4（b）所示，随着细颗粒减少，同一基质吸力下堆积土所能维持的含水量下降（持水性下降），与田湖南^[23]、Kim等^[24]的试验结论基本相符。

由式（3）可得 ${{\partial \phi } / {\partial t}} = {{\partial {\phi _{\rm{er}}}} /{\partial t}}$ ，即细颗粒的侵蚀损失将直接增大土体的孔隙率。图4（b）中也可解释为随着孔隙率增大，堆积土的持水性下降。此趋势与Gallipoli等^[25]所报道的孔隙率对水土特征曲线的影响规律较为一致。据此，本文采用孔隙率作为水土特征曲线中状态变量，间接反映细颗粒潜蚀对土体持水性的影响。采用van Genuchten模型为水土特征曲线基础模型：

${S_{\rm{L}}} = {\left[ {1 + {{\left( {{\alpha _{\rm{v}}}{p_{\rm{c}}}} \right)}^{{n_{\rm{v}}}}}} \right]^{ - \left( {1 - 1/{n_{\rm{v}}}} \right)}}$

(10)

式中， ${\alpha _{\rm{v}}}$ 和 ${n_{\rm{v}}}$ 为拟合参数， ${p_{\rm{c}}}$ 为基质吸力（ ${p_{\rm{c}}} = - {p_{\rm{L}}}$ ）。假定雨水入渗过程引发的土骨架变形很小可以忽略，孔隙率的演化将只受细颗粒含量演化控制；土体从初始状态经历不同的侵蚀阶段将对应不同的孔隙率。按照Gallipoli等^[25]提出的孔隙率相关模型计算非饱和土进气值参数 ${\alpha _{\rm{v}}}$ ，考虑渗流潜蚀过程中孔隙率演化对水土特征曲线的影响。

${\alpha _{\rm{v}}} = {\alpha _{\rm{e}}}{\left( {\frac{\phi }{{1 - \phi }}} \right)^\xi }$

(11)

式中， $\alpha _{\rm{e}} $ 和 $ \xi $ 为拟合参数。可见，非饱和土进气值将随孔隙率 $ \phi $ 增大（细颗粒侵蚀）而减小。根据式（10）和（11）对图4（b）中试验数据进行最小二乘法拟合（拟合参数 ${\alpha _{\rm{e}}} $ =0.95， $ \xi $ =4， $n_{\rm{v}} $ =1.8），其拟合结果也较为理想。

1.4.4 非饱和渗流法则

式（6）及（7）中液相流体的体积通量 $ {{J}}^{\rm{L}} $ 可通过达西定律予以计算^[20]：

${{{J}}^{\rm{L}}} = - \frac{{{K_{\rm{h}}}}}{{{\gamma _{\rm{L}}}}}\left( {\nabla {p_{\rm{L}}} - {\gamma _{\rm{L}}}} \right)$

(12)

式中， ${\gamma _{\rm{L}}}$ 为孔隙水的重度， ${K_{\rm{h}}}$ 为渗透系数。假定 ${K_{\rm{h}}} = {K_{{\rm{h}}0}}{K_{\rm{r}}}$ ，其中 ${K_{{\rm{h}}0}}$ 为土体初始饱和渗透系数， ${K_{\rm{r}}}$ 为渗透系数比率，用以考虑土体的孔隙率、饱和度、孔喉处细颗粒沉积等多因素影响，按如下形式计算：

${K_{\rm{r}}} = \exp \left( { - {\alpha _{\rm{v}}}{p_{\rm{c}}}} \right)\frac{{{{\left[ {(1 - {k_{\rm{r}}}{m_{\rm{de}}})\phi /{\phi _0}} \right]}^{{n_{\rm{r}}}}}}}{{1 + 2.5{c_{\rm{e}}}}}$

(13)

非饱和度的影响通过Gardner方程用基质吸力 ${p_{\rm{c}}}$ 考虑，其中 ${\alpha _{\rm{v}}}$ 由式（11）计算而得；孔隙率、细颗粒沉积堵塞及孔隙水粘度效应则参照Liu等^[18]，在Carman–Kozeny方程基础上进行了拓展，其中 ${n_{\rm{r}}}$ 为取值范围为2.5至3.5的经验指数， ${k_{\rm{r}}}$ 为有关细颗粒沉积堵塞效率的经验参数。

2 渗流潜蚀耦合过程模拟及边坡稳定性评价

第1节中液相质量守恒方程（6）、液相细颗粒组分质量守恒方程（7）、细颗粒侵蚀方程（8）及细颗粒沉积方程（9）将共同组成非饱和土渗流潜蚀耦合过程数学描述的控制性偏微分方程组（式（6）及式（7）中 ${{\partial {m_{{\rm{e}}'}}} / {\partial t}} = {{\partial {m_{\rm{er}}}} / {\partial t}} + {{\partial {m_{\rm{de}}}} / {\partial t}}$ ）；所选取的主要未知变量为：液体压强 $p_{\rm{L}} $ ，液相细颗粒体积浓度 $ c_{\rm{e}} $ ，孔壁处细颗粒质量分数 $m_{\rm{er}} $ 及孔喉处细颗粒质量分数 $ m_{\rm{de}} $ ；其中 $ \phi $ 根据初始孔隙率及固相细颗粒质量分数由式（3）计算； $ S_{\rm{L}} $ ， ${{\partial {S_{\rm{L}}}} / {\partial {p_{\rm{L}}}}}$ ， $ {{{J}}^{\rm{L}}}$ 分别根据 $ p_{\rm{L}} $ 由水土特征曲线函数式（10）、达西定律式（12）计算。这4个控制性偏微分方程被植入多场耦合隐式有限元开源程序Bil（ http://perso.lcpc.fr/dangla.patrick/bil/），以求解非饱和土渗流潜蚀过程的初始边界值问题。

降雨诱发滑坡一般为浅层滑坡，在均质土坡内滑裂面一般平行于坡面，可采用无限边坡稳定性分析方法计算边坡稳定性安全系数^[9]。文中，堆积土边坡中的土体被概化为均质的三相多孔介质。在此，将图2（a）所示的非饱和堆积土边坡简化为如图5所示的水位线固定的无限边坡，并采用Lu等^[26]提出的无限边坡模型对其进行降雨稳定性分析：

$\begin{aligned}[b] F_{\rm{s}}\left( {{\textit{z}},t} \right) =& \frac{{\tan \;\varphi }}{{\tan\; \beta }} + \frac{{2c}}{{\gamma {H_{\rm{s}}}\left( {\textit{z}} \right)\sin 2\beta }}- \\ & \frac{{S_{\rm{L}}^{}\left( {{\textit{z}},t} \right){p_{\rm{L}}}\left( {{\textit{z}},t} \right)}}{{\gamma {H_{\rm{s}}}\left( {\textit{z}} \right)}}\left( {\tan\; \beta + \cot \;\beta } \right)\tan\; \varphi \end{aligned} $

(14)

式中， $\;\beta $ 为边坡坡度， $c$ 为非饱和土体有效黏度， $\varphi $ 为其有效摩擦角， $\gamma $ 为水位线上土体的平均容重， ${H_{\rm{s}}}$ 为假定滑面距坡面的垂直深度（即稳定性评估位置）， $H_{\rm{w}}$ 为水位线深度。

堆积土中的渗流潜蚀过程也可能对其力学性质造成损伤。为此，在进行稳定性分析时，还将参考Wan等^[27]通过固相体积分数 $\left( {1 - \phi } \right)$ 的演化间接考虑细颗粒损失效应，以 ${{\left( {1 - \phi } \right)} / {\left( {1 - {\phi _0}} \right)}}c$ 及 ${{\left( {1 - \phi } \right)} / {\left( {1 - {\phi _0}} \right)}}\varphi $ 替换式(14)的土体有效黏度 $c$ 及有效摩擦角 $\varphi $ ，对比分析渗流潜蚀引发的力学性质损伤对边坡稳定性的影响。

模拟所需渗流潜蚀耦合模型及无限边坡模型的参数及赋值如下：控制细颗粒侵蚀参数（ ${k_{\rm{er}}}$ ， ${x_{{\rm{e}}'0}}$ ， $x_{{\rm{e}}'\infty }^ * $ ， ${{{v}}^ * }$ ， ${\alpha _{\rm{er}}}$ ），取为0.4、0.2、0.19、5.0×10^–6 m/s、0.1；沉积参数（ ${k_{\rm{de}}}$ ， ${\rho _{{\rm{e}}'}}$ ），取为0.2、2 600 kg/m³；水土特征曲线参数（ ${\alpha _{\rm{e}}}$ ， $\xi $ ， ${n_{\rm{v}}}$ ），取为0.95、4.00、1.80；渗透性参数（ ${\phi _0}$ ， ${K_{{\rm{h}}0}}$ ， ${k_{\rm{r}}}$ ， ${n_{\rm{r}}}$ ），取为0.34、7.5×10^–5 m/s、0.05、3；无限边坡模型参数（ $\;\beta $ ， $c$ ， $\varphi $ ， $\gamma $ ），取为25°、2 kPa、25°、24 kN/m³。除水土特征曲线参数外，其余参数参考了Lei等^[12]基于东川蒋家沟泥石流源区土体的参数取值，作为细颗粒潜蚀前（细颗粒含量为20%时）的土体材料属性参考值。

2.1 1维土柱非饱和入渗模拟

假定无限边坡模型中水位线位于6 m深处，并采用上述非饱和土渗流潜蚀有限元程序模拟6 m深一维非饱和堆积土柱在降雨入渗作用下的渗流潜蚀过程。本模拟假定土体骨架在降雨入渗过程中不发生变形，所采用的初始及边界条件为：土柱底部与地下水位线平齐，孔隙水压恒为0；土柱中初始孔隙水压 ${p_{\rm{L}}}$ 假定在重力作用下沿深度方向按 $ - {\gamma _{\rm{L}}}{\textit{z}}$ 线性分布；土柱中土体初始固相细颗粒质量分数 ${m_{\rm{er}}}$ 为340 kg/m³；初始沉积细颗粒质量分数 ${m_{\rm{de}}}$ 及初始液相细颗粒体积浓度 ${c_{\rm{e}}}$ 均为0；土柱顶部边界液相细颗粒体积浓度 ${c_{\rm{e}}}$ 保持为0，且施加强度为90 mm/h的降雨入渗通量 $q$ ，模拟清水入渗。

相较于Lei等^[12]，本模拟将着重于细颗粒迁移引发的渗透性及持水性演化对渗流潜蚀耦合过程的影响。在此，假定了3种工况：工况R（参照工况），降雨入渗不引发堆积土中细颗粒迁移（ $ {k_{\rm{er}}} = {k_{\rm{de}}} = 0$ ）；工况ED（只考虑渗透性演化），降雨入渗引发堆积土中细颗粒侵蚀，并在下游发生沉积，但水土特征曲线保持恒定；工况EDV（考虑渗透性及持水性演化），发生细颗粒侵蚀及沉积，且采用孔隙率相关的水土特征曲线。下面将对这些工况进行对比分析。

降雨入渗通量 $q$ 的持续作用下，非饱和堆积土柱内固相细颗粒质量分数 ${m_{{\rm{e}}'}}$ 、液相细颗粒体积浓度 ${c_{\rm{e}}}$ 在上述3种工况（R，ED，EDV）条件下不同时刻（0、1、3、5 h）的分布曲线如图6所示。相对于工况R（不发生潜蚀， ${m_{{\rm{e}}'}}$ 及 ${c_{\rm{e}}}$ 恒定为其初始值），工况ED及EDV中，雨水从土柱顶部的入渗触发了细颗粒迁移过程。随着降雨持时的增加，堆积土基质中的固相细颗粒逐渐在雨水渗流力作用下被侵蚀流失（图6（a）中 ${m_{{\rm{e}}'}}$ 逐渐减小）；与饱和土柱渗流不同，非饱和土柱中的水流渗透速率是随着浸润过程而不断增大的，从而在土柱不同高程处造成了不同程度的侵蚀。与此同时，侵蚀掉的细颗粒转化为液相悬浮细颗粒并伴随水流运移（图6（b）中 ${c_{\rm{e}}}$ 持续增加）；随着浸润峰的不断下移，雨水渗流路径上越来越多的细颗粒被雨水侵蚀携带，使得雨水中细颗粒浓度峰值不断增加。对比工况ED及EDV可见，由于考虑细颗粒侵蚀对土体持水性的减小，相同时刻工况EDV中更多的水份将入渗至深层土体，由此其深部土体潜蚀程度也更严重。

事实上，由于工况ED及EDV考虑了细颗粒的再沉积过程，液体中部分悬浮细颗粒会在孔喉处被拦截过滤，重新转化为固相土骨架的一部分（图7（a）中，工况R无沉积）。根据式（9），其沉积率正比于悬浮细颗粒通量；当液体中细颗粒浓度增大时，细颗粒在孔喉处被拦截的概率也逐渐增大。由于坡体表面的细颗粒浓度被固定为0以模拟清水入渗，由此坡体表层细颗粒的沉积率也恒为0；而土柱深处，随着细颗粒悬浮液的前移，固相沉积细颗粒含量也随时间逐渐增加。图6中，土柱内部沉积细颗粒的质量分数峰值在不同时刻主要集中于距入渗坡面0.5～1.0 m深处，即形成了图2（a）坡体概化模型中的细颗粒密集沉积层。此现象可根据本文模型得以较好地解释：在此峰值以上区域，由于清水的不断掺入，液相细颗粒浓度总是保持在较低值，越靠近坡面其细颗粒沉积率越低；在此峰值以下区域，夹杂着悬浮细颗粒的雨水浸润峰抵达此处的时间相对更迟，且流入该区域的液相细颗粒通量将随着沉积堵塞过程发展逐渐减小，其细颗粒沉积率也将随着深度增加而逐渐降低。两者综合作用在坡体浅层动态平衡，最终形成了细颗粒沉积量的峰值。

图7（b）给出了渗透系数比率 $ K_{\rm{r}} $ 的时空演化曲线。细颗粒迁移对渗透系数的影响是两方面的：一方面，细颗粒侵蚀过程同步增大了土体孔隙率（式（3）），从而增大土体渗透性；另一方面，细颗粒在孔喉处的沉积堵塞将减小土体渗透性。根据本文所采用的渗透系数演化方程（式（13）），在雨水入渗过程中，除饱和度增加非饱和渗透性外，这两种机制同时存在、相互竞争，三者共同决定土体的渗透系数比率演化。从图6（a）及图7（a）可看出，工况ED和EDV条件下上部区域（0～0.5 m深）中细颗粒侵蚀最严重，但沉积率相对较小，由此该区域孔隙率将显著增大；较之没发生渗流潜蚀的工况R，该区域将出现明显的土体粗化现象。由图7（b）可见，工况ED和EDV中细颗粒密集沉积层（图7（a）沉积细颗粒质量分数峰值）以上区域土体的渗透性显著增大，加快了坡面雨水入渗；较深处的细颗粒密集沉积层附近，悬浮细颗粒的密集沉积堵塞则使得该区域渗透性较之工况R急剧减小，形成一层相对不透水层，阻碍雨水的继续向下入渗。

比较图8（a）中饱和度演化曲线可见，由于相对不透水层的形成，工况ED及EDV中土柱上部的雨水下渗过程较之工况R明显减缓；当相对不透水层上面土体中雨水下渗速率小于其降雨入渗速率时，雨水便逐渐积累使得土体快速饱和。经过5 h的降雨，工况ED及EDV中相对不透水层上部土体已经完全饱和，而其下部土体还远未达到饱和状态。该数值模拟结果与王志兵等^[7]的水槽试验中观察到的现象完全一致。与之对应，工况ED及EDV中不透水层上部土体中的孔隙水压也快速增加（基质吸力快速减小），并在5 h连续降雨后产生了超孔隙水压（图8（b），说明地表径流已经形成）。当考虑细颗粒的侵蚀对土体持水性的减弱效应时（工况EDV），更多雨水不能被上部土体持有而下渗到下部土体，从而加速了不透水层的形成（相对于工况ED）：在5 h降雨后，图7（b）中工况EDV不透水层下部土体渗透性较之工况ED中渗透性显著减小；与此对应，图8（a）中工况EDV不透水层下部土体的饱和度由于更缺乏雨水的持续入渗补充，也明显比工况ED中低。且由于细颗粒的侵蚀降低土体的进气值，较之工况ED，相同的饱和度值在工况EDV中将对应更低的基质吸力（图8（b）中更高的孔隙水压）；由此，工况EDV较工况ED更容易形成超孔隙水压（图8（b））。

2.2 基于无限边坡模型的边坡稳定性分析

对该1维土柱进行降雨入渗有限元分析，并将所求解所得的“有效基质应力” $({S_{\rm{L}}}{p_{\rm{L}}}) $ 代入无限边坡模型式（14），计算出各工况降雨过程中无限边坡在不同时刻的稳定性安全系数。图9给出了各工况在降雨开始时及降雨持时5 h后的安全系数分布曲线。可见，在持续5 h降雨后，工况R中边坡浅层安全系数F_s远小于1，一直处于稳定状态。而考虑堆积土降雨–渗流潜蚀效应的工况ED及EDV中，边坡坡面浅层会形成相对不透水层，导致该层上部土体生成超孔隙水压（图8（b）），继而显著降低该区域的安全系数，诱发边坡浅层失稳。在持续5 h降雨后，图9中工况EDV在坡面以下0.25 m处F_s小于1，形成了初始失稳滑面，且其滑面深度与王志兵等^[7]及陈晓清^[5]针对东川蒋家沟泥石流源区相近坡角的堆积土边坡的试验结果较为一致。根据图8中工况ED及工况EDV相同的饱和度及孔隙水压演化形态可知，由式（14）计算出的两工况安全系数也将具有相同的演化形态；但较之工况ED，工况EDV中细颗粒潜蚀对土体持水性的减弱加速了其降雨–渗流潜蚀过程，进而加速了其边坡失稳过程。考虑力学性质损伤后各工况安全系数的分布曲线如图9中点线所示。可见，当考虑渗流潜蚀对土体力学性质损伤影响后，工况ED及EDV中由于固相细颗粒的流失（图6（a）），坡体上部土体强度随降雨入渗潜蚀过程不断减弱，加速坡体的浅层失稳。

3 结　论

对降雨作用下非饱和堆积土边坡中坡体及孔隙尺寸的细颗粒迁移现象进行了概化分析，基于多孔介质混合物理论，构建了可描述非饱和堆积土中细颗粒侵蚀–运移–沉积全过程的渗流潜蚀模型，并结合无限边坡模型提出了可定量考虑降雨–渗流潜蚀过程对土体渗透性、持水性及强度影响的堆积土边坡降雨稳定性分析方法。利用该模型对非饱和堆积土边坡降雨入渗过程的有限元模拟及稳定性分析，可得出如下结论：

1）结合考虑潜蚀效应的非饱和土渗透系数演化函数，模型能定量化地解释松散堆积土边坡中细颗粒迁移形成相对不透水层，诱发坡体浅层失稳现象的内在机理。

2）细颗粒侵蚀流失将减弱堆积土的持水性，加速浅层坡体雨水向下入渗，进而反馈加速坡体内渗流潜蚀进程。

3）细颗粒侵蚀流失减弱堆积土的抗剪强度，加速边坡浅层失稳。

4）降雨–渗流潜蚀过程引发的土体渗透性、持水性及强度演化均对堆积土边坡稳定性造成不利影响，由此对松散堆积体边坡进行降雨稳定性分析，应充分考虑其内在的细颗粒迁移效应。