强震区泥石流堵江-洪水灾害链式过程模拟研究——以汶川县下庄沟为例

罗玉婷,王立娟,马 松,唐 尧

(1.重大危险源测控四川省重点实验室,四川 成都 610045;2.四川安信科创科技有限公司,四川 成都 610045;3.四川省安全科学技术研究院,四川 成都 610045)

汶川地震是新中国成立以来最具破坏性的地震,产生了大量的崩塌和滑坡[1-2]。由于大量同震滑坡物质和极端天气条件,泥石流成为汶川震区常见的地质灾害之一[3-5]。在震后十余年,汶川强震区泥石流发生的频率和规模都急剧增加,对山区人民的生命财产、基础设施和重建工作构成了严重威胁。强震区的山地灾害一般以灾害链的形式发育,泥石流则是链式灾害效应中的重要环节[6],与单一灾害相比,链式灾害具有时间尺度长、范围广、破坏力大的特点,造成的人员伤亡与经济损失远大于各单一灾害的破坏效应之和[7]。汶川地震大幅降低了泥石流的降雨阈值,在地震后,汶川、北川、清平等地区多次因强降雨引发了大规模的泥石流灾害。泥石流堵江是其致灾过程中较为关键的环节,直接影响到其致灾的严重程度。汶川7.9级大地震后,次生地质灾害进入一个长达25年的活跃期,受地震影响区域将频繁发生泥石流灾害链[8]。

泥石流灾害链的形成过程极为复杂且具有突发性,人们往往难以对灾害发生及时作出反应。近年来,随着遥感观测和计算机应用技术的发展,研究人员开发了许多监测和模拟泥石流灾害链的技术手段,其中,数值模拟软件成为了重现并预测泥石流灾害链的有效工具[9-10]。Cao等[11]通过水流、泥沙输移和河床演变的耦合,提出了一个浅水流体动力学模型来描述溃坝水流在河床上的运动。江中舟[12]采用SPH的物理方法,提出一种基于粒子的固液两相泥石流侵蚀模型来模拟泥石流灾害过程。杨春阳[13]利用FLO-2D和RAMMS 2种软件,分别模拟不同暴雨条件下泥石流的活动特征和启动临界降雨量。然而,上述研究所涉及模型都只适用于单一灾害,难以模拟由几种具有不同物理机制的次生灾害组成的灾害链。目前,泥石流数值模拟方面的研究仍在不断进步、不断完善。虽然许多软件平台可以重现泥石流的运动过程,但大部分都忽视了泥石流所带来的链式效应,无法对灾害链所带来的危害进行有效的评估。

在目前的模型方法中,与泥石流有关的过程,如水文过程、滑坡和泥石流运动,大多是单独模拟的,忽略了不同过程之间所产生的相互作用或反馈,模型的预测能力受到极大的限制。为了更有效地评估灾害链产生的风险问题,应考虑采取一种综合集成的方法,同时耦合滑坡、泥石流运动、洪水运动各个过程。链生的灾害过程以及泥石流与洪水之间的相互作用很可能会造成新的致灾效应,会超过它们各自独立相加的影响,将2种过程分开考虑通常会严重低估其产生的危害。因此,在灾害链调查和风险评估中需要实现多灾害模拟,而数值模拟软件OpenLISEM在边坡破坏、泥石流冲出预测以及洪水运动模拟等方面都表现出优越的模拟效果。本文以汶川县2019年下庄沟泥石流为对象,通过利用OpenLISEM这一多灾害耦合模型来实现泥石流堵江灾害链的数值模拟研究。

1 研究区概况

下庄沟沟口地形沟道狭窄呈“V”型,岸坡陡峻,两岸坡度40°以上的面积占研究区总面积的49%,大量崩滑体松散固体物质堆积于坡脚及沟道内,在降雨作用下极易启动转化为泥石流。下庄沟最大高程为4 120 m,沟口地势最低处高程为1 450 m,相对高差达2 670 m,沟床平均纵坡降为262.4‰,高陡的地形有利于快速汇水,为堆积在沟谷两岸的崩滑体物源的侵蚀提供了充足的动力条件。

2019年8月20日,汶川县境内突降暴雨,导致岷江、杂谷脑河沿岸数十条沟暴发泥石流灾害。其中,下庄沟于8月20日凌晨暴发泥石流,泥石流发生前期累积降雨量达27.9 mm,诱发此次泥石流的小时雨强为12.7 mm/h,由于沟口处的杂谷脑主河宽度约90 m,此次暴发泥石流冲出的固体物质达10.05万m3,大量碎屑物质搬运出沟口后,在沟口形成大型泥石流堆积扇,堵塞杂谷脑河并形成堰塞湖,致使上游水位抬升近5 m,上游大量房屋与基础设施被淹没。此外,泥石流致使沟口G317国道被於埋近180 m,同时损坏下庄村房屋51间,掩埋蓉昌高速公路的多个桥墩,对上下游居民的生命安全构成巨大威胁。

2 模型理论

OpenLISEM模型中集成了滑坡、泥石流和山洪等多种灾害的理论过程,主要通过输入栅格图层的属性求解单元网格特定过程和控制流的微分方程,该模型可将水文过程和地表运动过程结合起来,在流域尺度上实现基于降雨事件的径流、洪水、滑坡和泥石流运动模拟。

a)地表流。在OpenLISEM模型中,地表流是地表径流和洪水的组合,其运输以及动力学特征可以用描述流体的微分方程来模拟,所有的流体运动计算都基于浅表流的圣维南二维非恒定流运动方程。

(1)

(2)

(3)

式中h——流深,m;u——流速,m/s;R——降雨量,m;I——降雨入渗量,m;g——重力加速度,m/s2;S——摩擦阻力项,m/s2;Sf——摩擦动力项,m/s2;n——曼宁摩擦系数。

b)边坡稳定性。边坡的稳定性评估是基于无限边坡模型,该方法计算了下滑力和抗滑力,其中下滑力的计算是基于破坏面平行于表面的假设。在OpenLISEM模型中可以提供动态的降雨输入选项,边坡破坏的触发机制为降雨事件中湿润锋下渗,湿润锋可通过增加土体总重量和提高地下水位来影响边坡稳定性。若湿润锋下渗未达到地下水位,则只增加土体总重量。湿润锋对边坡稳定性的影响可以见式(4):

(4)

式中 SF——边坡安全系数;c′——土壤有效黏聚力,kPa;c——根系的有效黏聚力,kPa;γ——土壤单位重度,kN/m3;γs——饱和土单位重度,kN/m3;γw——水的重度,kN/m3;β——边坡坡度,(°);φ′——有效内摩擦角,(°);Z——土壤厚度,m;Zw——初始地下水水位,m。

c)泥石流运动方程。为了模拟洪水和泥石流的动力学特征和相互作用,OpenLISEM模型中运用了一组使用广泛的两相泥石流方程,该方程基于物理的二相动量守恒定律,除包含压力和重力外,还包括黏滞力、非牛顿流体黏滞力、固相阻力和摩尔-库仑摩擦力[14]。这种方法可以在非黏性流、高浓度流和泥石流之间实现平稳过渡,解决了不同类型的流体之间的相互作用。

(5)

(6)

(7)

(8)

式中αs——固体体积分数;αf——液体体积分数;Pb——基面压强,kPa;NR——雷诺数;NRA——准雷诺数;CDG——阻力系数;ρs——固体密度,kg/m3;ρf——液体密度,kg/m3;γ——液体固体之间的密度比;x——流体垂直剪切的速度,m/s;ε——模型的高宽比;ξ——固体体积分数的垂直分布。

3 数据获取与参数厘定

3.1 数据获取

为探究下庄沟泥石流形成、运动以及堵江形成堰塞湖的整个灾害链过程,收集了下庄沟高精度的地形数据、无人机正射影像、卫星遥感影像以及下庄沟邻近雨量站的实时降雨监测数据,以此重现降雨事件下泥石流堵江灾害链的动态情境。在OpenLISEM数值模拟软件中,主要输入数据分为4类:高程数据(DEM)、土地植被数据、物源数据和降雨数据。

a)地形数据。本文的数值模拟基于1∶1 000等高线地形数据,经过与无人机航拍地形数据合成校正后,利用Arcgis平台生成研究区的数字高程模型(DEM),最终得到5 m×5 m的栅格数据(图1a)。

a)数字高程模型

b)物源数据。基于高分辨率遥感影像及无人机航拍影像,对研究区进行精细的遥感解译以此获取泥石流物源分布,并根据后期野外现场测量数据进行校准,物源厚度是泥石流数值模拟中极为重要的参数之一。地震滑坡为震后泥石流主要物源,其厚度计算依据前人研究提出的滑坡平均厚度与滑坡面积之间的函数关系,这里运用式(9)计算滑坡物源厚度[14],沟道物源厚度则根据野外侵蚀断面测量数据进行估算。

SD=1.105lnA-4.795

(9)

式中 SD——物源厚度,m;A——滑坡面积,m2。

计算后通过Arcgis平台对研究区物源的厚度属性进行赋值,并将物源矢量文件转换为栅格文件,导入QGIS软件后转换为模拟所需的Map文件(图1c)。

c)土地植被数据。基于遥感观测和野外调查,下庄沟流域内植被发育,主要以灌木为主,次为乔木以及草本植物,覆盖率达70%以上。为反映植被覆盖情况并准备基础数据,综合影像质量、成像时间和影像可使用率等方面选择了空间分辨率30 m的2020年7月27日Landsat7影像作为数据源,利用ENVI提取并估算了研究区的植被覆盖度(图1b)。不同的土地利用类型在灾害形成过程中,往往表现出不同的力学特征。基于对不同土地利用类型的遥感解译分类,根据OpenLISEM使用手册(2018版)对相应土地利用类型的多个属性进行赋值[16](表1)。

表1 土地利用类型相关参数取值(OpenLISEM使用手册,2018)

d)降雨数据。水文过程包括降雨、地表蓄水、截留、渗透、蒸发蒸腾和径流,模型中利用时间和空间的降雨强度值来模拟降雨。为了还原下庄沟“8·20”泥石流的灾害过程,需要输入触发泥石流的降雨事件的实时监测降雨过程,包括降雨历时与泥石流激发雨强。2019年8月20日凌晨3:00左右下庄沟暴发泥石流,本文收集了下庄沟邻近雨量站的降雨数据,并将泥石流暴发过程的雨量数据作为降雨事件输入软件进行模拟,见图2。

图2 降雨过程曲线

3.2 参数厘定

泥石流与洪水模型预测中有一个共同的难题,即输入数据和参数的准确性。只有在输入大量足够精确参数的情况下,才能使用综合方法来模拟边坡破坏、泥石流和洪水这一系列灾害过程。然而,由于近年来科学技术的发展,可获得的精密数据越来越多,数据获取问题已变得不那么困难。而本文为了进行更精准合理的模拟预测,针对研究区开展了多方面的现场灾害调查以及室内外试验(图3),并结合相关研究厘定模型中所需要的特征参数[15],见表2。

a)现场调查

表2 下庄沟OpenLISEM模拟主要参数

4 模拟结果分析与验证

4.1 结果分析

利用对应模型和相应参数对下庄沟“8·20”泥石流运动进行模拟,分析其影响范围及运动特征。此次模拟总时长为120 min,通过OpenLISEM模型计算,重现了下庄沟“8·20”泥石流启动—搬运—侵蚀—堆积整个过程的运动特征,图4分别是不同步长(1 step=45 s)的泥石流分布。从模拟结果可以看出,降雨事件的前期阶段t=step30时,随着降雨强度逐渐增大,松散土体达到饱和并受到径流的强烈冲刷,数个边坡逐渐开始失稳运动,此刻还未启动形成泥石流。t=step40时,大量靠近沟道的岸坡失稳破坏并快速向地势低洼处汇集,主要集中在下庄沟沟段下游部分,成为补给泥石流的主要来源之一。t=step50时,持续的降雨形成地表径流,松散固体物质不断在沟道汇集,在洪流强大的水动力作用下形成阵流,不断向下游推进。t=step100时,由于1号滑坡失稳,形成平均厚度约7 m、宽约35 m、长约80 m的滑坡坝,后由于上游汇聚的携沙水流的不断冲击,堰塞体发生溃决,流量及流速产生放大效应,高速向下游沟口推进。

a)step30

图5、6表明了监测点分布和泥石流沿程最大流深及流速的变化,距离沟口5 500 m的上游沟道处流深及流速较低,最大流速约为0.79 m/s,流深约2.62 m,这一阶段处于径流汇聚及能量积聚过程。由于重力势能不断转化为动能,流速随时间逐渐变大,但在局部弯道处,流速有所减缓,但整体上表现出增加的趋势,最大达3.39 m/s。随着动力作用的增强,泥石流流体的侵蚀作用及携沙能力也得到大幅提高,刮铲沟床及沟岸,使沟床普遍加宽,沟床下切深度加大,由于物源的不断汇入,泥石流规模变大。当泥石流流体达到距沟口1 300 m处的1号滑坡堵溃点,由于堵塞作用流速急剧下降,约为0.58 m/s,泥石流深度约为8.32 m,在泥石流流体持续的冲击作用下,滑坡坝失稳,产生溃决放大效应,导致流速急剧增加至4.73 m/s,冲出沟口,在沟口的最大堆积厚度为7.39 m。

图5 泥石流动力过程监测点位置分布

图6 泥石流沿程最大流深及流速变化曲线

t=step120时,泥石流继续向前推进至G317国道距离杂谷脑河道约10 m处(图7),到达沟口堆积扇顶部,此时泥石流流深普遍在2.7～3.4 m范围内,流速在4.2～4.7 m/s。t=step130时,泥石流进入杂谷脑河,由于杂谷脑河河床纵坡降较缓,泥石流整体推进速度降低至1.5～2.7 m/s,冲出的固体物质基本堵塞河道,形成长约146 m,宽约51 m,平均厚5.5 m的堰塞体。t=step140及t=step150时,泥石流完全堵塞主河,分别形成了长、宽、厚为236.0、74.0、6.2 m和312.0、79.0、7.1 m的堰塞体(图7),因杂谷脑河对堆积扇的扰动以及冲刷作用,泥石流堆积扇平面形态呈偏向下游的长条形。主河完全被堵塞后,造成上游水位抬升,居民区房屋及公路被淹没。坝后及监测点A、B、C的洪水水深时间曲线见图8,上游居民区淹没深度为3～5 m,与实际调查情况较为吻合。

图8 监测点洪水过程曲线

a)step120

4.2 模拟结果验证

前期野外实地测量得知,下庄沟“8·20”泥石流堆积范围为2.17×104m2,堆积扇平均厚度5.2 m左右,冲出方量约10.5×104m3。为验证模拟结果的准确性,以野外实际调查情况与模拟计算堆积结果对比验证(表3),经过对比计算,模拟堆积范围与实测范围重叠区域面积2.04×104m2,形态大小基本一致。同时,由式(10)计算得到模拟精度,经过计算,模拟精度达到84.9%,符合数值模拟精度要求。

表3 模拟结果与野外实测验证对比

(10)

式中A——泥石流数值模拟精度,%;S0——模拟与实测堆积范围重叠区域,m2;SM——野外实测堆积范围,m2;SN——模拟堆积范围,m2。

5 结论

震后泥石流的形成是地震与降雨共同作用的结果,研究其运动过程和灾害链效应,有助于实现泥石流灾害链风险评价的精准性和有效性。通过对2019年8月20日下庄沟泥石流堵江灾害链的调查,首次利用OpenLISEM重现了下庄沟“8·20”泥石流从启动、搬运、侵蚀、冲出堵江到上游洪水淹没整条灾害链的动态过程,揭示了泥石流堵江灾害链的形成过程与运动特征,得出以下认识。

a)作为震后泥石流堵江-洪水灾害链的典型案例,下庄沟泥石流堵江灾害链的形成主要是流域内滑坡体失稳造成沟道堵塞,滑坡坝失稳,产生溃决放大效应,导致流速急剧增加至4.73 m/s,冲出沟口物质达10.05万m3。随着泥石流形成的堆积扇渐进向河道推进,主河过流断面逐渐缩减,造成上游的洪水淹没危险。

b)OpenLISEM模型应用于下庄沟“8·20”泥石流灾害,模拟重现了泥石流的形成—运动—冲出—堵江—洪水动态过程,不仅确定了泥石流的危险范围,同时获得了堵江后洪水淹没深度动态图,可以有效评估后期该区域可能存在的灾害风险,为当地风险管控提供依据。通过将模拟结果与野外实测数据对比,模拟的准确度达84.9%,表明该模型对滑坡-泥石流-洪水的灾害链过程模拟具有较好的适用性,尤其在同震滑坡物源充分发育的汶川震区具有极大的意义。

c)灾难性的地震对环境的影响是一个长期的过程,汶川震后的十余年来,群发性泥石流事件一次又一次对山区人民的生存环境造成毁灭性的打击。在对泥石流灾害事件的模拟中,如果只单独考虑单一灾害过程的模拟,忽略其他次级过程相互作用机制,这会对模拟结果的预测精度和应用造成较大影响。另外,泥石流的发生往往会伴随洪水,将这2个过程分开考虑通常会严重低估其产生的风险及危害。