泥（石）流是植被发育不良山区常见的自然灾害，是土体在水力侵蚀作用下十分活跃的现代动力学地貌过程^［1-2］。沙漠流域的泥流灾害则大多为风水复合侵蚀作用的产物。季节交替的风蚀、水蚀过程使得大量风沙土堆积物充填于沟道，并由上游暴雨洪水裹挟起动，形成典型的沟谷型沙漠泥流，乃至演变为泥流灾害，威胁着流域生态系统稳定，一直以来都是干旱、半干旱区土壤侵蚀及水土流失防控研究的热点^［3-6］。

近年来，以泥沙搬运堆积历史过程反演^［7-9］、野外观测与遥感分析^{［3-6，10-12］}等为主要手段的相关研究，针对以中国北方干旱区和黄河上游“十大孔兑”区域等为代表的沙漠流域在风水交互作用下的泥流灾害过程等方面，取得了诸多有益的成果。Ta等^［5］基于黄河上游“十大孔兑”高含沙洪水（泥流）过程原型观测数据，通过理论分析，首次明确了该区域风水复合侵蚀诱发泥流灾害的环境因子等，进而从气候、地形等宏观方面揭示了其发生条件。然而，有关其微观角度的发育条件、形成机制及灾害防控等研究，却鲜有报道。

部分学者以山区泥石流为研究对象，采用模型试验等手段^［13-15］，通过改变试验土体物质组成、沟床坡度、模拟降雨强度等，分析降雨量变化下泥沙颗粒体的稳定性以及初始运动^［16-19］等，一定程度上揭示了其发生发展的动态机制^［20-21］，对于研究沙漠泥流微观发生机理而言，具有方法上的可借鉴性。但山区泥石流体中的细颗粒成分以黄土为主^［22］，且常夹杂有大石块；而干旱区沙漠沟道中的地表物质组成则为黄土与沙漠粗沙掺混后的堆积物，这种堆积体结构上不但具备黄土颗粒的细观连续性，又具备沙漠粗沙颗粒的突变性。这种非均质且不连续的特征使得流域沟道内的泥流堆积体在水流作用下的失稳破坏模式和运动机理较普通的堆积体有所差异，即常规的堆积体破坏模式很难直接解释沟谷型沙漠泥流堆积体的运动。

鉴于此，本文基于团队前期研究基础，以黄河上游宁蒙河段典型风水复合侵蚀流域为研究区域（图1），重点针对堆积体堵塞沟道极端条件下的沟谷型沙漠泥流，从堆积体起动这一成灾的首要环节入手，通过模型试验和数值模拟等方法，概化模拟上游暴雨洪水入渗条件下堆积体内孔隙水压力变化特性及其起动过程，以期从微观角度指征其成灾机理，进而为灾害预警、防控奠定理论基础。该区域风蚀水蚀交替，风沙水沙过程剧烈，发育典型的沙漠宽谷，分布着中国乃至世界上地表风水复合侵蚀强度最大的区域^{［6，23-24］}，对本研究而言，具有不可替代的地域优势。

黄河上游宁蒙段上起下河沿，下至头道拐，全长1 080 km，是沙漠、河流交互演化过程的典型区域。该河段沿岸山洪沟众多，均为狭长型，且大多与黄河垂直相接，沟道内固体物质组成自西向东由黄土向沙漠粗沙过渡，风蚀、水蚀作用的主导地位亦有所差异。风水复合侵蚀诱发的高含沙洪水乃至泥流灾害入黄，严重影响着黄河健康，尤其以“十大孔兑”最为显著，风水复合侵蚀强度比例可达1∶1^{［6，25-26］}。

本文选取宁蒙河段最西端的中卫段作为主要研究区域。该区域黄河北岸为腾格里沙漠，黄河南岸山洪沟道地势南高北低，固体物质组成以黄土为主，掺混有部分沙漠粗沙。采样分析发现，沙漠粗沙颗粒含量为0—20%，泥流灾害以水蚀作用为主导。冬春两季风力作用搬运沙漠粗沙至沟道底部，与沟道两侧岸坡重力侵蚀滑塌产生的黄土混为一体，或形成新的堆积体，或将历史洪水残留的泥流堆积体加高培厚，形成足以堵塞沟道的沙漠泥流堆积体。雨季来临时，上游暴雨洪水裹挟泥流堆积体起动，逐渐演变为高含沙洪水乃至泥流灾害。

研究区域中黄河南岸四条山洪沟道相对比较顺直。经实地量测，沟道内残留的泥流堆积体形状虽不规则，但总体为梯形；故结合前人工作经验^［27-28］，最终概化为梯形堆积体开展模型试验；在四条山洪沟道中采集堆积体土样，取样点位置分别为37°27′00″N、105°00′40″E，37°28′25″N、105°09′40″E，37°28′49″N、105°13′43″E，37°29′44″N、105°16′47″E，主要在中游及沟口处；在北岸腾格里沙漠中（37°28′55″N、104°57′28″E）采集沙漠粗沙试样；分析样本主要物理性能指标，作为本文试验和数值模拟中5种工况的设计依据。

为精准控制边界条件，便于试验量测起见，将研究区域的顺直沟道概化为直槽，设计制作可变坡度的有机玻璃直槽概化模型试验装置（长2 000 mm、宽300 mm、高300 mm，图2）；自制梯形模具，作为堆积体成型的侧限，用以制备顶宽50 mm、高210 mm、上下游坡比均为1∶1.4、顺水流底部长度650 mm的试验堆积体。堆积体初始质量含水率经由前期测验设定为满足堆积体成型要求的最小值1.3%，用以还原研究区域气候干旱、沟道堆积体初始含水量小的真实物理图景。试验前将土体分3次填入模具，每层厚70 mm，使用木槌夯实，直至达到设计干表观密度（环刀法测定）后停止，拆模并整平处理后开始试验。前期预试验发现，上游来水条件下，堆积体上游坡脚处首先发生入渗，但最终出现位移变形的区域位于下游坡脚处，因此，在试验堆积体内上中下游分别布设了5个DMKY系列应变式孔隙水压力传感器，并配套相应的DM-YB1808动静态电阻应变仪进行孔隙水压力的实时监测。

根据研究区域的土体物质组成将试验工况设为不同质量比例下黄土与沙漠粗沙的混掺物，共5组：沙漠粗沙0掺混、5%掺混、10%掺混、15%掺混、20%掺混。据此开展上游来水条件下的堆积体内孔压特性及起动过程研究。

本文采用基于有限元分析软件ABAQUS 2020中Mohr-Coulomb本构模型^［29］和有效应力原理的堆积体流固耦合模型^［30］，其中，反映土体应力场变化导致破坏的剪切屈服面函数为：

式中：φ为q-p应力面上Mohr-coulomb屈服面的倾斜角，称为材料的摩擦角，0°≤φ≤90°；c为材料黏聚力，kPa；R_mc（Θ，φ）控制屈服面在π平面上的形状，由式（2）进行计算。

式中：Θ为极偏角；cos(3Q)=r3/q3，r是第三偏应力不变量J₃。

为避免屈服面尖角导致塑性流动方向不唯一的问题，采用连续光滑的椭圆作为塑性势面，计算式如下：

式中：ψ为剪胀角；c|₀为初始黏聚力；ε为子午面上的偏心率。

式中：e为π平面上的偏心率，控制π平面Θ=0°—60°的塑性势面形状。

反映土体中流固耦合效应的有效应力原理可表达为：

式中：σijs为固体颗粒之间的应力；p_w为孔隙水压力；ϕ为孔隙度。

网格划分方案采取关键部位加密划分方案，并为了避免积分运算时的“沙漏现象”，选取流固耦合分析中有着较强收敛性的孔压单元CPE6MP。模拟工况分别为：①改变相应于模型试验中5种混掺比例堆积体的密度和渗透性等条件，提取特征点A处的孔隙水压力变化，结合模型试验对所得结果进行分析验证；②设置5种水位上升速率，以堆积体上游侧水位蓄至堆积体高度的95%为准，分别为10、20、50、80、110 min，用来间接表征降雨强度的变化，通过特征点A处的孔隙水压力变化反映堆积体起动中的孔隙水压力变化特性。

为简便计算，本文将水位上升过程中的孔隙水压力近似地看作不同水位作用下的静水压力来计算，可将水位与时间的关系转换为孔隙水压力与时间的关系，再结合试验中的孔隙水压力实测数据进行拟合分析，以验证其正确性。

就水位与时间的关系而言，首先假设堆积体迎水面封闭（不发生渗流），结合上游来水量等试验工况，对得到的数据进行拟合，拟合结果为：

式中：h为上游的蓄水位高度，m；t为时间，s。

其次，真实试验现象表明，堆积体前的槽蓄水流近似于三角形变化区域S（图3）。根据堆积体前的汇水面积的几何关系可得：

式中：S为上游蓄水截面面积，m²；b为上游蓄水区域长度，m；Q为上游蓄水量，m³；q为水流流量，m³·s^-1；L为试槽宽度，即堆积体垂直于水流方向上的尺寸，m；α为试槽与水平面的夹角，即沟道的坡度；θ为堆积体坡脚的角度。

结合式（7）—（9）可得堆积体前水位与时间变化关系为：

将试验的各项物理参数代入式（10）可得：

对比式（6）和式（11）发现，堆积体前所积水位与时间呈幂函数关系。但考虑到实际堆积体入渗时需要消耗部分水量，时间一定的情况下，式（11）得出的水位较小，因此，本文中水位与时间的变化关系可按式（11）表达。

结合式（11）可得孔隙水压力与时间变化关系为：

式中：P为孔隙水压力，kPa；γ_ω为水的容重，kN·m^-3。

根据基本假设与试验现象建立堆积体起动中孔隙水压力随时间变化模型（图3）。前期预试验表明，上游蓄水位到达限定水位之前，即t₀—t_i 的蓄水阶段，湿润锋从上游坡脚处（1^#传感器附近）开始，在堆积体内呈圆弧状推进至堆积体中部（2^#传感器附近），该阶段孔隙水压力变化最为显著，因此，对5种试验工况下渗流过程中1^#、2^#孔隙水压力传感器实测变化曲线进行了拟合分析，并对流固耦合数学模型进行验证。

两个测点处在5种工况下的孔隙水压力均随时间呈现出先增大，后趋于稳定的变化特性（图4）。拟合结果表明，其变化规律与时间呈指数为0.5的幂函数模型关系，模型系数0.041—0.088，且决定系数R²大于0.709，吻合良好。

1^#测点处的孔隙水压力量测值普遍大于2^#测点处，这是由于两个测点处的孔隙水压力传感器埋设高度不同，静水压力、水土压力、内部应力等存在差异。

随着堆积体中粗沙颗粒的增加，渗流过程中孔隙连通性增强，通气孔隙度增大，故较之纯黄土堆积体而言，孔隙水压力达到峰值后率先开始产生下降趋势，而黄土颗粒间因为黏聚力强，孔隙水压力达到峰值后尚有一定的维持稳定的能力。

为验证本文所采用的流固耦合模型，开展与模型试验堆积体同尺寸、同工况的数值模拟，并选取上游侧1^#传感器附近的A点作为特征点进行孔隙水压力变化特性分析（图5）。

孔隙水压力变化特性仍与时间呈现为指数为0.5的幂函数关系，决定系数R²均大于0.843，拟合程度良好，进一步验证了堆积体起动过程中孔隙水压力随时间的幂函数关系模型。再经与模型试验结果对比，模拟值与试验值的拟合结果定量吻合，进而验证了本文所建立的堆积体流固耦合模型的正确性和适用性。

由于受野外泥流灾害突发性强、泥流过程危险度高、观测不便、边界条件复杂、参数量测困难等因素限制，尚无法开展真实沟道中的泥流堆积体起动试验，故根据野外量测结果，将本文试验条件下的泥流堆积体尺寸放大10倍，并基于本文采用的流固耦合模型，分析起动过程中孔隙水压力变化特性。为进一步还原降雨因素对堆积体起动的影响，模拟以堆积体前水位上升速率间接表征上游降雨强度的变化，5种工况如上文所述。提取各工况下A特征点处的孔隙水压力变化曲线，并对曲线进行指数为0.5的幂函数拟合，拟合结果如图6所示。

由图6中的结果可知，孔隙水压力随时间仍呈现出指数为0.5的幂函数递增关系，决定系数R²均在0.929以上，拟合程度良好。对比5种工况下的孔隙水压力变化过程发现，上游降雨强度大、产汇流过程迅速，堆积体前蓄水位上升速率快，初期水流渗透不够充分，因而孔隙水压力值较小；随着堆积体前水位壅高时间的持续，堆积体前静水压力持续增大，孔隙水压力有增大的趋势，堆积体发生起动破坏的概率变大。

沟谷型沙漠泥流的起动是堆积体内渗流作用下孔隙水压力变化导致土体饱和液化，进而破坏失效的复杂非线性过程。王如宾等^［31］研究了降雨入渗对堆积体变形失稳破坏的影响，认为孔隙水压力是重要的诱发因素，且通过模型试验发现孔隙水压力随降雨强度和时间的增长而增大。除此之外，部分学者^［32-34］对堆积体的失稳破坏做了大量的研究，均指出孔隙水压力随湿润锋到达后会首先呈现出增大的趋势。本文研究结果与上述学者观点一致，但进一步建立了孔隙水压力增大过程中随时间变化的幂函数模型。

由传感器实测的实时孔隙水压力变化曲线波动剧烈，真实反映了土体应力场与水流渗流场之间复杂的动态耦合响应过程（图4）。图4—6中拟合的幂函数模型曲线系数变化范围较大（0.041—0.832），笔者初步分析认为，该系数或许综合反映了野外真实沟道中的堆积体型态尺寸、混掺均匀度、复杂地形、水流流态等复杂条件。下一步通过系数拟合优化后的幂函数模型将更加精确地指示堆积体中孔隙水压力的变化特性，本文不再赘述。

左自波等^［35］通过对多种堆积体土样进行模型试验探究了堆积体的滑坡失稳规律，认为土体颗粒级配对土体破坏模式有显著影响，并指出堆积体湿润锋到达后其孔隙水压力随之增加，达到峰值后稳定，失稳后立即降低。

图7为本文试验过程中1^#、2^#孔隙水压力传感器记录的全过程变化曲线，结果同样表明，孔隙水压力随时间先增大，直至达到渗流饱和，孔隙水压力达到峰值后趋于稳定，继而迅速减小，反映了堆积体屈服失效后的瞬间起动过程。

分析其原因，笔者认为，随着孔隙水压力的增大，堆积体内孔隙连通性增强、孔隙度动态变化、土体应力场相应变化，经历了应力场→体积应变→孔隙率→孔隙水压力→渗流场→塑性应变→屈服破坏的系列复杂变化过程。

团队前期对黄河上游宁蒙段泥流过程的野外观测数据^{［5-6，26］}表明，洪峰涨退最剧烈的时段平均持续约50 min，该时段也是堆积体前蓄水的关键时期，因此本文选取了50 min蓄满和20%沙漠粗沙混掺堆积体，深入开展上游水位上升过程中堆积体的塑性变形及位移数值模拟研究，以期分析其起动机制。

计算终止时，孔隙水压力的最大值出现在堆积体上游坡脚处，最小值则位于坡顶处（图8）。由于水流与土体之间的耦合作用，形成了如图8B所示的堆积体塑性应变区，其最危险区域位于下游坡脚处，并不断向上延伸，直至变形发展到一定程度后，形成塑性贯通区，最终表现为如图8C所示的堆积体下游侧边坡首先出现脱离整体向下游滑动变形的趋势，继而堆积体整体起动，这与试验现象高度一致，进而揭示了流固耦合作用下沟谷型沙漠泥流堆积体的起动模式与机制。

本文通过对沟谷型沙漠泥流堆积体起动中的孔隙水压力变化特性试验分析，并采用基于Mohr-Coulomb本构模型和有效应力原理的流固耦合模型，揭示了孔隙水压力变化作用下堆积体起动过程的细观机理，建立了水位上升时堆积体起动过程中孔隙水压力与时间二者之间指数为0.5的幂函数模型，揭示了流固耦合作用下堆积体内应力场→体积应变→孔隙率→孔隙水压力→渗流场→塑性应变→屈服破坏的复杂非线性变化过程和起动机理。

本文从起动这一成灾的首要环节入手，为沟谷型沙漠泥流预警防控提供了基于孔隙水压力变化监测的新方法和新途径，为未来结合信息化手段进行泥流灾害的防治奠定了坚实的理论基础。