初始土体含水率对泥石流方量影响的数值模拟

谢鑫鹏,杨 龙,张 明,焦伟之,李成文,谢骏锦

(中国地质大学(武汉)工程学院，湖北 武汉 430074)

我国特殊的地质构造和复杂的地理环境，导致各种地质灾害频发，其中以泥石流地质灾害最为严重。泥石流地质灾害不仅会对生态环境造成破坏，而且严重威胁着人民的生命和财产安全[1]。泥石流是在降雨条件下，于沟谷中形成的携带大量泥砂、石块等固体物质的复杂混合流体[2-4]。目前，关于降雨触发泥石流地质灾害的研究较多。如：Wilson等[5]和Miller等[6]将降雨诱发的泥石流地质灾害归纳为两种类型，一种是持续时间较短且短时间内降雨强度较大诱发的泥石流，另一种是持续时间较长且长时间内降雨强度较小诱发的泥石流；Rudolph等[7]分析了前期土体含水率和降雨类型对微地貌变化的响应；Coe等[8]分析了泥石流启动前后土体的含水量变化特征，认为径流对冲沟中堆积物的冲蚀作用是泥石流启动的主要原因;Berti等[9]对泥石流沟道堆积物的渗透性进行了研究，发现降雨条件下沟道堆积物的渗透性越差土体越易达到饱和，并得出水动力条件是泥石流启动的主要原因之一;宋兵等[10]通过对鲁家沟泥石流进行详细的现场调查，指出鲁家沟泥石流是在暴雨作用下顺层斜坡失稳破坏，堰塞沟道并溃决形成的；毛硕等[11]研究认为野牛沟泥石流的启动主要是因为强降雨条件下导致的雨水快速汇集效应造成的；苏远东等[12]通过对热河沟泥石流的形成过程进行详细分析，得出热河沟泥石流的形成机理为“暴雨-溃决”型；王毅等[13]以毛家湾隧道口沟泥石流为研究对象，对该泥石流的运动形式与强降雨之间的关系进行了深入研究；王征强等[14]以后沟泥石流为研究对象，通过对其形成机制进行分析得出充足的水源是该泥石流形成的重要条件之一；杨敏[15]以小秦岭金矿区矿渣泥石流为例，通过对其形成的主要影响因素进行分析，得出该泥石流为暴雨型水石流，其启动方式主要包含“河谷起动型”和“沟道坡面滑塌—沟道堵塞—形成溃决”两种。

近年来，随着计算机技术的革新与发展，越来越多的学者采用数值模拟的方法来研究泥石流的启动过程和危险性。如：O’Brien等[16]通过对不同类型的泥石流进行分析研究，建立了FLO2D模型，为泥石流启动过程的数值模拟提供了依据；胡明鉴等[17]采用二维颗粒流程序(PFC2D)分析了降雨作用下松散碎屑物质启动形成泥石流的过程及其与土体含水率的关系；刘鑫磊等[18]为了预测溃决型泥石流冲出量特征值与威胁范围，基于FLO2D软件对不同降雨频率下牛圈沟泥石流的危险性进行了数值模拟研究；胡卸文等[19]利用基于有限体积法的CFX软件，对江口沟泥石流危险区范围和速度场分布情况进行了数值模拟研究；左正金等[20]利用PFC2D软件对栾川县柿树沟泥石流的失稳启动过程进行了数值模拟研究，结果表明强降雨是柿树沟沟内松散堆积体启动形成泥石流的触发因素。

目前的研究主要关注降雨对泥石流启动及运移过程的影响，而针对泥石流失稳方量与土体含水率之间关系的研究较少。因此，本文以新疆维吾尔自治区苏约克河流域内典型的泥石流沟——乔库而泥石流为研究对象，通过分析其工程地质条件，建立了离散元数值模型，并利用PFC2D软件对该泥石流物源区堆积物的失稳方量与土体含水率之间的关系进行了数值模拟研究，以为泥石流地质灾害的预警预报提供科学依据[21-22]。

1 乔库而泥石流概况

如图1所示，在新疆维吾尔自治区苏约克河流域内19 km长的河谷中，野外调查发现发育小型(V<2万方)、中型(2万方≤V<小于20万方)、大型(20万方≤V<50万方)和巨型(V≥50万方)泥石流沟共计39处，且该区段发育3条断层，岩性以砂岩、泥岩、玄武岩为主，沿断层两侧岩体较为破碎，经风化作用坡脚处堆积大量的风化物，为泥石流提供了丰富的物源，每逢雨季该流域内泥石流地质灾害频发。为了保证当地人民的生命和财产安全，亟待对该流域内泥石流失稳方量与土体含水率的关系进行研究。

本文选取苏约克河流域内典型的泥石流沟乔库而泥石流为研究对象，该泥石流沟位于苏约克河的西岸(见图1中的102号泥石流)，主要威胁对象为河流漫滩上的牧场。乔库而泥石流的岩性主要为第三系的红色砂泥岩和白垩系的黑色砂岩、砾岩，该泥石流区被一逆冲断层横切(见图1)，断层下盘为红色的砂泥岩，岩体表面较为破碎，受风化作用在其横穿的斜坡坡脚处往往堆积有小规模的细粒崩坡积层；上盘为黑色的砂岩、砾岩，在断层面附近岩体极为破碎，形成的崩坡积物是乔库而泥石流的主要物源。

乔库而泥石流沟的整个流域面积为4.42 km2，其物源区、流通区和堆积区如图2所示。其中，物源区主要为临近断层面的上盘、沟谷南侧破碎山体失稳堆积在坡脚处的滑坡和崩塌堆积物，堆积物的方量约为120万方；流通区分布在断层下盘的沟谷底部，长度约为1 500 m，两侧为红色砂泥岩组成的沟壁；堆积区面积约为0.11 km2，堆积物呈扇形堆积，平均厚度约为4 m，体积约为44万方。根据这三个区域得到的乔库而泥石流典型的剖面如图3所示。

图1 新疆维吾尔自治区苏约克河流域内泥石流分布图

图2 乔库而泥石流的平面图

图3 乔库而泥石流典型的剖面图

2 泥石流物源区堆积物失稳方量与土体含水率之间关系的数值模拟研究

本次针对乔库而泥石流物源区堆积物失稳方量与土体含水率之间关系的数值模拟研究的具体思路为：首先，根据遥感影像和现场调查，估算乔库尔泥石流的全部物源量以及威胁牧场和牧民的最小泥石流冲出量；然后，利用PFC离散元数值方法反演乔库而泥石流在相当于最小泥石流冲出量的土体失稳启动和物源区堆积物土体全部启动并正常运行至河谷两种工况下，土体所对应数值模型中颗粒的细观参数；随后，利用双轴压缩试验获取上述细观参数对应的土体宏观剪切强度参数；最后，结合直剪试验所得到的泥石流物源区堆积物土体剪切强度参数随土体含水率变化的曲线，得到土体相应的含水率。

乔库而泥石流物源区堆积物土体为黄褐色的粉质黏土，硬塑，结构较为密实，土体粒径较小，土体表面分布块石，块石粒径约5～10 cm。如图2所示，在乔库而泥石流物源区堆积物土体(B01)中取了1个样品，利用样品在含水率分别为13%(天然)、20%、27%和35%(饱和)条件下进行了直接剪切试验，样品每种含水率对应100 kPa、200 kPa、300 kPa和400 kPa 4种不同的法向荷载,最终获取的泥石流物源区土体(B01)剪切强度参数随含水率的变化曲线，见图4。

图4 乔库而泥石流物源区土体剪切强度参数随含水率 的变化曲线

2. 1 数值模型建立

根据如图3所示的乔库而泥石流典型剖面建立其数值模型。在数值模拟过程中，为了清晰地观察泥石流的启动和运移过程，将物源区的粉质黏土划分为4个部分，分别用颗粒填充。考虑到计算速度与结果可靠性之间的平衡，本次模拟采用的颗粒最小半径为0.25 m、最大半径为1 m，颗粒总数为16 785个，采用高斯分布。同时，为了更好地模拟乔库而泥石流物源区散体物质的特征，颗粒即球体之间的接触模型选用接触黏结模型。除了物源区的土体，其余的边界都用墙(Wall)来代替。在不影响结果可靠性的前提下为了提高数值计算效率，最终建立的乔库而泥石流数值模型，见图5。

图5 乔库而泥石流数值模型

2. 2 土体物理力学参数确定

本次数值模拟采用墙(Wall)和球(Ball)的计算参数，如表1所示。墙的参数包括法向刚度、切向刚度、摩擦系数，均取常用值，在数值模拟中采用固定参数。球体用接触黏结模型，模型中需要8个参数，但是多参数反演不仅会大大增加计算量而且会导致多解情况的加剧。为了尽量避免这一情况，首先利用工程地质类比法[23]确定密度、法向刚度、切向刚度、半径系数等对泥石流动力学过程影响较小的参数；然后，通过数值反演获取法向抗拉强度、切向抗拉强度、黏聚力和内摩擦角等对泥石流动力学过程影响较大的参数。

表1 乔库而泥石流物源区土体的细观参数

2. 3 模拟工况设置

根据遥感影像和现场调查，假设当泥石流堆积物约为1 m厚时，就会对牧场和牧民产生威胁；而乔库而泥石流堆积区的面积为11万方，因此估算当泥石流一次性冲出11万方的冲积物时，就会对泥石流扇前方和侧缘的牧场和牧民产生威胁。同时根据估算，乔库而泥石流物源区的全部物源为53万方。

本次数值模拟反演了两种工况：①工况1下，乔库而泥石流物源区至少有11万方物质失稳，这种工况下反演出的土体剪切强度参数和含水率是泥石流启动所需要的土体最低剪切强度参数和含水率；②工况2下，乔库而泥石流物源区的物质全部启动，形成的泥石流到达河谷，这种工况下反演出的土体剪切强度参数和含水率是泥石流全部失稳启动所需要的土体剪切强度参数和含水率。

2. 4 模拟结果分析

2.4.1 乔库而泥石流启动和运移过程反演结果分析

经过多次反演计算，当球体的法向抗拉强度、切向抗拉强度、黏聚力和内摩擦角分别等于20.4 kPa、22.7 kPa、20 kPa和10°时，模拟所得到的乔库而泥石流启动和运移过程的结果与工况1要求最为符合。图6为工况1下乔库而泥石流启动和运移过程的数值模拟结果，数值反演计算得到的该泥石流堆积体的分布特征如图6中150×104steps所示。

由图6可见，虽然该泥石流物源区整体失稳，但是最终大部分土体(深蓝色、绿色和红色部分)稳定下来，几何形态没有发生变化，泥石流前缘蓝色部分土体发生解体并有较大的位移，这部分土体的方量约为10.1万方，与估算的最小威胁泥石流的土体方量接近。

经过多次反演计算，当球体的法向抗拉强度、切向抗拉强度、黏聚力和内摩擦角分别为12.8 kPa、15.4 kPa、15 kPa 和8°时，模拟得到的乔库而泥石流启动过程的结果与工况2要求最为接近，即此时乔库而泥石流的流动距离能达到河谷。图7为工况2下乔库而泥石流启动和运移过程的数值模拟结果，直观地显示了该泥石流物源区土体从失稳、运移到动能耗尽并最终停积下来的全过程，数值反演计算得到的该泥石流堆积体的分布特征如图7中150×104steps所示。

图6 工况1下乔库而泥石流启动和运移过程的数值 模拟结果

图7 工况2下乔库而泥石流启动和运移过程的数值 模拟结果

由图7可见，由物源区运移到流通区和堆积区的土体方量约为48万方，该泥石流的最大运动距离约为2 450 m，且部分物源运移至苏约克河道内。

2.4.2 乔库尔泥石流物源区土体双轴压缩试验结果分析

继续利用上述工况1和工况2反演所得到的球体细观参数包括法向抗拉强度、切向抗拉强度、黏聚力、内摩擦角和表1中其他参数，通过建立双轴压缩试验数值模型(见图8)，对乔库而泥石流物源区土体进行双轴压缩试验，可得到两种工况下乔库而泥石流启动时物源区土体的剪切强度参数[24]，其试验结果见图9和图10。

图8 双轴压缩试验数值模型

图9 工况1下乔库而泥石流物源区土体双轴压缩试验结果

本次双轴压缩试验选用的颗粒尺寸和黏结模型(接触黏结模型)与上述泥石流数值模拟相同。在进行双轴压缩试验时，为了保证试样具有足够的颗粒并消除试样的尺寸效应，将双轴压缩试验的数值模型尺寸设定为240 m×480 m，颗粒的数量为8 920个(见图8)。根据泥石流数值模拟得到的细观参数，设置双轴压缩试验颗粒的细观参数，并分别在围压为50 kPa、100 kPa、150 kPa条件下对乔库尔泥石流物源区土体进行双轴压缩试验。

由图9(a)、(b)、(c)可知，通过采用工况1中球体的细观参数在围压分别为50 kPa、100 kPa、150 kPa条件下对乔库而泥石流物源区土体进行双轴压缩试验获取的土体应力-应变曲线，可得到该土体的峰值剪切强度分别为164 kPa、246 kPa、304 kPa；由图9(d)可知，根据采用3种围压及对应的峰值剪切强度绘制的土体摩尔破坏包络线，可得到乔库而泥石流物源区土体的黏聚力为19.2 kPa、内摩擦角为21.5°。

由室内直剪试验得到的乔库而泥石流物源区土体黏聚力、内摩擦角与土体含水率的关系曲线(见图4)可知：当土体黏聚力为19.2 kPa时，对应的土体含水率为33%时，当土体内摩擦角为21.5°时，对应的土体含水率约为31%，说明当物源区土体含水率为31%～33%时，即土体快达到饱和状态时(土体饱和含水率为35%)，乔库而泥石流物源区才约有10.1万方的土体失稳启动并有一定的运动位移。

图10 工况2下乔库而泥石流物源区土体双轴压缩试验结果

由图10(a)、(b)、(c)可知，根据采用工况2中球体的细观参数在围压分别为50 kPa、100 kPa、150 kPa条件下对乔库而泥石流物源区土体进行双轴压缩试验获取的土体应力-应变曲线，可得到该土体的峰值剪切强度分别为115 kPa、170 kPa、230 kPa；由图9(d)可知，根据采用3种围压及对应的峰值剪切强度绘制的土体摩尔破坏包络线，可得到乔库而泥石流物源区土体的黏聚力为12 kPa、内摩擦角为9°。

由室内直剪试验得到的乔库而泥石流物源区土体黏聚力、内摩擦角与土体含水率的关系曲线(见图4)可知：当土体黏聚力为12 kPa、内摩擦角为9°时，该泥石流物源区土体含水率已经超过其饱和含水率，说明乔库而泥石流物源区的堆积物要全部转化为泥石流并达到苏约克河河谷，土体含水率要超过其饱和含水率。

由于乔库而泥石流物源区的地形较高较陡，除了物源区本身接受的雨水外，其余汇水区汇集的雨水很难渗入到物源区土体中，加上研究区年降雨量少，可以推测乔库而泥石流物源区土体的含水率很难高于其饱和含水率。因此，本研究认为乔库而泥石流物源区土体在降雨作用下仅能失稳，如果没有其他因素导致这些土体在启动和运移过程中剪切强度降低，它们将很难形成远程运动的泥石流。故乔库而泥石流的运移机理是值得进一步深入研究的问题。

3 结 论

本文通过对新疆乌恰县托云乡乔库而泥石流分别处于最小泥石流冲出量和物源区土体全部失稳启动两种工况下土体的含水率特征进行了数值模拟分析，得到如下结论：

(1) 乔库而泥石流达到最小冲出量时(10.1万方)，物源区土体含水率为31%～33%，即土体快达到饱和状态(土体饱和含水率为35%)，此时物源区土体的黏聚力为19.2 kPa、内摩擦角为21.5°。

(2) 乔库而泥石流物源区的土体全部失稳启动时，物源区土体含水率需要超过其饱和含水率，此时物源区土体的黏聚力和内摩擦角分别为12 kPa和9°。

(3) 乔库而泥石流物源区土体在降雨作用下仅能失稳，如果没有其他因素导致这些土体在启动和运移过程中剪切强度降低，它们将很难形成远程运动的泥石流。因此，有必要进一步深入研究乔库而泥石流的运移机理。