基于FLO-2D数值模拟的工布江达县城泥石流灾害危险性评价

张浩韦,刘福臻,王军朝,张佳佳

1.西南石油大学土木工程与测绘学院,四川 成都 610500;

2.中国地质调查局地质灾害防治技术中心, 四川 成都 611734

0 引言

泥石流危险性评价被认为是减少泥石流危害的有效途径(张春山等,2004),国内外学者对泥石流灾害的危险性已进行了一定程度的研究,方法众多。较为传统的方法主要包括层次分析法(铁永波和唐川,2006;Yalcin,2008)、模糊数学法(王欣宝等,2000)、粗糙集理论(匡乐红等,2006)、灰色系统评价法 (罗冠枝和徐林荣,2008)以及神经网络评价法(谷秀芝等,2010)等。这些传统方法多是基于数学方法,通过采用半定量或定性的分析方法建立评价模型来评价泥石流的危险性(李阔和唐川,2007)。随着计算机技术的发展,GIS技术与传统评价方法开始相互结合(王骏等,2014),使得传统评价方法的科学性日益增强,但这些评价方法依然无法精确量化泥石流对居民聚集区的影响,难以对泥石流的运动和堆积进行刻画。

与传统评价方法不同,数值模拟评价的重点在于模拟泥石流在复杂地形下的运动与堆积情况(龚柯等,2017;陈明等,2018),量化泥石流泥深、流速和堆积范围等特征(梁鸿熙等,2016;罗玉婷等,2020;刘晶晶等,2020;张建石,2020;李明威等,2021),重现泥石流对人口聚集区的影响。近年来,FLO-2D作为二维动力模拟模型,已广泛应用于泥石流灾害的危险性量化,并表现出了令人满意的效果(Lin et al., 2011;杨涛等,2018;侯圣山等,2021)。在区域泥石流危险性的研究中,FLO-2D模拟程序能够确定泥石流流体深度与流速的时空分布(张鹏等,2014;徐惠梁,2018),并以此为基础划定泥石流危险性分区(常鸣等,2019;王高峰等,2020;曹鹏等,2021),圈定潜在威胁较大的人口聚集区,实现对泥石流灾害风险的定量分析。

西藏林芝地区因其活跃的构造运动,形成了山势陡峻和沟谷深切的高山峡谷地貌,相对高差大(蒋清明等,2017;张佳佳等,2018);地层岩性复杂多变,岩石风化剧烈,岩土体完整性差,物源丰富(张奋翔,2019);降雨在时空分布上较为集中。历史上,西藏林芝地区的古乡沟(鲁安新等,2006)、天摩沟(高波等,2019)、卡达沟(钟鑫等,2018)、扎木弄沟(李俊等,2017)均暴发过破坏程度严重的泥石流灾害。近年来,林芝市工布江达县境内受频繁的地震活动影响,地质灾害数量上升,灾害点密度增大,物源储量陡增,在遇降雨和冻融等极端工况条件下,泥石流频发。该地区泥石流普遍具有雨季活动频率高、多点暴发、暴发以中—小型规模居多的特性(陈宁生等,2011),其暴发时冲击力强且破坏力巨大,对居民聚集区造成了巨大威胁。

工布江达县城目前已开展了重点场镇地质灾害调查工作,但并未针对周边泥石流灾害进行定量的危险性评价。为更好探求工布江达县城泥石流危险性的时空分布特征,确定其具体危险范围和危险性大小,文章以野外现场调查及光学遥感影像为数据基础,选取4种不同暴发频率(100年一遇、50年一遇、20年一遇及10年一遇),通过模拟泥石流冲出特征(赵东亮等,2021),确定泥石流强度并划分危险性分区,以期为当地泥石流的预警与治理提供科学参考。

1 区域地质环境

工布江达县位于西藏自治区东南部,林芝市西北部,县域内山势陡峻,沟谷深切,属典型的深切割高山河谷地貌。研究区为工布江达县城所在区域,其整体坐落于尼洋曲河谷(图1),是西藏地区典型降雨型泥石流活跃区域之一。

区域地质资料显示研究区位于冈底斯陆块、雅鲁藏布江结合带内强烈挤压和碰撞的嘉黎-然乌、多其木-东久强震带附近,构造活动强烈,历史上地震活动频繁,地震基本烈度为Ⅶ度。县城后侧斜坡主要出露晚白垩世花岗岩和早奥陶世松多岩群石英片岩,表层主要被崩坡积物、冲洪积物和泥石流堆积物覆盖,研究区剖面图如图2所示,剖面位置见图1。

图1 工布江达县城地理位置图Fig.1 Location of Kongpo Gyamda(a) Location map of Nyingchi; (b) Location map of Kongpo Gyamda;(c) Topographic map of the study area

图2 工布江达县城A-A’地质剖面图Fig.2 A-A′ geological section of Kongpo Gyamda

工布江达县城泥石流形成区整体地形陡峻,植被稀少,在构造、风化及冻融作用下,区内广泛发育板岩和灰岩地层松散堆积物;流通区地表径流对沟道两侧松散堆积物的入渗使坡体稳定性受到破坏,成为泥石流的固体物质来源。

工布江达县城处于半湿润高原季风气候带,多年平均降水量640.1 mm,日最大降水量45.2 mm,小时最大降水量16.5 mm(工布江达县地方志编纂委员会,2008)。受海拔与季风气候影响,降雨时空分布差异大,5月至9月为工布江达县的雨季,降水量占全年降水量的85%。沿沟口至山顶,降雨量随海拔升高增大,汛期大强度的集中降水是引发区内泥石流的重要因素。

工布江达县的地形、物源以及水源条件极易引发泥石流灾害。县城周边主要分布有5处泥石流灾害点(图3),左右对称分布于尼洋曲两岸的人口聚集区,灾害点密度大,严重威胁着居民聚集点安全。

图3 工布江达县城泥石流沟平面图Fig.3 Plan view of the debris flow gullies in Kongpo Gyamda(a) Jiedigang debris flow gully and Dicun debris flow gully; (b) Genianglang debris flow gully; (c) Zhexizhefengpu debris flow gully; (d) Karepu debris flow gully; (e) Distribution map of debris flow gullies

2 研究区泥石流发育特征

县城周边泥石流沟沟谷多呈“V”型,主沟沟床多顺直,平均坡降102‰~282‰,局部区域坡降可达500‰;高差较大,可达979~1970 m;平均坡度为40°~50°,补给段长度占比58%~67%。流通区植被多以低矮灌木和草地为主,形成区植被不发育、地表多以残坡积堆积体和崩滑堆积体为主的松散物质,潜在松散物源量高达60×104~300×104m3,大型潜在松散物源具体位置见图4,泥石流沟具体参数见表1。

表1 泥石流沟特征值参数Table 1 Characteristic parameters of the five debris flow gullies

该区域泥石流属于典型的降雨型泥石流(高波等,2019),在降雨量达到泥石流启动临界值时,斜坡表层径流带动滑坡体表面和前缘松散物质向下输移,形成区松散固体物质启动形成并快速涌入沟道。泥石流在沟道内快速运动过程中,高势能将大部分能量转换成泥石流流动速度和对沟床及沟道两侧斜坡的侵蚀,使得泥石流在流动过程中动储量聚集性增大,以“滚雪球”的方式不断为泥石流补充固体物质,加大泥石流的冲出规模。

3 泥石流危险性评价

FLO-2D软件模型是1988年由O′Brien提出的,它将数字高程模型(DEM)划分为相同大小、规则的地形格网,利用非牛顿流体与中央有限差分法求解泥石流运动的控制方程,以数值定量的方法来模拟泥石流冲出的流动过程、堆积范围并评价危险区域。受理论模型限制,计算过程中需要满足以下限制和假设条件 (O′Brien, 2009;Castelli et al, 2017;丛凯等, 2019):①假定流体为静水压力分布;②假定格网内各参数(高程值、粗糙系数)保持一致;③假定为浅水波模式;④假定差分时间间隔内为固定恒流;⑤不考虑沟道侵蚀现象;⑥不考虑流动过程中的跳跃和震荡现象;⑦不考虑泥石流对于工程结构的损毁现象。

通过运动方程以及连续方程,可以计算出每个网格中x方向上和y方向上的流体流速和堆积深度,进而得知流体的运动范围。通过连续方程控制泥石流的质量守恒,利用运动方程计算出相邻网格间流体的速度变化。模型方程如下:

(1)运动方程

公式中:i为降雨强度,mm/h;h为流体流深,m;t为泥石流流域降雨历时,h;λ、θ分别为x轴、y轴方向上的平均流速,m/s。

(2)连续方程

公式中:Sox、Soy分别为x、y方向的泥石流沟床坡降,‰;Sfx、Sfy分别为泥石流沟x、y方向上的摩擦坡降,‰。

若流体内部固体物质浓度较大,则流体运动过程中,固体物质的颗粒间碰撞增多,流体的扩散应力变大,需要考虑泥石流运动时颗粒之间的碰撞对泥石流流动阻力的影响。

(3)流变方程

公式中:Sf为摩擦坡降,‰;Sy为屈服坡降,‰;Sv为粘性坡降,‰;Std为紊流分散坡降,‰;τy为屈服应力,MPa;m为泥石流流体重,kg;K为层流阻力系数;h为流体流深,m;η为流体粘滞系数;n为曼宁系数(贾涛等,2015;王纳纳等,2014)。

3.1 数值模拟

3.1.1 数据来源与数据处理

数据模拟主要利用的数据有:DEM、遥感影像、地质资料、降雨数据和地质灾害数据资料。高程DEM数据来源于日本宇宙航空研究所(JAXA)的12.5 m ALOS PALSAR传感器获取的数字高程模型,高精度遥感影像图来源于Google Earth 1m分辨率遥感影像图,地质资料来源于中国地质调查局的1∶25万区域地质图,降雨数据来源于工布江达县气象局,地质灾害调查数据来源于工布江达县地质灾害调查数据库。所需主要数据说明如表2所示。

表2 数值模拟涉及数据说明Table 2 Data description of the numerical simulation

在ArcGIS里将DEM转化为FLO-2D数值模拟软件能够识别的ASCII文件,根据所划定泥石流流域的高分辨率DEM作为基础数据,依据泥石流沟的流域面积、泥石流峰值流量以及电脑处理性能建立10 m×10 m的计算格网用以划分泥石流流域边界,并依据沟道位置确定计算域,随后对格网进行高程赋值,完成对地形数据的预处理。

3.1.2 基础参数设定

泥石流是一种固-液两相体,高浓度的粘性泥石流流体夹杂着大量的固体(王骏等,2016),且沿着深度方向和前进方向呈不均匀分布,固相和液相所占的比例很大程度上决定着泥石流流体的流变特征,同时也对泥石流堆积体的前缘形态有一定程度上的影响(胡宏伟,2019)。在模型中需要确定泥石流固体物质比重值以及体积浓度,由于现场配制泥石流浆体难度较大,此次泥石流流体重度采用查表法求得。按照《泥石流灾害防治工程勘查规范》 (中华人民共和国国土资源部,2006;DZ/T0220—2006)附录H填写泥石流调查表并按附录G进行易发程度评分,根据表附录G.2采用查表法分别确定县城周围五条泥石流沟的泥石流重度和泥沙修正系数。根据泥沙修正系数公式反求得到泥石流固体物质比重,从而确定泥石流体积浓度(表3;余斌,2010)。

公式中:γc为泥石流重度,t/m3;γw为水的重度,取1.0 t/m3;γH为泥石流固体物质的重度,t/m3。

曼宁系数(n)是表征流体流经不同摩擦阻力介质时的流变特征的参数,在FLO-2D数值模拟软件的计算过程中,需要确定不同微地貌特征区的曼宁系数。此次计算中依据FLO-2D使用手册所提供的建议取值表,根据泥石流流域范围内地面条件确定整体流域曼宁系数(表3),再根据实际调查以及高分辨影像数据对典型的微地貌特征重新赋值。

通常情况下,在泥石流模拟的初始设置中,默认泥石流与暴雨同时发生,即洪峰经历的时间等于泥石流历经时间。经过调查发现工布江达县暴雨引起的泥石流多为单峰型涨落曲线,单峰洪水历时为汇水面积,km2),泥石流历时取值如表3所示。

表3 工布江达县城泥石流基础参数拟定Table 3 Basic parameters of the five debris flow gullies in Kongpo Gyamda

屈服应力与粘滞系数(α1、β1、α2、β2)与体积浓度相关性最大,呈指数相关;在泥石流流体运动过程中泥石流分层流动现象产生层流阻力,阻力系数(K)主要受到地面条件影响。实际现场测量难以获取这些参数,文章根据FLO-2D软件使用手册所给的建议值结合Chang et al.(2017)模拟小流域泥石流参考取值进行赋值。工布江达县城5条泥石流沟的体积浓度以及地面条件大致相同,通过泥石流冲出模拟实验对赋值参数进行修正,最终参数取值见表4。

表4 工布江达县城泥石流基础参数取值Table 4 Values of basic parameters of debris flows in Kongpo Gyamda

在运用FLO-2D进行模拟泥石流时,集水点的选取至关重要,大部分的泥石流起动点在松散物源聚集处、在沟道物源堆积处和水动力条件充足的位置(杨涛等,2017)。综合考虑工布江达县城泥石流水流条件和物源条件,将泥石流启动点大致选在形成区和流通区交界点附近。再根据遥感数据解译的大规模崩塌滑坡体,结合实地踏勘泥石流沟的灾害发育特点,最后确定实验的集水点位置。

藏东地区工布江达县的中小流域暴雨洪水特征与四川西南山区相似,故可以参照《四川省中小流域暴雨洪水计算手册》 (四川省水利电力厅,1984),以工布江达县1990年至2019年的降雨数据作为基础,按手册给出的建议公式计算出不同降雨频率下的暴雨洪峰流量。通过Fang et al.(2019)介绍的五边形法确定流域清水流量过程曲线,在实验中假设泥石流与暴雨同频率、同步发生且暴雨洪水设计流量全部转化为泥石流流量,泥石流流量过程线为单峰型涨落曲线。泥石流运动过程中会产生一定的放大效应,故最终输入的泥石流流量过程数值为泥石流清水流量乘以体积膨胀系数BF=1/(1-Cv)(图5)。

图5 工布江达县城泥石流沟不同频率下流量过程曲线Fig.5 The flow-process curves of the debris flow gullies in Kongpo Gyamda under different frequencies(a) Dicun debris flow gully; (b) Jiedigang debris flow gully; (c) Genianglang debris flow gully; (d) Zhexizhefengpu debris flow gully; (e) Karepu debris flow gully

3.1.3 模拟验证

首先模拟工布江达县城2013年8月暴发的20年一遇规模的泥石流,得出县城20年一遇的降雨频率下周边泥石流泥深、流速和堆积特征。

模拟结果显示,20年一遇的泥石流冲出过程中,受沟道影响明显,最大堆积厚度主要集中在沟道内转角以及沟口处,最大堆积厚度可达3.7 m,堆积区堆积厚度沿沟口至堆积前缘逐渐下降;最大流速位于沟道内地形急剧变化处,流体在狭长顺直的流通区运动速度最快,大于3 m/s,最大速度可达7.2 m/s,泥石流流速在沟口处迅速降低,沿堆积区逐渐歇止。部分泥石流未形成有效冲出,对下方聚集区无威胁;地村泥石流、结底岗泥石流下方居民聚集区遭受小规模轻度淤积,与实际调查情况一致(图6)。

图6 20年一遇的暴雨条件下模拟结果Fig.6 Simulation results under the heaviest rain conditions in recent 20 years(a) Velocity; (b) Depth

为了检验FLO-2D模型对于工布江达地区高寒高海拔泥石流沟的适用性,以2013年8月暴发的20年一遇降雨工况下的泥石流模拟计算结果为例。将模拟泥石流堆积范围和冲出方量与实际野外调查、访问和遥感解译所得结果进行统计分析与比较,再对其误差分析,结果如表5所示。

表5 模拟结果与实际情况对比表Table 5 Comparison between simulation results and actual conditions

结果显示:因FLO-2D无法模拟泥石流冲出过程中的沟道侵蚀作用,堆积范围模拟结果大多小于实际结果;从2013年8月大规模暴发泥石流之后每年雨季均有小规模物源冲出,与实际调查值相比,模拟泥石流冲出的方量较小。近年来结底岗沟泥石流堆积扇上方人为改造作用明显,多开垦土地改造农田和修建房屋,模拟结果大于调查结果;卡热普沟口处修建有拦沙池,堆积范围和冲出量模拟与实际差异较大。其余泥石流沟堆积范围误差率为-19.3%~9.0%,冲出方量误差率为-15.35%~18.2%,模拟结果较好。

采用同样的方法模拟工布江达县城周边5条泥石流沟在不同暴发频率下 (100年一遇、50年一遇、10年一遇)的堆积特征,结果如表6所示。

表6 不同频率下的泥石流模拟结果表Table 6 Table of simulation results of the debris flow gullies at different frequencies

3.2 工布江达县城泥石流危险性评价

目前泥石流危险性分级标准的影响因子主要有泥深、流速以及暴雨暴发的频率和强度等。通过泥深(h)、流速与泥深乘积结果(vh)划分泥石流强度(表7),结合不同暴雨重现周期按泥石流冲出的区域分为高危险区、中危险区、低危险区三种类型,运用ArcGIS对各危险区进行赋值,构建泥石流危险性逻辑关系(图7),并绘制工布江达县城泥石流不同分级的危险区划图(图8;唐川等,1994;Chang et al, 2020)。

图7 工布江达县城泥石流危险性评定标准Fig.7 The risk assessment standard for the debris flows in Kongpo Gyamda

图8 泥石流危险性分区图Fig.8 Zoning map of the debris flows

表 7 泥石流强度划分表Table 7 Table of debris flow intensity

4 结果分析

在不同频率的降雨条件下,对工布江达县城泥石流堆积情况进行模拟(表6),可知:县城5条泥石流沟均有不同规模冲出,各娘朗沟和卡热普沟在100年一遇的降雨条件下预估堆积范围分别为4.19 km2和7.62 km2,虽有一定方量物源冲出,但受到沟内地形条件影响和工程治理等因素,未形成有效冲出,对下方居民聚集区造成危险可能性较小;地村沟、结底岗沟及折西折凤普沟在10年一遇的降雨条件下,堆积范围分别为7.97 km2、3.84 km2、11.37 km2,而在100年一遇的暴发条件下预测堆积范围可达12.99 km2、16.87 km2和43.05 km2,且折西折凤普沟堆积区平均泥深可达0.42 m,对堆积扇上方居民聚集区将造成影响;通过对5条泥石流模拟结果发现,随着降雨量的增加,泥石流平均泥深、流速、堆积范围、冲出方量显著增加。

对危险性较高的3条沟进行统计(表8),结底岗沟泥石流堆积扇高危险区占比较小,仅占比5.79%,地村和折西折凤普高危险区占比分别为17.48%和33.14%。结底岗泥石流出口处有一弯道,泥石流动能被减弱,物源淤积,虽有一定物源冲出沟口,但高危险性占比较少;地村泥石流由于泥石流流域面积小、形成区地质灾害较少,物源储量一般,且流通区较短,泥石流冲出规模有限,低危险性区域占比较多;而折西折凤普流域面积大,集水性好,物源丰富,且流通区沟道顺直无堵塞现象,导致大量物源冲出。

表8 不同降雨频率下泥石流堆积区危险性分区统计表Table 8 Risk zoning table of debris flow accumulation areas under different rainfall frequencies

通过危险性分区图对3条高危泥石流沟堆积扇上方人口聚集地淤埋区域进行对比(图9—图11),统计分析泥石流对堆积扇上方居民区影响程度(表9)。结果可知:地村沟受威胁性最大,泥石流高危险区占人口聚集区比例高达80.6%,且有被全部淤埋的可能,结底岗沟和折西折凤普沟受影响区域分别为33.7%和62.7%,结底岗沟和折西折凤普沟高危险区占人口聚集区比例分别为为0.5%和18.4%。

图9 地村沟堆积扇人口聚集区灾害发生前后对比图Fig.9 Comparison of before and after the disaster in the densely populated area of the Dicun accumulation fan(a) Before the disaster; (b) After the disaster

图11 折西折凤普沟堆积扇人口聚集区灾害发生前后对比图Fig.11 Comparison of before and after the disaster in the densely populated area of Zhexizhefengpu accumulation fan(a) Before the disaster; (b) After the disaster

表9 人口聚集区受影响面积占比Table 9 Proportion of the affected areas in densely populated area

5 结论

文章结合工布江达县周边泥石流特征,利用FLO-2D数值模型软件重现4种不同暴发频率的泥石流运动和堆积特征,以此评价工布江达县城周边5条典型泥石流沟对于人口聚集区的危险性,取得以下认识。

图10 结底岗沟堆积扇人口聚集区灾害发生前后对比图Fig.10 Comparison of before and after the disaster in the densely populated area of the Jiedigang accumulation fan(a) Before the disaster; (b) After the disaster

(1)县城周边泥石流荒漠地带风化严重,岩体结构面发育,坡面堆积有风化剥蚀碎屑物以及薄层坡积泥质碎石土,物源丰富且参与泥石流程度高;泥石流沟多为“V”型谷,沟道顺直,在泥石流活动过程中动能将聚集式的增加;相对高差大,泥石流势能条件较好,在汛期短时强降雨下,物源冲出可能性极大。

(2)数值模拟的冲淤特征与实际相符,确定了该评价方法对于藏东林芝地区人口聚集区周边泥石流危险性的可行性。

(3)在4种不同降雨条件下,卡热普沟和各娘朗沟未形成有效冲出,对下方居民聚集区造成危险可能性较小;地村沟、结底岗沟及折西折凤普沟堆积范围较广,对居民聚集区影响较大,建议对泥石流工程治理或居民地进行搬迁避让。

(4)绘制危险分布图可以确定威胁聚集区具体位置,能够为人口聚集区防灾预案以及紧急避难点提供有效的参考。