基于Visual Modflow的关门山沟堰塞坝渗流稳定模拟研究

李明忠,林 妮

(长江勘测规划设计研究有限责任公司,湖北 武汉 430010)

0 引 言

滑坡、崩塌、泥石流堵塞天然河道,形成壅塞体,上游来水在库内蓄积,形成堰塞湖[1]。堰塞湖水位不断上涨,将首先造成淹没灾害,水位继续上涨,渗流作用加强,可能出现管涌,进一步发展可能导致堰塞湖溃决。溃决后,库水宣泄,下游洪水可相当于千年一遇至万年一遇[2],将毁坏沿途城镇、耕地,造成大量生命财产损失,同时,溃坝洪水将可能引起新的滑坡、崩塌等次生灾害,进一步加剧损失[3]。

“5·12”汶川地震以来,国内学者开始关注堰塞湖的渗流问题,徐文杰等研究了肖家桥堰塞湖渗透稳定性,通过出逸流量及高程差,反算流量系数,推求了肖家桥堰塞湖水位在743 m和761 m时的渗透坡降,二者均小于材料允许比降,在缺少地质勘察的情况下,初步判断肖家桥堰塞坝不会出现渗透破坏[4]。严祖文等[5]计算了唐家山堰塞湖非稳定渗流和稳定渗流两种工况的渗流稳定性,得出了唐家山堰塞坝渗透破坏可能性较小的结论。石振明等[6]从渗流的基本规律、数值方法、模型试验等方面总结了堰塞坝组成材料渗流的机理和研究方法,指出土体渗透变形和渗透稳定性由几何条件和水力条件共同决定。胡卸文等[7]对唐家山堰塞湖渗流问题进行了数值模拟,得出结论:下游坝脚处某些部位会发生零散渗透破坏,但对于堰塞坝的整体稳定性影响不大。徐轶等[8]对堰塞湖应急处置工程措施及典型案例进行分析,较全面地总结梳理了堰塞湖应急处置工程措施的经验及技术发展趋势。

本文基于关门山沟堰塞坝实测土层物理参数,采用Visual Modflow可视化三维地下水模型软件模拟渗流场,分析渗透坡降变化规律,判断在各模拟水位条件下,堰塞坝是否发生渗透破坏及渗透破坏形式。根据渗透坡降变化规律,提出发生渗透破坏时堰塞湖临界水位,对堰塞湖渗流稳定性研究具有理论和实际参考价值。

1 关门山沟堰塞湖概况

关门山沟堰塞湖位于都江堰市白沙河上游关门山沟,“5·12”汶川特大地震时右岸百余米高的山体在强震作用下形成滑坡,滑坡体快速冲击左岸基岩,淤积河道形成堰塞湖。堰塞坝下游河底高程1 648.0 m,坝顶高程1 728.0 m,上游水位1 705.0 m,下游坝脚水位1 649.0 m。堰塞坝高80 m,宽230 m,顺河长约500 m,不计算右岸滑坡体部分,估算土石方量2.7×107m3,蓄水量约370万m3。堰塞坝主要由碎石土组成,结构松散,滑坡陡壁高150 m,坡度50°～60°,右岸边坡尚存大量物源,有再次下滑的可能。堰塞坝以上主河道长约11 km,河道坡降13.5%,堰塞湖以上集雨面积约56 km2。根据SL 450-2009《堰塞湖风险等级划分标准》,关门山沟堰塞湖风险等级为Ⅱ级。

关门山沟堰塞湖所在区域海拔高程在1 200～3 750 m之间,地貌特征主要为构造剥蚀型,其次为侵蚀堆积地貌。堰塞湖处于扬子地台与松潘-甘孜地槽区之间的构造过渡带上,位于龙门山构造带北段,龙门山主中央断裂带附近。根据四川防震减灾信息网公布的地震烈度分布图,“5·12”汶川地震对工程场地的影响烈度为Ⅹ度。堰塞坝主要由碎石土构成,通过钻孔及物理指标试验,获取碎石土物理参数见表1。

表1 堰塞坝物理参数

对于原河道覆盖层,由于没有进行物理指标试验,参考白沙河上游地区相关文献确定其参数。根据太沙基公式,综合考虑计算成果及实验测量,确定关门山沟堰塞坝允许坡降取值0.40～0.55,河道覆盖层取值0.47～0.60。

2 三维渗流模型

Visual Modflow由加拿大Waterloo水文地质公司开发研制,基于有限差分法理论进行计算,作为三维地下水水流和溶质数值模拟评价的标准可视化专业软件系统,运用广泛[9]。基于常密度的地下水三维流动基本方程:

(1)

式中:Kxx,Kyy,Kzz分别为渗透系数在x,y,z方向的分量;H为作用水头;w为单位时间内单位体积流进或流出的水的体积通量;Ss为储水系数,表示含水层地下水头变化时,由于含水层垂向压缩和地下水弹性膨胀,单位体积含水层变化的水的体积;t为时间。

(1) 第一类边界条件。在边界上,各时刻每一点水头是给定的:

(2)

式中:H(x,y,z,t)为边界S1上点(x,y,z)在t时刻的水头,f1(x,y,z,t)为S1上的已知函数。

(2) 第二类边界条件。已知某一部分边界单位面积上流入(流出)的流量,表示为

(3)

式中:n为边界S2的外法相方向,q1为边界z单位面积的补给量,为已知函数。

Modflow采用迭代法求解渗流差分方程,在计算之前,需要给定一个初始水位值,并求解下一个时间段的水位值;所求得的水位值作为下一步初始水位,再次迭代,直到两次迭代结果相差足够小为止[10]。

采用Visul Modflow,建立关门山沟堰塞湖三维渗流模型。取两岸山体及河床下伏基岩为相对不透水层,上下游作为自由流动面,根据渗透系数对堰塞坝进行分区,分区内考虑土体性质均一,在渗流模拟过程中,土体物理力学性质不发生变化,假定堰塞湖上游水位和下游出水口处水位不变,并作为定水头边界处理。

2.1 网格划分及参数给定

本次模拟选用稳定渗流,堰塞湖顺河长500 m,宽230 m,堰塞坝三维模型见图1。模型平面上划分为200×100的网格,将堰塞坝概化为均质,纵向划分为堰塞坝及原河道覆盖层两层[11]。第一层渗透系数赋值为6.5×10-5m/s,第二层渗透系数赋值为1.3×10-5m/s。第一层贮水率为1×10-4,重力给水度、有效孔隙度、总孔隙度分别为0.2,0.18和0.36;第二层贮水率取为0.5×10-5,重力给水度、有效孔隙度、总孔隙度分别取为0.03,0.085和0.12。

图1 堰塞坝三维模型Fig.1 Three-dimensional model of barrier dam

2.2 边界条件给定

堰塞坝地下水运动形式主要是上游湖水补给,下游出口排泄。假设在模拟时间段内,水位变化微小,上下游均采用定水头边界[12],堰塞坝与两岸、下伏基岩接触部位不设定边界条件,系统默认其为隔水边界,顺河方向上,上下游坝坡作为自由临空面[13]。上游水位分别取1 705,1 710,1 715 m和1 720 m,下游水位取1 649 m。

3 渗流模拟分析

对上述4种水位条件下堰塞湖渗流场进行模拟分析,见图2～5。

图2 水位1 705 m渗流场Fig.2 Seepage field diagram at 1 705 m

图3 水位1 710 m渗流场Fig.3 Seepage field diagram at 1 710 m

图4 水位1 715 m渗流场Fig.4 Seepage field diagram at 1 715 m

图5 水位1 720 m渗流场Fig.5 Seepage field diagram at 1 720 m

3.1 渗流特性分析

对4种水位条件下堰塞坝渗流场进行分析,当水位为1 705 m时,堰塞坝浸润线较低,渗水主要由堰塞坝与覆盖层分界处溢出,与现场踏勘情况相符。坝轴线以上,渗透坡降较小且变化也较小,坝轴线以下,渗透坡降变化加大,在坝脚附近达到最大。当水位上升到1 710 m时,堰塞坝浸润线提升,溢出点向上游移动,最大渗透坡降增大。随着水位继续上升、浸润线将进一步升高,溢出点进一步向上游移动,最大渗透坡降持续增大。

3.2 渗透破坏临界水位确定

堰塞坝发生渗透破坏的条件是渗透坡降达到允许坡降。堰塞坝在不同水位条件下出口处渗透坡降见表2。从表2可知,渗透坡降与水位呈正相关性,最大渗透坡降出现在下游出水口附近,堰塞坝平均坡降及最大坡降均大于覆盖层,且堰塞坝层的渗透坡降变化速度大于覆盖层。

表2 堰塞坝出口处渗透坡降

当堰塞湖水位在1 705 m和1 710 m时,两土层的最大渗透坡降均小于允许坡降,堰塞坝处于稳定状态。当堰塞湖水位达到1 715 m时,原河道覆盖层的渗透坡降仍小于允许坡降,但堰塞坝层最大渗透坡降大于允许坡降下限,一些较细的颗粒将被带走,堰塞坝局部出现渗透破坏[14]。当水位达到1 720 m时,堰塞坝层最大渗透坡降超过允许坡降上限,覆盖层渗透坡降达到允许渗透坡降上限,在坝脚处将出现较大范围的渗透破坏。如果水位进一步上升,将引起更大范围的渗透破坏,甚至形成贯穿性通道,危及堰塞坝整体安全。

最大渗透坡降与允许坡降对比见图6,根据渗透趋势,随着堰塞湖水位的上升,两土层的渗透坡降均逐渐增大。在一定范围内,水位上升主要影响堰塞坝层,当堰塞湖水位达1 712 m时,堰塞坝层最大渗透坡降达到其允许渗透坡降下限,堰塞坝局部细颗粒将被水流带走;当水位达到1 714 m时,覆盖层渗透坡降达到其允许渗透坡降下限,覆盖层局部细颗粒也将移动,水位继续上升;当水位达到1 719 m时,堰塞坝层最大渗透坡降达到其允许渗透坡降上限,堰塞坝层将出现更大范围的渗透破坏;当水位达到1 722 m时,覆盖层渗透坡降达到其允许渗透坡降上限,覆盖层也将出现更大范围的渗透破坏,严重影响堰塞坝安全。因此,将1 719 m确定为关门山沟堰塞湖渗透破坏临界水位。

图6 最大渗透坡降与允许坡降对比Fig.6 Comparative analysis of the maximum seepage gradient and the allowed gradient

随着堰塞湖水位的抬升,关门山沟堰塞坝表层碎石土发生渗透变形,稳定性降低,但坝体整体稳定。预测堰塞坝溃决模式为:下游侧碎石土层因渗透发生破坏,出水口甚至发生局部坍塌,但由于大块石骨架的存在,由渗流引起整体溃决可能性不大。随着堰塞湖蓄满,漫顶导致碎石土被侵蚀、淘刷,水流速度加大,带动进一步冲刷下切,形成以漫顶冲刷为主的大范围溃决。

4 结 论

(1) 本文采用Visual Modflow三维地下水可视化软件,建立关门山沟堰塞坝三维模型。将渗流区分成堰塞坝层和原河道覆盖层,取两岸山体及下伏基岩作为相对不透水层,假设堰塞湖水位和下游出水口水位不变并作为定水头边界处理。模拟了堰塞湖在1 705,1 710,1 715 m和1 720 m四种水位条件下的渗流场。随着水位的上升,堰塞坝渗透坡降逐渐增大,且最大渗透坡降出现在下游出水口附近。

(2) 基于渗流模拟成果,得出4种水位条件下出水口平均渗透坡降和最大渗透坡降,并与筑坝材料允许坡降进行比对。通过渗透坡降变化趋势,确定1 719 m为堰塞坝渗透破坏临界水位。