易贡滑坡-碎屑流-堰塞坝溃坝链生灾害全过程模拟与动态特征分析

戴兴建，殷跃平，邢爱国

(1.上海交通大学船舶海洋与建筑工程学院，上海 200240；2.中国地质环境监测院，北京 100081)

0 引言

高山峡谷地区地质灾害频发，其中高速远程滑坡占比突出。高速远程滑坡具有高位、高速、远程等特点，易诱发次生灾害形成灾害链，造成大规模的人员伤亡和财产损失[1-3]。当滑坡运动方向有水体分布时，滑坡碎屑堵江易形成堰塞湖，一旦坝体溃决将对下游区域造成巨大破坏。世界范围内发生过多次灾难性的滑坡堵江事件[4-6]，我国也是此类灾害的多发国家[1,7]。2000年4月9日西藏波密县易贡乡扎木弄沟发生巨型滑坡，滑坡碎屑堵塞易贡藏布河形成堰塞湖，堰塞坝最终溃决冲击下游区域，对下游区域造成了巨大的破坏。

国内外对于易贡滑坡的研究主要集中在现场调查、室内试验和滑坡机理研究[8-11]，对于易贡滑坡的动力学研究比较少，且仅停留在灾害链部分阶段的模拟研究[12-14]。对滑坡-堰塞坝溃坝-洪水演进完整灾害链的动力学研究几乎为空白，因此对于易贡滑坡完整灾害链动力致灾过程的模拟研究具有重要的学术意义和应用价值。

在现场调查资料和上述研究成果的基础上，本文依据遥感影像数据建立三维数值模型，结合DAN3D和FLOW3D模拟易贡滑坡-碎屑流-堰塞坝溃坝链生灾害全过程分析，以期揭示高山峡谷区滑坡-堰塞坝溃坝链生灾害的动力学特征。

1 地质环境条件

2000年4月9日西藏自治区波密县易贡乡扎木弄沟(E94°55′～95°，N30°10′～30°15′)发生巨型山体滑坡，滑坡点位于雅鲁藏布江支流易贡藏布河左侧(图1a、1b)。易贡地区位于西藏东南部，属于高山峡谷地貌，该地区地势陡峻，河谷纵横，气候复杂。易贡藏布河是雅鲁藏布江下游北侧的一条支流，长286 km，流域面积13 533 km2，日平均径流量3.27×107m3。

由于气温转暖、冰雪融化、持续降雨等原因，约9.0×107m3滑体高速下滑，冲入易贡藏布河，滑坡后缘高程5 520 m，河道高程2 190 m，垂直高差3 330 m，最大水平移动距离10 000 m(图1c)。滑坡碎屑最终堵塞易贡藏布河形成长宽2 500 m，平均高度60 m，体积约3.0×108m3的堰塞坝[8]。2000年6月10日，坝体最终溃决，对下游区域造成了严重的破坏。

图1 易贡滑坡位置图Fig.1 Location of the Yigong Landslide

滑坡区出露的地层，由老至新为前震旦系冈底斯变质岩群(AnZgd)、喜山期花岗岩、第四系堆积物。海拔4 100 m以下为前震旦系冈底斯岩群，总体倾向SSW，地层倾角在30°以上。4 100 m以上为喜山期花岗岩，在扎木弄沟内与前震旦系冈底斯变质岩群呈侵入接触关系(图2)。

2 滑坡-碎屑流分区

滑坡高差3 330 m，水平距离10 000 m。根据滑坡前后遥感影像和现场调查情况，将滑坡区分为滑源区Ⅰ、流通区Ⅱ、碎屑流堆积区Ⅲ(图1c，图2)。

(1)滑源区：位于扎木弄沟沟谷源区，高程范围4 100～5 520 m，水平距离2 000 m。顶部尖峭，终年为冰雪覆盖，崩滑体呈楔型，岩体主要为喜山期花岗岩。滑坡发生后裸露出其底部形态，形成“V”型悬谷(图3a)。

(2)流通区：位于扎木弄沟沟谷，高程范围2 600～4 100 m，水平距离3 500 m。沟谷上宽下窄，呈 NE-SW 走向，沟内原堆积有50～80 m厚碎屑堆积物，滑坡发生后沟谷内岩石出露，沟谷两侧陡壁上有明显的擦痕(图3b)，说明区内早期堆积物已全被裹挟铲刮。

(3)碎屑流堆积区：位于扎木弄沟口至藏布河南岸，高程范围2 190～2 600 m，水平距离4 500 m。约3.0×108m3滑坡物质呈喇叭状堆积于此(图3c)，形成长宽各2 500 m，高60 m堰塞坝。

图2 易贡滑坡纵断面Fig.2 Longitudinal profile of the Yigong landslide

图3 滑坡-碎屑流灾后地貌Fig.3 Photograph showing the study area after Yigong rockslide-debris avalanche

3 灾害链全程动力学分析

3.1 滑坡-碎屑流阶段模拟

3.1.1DAN3D简介

DAN(Dynamic Analysis)是加拿大学者HUNGR[15]开发的等效流体动力分析软件。DAN3D对于高速远程滑坡的模拟已经非常成熟[14-17]。DAN3D模拟过程中，滑体上边界为自由边界，滑体下边界的应力包括重力作用产生的正应力和运动时产生的剪应力，其中剪应力由流变模型决定。模拟选用Frictional模型及Voellmy模型对易贡滑坡-碎屑流运动过程进行模拟分析。

Frictional模型假设抗剪强度(τ)仅与正应力(σ)有关，即：

τ=σ(1-ru)tanφ

(1)

式中：ru——孔隙水压力系数(即孔隙水压力与基底正应力之比)；

φ——动摩擦角/(°)。

其中，孔隙水压力系数和动摩擦角之间有如下关系：

φb=arctan(1-ru)tanφ

(2)

式中：φb——有效摩擦角/(°)。

Voellmy模型假设滑体受到的阻力来源于摩擦力和湍流造成的阻力之和，即：

τ=σf+ρgv2/ξ

(3)

式中：f——摩擦项，等同于tanφb；

ξ——紊流项，与碎屑流的速度和密度有关。

3.1.2模型及参数选取

参考易贡滑坡已有研究成果[18]，建立体积为9.0×107m3滑体模型。依据SRTM-3遥感影像建立滑动路径数字高程模型(图4)。运用DAN3D对易贡滑坡-碎屑流阶段进行数值模拟，通过试错反演，滑源区采用Frictional模型，流通区和碎屑流堆积区采用Voellmy模型，反演得到计算参数见表1。

图4 易贡滑坡数字高程模型Fig.4 Digital elevation model of the Yigong Landslide

表1 流变模型和计算参数

3.1.3模拟结果分析

通过DAN3D模拟得到滑坡-碎屑流阶段不同时刻滑坡体堆积形态如图5所示。

模拟结果显示滑坡-碎屑流过程持续300 s，最大水平运动距离10 000 m，垂直高差3 330 m，平均速度35 m/s。滑体自4 000 m高程剪出后，沿NE-SW方向运动，30 s时滑体前缘进入沟谷，经过铲刮体积快速增大，90 s时前缘冲出沟口(高程2 600 m)，在沟口外开始呈扇形扩散。至150 s，滑体运动至藏布河道，210 s时，滑体前缘运动至河流对岸并停止。至300 s，滑体后缘停止运动，堆积形态呈现为喇叭状不再变化。

图5 滑坡-碎屑流过程不同时刻滑坡堆积形态Fig.5 Time-lapse image of deposit distribution during the rockslide-debris avalanche process

最终堆积形态如图6所示，滑坡堆积物呈喇叭状分布于藏布河道，最大堆积厚度100 m，平均厚度60 m，前沿最宽3.2 km。模拟得到滑坡堆积范围与实际接近，平均厚度与实测堰塞坝平均厚度(60 m)相同。

铲刮分布结果如图7所示，沟谷内铲刮较深位置在高程2 600～4 000 m，铲刮深度超过60 m，最大铲刮深度100 m，与早期碎屑堆积物厚度(50～80 m)接近。铲刮方向呈NE-SW，与谷内老滑坡堆积区形态大致相符。

速度分布模拟结果如图8所示，滑坡从海拔5 520 m加速下滑，在高程范围3 200～3 800 m时速度超过70 m/s，最大速度达90 m/s，接下来速度逐渐减小，到沟口位置速度接近于0，运动停止。从滑坡纵断面图(图2)可以看出，在滑动路径中高程范围3 200～3 800 m坡度最大，高程3 200 m以下地势逐渐平缓，模拟结果符合实际情况。

图6 滑坡-碎屑流最终堆积分布图Fig.6 Final deposit distribution of rockslide-debris avalanche

图7 滑坡-碎屑流铲刮分布图Fig.7 Scraping distribution of rockslide-debris avalanche

图8 滑坡-碎屑流最大速度分布图Fig.8 Maximum velocity of rockslide-debris avalanche

3.2 溃坝后洪水演进过程模拟

3.2.1FLOW3D简介

FLOW3D软件是基于CFD (Computational Fluid Dynamics)解算技术的计算仿真软件。其VOF算法对流体自由表面进行追踪计算，可对多种复杂流体运动现象进行精准模拟。选用溃坝模拟常用模型RNGk-ε模型进行模拟。

3.2.2三维模型建立及参数选取

依据SRTM-3遥感影像数据建立山体模型，根据DAN3D堆积区形态模拟结果(图6)，建立长宽2 500 m，高60 m喇叭状堰塞坝模型(图9)，堰塞坝材料设置为泥沙冲刷模式，以模拟坝体逐渐溃坝过程，堰塞坝材料参数见表2。并根据真实泄洪情况建立截面尺寸150 m×30 m×24 m，长1 500 m引流渠模型。设置水体模型水流深度为实测堰塞湖水深60 m。组合模型形成溃坝模拟三维模型(图10)。

表2 堰塞坝材料模型参数

图9 堰塞坝模型Fig.9 Landslide dam model

图10 研究区Flow3D三维模型Fig.10 Three-dimensional model of study area

计算区域内设置3.2×107个尺寸为25 m×20 m×10 m的网格进行模拟运算。监测点设置如图10所示，在堰塞坝坝口，坝尾分别设置监测点1、2，河道曲折度较大处设置监测点3、4，距离堰塞坝17 km通麦大桥处设置监测点5。

3.2.3模拟结果分析

通过FLOW3D模拟溃坝后洪水演进过程得到不同时刻水流深度、流量等水动力特征参数。模拟时间12 000 s。

洪水演进过程中水流深度变化如图11所示，初始水深60 m，300 s时，堰塞坝开始过流，1 500 s洪水到达监测点3，湖水深度降至46 m，2 800 s达到监测点4, 湖水深度降至45 m。3 200 s时达到通麦大桥。在8 000 s的时候，堰塞湖水深降至30 m左右不再变化，达到出入库平衡。

图11 洪水演进过程不同时刻水流深度变化图(单位: m)Fig.11 Time-lapse image of variation of water depth during the outburst flood process

洪水演进过程中各监测点水流深度变化如图12所示，泄洪过程中堰塞湖水深在8 000 s时降至30 m后逐渐达到稳定。其他监测点水流深度随洪水到达上升至一定高度后小幅波动然后急剧增大，该趋势表现出堰塞坝逐渐溃决过程，坝体完全溃决后，各监测点水流深度快速达到峰值然后逐渐降低。下游河道狭窄位置水流深度大，通麦大桥处在7 000 s前后达到峰值约100 m。

各监测点流量监测结果如图13所示，模拟得到通麦大桥处在5 100 s达到洪峰流量约130 000 m3/s，与实测洪峰流量120 000 m3/s[8]较为接近。整个洪水演进过程中，流量峰值出现在监测点4即河道最急处，此处在4 950 s时达到峰值流量240 000 m3/s。10 000 s时，各位置流量基本达到稳定，通麦大桥上游区域稳定在70 000 m3/s，通麦大桥处稳定流量约50 000 m3/s。

图12 不同时刻监测点水流深度变化曲线Fig.12 Water depth at different time steps of monitoring points

图13 不同时刻监测点流量变化曲线Fig.13 Flood discharge at different time steps of monitoring points

4 结论

基于遥感影像建立三维数值模型对易贡滑坡-堰塞坝溃坝-洪水演进完整灾害链进行反演分析。得到如下结论：

(1)运用DAN3D的F-V-V模型能够反演易贡滑坡-碎屑流运动过程，滑坡-碎屑流过程持续300 s，平均速度35 m/s，最大滑动速度90 m/s，平均堆积厚度60 m，碎屑流最终堆积形态与实际接近。

(2)基于DAN3D的碎屑流最终堆积分布结果建立堰塞坝模型，运用FLOW3D能够反演溃坝后洪水演进过程，模拟结果显示通麦大桥处洪峰流量130 000 m3/s，与实测值接近。

(3)DAN3D和FLOW3D的结合使用能够较好的反演滑坡-堰塞坝溃坝-洪水演进完整灾害链动力致灾过程，并得到相关动力学参数。对于此类链生灾害的系统性模拟分析及防治工程提供了研究思路。