基于Massflow的金沙江白格滑坡数值模拟及堵江风险预测

曹水合，吴新明，钟东

(四川省地质矿产勘查开发局九一五水文地质工程地质队，眉山 620020)

白格滑坡位于西藏自治区江达县波罗乡白格行政村北面，金沙江右岸。该滑坡分别于2018年10月10日(简称“10.10”滑坡)和2018年11月3日(简称“11.03”滑坡)间隔24 d先后发生两次特大型滑坡，并两次形成堰塞体阻断金沙江，形成巨大的回水库容。其溃坝式泄流影响范围之大、波及范围之广，造成的经济损失巨大、社会影响深远，以及整个过程中应急响应，已经超出一般意义的地质灾害威胁和危害对象统计范畴，在国家地质灾害防治史上具有里程碑的意义。

白格滑坡经两次滑动后，松散物质已大为减少，但滑坡体对山体挠动作用巨大(主要表现为拉拽、刮、铲、刷)，滑坡三面由于临空卸荷作用不断加强，变形迹象十分明显，后缘及两侧发育多条深大裂缝，仍存在大量潜在不稳定块体(残留体)。残留体状态为滑坡，周界明显，发生了明显的下滑、错动，仍具备进一步下滑的趋势。因此，查明残留体基本特征、进行稳定性分析及堵江风险预测意义重大。

根据相关文件，四川省地矿局九一五队承担并完成了金沙江白格滑坡应急勘查任务。本文作者全面参与了该滑坡后期应急勘查及残留体应急处置，完成了滑坡区无人机航空摄影测量、工程地质测绘、钻探、物探等详细勘查工作，基本查明了滑坡残留体特征，并进行了稳定性分析评价，为本次滑坡残留体失稳堵江过程数值模拟提供了重要依据。

目前，利用数值模拟软件对高位滑坡堵江风险定量评价取得了一定进展。欧阳朝军[1]、范宣梅[2]等采用Massflow软件对白格“10.10”和“11.03”两次滑坡运动过程进行了模拟分析，并与两次滑坡实际堆积特征进行了对比分析，取得了较好的效果；周礼等[3]利用PFC软件进行了白格滑坡运动过程特征数值模拟与危险性预测研究，揭示了两次滑坡失稳、运动、铲刮、碰撞、停积等动力学特征。本文是在应急勘查的基础上，采用Massflow软件完成了滑坡残留体失稳堵江过程数值分析。

1 滑坡残留体特征及发展变化趋势

1.1 滑坡残留体特征

许强等[5]通过应急调查、多期历史卫星影像遥感解译、Insar、无人机航拍等技术手段对“10.10”和“11.03”白格两次滑坡-堰塞堵江基本特征进行了详细的分析;冯文凯等[6]进一步分析“10.10”滑坡形成机制和发展演化过程。两次滑动产生的卸荷牵引作用强烈，在滑坡区后缘及中后部两侧产生强烈变形，形成一系列潜在不稳定块体(残留体)。为尽量减缓残留不稳定块体的变形发展并提高其稳定性，降低残留体大规模滑动再次堵江的风险，相关部门及单位对残留体采取了以削方减载为主的应急处置措施。根据削方后残留体变形特征、物质组成及空间位置关系，将残留体分为3个区块：滑坡后缘为K1残留体;滑坡下游侧边界为K2残留体;滑坡上游侧边界为K3残留体(见图1)。其中，K1区可分为K1-1、K1-2。此外，削方弃渣和局部滑塌残留于滑床部位的堆积体QZ(图1)形成新的堵江风险源[7]。

图1 滑坡残留体分布图

根据勘查资料估算得到各分区残留体方量为：K1-1区约209.2×104m3，K1-2区约24.6×104m3，K2区约62.7×104m3，K3区约93.6×104m3，QZ区约162.6×104m3，总计约552.7×104m3。

1.2 残留体稳定性及发展变化趋势

根据1 a的调查、勘察和空天地监测，滑坡残留体仍在继续调整应力及变形，其规模和稳定性受内在整体滑面和外在地下水、降雨、冻土等影响[8]，变形速率时快时慢，但变形仍未收敛。由中国地质调查局成都地调中心牵头、四川大学配合开展了滑坡残留体深部位移监测工作，共计布设7个，其中4个孔(ZK1、ZK8、ZK17、ZK18)监测到变形，见图2。K1-1区的平均变形速率最小，为0.22 mm/d(ZK1)，总体上呈收敛趋势；K3区的平均变形速率为 3.33 mm/d；K2区变形速率最大，参照ZK18剪断时的变形量，推测速率应在20 mm/d左右。

图2 深部位移监测条形图

结合宏观变形及变形监测成果综合分析，K1-1区后缘裂缝发育平缓趋于稳定，K1-2蛇纹岩山体两侧及后缘下错明显且缝持续发展，存在失稳趋势；左侧K3残留体形成整体变形剪切、拉张，但滑坡整体解体，尚难以整体下滑，但在上部岩体失稳冲击作用下可能整体滑动；右侧K2残留体前期与主滑体斜交，受主滑坡挤压，两次滑坡后迅速卸荷，垂直挤压方向(临空面)产生大量拉张裂缝。由于K2残留体前缘临空条件极佳，且岩性破碎、地下水在其前缘已形成了溢出带，后期持续变形明显，目前尚难以形成稳定老滑坡，存在再次滑坡的风险。滑坡槽弃渣堆积体QZ主要以小规模溜滑的方式堆积于滑槽缓坡平台(10°～17°)，并形成自然休止角，天然状态基本稳定，但在上部滑坡冲击及刮铲作用下形成新的滑坡物源，增加堵江风险，见图4。

2 Massflow数值模拟预测

2.1 技术原理

Massflow软件是基于深度积分的连续介质力学理论，利用改进的MacCormack-TVD有限差分方法，满足质量守恒方程式(1)以及x、y、z方向上的动量守恒方程式(2)、(3)。

(1)

(2)

式中，ρ为流体密度(kg·m-3)；h为流体高度(m)；t为时间(s)；u,v为分别为x方向和y方向上的流体速度(m/s)；kap为土压力系数；g为重力加速度(m·s-2)；(τzx)b,(τzy)b分别为底部zx方向和zy方向的剪应力(Pa)。

其中，底部剪应力τb遵循库伦破裂准则式(4)。

τb=c+ρghtan(φ)

(4)

式中，c为粘聚力(kPa)；φ为内摩擦角(°)。

2.2 数据处理流程

通过滑坡勘查资料收集，利用ArcGIS软件平台3D Analyst 扩展模块里的空间分析功能建立地形、滑体等数字高程模型，导出Massflow基础数据格式。在Massflow 软件中进行相关参数设置及模拟计算过程，并将结果导出至ArcGIS软件平台完成成果分析及可视化。具体流程见图3。

图3 数值模拟处理流程

2.3 堆积模拟结果

由于滑坡残留体空间位置及地质结构差异，其失稳和滑动模式有所不同。其中残坡积物质和片麻岩碎裂体以局部崩解、小规模滑塌崩落为主，滑动过程中逐步解体，并散落在滑坡槽，其影响较小；原主滑源区的蛇纹岩被石棉和滑石面切穿后，形成风化面和软弱面，形成块状岩块崩滑体，从而发展为整体、高位、崩落型滑坡。对各残留体失稳精确预判难度加大，为此我们对各种失稳进行组合。根据上述金沙江白格滑坡残留体(K1-1、K1-2、K2、K3)和滑槽弃渣堆积体(QZ)不同失稳块体组合，建立共计8种不同模型。按照上述计算原理和处理流程，运算得到不同组合模型下滑坡堵江预测。

由图5可见，模拟预测的5块潜在不稳定体，从失稳后基本沿着原滑坡槽向下运动到金沙江，整个运动与堆积过程持续约80 s。在80 s时大部分颗粒运动到滑坡下部，堵塞金沙江，形成长条形堰塞坝体，少部分颗粒残留在滑坡表面。

图4 典型工程地质剖面图(A-A′)

图5 t=80 s时刻不同工况滑坡堆积厚度

根据表1，预测未来再次堵江形成堰塞坝，最大厚度为58 m，坝顶高程2 959 m。监测数据表明K2块体变形强烈，稳定性差，下滑后将刮铲滑坡槽弃渣堆积体QZ，为最不利块体组合。由图4和表1可知，最不利组合2形成堰塞坝坝高38 m，坝顶高程2 938 m。金沙江白格“10.10”滑坡堰塞体主堆积高程2 998～3 005 m，“11.03”滑坡比第一次滑坡形成的堰塞坝平均高出30 m。两次滑坡均形成堰塞体，其中“10.10”滑坡后，10月12日堰塞坝自然泄流，最大溃决流量10 000 m3/s；“11.03”滑坡后，经人工开挖引流槽干预，堰塞体于11月12日溃决，相应的最大溃决流量31 000 m3/s[9]。鉴于前两次滑坡的堆积厚度来看，本次预测模拟河道形成的堆积体极有可能再次堵塞金沙江。

表1 8种失稳组合滑坡堰塞坝几何特征统计表

3 结论

(1) 根据1 a的调查、勘察和空天地监测，滑坡残留体仍处在应力调整及变形状态，K2、K3及K1-2变形速率较大。K2块体最为变形强烈，稳定性差，下滑后将刮铲滑坡槽弃渣堆积体(QZ)，为最不利块体组合。

(2) 本文采用Massflow软件平台完成了不同组合计算模型下滑坡残留体堵江过程数值模拟，并结合ArcGIS平台进行数据处理和空间可视化分析，较为准确地揭示了滑坡随时空运动及堵江过程。

(3) 模拟预测表明，不同残留体组合形成堰塞坝最大高约58 m。K2+QZ最不利组合形成的堰塞坝高为38 m，存在较高的溃坝洪水风险。因此，继续加强白格滑坡残留体的监测预警，并科学开展综合工程治理。

(4) 本文通过白格滑坡残留体危险性预测模拟分析表明，该研究方法对高位滑坡动力学特征及风险评价研究具有一定的参考和借鉴意义。

致谢:感谢四川大学邓建辉教授团队提供的滑坡残留体深部位移监测数据。