基于Massflow模型的青龙沟台风暴雨型泥石流运动特征研究

熊朝正,吉 锋,石豫川

(1.四川藏区高速公路有限责任公司,四川 成都 610041;2.成都理工大学 地质灾害防治与地质环境保护国家重点实验室,四川 成都 610059)

近年来随着全球气候变化,台风对中国东南沿海地区造成了严重影响[1-2]。伴随台风而产生的暴雨是一种复杂多变的降雨类型,具有雨量大、降雨集中,极端雨强大、雨型单峰变化、持续时间短等特点[3-7]。由台风暴雨直接诱发的台风暴雨型泥石流是一类特殊的灾害类型,此类泥石流凭借其显著的突发性、群发性和破坏性,愈发引起人们的重视。开展台风暴雨型泥石流运动特征研究,具有十分重要的理论和现实意义[8]。

泥石流作为一种非均匀混合介质,具有复杂的物理过程和动力学特征[9],随着计算机技术、数值算法、本构模型等方面的发展成熟,数值模拟已成为现在泥石流运动特征研究的主要方法,其中应用较多的软件主要有FLO-2D、DAN、MassMov2D、RAMMS、EDDA、Massflow等。唐得胜[10]将FLO-2D应用于都江堰龙溪河流域泥石流沟数值模拟,结果显示重现期为20年一遇时的模拟结果与“8.13特大泥石流”对比精确度在60%～87%。Cheon等[11]将DAN-3D软件与机器学习方法相结合,对韩国Umyeon山泥石流进行泥石流宽度、体积及冲击力的计算,在生成危险性评估图的基础上选择最佳拦挡设施放置位置。Beguería[12]应用MassMov2D模型分别对奥地利蒂罗尔和法国阿尔卑斯山2个泥石流进行反分析,模拟泥石流堆积物的最终分布和厚度。宋兵等[13]探讨了RAMMS在泥石流模拟中的运用方法,通过模拟北川县乡白沙沟泥石流,得出泥石流冲出量与雨洪法计算量(20年一遇)误差不超过5%,可信度较高。Chen等[14]提出一种三维综合数值模型EDDA,该模型考虑了泥石流过程中密度、屈服应力和动力黏度的变化,验证了该模型在流域尺度模拟中的性能。

Massflow基于深度积分的连续介质力学方法[15],能揭示滑坡、泥石流、碎屑流、山洪等山地灾害的时空演化全过程,把三维计算问题简化为二维,有效地提高了计算效率。丁邦政[16]运用Massflow对尾矿库溃坝泥石流进行了数值模拟,计算泥石流最大冲起高度、冲击力、流深、流速等数据,在此基础上对桥梁进行泥石流冲击作用响应分析和风险评估。段学良等[17]利用无人机航拍三维建模生成DEM,结合Massflow软件,比选后采用更精准的Manning摩擦模型对西藏仁布杰仲沟泥石流进行灾害模拟评估,将研究区划分为4个危险区,为泥石流防治提供了重要依据。

本文运用Massflow软件进行台风暴雨型泥石流运动特征数值模拟,对比分析泥石流的运动堆积过程、泛滥范围及堆积厚度等动力学特征。研究成果不仅能给台风影响地区的地质灾害防治提供参考,同时也能为今后台风诱发泥石流的深入研究提供计算方法及手段参考。

1 研究区概况

研究区拟建抽水蓄能电站由下水库(坝)、上水库(坝)、地下厂房、输水系统等主要建筑物组成。青龙沟位于上水库区,其中部为上水库蓄水区,拟建水工大坝横跨青龙沟将其截断,为更好地评价青龙沟对库区的影响,以坝址为界,将青龙沟分为上下两段。青龙沟流域面积约1.84 km2,长度约2.86 km,高差为580 m,整体纵坡降为203‰,流域内植被发育,植被覆盖率大于60 %,其中以乔木、竹为主。青龙沟地形总体上属于构造剥蚀低山地貌,两侧支沟和小型的冲沟很发育。据野外地质调查,调查段共发育有13条主要的一级支沟(Qs1—13)及若干的小型冲沟(图1),两岸地形总体较陡,坡度多为30～50°。通过调查青龙沟物源点有13处,共有松散物源总量约185.34万m3,不稳定物源量约17.56万m3。因此青龙沟具有对泥石流暴发有利的地形和物源条件。

图1 青龙沟及支沟分布

2019年8月10日在“利奇马”引发的台风暴雨作用下,当日下午16时左右(即在台风暴雨短历时超强降雨发生时刻),启动青龙沟暴发了大规模泥石流。通过访问居住在青龙沟附近的居民,得知青龙沟泥石流暴发持续时长约为15 min。由于泥石流流量过大,无法沿原始沟道正常运动排泄,最终冲出原始沟道,致使沟道改变,在左侧形成新的泥石流沟道,并大量扩散堆积(图2),青龙沟泥石流在上水库区造成了严重的破坏,摧毁部分居民房屋,中断工程进度。通过收集研究区泥石流暴发过程实际台风暴雨资料,结合文献分析[18-20],台风暴雨具有突出的短历时强降雨特征,能在短强降雨时刻引发泥石流,青龙沟泥石流的暴发也符合这一特征。台风对泥石流的影响表现于通过树木松动表层物源,由于台风作用难以定量,该耦合过程较为复杂,且主要以台风暴雨作用为主,本文主要考虑台风暴雨作用下泥石流运动特征。

图2 青龙沟沟道改变

2 Massflow数值模拟工作

针对青龙沟泥石流的运动堆积特征,采用中科院山地所研究开发的一款地表动力过程数值模拟软件Massflow[21-22],来进行台风暴雨型泥石流运动特征研究,能取得良好的效果。该软件以深度积分连续介质力学方法为理论基础,运用MacCormack-TVD有限差分法来求解以上方程[23-24],采用Fortran和C#编程语言,可跨越Windows和Linux平台,并结合MPICH分布式和OpenMP共享内存并行计算技术,有效地提高了计算效率,减少了计算时间。Massflow通过对大区域切割和重点区域网格细化的方式提高了计算规模;此外还提供了修改物理模型及微分方程的源代码,支持用户通过二次开发编入自定义模型[25]。

2.1 地形数据处理

本次研究收集到青龙沟实测1∶2 000地形等高线和高精度遥感影像图。通过ArcGIS软件将原始高程数据进行转化,获取青龙沟DEM,本文根据实际情况建立4 m×4 m的计算网格数据。进一步运用ArcGIS软件转换工具中栅格转ASCII功能,将上一步获得的DEM数据转换为ASCII格式(txt文件)的高程文件,即可导入Massflow软件进行模拟计算。计算模型见图3。

图3 青龙沟程序计算模型

2.2 泥石流特征值设定

本文采用基于泥石流沟易发程度数量化评分的泥石流容重确定方法,确定青龙沟泥石流重度为16.34 kN/m3,相应1+φ取值1.637(φ为泥石流泥沙修正系数)。拟设计非台风暴雨(工况1)及台风暴雨(工况2)两种工况进行对比模拟分析。如前文所述,研究区实际台风降雨过程具有突出的台风暴雨特征(即短历时强降雨特征),是青龙沟泥石流的暴发的关键因素,因此,工况2采用实际降雨资料进行计算模拟。同时收集研究区往年雨季非台风暴雨资料与台风暴雨相比较,峰值特征不突出,过程雨量分布较为均匀,因此工况1将降雨过程概化为雨强均匀分布,根据青龙沟泥石流防治工程体系50年一遇的设计标准,计算工况1的24 h暴雨量并24 h平均,得到设计雨强为10.5 mm/h。2种工况的降雨过程见图4。

图4 模拟降雨工况设计

根据研究区泥石流沟谷实际汇流情况,运用实际降雨资料及当地洪水计算手册,计算概化泥石流流量曲线,结合已有学者的研究,本文采取“概化五边形模型”来进行相关计算[26]。根据青龙沟泥石流汇流特征,分别设定泥石流启动点1(Q1)和启动点2(Q2),流量过程曲线的持续时长为900 s。各启动点在不同工况下的泥石流“概化五边形”流量曲线计算结果见图5、6。

图5 启动点1泥石流流量曲线

图6 启动点2泥石流流量曲线

2.3 基底摩擦模型及参数

结合本文泥石流灾害特点,参考已有学者的研究成果,选取Voellmy模型作为模拟的基底摩擦模型。结合青龙沟8月10日泥石流暴发后的实际泥位线及堆积范围,进行大量反演工作,最终确定基底摩擦系数取值0.11,湍流系数取值260 m/s2。

3 泥石流运动特征对比分析

工况1条件下模拟结果显示,库区内的模拟堆积范围与实际堆积范围基本吻合,仅在北西侧边缘略大于实际堆积范围(图7),分析认为是由于该处居民修建住房时对局部地形进行了一定的堆积抬升,导致实际泥石流的泛滥受到一定限制。在青龙沟下段沟道调查处模拟最大泥深为173 cm,与实际泥位线高度170 cm仅相差3 cm(图8),整体上模拟结果准确性较高。

图7 堆积范围对比

图8 调查处泥位线高度

根据青龙沟泥石流在非台风暴雨型及台风暴雨型降雨工况下运动特征模拟结果,泥石流运动方向及整体运动过程基本一致。启动点1处泥石流首先向库区中心运动,随后与启动点2处泥石流汇合。交汇后的泥石流都表现出叠加增强效果,加快了泥石流在青龙沟下段沟道中的前进速度。当运动至沟口后开阔的地形使得泥石流开始快速停留堆积,均形成了典型的扇状堆积区。2种工况条件下,泥石流的最终堆积区也基本一致,主要在上库库区、青龙沟下游沟道、青龙沟沟口3处,分别为堆积区1、堆积区2和堆积区3(图9、10)。

图9 工况1条件下泥石流流深分布

图10 工况2条件下泥石流流深分布

在非台风暴雨工况下,泥石流分别在700、200 s到达监测点j1、j2,约在1 700 s时两处泥石流汇流,汇流后的泥石流在1 800 s运动至j3,2 000 s时运动至j4,2 200 s时到达j5并在沟口停留堆积。在台风暴雨型降雨工况下,泥石流分别在600、100 s到达监测点j1、j2,约在1 400 s时两处泥石流汇流,汇流后的泥石流在1 500 s运动至j3,1 600 s时运动至j4,1 900 s时到达j5并在沟口停留堆积。泥石流到达各监测点的时间分布见图11,由此可见台风暴雨型泥石流暴发后,在更少的时间内到达沿途监测点位置,能更快对沿途沟道及建筑物造成破坏,居民的逃生时机更短,这也是台风暴雨型泥石流通常会造成更大生命财产损失的重要原因之一。

图11 到达监测点时间对比

取2种工况下流深区别较大的j1、j2、j5进行统计,见图12,各监测点在相同时刻泥石流流深以及最终流深,台风暴雨型泥石流均大于非台风暴雨型泥石流。在监测点j5表现最为明显,2 500 s时工况1最终流深5.86 m,工况2最终流深7.47 m。

图12 监测点流深对比

图13为两工况条件下各堆积区最大流深统计,可见工况2条件下,3个主要堆积区的泥石流流深均大于工况1,在堆积区3最为明显,工况1泥石流流深为6.06 m,工况2为7.75 m,约为前者1.3倍,表明台风暴雨型泥石流具有更大垂直影响高度的特征。图14为两工况条件下各堆积区泛滥面积统计,在堆积区1,工况2略大于工况1,堆积区2泛滥范围基本一致,最主要差别在沟口堆积区3,工况1泛滥范围为18 824.9 m2,工况2为32 279.5 m2,约为前者1.7倍。结合两工况条件下最终泛滥范围对比(图15)不难看出,台风暴雨型泥石流在沟口堆积区的影响范围明显大于非台风暴雨型泥石流,具有更大平面影响范围的特征。

图13 堆积区最大流深

图14 堆积区面积

综上分析,台风暴雨型泥石流在整个运动过程中具有更快的前进速度,龙头能在更短的时间内对沿途沟道及建筑物造成影响和破坏,此外,其泛滥过程中的垂直影响高度和平面影响范围也都大于非台风暴雨型泥石流,说明台风暴雨型泥石流具备更大的泥石流物质量,其携带能力及破坏能力都远强于非台风暴雨型泥石流,在没有防治措施的情况下对工程建筑安全及居民的生命财产安全具有很大的威胁能力。

4 结论

a)台风暴雨型泥石流及非台风暴雨型泥石流运动方向及整体运动过程基本一致。启动点1与2处泥石流汇合,表现出叠加增强效果,运动至沟口后快速停留堆积,形成上库库区、青龙沟下游沟道、青龙沟沟口3处主要堆积区。

b)按50年一遇降雨强度设计的非台风暴雨条件下,泥石流分别在700、200、1 800、2 000、2 200 s运动至监测点j1—j5,堆积区最大堆积厚度分别为2.65、3.77 、6.06 m,堆积面积分别为33 378.4、3 349.8、18 824.9 m2。

c)在台风暴雨型降雨条件下,泥石流分别在600、100、1 500、1 600、1 900 s运动至监测点j1—j5,堆积区最大堆积厚度分别为3.08、3.80、7.75 m,堆积面积分别为36 324.4、3 700.4、32 279.5 m2。

d)台风暴雨型泥石流具有很快的前进速度,能在短时间破坏沿途沟道及建筑物,居民的逃生时机较短,同时还具有高度影响范围和平面影响范围较大的特征,破坏范围广,因此对于台风暴雨型泥石流更需注重防治措施的修建。