强震区泥石流防治工程效果数值模拟分析

唐 海, 唐 川, 陈 明, 王晓迪, 罗玉婷

(成都理工大学 地质灾害防治与地质环境保护国家重点实验室, 四川 成都 610059)

泥石流是山区常见的自然灾害之一，水源和丰富的松散物源是引发泥石流灾害的重要因素[1]。汶川地震后，震区地质环境更加脆弱，崩塌滑坡大量发育，植被破坏严重，水土加剧流失，为泥石流发活动供了丰富的物源，震后暴雨诱发的群发性及特大泥石流灾害的频繁发生[2-4]。根据据黄润秋等[5]的研究成果，汶川地震后泥石流将在20 a内处于活跃状态。随着大量泥石流治理工程的兴建，地震灾区的大多数泥石流沟已得到了相应的工程治理，但有的防治工程没能达到预期的防治效果，最典型的当属2010年8月13日文家沟泥石流。“5·12”汶川地震后文家沟多次发生泥石流灾害，为预防泥石流的发生，2010年7月，在沟内修建了19座谷坊坝群，41座潜坝和一座长215 m高8 m拦挡坝[6]。2010年8月13日，文家沟流域遇大暴雨引发文家沟特大山洪泥石流，已修建的防治措施在在泥石流的冲击下直接被摧毁，造成清平乡学校、医院等400余户房屋受损，桥梁、公路被毁，6000多人受灾，直接经济损失达6亿元[7]，损失惨重。为了更好地预防泥石流灾害的发生，进一步研究泥石流的防治效果很有必要。考虑到采用数值模拟对防治工程减灾效果测验的经济性与简便性，可以对治理工程的防治效果进行数值模拟研究，以确保治理工程的安全性与可靠性。

FLO-2D是一款适用于模拟泥石流灾害的数值模型，能较完整的分析泥石流的运动状态特征，尤其在小流域泥石流的应用上。在运用FLO-2D模拟泥石流方面，杜雪剑[8]采用FLO-2D对红椿沟在两种重现周期下的爆发过程进行模拟研究，得到了泥石流的运动和堆积特征，杨海龙[9]应用FLO-2D对都江堰银洞子沟泥石流模拟和工程治理效果评价，王俊[10]利用FLO-2D软件对联合溪沟泥石流的运动、淹没过程进行模拟，在此基础上提出综合防治建议，杨涛[11]、刘鑫磊[12]运用FLO-2D对不同降雨频率下的溃决型与非溃决型泥石流进行模拟研究。在国内FLO-2D 模型已被广泛的运用于泥石流的模拟并取得了较好的效果，因此，本文采用FLO-2D软件，以2010年龙池镇黄央沟“8·13”泥石流为原型，结合雨量数据、数字地形数据及泥石流的各项流变参数，模拟在不同的降雨频率下泥石流的运动过程和堆积范围，对泥石流的防治工程效果进行分析评价，为今后泥石流的研究和防治提供借鉴。

1 研究区概况

黄央沟位于都江堰市龙池镇龙溪河流域，沟口地理坐标为103°33′39.0″E，31°04′16.5″N。黄央沟由一条主沟和一条小次级冲沟组成，流域面积仅0.68 km2，沟长约1.95 km，纵坡平均坡降435‰。黄央沟总体而言较为顺直，坡度大，中上游为宽谷地貌，呈U字形，两侧谷坡坡度45°～75°，下游沟谷较狭窄，呈V字形，两侧谷坡坡度40°～60°(图1)。“5·12”地震后沟内山体破碎，斜坡稳定性变差，大量崩塌滑坡堆积物堆积在沟谷内，泥石流物源量剧增，致使黄央沟在暴雨后已发生多次不同程度的泥石流。根据气象资料显示，2010年8月13日龙池镇遭遇特大暴雨，最大1 h降雨量达75 mm，为20 a一遇暴雨，直接造成龙溪河流域暴发群发性泥石流，黄央沟也在其中，仅数天之后，8月18日龙池再次遇强降雨，黄央沟再次暴发泥石流，这两次泥石流导致龙溪河河道被严重挤压，河床被压迫迁移30～40 m。2013年7月9日龙池镇降雨量达185.4 mm，9日晚间黄央沟再次暴发泥石流。由此可见，有利的地形和丰富的物源为泥石流提供了良好的发育条件，在强降雨作用下，黄央沟易暴发泥石流。

图1 龙池镇黄央沟流域示意图

2 FLO-2D理论基础

2.1 适用条件

FLO-2D为O’Brien(1988)所设计的用于模拟洪水与泥石流的软件，其运动控制方程的求解是利用非牛顿流体模型与中央有限差分数值方法。受理论模型的限制，FLO-2D在模拟有以下假设和限制[13-14]： ①差分时间间隔内为稳定流； ②水道断面和粗糙度是规则的； ③净水压力分布； ④每一个网格有单一的高程及N值； ⑤无法模拟水道的刷深现象； ⑥无法模拟上游崩塌可能对泥石流造成的影响； ⑦假定泥石流流通区为定沟床模式，故无法模拟沟道的侵蚀现象。

FLO-2D适用于泥石流运动状态的模拟，模拟结果可以得到泥石流的流速及流深，但不能模拟泥石流的侵蚀过程，也不能模拟泥石流对工程结构造成的损毁情况，因此本文对已修建工程防治措施后的模拟中未考虑工程措施的损毁情况。

2.2 控制方程

FLO-2D的运动控制方程是利用非牛顿流体模型与中央有限差分数值方法求解，从而得到泥石流在水平方向和垂直方向的流速、流深、冲於深度，其连续方程和运动方程为[15]：

(1) 连续方程：

(1)

式中：t——泥石流持续时间；h——流动深度；I——水力坡降；u，v——水平方向和垂直方向上的平均流速。

(2) 运动方程：

(2)

(3)

式中：Sfx，Sfy——摩擦坡降；Sox,Soy——沟床坡降。

(3) 根据O’ Brien提出二项流变模型，将泥石流的流体总应力τ分为5项应力之和：

τ=τc+τmc+τv+τt+τd

(4)

式中：τc——凝聚型屈服应力；τmc——莫尔—库伦剪应力；τv——黏滞剪应力；τt——紊流剪应力；τd——扩散剪应力。

3 模拟参数选取

3.1 地形参数

在模拟之前收集了黄央沟地区的高精度遥感影像图和1∶1万等高线，利用ArcGis软件将修正后的等高线转化为DEM，然后将DEM的文件格式转换生成FLO-2D模型中能够识别的ASCⅡ高程文件。将高程文件导入FLO-2D后，然后对黄央沟的地形进行网格划分。网格的大小决定了模拟时间和模拟结果的准确性，经过多次的模拟调整，最终将网格划分为5 m×5 m大小，这既能兼顾模拟效率，又能满足模拟精度需要。

3.2 体积浓度和重度

泥石流的体积浓度与流体性质、堆积规模和堆积形态密切相关，通常要确定泥石流的体积浓度要到泥石流发生地采集样本，但泥石流发生时的样本很难取到，因此参考了王纳纳[16]、程思[17]在模拟时黄央沟时体积浓度的取值，结合黄央沟泥石流勘察报告[18]和FLO-2D使用手册，取降雨频率P=5%,时体积浓度为Cv=0.57，P=2%时Cv=0.6。根据泥石流特征及野外数据，确定黄央沟泥石流重度为1.67 g/cm3。

3.3 流量过程线

在应用FLO-2D进行模拟时，需要输入的是清水流量过程线，参照《四川省中小流域暴雨洪水计算手册》，按水文方法计算出不同降雨频率下的暴雨洪峰流量，然后运用简化的五边形法则绘制出黄央沟泥石流清水流量过程线。但清水流量并不能代表泥石流流量，泥石流流量应是清水流量乘以膨胀因子BF，BF=1/(1-Cv)，其中，Cv表示体积浓度。不同降雨频率下的暴雨洪峰流量和泥石流洪峰流量详见表1。

在运用FLO-2D进行泥石流的数值模拟时，集水点的选取至关重要。黄央沟流域在“5·12”大地震后不良地质现象发育，沟内多处发生规模较大的崩塌滑坡，致使大量的碎块石滞留于沟谷内，为泥石流的发生提供了大量的松散物源，是泥石流起动的最佳场所，所以本文选择松散物源聚集处作为泥石流的集水点。

表1 不同降雨频率下峰值流量

3.4 防治工程的设置

黄央沟在2010年“8·13”，“8·18”泥石流之后进行了工程治理，主沟道内修建了2座拦砂坝，工程平面分布图和拦砂坝工程数据见图1和表2。在运用FLO-2D进行数值模拟时，可以将坝体插入到数值模拟的网格中，然后对坝体高度、方向和位置进行赋值，再导入其他各项参数之后便可进行模拟。

表2 黄央沟拦砂坝工程数据

3.5 其他参数

在前人的模拟研究的基础上[19]，本文模拟参数选取α1=0.811，α2=0.00462。β1=13.72，β2=11.24。曼宁系数n、层流阻滞系数K根据FLO-2D使用手册建议选取n=0.1，K=2285；本文泥石流的模拟时间1.3 h，与黄央沟“8·13”泥石流运动的总时间一致。

4 模拟结果及分析

本文先研究了黄央沟泥石流在自然条件下的运动过程，然后针对黄央沟存在防治工程的情况展开模拟研究，分析黄央沟在有工程防治措施的情况下泥石流的冲出范围、淤积厚度及坝后的回於情况，以此来分析评价工程防治措施的效果。黄央沟在自然条件和有工程防治措施两种工况下的数值模拟结果如图2所示。

图2 黄央沟数值模拟结果

在自然条件下，当降雨频率P=5%时，沟口堆积区面积为1.29×104m2，最大堆积深度3.45 m，平均堆积深度2.03 m，冲出固体物质总量2.62×104m3；当降雨频率P=2%时，沟口堆积区面积为1.54×104m2，最大堆积深度3.97 m，平均堆积深度2.3 m，冲出固体物质总量达3.60×104m3。由以上模拟结果可以看出，泥石流的堆积区面积、堆积深度及冲出量与降雨频率呈正相关，P=2%的冲出量是P=5%的1.37倍，且当降雨频率为20 a一遇时，泥石流已造成龙池景区旅游公路中断，并且大量固体物质冲入龙溪河造成堵塞。

在有治理工程的情况下，因治理工程的拦挡作用，泥石流所裹带的大部分的固体物质停留在坝后，拦挡坝后出现明显的回於现象，泥石流在大坝后的淤积深度明显大于沟道的其他位置。当流体深度接近坝高时，泥石流漫过坝顶继续向下运动，到达沟口后速度逐渐降低，最终冲出沟口形成堆积扇。与自然条件下相比，其冲出量、堆积范围、堆积深度都明显降低。当降雨频率P=5%，沟口堆积范围2 100 m2，最大堆积深度2.13 m，平均堆积深度1.1 m，冲出量2 300 m3，与同周期无工程措施相比，堆积范围缩小6.14倍，固体物质冲出量减小11.39倍。当降雨频率P=2%时，堆积区面积为5 700 m2，最大堆积深度2.87 m，平均堆积深度1.5 m，固体物质冲出量体积8 600 m3，相比自然条件下减小了4.18倍，冲出范围减小了2.7倍。两种工况模拟结果详见表3。总体来说，拦挡坝充分发挥了拦蓄泥石流、回於和稳沟固坡的作用，减少了物源向下游的启动量，泥石流的冲出量、冲出规模也到了有效的控制，避免了对龙溪河的堵塞，保证了龙溪河下游龙池场镇的安全。

为验证模拟结果的准确性，以黄央沟已暴发过“8·13”泥石流为例，引用模拟堆积区与实际堆积区的重叠面积对模拟结果的准确性进行评估[20]，相应的计算公式为：

式中：F——模拟精度；So——模拟堆积区与实际堆积区的重叠面积；Sm——模拟堆积区面积；Sr——实际堆积区面积。

对比验证2010年“8·13”泥石流暴发后的实际调查情况与数值模拟结果(表4)。模拟精度78.1%，模拟结果与实际基本吻合，说明了模拟参数的选取是合适的，这为不同降雨频率及有工程措施条件下模拟结果的准确性提供了保障。

表3 不同工况条件下模拟结果统计结果

表4 模拟堆积区与实际堆积区对比

5 讨 论

为了研究黄央沟的防治工程对泥石流的拦挡效果，分别统计了两种降雨频率下1号坝和2号坝后泥石流的回於范围、回於总量和平均回於厚度(表5)。根据表5的统计结果可以看出，拦砂坝后的淤积总量远远大于设计库容，超出的部分固体物质淤积在沟口处，这对沟口处龙池景区的唯一1条旅游公路和沟口居民产生了威胁，一旦爆发泥石流，冲出的松散固体物质将堵断公路，给过往行人车辆造成不便，切断了人员和物资的出入，给当地旅游经济发展带来重大影响。事实上，黄央沟的防治工程修建完成后，2013年7月9日遇强降雨再次暴发泥石流，拦挡坝并没能将泥石流彻底拦截在坝后(图3)，仍有固体物质漫过坝顶后冲出沟口，於埋了沟口的公路(图4)。而且由于黄央沟沟道狭窄，两侧斜坡陡峭，拦砂坝在於满后大型机械无法大量使用，这给清淤工作带来了很大的难度。但若不清淤，大坝则失去了防灾减灾的效果，为后期泥石流的活动埋下了严重的隐患。另一方面，虽拦砂坝的回於效果有利于减缓沟道纵坡比降、降低沟道侵蚀速率的作用[21]，但黄央沟2座拦砂坝的回於长度约100 m，而黄央沟的主沟长1.95 km，平均纵坡降435‰，所能减缓的沟床纵坡降的长度有限，相对于整个沟长而言微不足道。

表5 黄央沟数值模拟坝后堆积特征统计结果

图3 泥石流冲出2号拦砂坝

图4 泥石流於埋沟口公路

6 结 论

(1) 在自然条件下，当降雨频率P=5%时，黄央沟就会堵塞龙溪河，增加龙溪河的泄洪压力，降雨频率减小，泥石流冲出量、冲出范围、堆积深度随之变大。对比验证了黄央沟“8·13”泥石流实际结果，模拟精度为78.1%，模拟结果准确性较高。

(2) 对两种重现周期下泥石流的防治工程效果分析，结果表明：当降雨频率P=5%时，泥石流的冲出量减小11.39倍，堆积范围减小6.14倍；当降雨频率P=2%时，泥石流的冲出量减小4.18倍，堆积范围减小2.7倍，相比于自然条件泥石流的冲出量和冲出规模得到了一定的控制，不会堵塞龙溪河，保证了下游龙池场镇的安全，已有的治理工程有利于整个沟域泥石流的防治。

(3) 黄央沟的防治工程虽取得了一定的防灾减灾效果，但拦挡坝的设计库容偏小。对坝后泥石流堆积量的统计表明，当降雨频率P≤5%，坝后泥石流堆积量远大于大坝的设计库容，2座拦砂坝并不能将泥石流彻底拦截在坝后，超出的固体物质会影响沟口旅游公路的正常通行。