基于RAMMS的冰湖溃决型泥石流演进模拟及危害性

刘府生, 周瑞宸, 孙红林, 胡卸文, 黎尤, 徐玉龙

(1. 中铁第四勘察设计院集团有限公司, 武汉 430063; 2.西南交通大学地球科学与环境工程学院, 成都 610031)

由于温室效应不断积累导致全球气候变暖,冰川消融速率呈明显上升趋势,高海拔或高纬度山区冰湖规模及数量急速增加,冰湖溃决事件高发频发[1-2]。近50年喜马拉雅地区至少记载有20次大型冰碛湖溃决灾害事件,其中3/4发生在西藏[3]。在高原铁路规划建设的背景下,冰湖溃决型泥石流作为目前西藏境内最主要的地质灾害之一,研究其形成演化特征具有重要意义[4-5]。西藏冰湖溃决泥石流主要分布于藏东南地区,其造成的后果往往比常见的暴雨泥石流更为严重,不仅直接对下游的房屋建筑、公路、桥梁等工程及人员安全造成巨大危害,还极易诱发沿岸崩塌、滑坡等次生灾害[6-7]。例如,1981年7月11日聂拉木县樟藏布沟发生特大规模冰湖溃决型泥石流,致使中尼公路聂友段多数公路及桥梁被毁,200余人罹难[8]。2013年嘉黎县然则日阿错发生溃决,导致周边居民和基础设施不同程度受灾,直接经济损失达2.7亿元[9]。2020年6月嘉黎县吉翁错冰湖发生溃决,泥石流淤积沟口、堵塞主河道,对下游村镇构成严重威胁[10]。

冰湖溃决型泥石流通常暴发突然、能量巨大,具有持续时间短、洪峰流量大、波及范围广等特点。国内外许多学者对冰湖溃决的致灾诱因及机理进行了研究,总的来说分为两类:第一类是冰崩(冰滑坡)、强降水、地震等外部诱因造成水位上涨、激发涌浪,引起冰碛湖自身失衡而导致的冰湖溃决;第二类是冰碛坝内死冰消融、堤坝管涌扩大等内部诱因,导致坝体失效、引发洪水,沿途冲刷侵蚀形成冰湖溃决型泥石流[11-12]。冰湖溃决诱因众多而复杂,且灾害一旦发生后果严重,因此对冰湖溃决进行风险评价十分重要,这也是冰湖溃决型泥石流演进模拟的必要前提。冰湖溃决风险评价方法大致分为定性、半定量和定量三类,随着研究的深入与技术的进步,逐步向评价指标定量化的方向发展。根据历史上已溃决冰碛湖的不同特点,国内外学者提出了相应的冰湖溃决风险评价指标体系。Mckillop等[13]以遥感数据为基础,根据冰碛湖、冰碛坝、母冰川和冰湖盆的特点及其相互关系归纳出18项评价指标,并筛选出4个参数、基于逻辑回归方法建立了冰湖溃决风险概率方程。Ermini等[14]基于对84座堰塞坝的统计分析(阿尔卑斯山脉和亚平宁山脉36座,日本17座,美国和加拿大20座,新西兰和印度等国家11座),选取坝体体积、流域面积和坝高作为评价因素,提出采用无量纲堵塞指数(dimensionless blockage index,DBI)来评估堰塞坝稳定性。位宏等[15]选取母冰川、冰湖和终碛坝3类指标共8个因子构建了新疆冰湖溃决风险模型。余斌等[16]基于冰碛湖溃决机理,提出采用6个无量纲影响因子的冰湖溃决易发性评价模型。

冰湖溃决数值模拟旨在通过数学公式或物理模型估算溃决洪水灾害动力学特征参数、重现或预测冰湖溃决演进过程及下游危险性范围。Wang等[17]利用BREACH模型对西藏龙巴萨巴湖和皮达湖进行溃坝风险评价及两湖全溃模拟分析,并验证了该模型的可靠性。韩松林等[18]基于HEC-RAS一维非恒定流水动力模型,建立了上游冰湖至下游河道的洪水演进数学模型,模拟了桑旺错溃坝、洪水下泄的演进过程。刘建康等[19]通过RAMMS数值模拟分析了西藏波密丹卡弄巴冰湖在全溃条件下的演进过程及致灾风险。刘波等[20]利用RAMMS内置泥石流模块的Vollemy-Salm单相流模型,模拟了西藏洛隆巴曲冰湖溃决型泥石流在4种极端工况下的演进过程。王翔等[21]采用FLO-2D水文动力学模型对藏东南地区典型冰湖强宗克措和吉莱普措进行溃决泥石流模拟,为下游地区灾害防治措施提供依据。

总的来说,国内外学者有关冰湖溃决洪水及泥石流研究取得了较好的进展,但由于基础资料的缺失及调查方法的局限,对冰湖溃决风险评价仍缺乏完善的体系[22-23]。现以西藏洛隆冻错曲冰湖为例,在查明泥石流孕灾条件的基础上,对泥石流的动力学特征参数及进行计算;采用无量纲堵塞指数(DBI)对冰湖堰塞体进行溃决风险判别;采用经验公式分别计算2种溃决模式下的溃口峰值流量;基于三维动态模拟软件RAMMS分析冰湖泥石流溃决演进过程,定量评价潜在威胁区危害性,以期为研究区公路、铁路等重大工程的风险管理与防治设计提供参考。

1 流域概况

冻错曲冰湖泥石流沟位于西藏昌都市洛隆县腊久乡,拟建工程位于该泥石流沟沟口、冻错曲主河下游3 km处(图1)。冻错曲冰湖泥石流沟流域面积226.27 km2,主沟长度31.44 km,流域范围内最高海拔约5 180 m,最低海拔约3 815 m,相对高差1 365 m,主沟平均纵坡降27.27‰。

图1 冻错曲流域全貌图Fig.1 The view of the Dongcuoqu watershed

1.1 地质环境条件

冻错曲冰湖泥石流沟位于伯舒拉岭山脉东北侧,为高山峡谷地貌,受区域环境因素和地势影响,冰湖发育,该区域内泥石流以冰川-暴雨混合型泥石流为主[19]。该流域地势东高西低,沟道呈宽浅“V”形,支沟较为发育,谷坡25°～70°,谷宽20～100 m。沟口与巴曲、义俄4号沟交接,附近堆积扇由于人工改造较少,基本保持原地貌形态。研究区在区域构造上属隆格尔至纳木错断裂带,区内延伸60 km,产状203°～210°∠45°～50°,表现为逆推覆构造,斜切冻错曲冰湖上游汇水区。地表出露基岩主要为石炭-二叠系来姑组石英砂岩、白垩纪花岗岩、侏罗纪二长花岗岩(图2)。流域内支沟发育,主沟道纵坡降存在明显的分段特性,可分为上游汇水区、形成流通区和沟口堆积区(图3)。

图2 冻错曲冰湖泥石流沟物源分布图Fig.2 Distribution of material sources in the Dongcuoqu glacier-lake debris flow watershed

图3 冻错曲冰湖泥石流沟主沟道纵断面图Fig.3 Longitudinal section of the main channel of the Dongcuoqu glacier-lake debris flow

1.2 物源条件

冻错曲冰湖泥石流沟受常年积雪影响,流域内松散物源以沟道物源、坡面物源、崩滑物源、冻融物源为主(图2)。沟域内可启动物源共264处,总静储量约9.5 × 107m3,总动储量约1.7 × 106m3。冻融物源及冰碛物稳定性较好,只有在大型泥石流铲刮下启动,相比之下崩滑物源和沟道物源在暴雨条件下更易于启动[24]。

1.3 水文条件

研究区所在区域属温带半干旱季风气候区,日照充足,干湿季节分明,多年平均降雨量490.23 mm,主要集中在6—9月,最大日降雨量39.2 mm。该流域年平均气温4.5～6.4 ℃,最高气温30.1 ℃,多出现在6—8月,多年平均无霜期128 d。冻错曲流域属怒江水系,水源主要来自大气降水、冰川融水和地下水,沟内流水常年汇入冻错曲,流量随季节变化较大。

2 泥石流动力学特征参数

以拟建工程位置作为计算断面,根据规范公式计算冻错曲冰湖泥石流在不同降雨频率下的动力学特征参数。

2.1 泥石流流速

根据《泥石流防治工程勘察规范(试行)》(T/CAGHP 006—2018),冻错曲泥石流的易发性评分为83分,根据查表法得泥石流容重为1.572 t/m3,属稀性泥石流,故采用式(1)计算泥石流流速。

(1)

式(1)中:VC为泥石流断面平均流速,m/s;γH为泥石流固体物质容重,t/m;φ为泥沙修正系数,据查表法取0.54;1/n为清水河床糙率系数,取1.2;R为水力半径,一般可用平均泥深H代替,m;I为泥石流水力坡度,取全流域平均纵坡降26.12‰。

计算结果如表1所示。

表1 冻错曲泥石流拟建工程处流速Table 1 Flow velocity of the Dongcuoqu debris flow at the location of the proposed construction site

2.2 泥石流峰值流量

假设泥石流与暴雨同步发生,且暴雨流量完全转化为泥石流流量。研究表明,冰湖溃决型泥石流对比暴雨型泥石流具有流量放大效应,可将冰湖溃决考虑为沟道堵塞溃决的一部分,采用放大堵塞系数的雨洪法[25-26]。首先获取不同降雨频率下计算断面的暴雨洪水设计流量,再通过采用合适的堵塞系数对泥石流流量进行快速评估。根据式(2)计算泥石流峰值流量。

QC=(1+φ)QDDC

(2)

式(2)中:QC为泥石流峰值流量,m3/s;QD为暴雨洪水设计流量,m3/s;DC为泥石流堵塞系数,取5.5。

计算结果如表2所示。

表2 冻错曲泥石流拟建工程处峰值流量Table 2 Peak flow of the Dongcuoqu debris flow at the location of the proposed construction site

2.3 一次泥石流过程总量

由于泥石流比一般洪水更具暴涨暴落的特点,一次泥石流过程一般均比较短,泥石流过程线可以概化成五角形,故一次泥石流过程总量按照式(3)进行计算。

Q=0.264TQC

(3)

式(3)中:Q为一次泥石流的总量,m3;T为泥石流历时,s。

一次泥石流冲出的固体物质总量按式(4)计算。

QH=Q(γC-γw)/(γH-γw)

(4)

式(4)中:QH为一次泥石流固体物质冲出量,m3;γc为泥石流容重,t/m3;γw为清水容重,t/m3。

计算结果如表3所示。

表3 冻错曲泥石流一次冲出总量Table 3 Total volume of run-out materials of the Dongcuoqu debris flow

2.4 泥石流冲击力

泥石流冲击力是泥石流防治工程设计的重要参数,包括(流体)整体冲压力和(石块)单块最大冲击力。采用式(5)计算泥石流整体冲压力。

(5)

式(5)中:δ为泥石流整体冲压力,kPa;λ为受力物形状系数,取1.33;g为重力加速度,m/s2;α为受力面与泥石流冲压力方向的夹角。

考虑石块正面撞击下游桥墩,单块最大冲击力按式(6)计算。

(6)

式(6)中:F为石块冲击力,kN;γ为动能折减系数,取0.3;W为石块质量,t;α为石块运动方向与受力面的夹角;C1、C2为石块、桥墩的弹性变形系数,取C1+C2=0.005。

计算结果如表4所示。

表4 冻错曲泥石流整体冲击力及单块最大冲击力Table 4 Impact pressure and maximum impact force of rock mass of the Dongcuoqu debris flow

3 冻错曲冰湖发育特征

3.1 冰湖基本特征

冻错曲冰湖距泥石流沟沟口约8.3 km,东侧湖口地理坐标:30°18′18.26″N,96°4′0.41″E。湖面呈长条状,轴向长度7 100 m,平均宽度 290 m。非汛期湖面高程为 4 055 m,汛期湖面高程为 4 060 m,常年处于溢流状态。经过对Landsat-8卫星影像解译,冻错曲冰湖面积约2.05 km2。

冰湖库容是反映冰湖蓄能大小与冰湖潜在危害性的重要指标,既能为冰湖危险程度初步评价提供参考,也是峰值流量及溃决洪水演进的估算基础。目前尚不能通过遥感等手段直接测量冰湖库容,一般采用经验公式间接估算。Huggel经验公式样本取自与青藏高原地质条件相近的阿尔卑斯山脉,且公式误差波动较小[27-28],因此采用式(7)估算冻错曲冰湖平均深度,进而计算冰湖库容,计算结果见表5。

表5 冻错曲冰湖基本参数Table 5 Basic parameters of the Dongcuoqu glacier lake

(7)

式(7)中:V为冰湖库容,m3;A为冰湖面积,m2;D为冰湖平均深度,m。

3.2 冰湖堰塞体稳定性分析

2018年11月及2019年8月两次现场踏勘结果表明,冻错曲冰湖堰塞体为两岸大型崩滑堆积及右岸支沟泥石流堆积综合形成,坝体主要以孤、块石为主,其左岸大型崩塌堆积体结构松散、植被稀少或无,而右侧崩塌及局部泥石流堆积体上有低矮灌木分布。野外调查结果及遥感影像显示,堰塞体高度约45 m,顺沟长度约458 m,坝体顶宽约320 m,估算其总体积为6.59×106m3(图4)。

图4 冻错曲冰湖堰塞体概况Fig.4 Overview of the Dongcuoqu glacier-lake barrier dam

冻错曲冰湖常年处于稳定渗流状态,在坝顶处无水漫流,但在坝体内存在一个较为稳定的渗流通道,沟内可以观察到在坝顶下游约400 m开始出现稳定水流沟道,说明该坝体具有稳定的渗透浸润线,因此对冻错曲冰湖的溃决风险评价尤其重要。冰湖溃决风险评价的主要内容是对冰湖堰塞体的稳定性评价[12],在缺乏坝体材料土石强度参数以及土力学参数的情况下,常采用堰塞坝稳定性快速评价方法[29-30]。钟启明等[31]总结了国内外学者基于形态学指标,提出的5种快速评价堰塞坝稳定性的经验公式,并对世界范围内421 个堰塞坝实例进行了评价,结果表明DBI和逻辑回归模型指标Ls(AHV)的准确性较好。DBI较Ls(AHV)错判率更低,综合考虑了包括漫顶和管涌等各种形式坝体溃决,且样本具有广泛性和代表性,因此采用DBI来判别冻错曲冰湖堰塞体稳定性。DBI表达式为

(8)

式(8)中:Ab为堰塞坝控制流域面积,km2;Hd为坝体高度,m;Vd为堰塞坝体积m3。

通常认为当DBI<2.75时,坝体稳定;当2.753.08时,坝体不稳定。

基于遥感影像划定堰塞体控制流域面积,结合现场调查结果对其高度和体积进行估算。计算得到冻错曲冰湖堰塞体的DBI为3.01,位于非稳定区与稳定区之间,从图5中也可以发现有为数不少的不稳定堰塞坝在缓冲区甚至稳定区,同样对于稳定堰塞坝也有此类情况出现,因此不排除冻错曲冰湖有溃决的可能性[14]。

图5 冻错曲冰湖堰塞体DBI判别图[14]Fig.5 DBI discriminant diagram of the Dongcuoqu glacier-lake barrier dam[14]

3.3 溃口峰值流量计算

冰湖溃决溃口峰值流量的计算结果通常与堰塞坝的形态特征、坝前水深等因素有关,计算方法则主要根据溃决模式有所不同。对于冻错曲冰湖,分别在瞬时部分溃决和瞬时全部溃决2种溃决模式下,对其溃口峰值流量进行计算。

刘建康等[32]统计早期西藏地区冰湖溃决记录,认为洪水一般在下泄1/3库容时出现峰值流量。因此亦采用1/3溃决模式对冻错曲冰湖进行瞬时部分溃决计算,并且考虑极端工况下冻错曲冰湖全部溃决模式作为参考。对于峰值流量的计算,国内外不少学者对此进行研究,并提出了许多计算公式,但是这些公式中绝大部分都是基于冰湖库容这一单一因子对峰值流量估算[33]。这里分别采用Froehlich针对溃口高度提出的式(9)计算冻错曲冰湖瞬时部分溃决溃口峰值流量qp[34],以及谢任之[35]利用特征线法归纳出的式(10)计算瞬时全部溃决溃口峰值流量Qp,表达式为

qp=0.607V0.295h1.24

(9)

(10)

式中:V为冰湖库容,m3;h为溃决口高度,m;λ为流量参数,λ=m2m-0.5[2/(2m+1)]2m+1,m为沟道断面指数,取1.25;B为堰塞体宽度,m;g为重力加速度,m/s2;H1为溃坝前最大水深,取冰湖平均深度,m。

表6 不同溃决模式溃口峰值流量Table 6 Peak flow of different dam failure modes

(11)

(12)

式中:k为洪峰流量系数,k=1+(γC-γw)/(γH-γc);γc为泥石流容重,t/m3;γw为清水容重,t/m3;γH为泥石流固体容重,t/m3。

4 冰湖溃决泥石流数值模拟

4.1 RAMMS

RAMMS采用改进的VS(Voellmy-Salm)单相流连续介质模型,泥石流流体被假设为一种非稳定以及非均质流体,通过双参数Voellmy摩擦模型来描述流动岩屑之间的摩擦行为[37]。研究表明,VS模型能够准确模拟泥石流的运动,并已被广泛运用于泥石流的动力学研究[38-39]。此外,RAMMS还能够将结果导出到GIS、修改地形数据(例如设置构筑物)、增加附加参数(屈服应力等),并增强泥石流的预测效果[40]。

VS模型假设剪切变形集中在流动基面附近[37],将摩擦阻力分为两部分:一是静态摩擦阻力,用与法向应力成比例的干库仑型摩擦因数μ表示,不依赖于速度;二是运动阻力,与速度平方阻力或黏性湍流摩擦因数ξ相关。这两项都提供了流动阻力,可以综合模拟高速区域和速度较小的尾部的流动行为。改进后的总摩擦阻力S表示为

(13)

式(13)中:N为法向应力,由ρhgcosφ表示;ρ为密度,kg/m3;g为重力加速度,kg/m2;φ为流体材料内摩擦角,(°);h为流体高度,m;u为流体速度,m/s;N0为流体材料屈服应力。

当法向应力N较小(低流体高度h)时,剪应力S的数值将迅速从0变为N0;当法向应力N较大时,屈服应力N0会使剪应力S增大,从而使泥石流提前停止。

4.2 模型设置

RAMMS为泥石流模拟提供物源释放和水力释放两种方式。水力释放基于流量和时间的关系,通过在给定位置处释放一个流量对泥石流进行模拟,计算准确度较好,该方法可以显著减少模拟时间(只需计算冰湖下游而不是全流域),且溃决洪水流量曲线恰好可以作为数值模型的输入以进行下一步模拟,相比物源释放更具优势。因此,基于冻错曲冰湖堰塞体稳定性评价结果,将不同溃决模式下冰湖溃口峰值流量输入RAMMS,分别得到2种场景下的冰湖泥石流溃决演进过程:瞬时部分溃决和瞬时全部溃决。

考虑计算机的最佳计算能力,将ASTER GDEMV3获取的30 m分辨率DEM处理后作为地形数据,基于网格进行数值分析。根据RAMMS水文手册,可以将泥石流流量-时间曲线概化为三点式水文过程图,即通过定义泥石流峰值流量Qmax及其对应时刻t1、结束时间t2进行水力释放,具体参数见表7。根据野外调查结果和物理原理计算进行参数率定,最终模型基本参数设置为:密度为1 572 kg/m3,重力加速度取9.8 m/s2,干库仑型摩擦因数为0.18,黏性湍流摩擦因数为192 m/s2。

表7 RAMMS水力释放参数Table 7 Parameters of hydrology in RAMMS

4.3 溃决演进过程分析

根据RAMMS模拟结果,冻错曲冰湖泥石流在2种溃决模式下的演进过程相似,根据时序和流体运动特点,可归纳为4个阶段:初始溃决阶段、加速运动阶段、减速运动阶段、沟口停淤阶段。

在初始溃决阶段,冰湖堰塞体在极端地震或强降雨等外动力地质作用下发生破坏,冻错曲冰湖开始发生溃决,湖水倾泻而下,此时瞬时部分溃决和瞬时全部溃决溃口处最大流深可分别达到34.62 m和50.96 m[图6(a),图7(a)]。在加速运动阶段,泥石流铲刮沿程沟道,夹带沟道内泥沙、碎石等松散固体物质向下游运移,泥石流在下游沟道狭窄段流速相对更大;对比到达沟口的时刻,部分溃决相对于全部溃决具有一定的滞后性[图6(b),图7(b)]。冻错曲冰湖下游主沟纵坡降仅有21.9‰,冲出沟口后,泥石流流速逐渐变小,进入减速运动阶段,由于溃决洪水总量较大,存在一个持续冲出沟口的过程[图6(c),图7(c)]。 随着溃决过程的结束,冲出的泥石流在沟口漫流后自然减速停淤,RAMMS会通过泥石流流体低通量平衡自动停止,冻错曲冰湖瞬时部分溃决和瞬时全部溃决演进过程分别历时5 h和3.5 h[图6(d),图7(d)]。

图6 冻错曲冰湖泥石流瞬时部分溃决演进过程Fig.6 Evolution process of instantaneous partial outburst of the Dongcuoqu glacier lake debris flow

图7 冻错曲冰湖瞬时全部溃决演进过程Fig.7 Evolution process of instantaneous complete outburst of the Dongcuoqu glacier lake debris flow

4.4 危害性评价

根据2种溃决模式下泥石流的影响范围及特征参数可以评价其对下游工程的影响。如图6(e)、图7(e)所示,冻错曲冰湖瞬时部分溃决和瞬时全部溃决堆积范围均经过工程位置,其中全部溃决的影响范围更大;潜在威胁区泥石流平均深度分别达4.87 m和8.26 m,最大流速分别为5.29 m/s和5.74 m/s[图6(f)、图7(f)],将对拟建工程构成一定影响。考虑到拟建工程仅从泥石流堆积区边缘通过,可通过提高路基和桥梁标高以保留一定的安全高度,并采取排导槽或导流堤等防护措施对泥石流进行引流和拦截,避免其直接对拟建工程直接冲击。由于冰湖溃决泥石流致灾的复杂性与随机性,在进行防治工程设计时,需特别注意对模拟结果的解读和使用。

5 结论

(1)冻错曲冰湖泥石流沟流域面积226.27 km2,主沟长度31.44 km,相对高差1 365 m,平均纵坡降27.27‰。流域内松散物源丰富,沟内流水常年汇入冻错曲冰湖。考虑冰湖溃决流量放大效应,100年一遇冻错曲泥石流在拟建工程处洪峰流量为2 509.10 m3/s,流速为4.99 m/s,整体冲击力为53.18 kPa。

(2)冻错曲冰湖面积2.05 km2,平均深度46.56 m,库容9.54 × 107m3。无量纲堵塞指数(DBI)分析结果显示冻错曲冰湖堰塞体位于非稳定区与稳定区之间,存在发生溃决的可能性。运用RAMMS的Voellmy-Salm单相流模型对冻错曲冰湖泥石流进行模拟,结果表明其溃决演进过程可归纳为初始溃决、加速运动、减速运动、沟口停淤4个阶段。

(3)在2种溃决模式下,冻错曲冰湖泥石流在潜在威胁区平均深度分别为4.87 m和8.26 m;在瞬时全部溃决场景下,拟建工程处最大流速为5.74 m/s,峰值流量为2 843.38 m3/s,与规范公式计算结果接近。通过提高路基和桥梁标高以保留一定的安全高度,并采取排导槽或导流堤等防护措施对泥石流进行引流和拦截,避免其直接对拟建工程直接冲击。