滲流和溢流導致冰磧壩潰決的數(shù)值分析

王亞妮，劉德高，董千里，申國強

（1.淮安市水利勘測設計研究院有限公司，江蘇 淮安 223005；2.水發(fā)規(guī)劃設計有限公司，山東 濟南 250014；3.承德水創(chuàng)工程項目咨詢有限公司，河北 承德 067000）

1 冰磧壩潰決成因

受氣候變化影響，新疆的冰川區(qū)常因冰湖潰決洪水（GLOF）引起洪水和泥沙災害[1]。冰湖潰決通常是由于滲漏和湖水溢出侵蝕大壩，引起冰磧壩破壞而發(fā)生的。冰磧壩相對較弱，可能會突然破裂，導致大量水和沉積物突然排放，造成下游地區(qū)的災難性洪水，對群眾的生命、財產(chǎn)、土地和基礎設施造成嚴重破壞[2]。為了避免和最大程度地減少群眾生命財產(chǎn)損失，迫切需要一種機制方法來研究冰湖潰決事件及洪水對下游的影響機制。

受到客觀條件的限制，對冰湖潰決事件及其下游洪水的研究非常有限。目前冰湖潰決的研究主要集中在衛(wèi)星觀測、歷史記錄和預防措施等方面，而對冰磧壩的破壞機理和侵蝕過程仍缺乏深入的認識[3]。學者們通過收集已有的冰湖潰決事件中的相關資料建立了經(jīng)驗模型，以預測冰湖潰決的洪峰流量與湖泊體積、深度和面積的關系，并建立了物理模型來計算冰湖潰決的特征[4]。此外，洪峰流量取決于壩體特性、破壞機制、泥沙特性等多種因素，而這些因素并不包含在經(jīng)驗關系式中。NWS-BREACH、HEC-RAS 潰壩模型等物理模型通常用于預測冰湖潰決的出流過程線，但這些模型的局限性在于沒有考慮壩面邊坡破壞、壩面切頭和側(cè)向侵蝕等實際機制[5-6]。此外，一些冰湖潰決事件的下游淹沒或影響可通過洪水演進模型進行研究，這些模型均未考慮河床泥沙顆粒的沖淤過程。

冰磧壩和滑坡壩是多種顆粒大小不等的非均質(zhì)混合物。冰磧壩與通常的天然壩并無實質(zhì)區(qū)別，冰磧壩被視為天然壩的一種。因此，即使在冰磧壩的情況下，也可以采用與滑坡壩相同的試驗和計算處理[7]。本研究建立了一個數(shù)值模型來計算由于滲流和溢流破壞冰磧壩而導致的冰湖潰決的特征，并對滲流和冰磧壩破壞進行了數(shù)值分析。為了計算冰磧壩的孔隙水壓力和大壩的邊坡穩(wěn)定性，將滲流模型和邊坡穩(wěn)定性模型合并到流量和壩面侵蝕的數(shù)值模型中，并將模擬結(jié)果與實驗結(jié)果進行比較。

2 數(shù)值模型

2.1 滲流模型

滲流是引起冰湖潰決的主要因素，Richards 方程是描述滲流過程的基本控制方程，可用于模擬非均質(zhì)土體飽和-非飽和滲流過程。因此，本文通過使用混合格式的Richards 方程計算冰磧壩非飽和土壤的孔隙水壓力變化，冰磧壩非飽和土壤孔隙水壓力變化Richards方程計算公式如下：

式中：ψ——壓力水頭，m；

Kx和Kz——x 和z 方向的水力傳導率，m/s；

C——比濕容量；

θ——土壤的體積含水量，%

x——水平空間坐標，m；

z——向上為正的垂直空間坐標，m；

t——時間，s。

蓄水系數(shù)與滲透系數(shù)的關系對于使用公式計算非飽和土壤系統(tǒng)中的水分運動非常重要。本文給出了廣泛使用的van Genuchten 本構(gòu)關系進行了冰磧壩中的滲流分析，計算了儲水系數(shù)、滲透系數(shù)和比濕容量，如：

式中：Se——有效飽和度；

θs、θr——分別為泥沙混合物的飽和含水率和殘余含水率，%；

α、η——為與土壤基質(zhì)勢有關的參數(shù)，通過使用土壤-持水曲線的曲線擬合確定；

ψ——壓力水頭；

Ks——飽和導水率，m/s；

m——1-1/η。

此外，Brooks 和Corey 提出的本構(gòu)關系給出了儲水系數(shù)，滲透系數(shù)和比濕容量的描述：

式中：ψb——空氣進入值；

λ——孔徑分布指數(shù)；

δ——經(jīng)驗常數(shù)。

通過滲流模型計算的不同案例含水率模擬值見表1。

表1 含水率模擬值匯總

2.2 邊坡穩(wěn)定性模型

許多研究人員使用簡化的Janbu 方法來計算滑移表面的安全系數(shù)Fs，如下：

式中：c——壩體材料的內(nèi)聚力，N/m2；

li——每個切片底部的長度，m；

α——每個切片底部的斜率；

Wi——每個切片包括地表水的重量，N；

uwi——每個切片底部的平均孔隙水壓力，N/m2；

φ——內(nèi)摩擦角，°。

簡化的Janbu 方法未考慮非飽和土壤中產(chǎn)生的負孔隙水壓力。負孔隙水壓力導致剪切強度增加。因此，非飽和土的抗剪強度（τ）定義為：

式中：σn——總法向應力，N/m2；

ua——孔隙氣壓，N/m2；

uw——孔隙水壓力，N/m2。

通過使用剪切強度方程方程，滑移面的安全系數(shù)可以定義為：

式中：Pi——第i 個切片的大氣壓強；

uαi——第i 個切片的孔隙氣壓；

mαi——第i 個切片底部的重量。

考慮吸力和不考慮吸力的邊坡穩(wěn)定性模型模擬的滑動面模擬結(jié)果如表2 所示。

表2 滑動面模擬結(jié)果

2.3 流動和侵蝕模型

泥沙-混合流的二維動量方程、水流的連續(xù)性方程和泥沙顆粒的連續(xù)性方程可以表示為：

式中：M——x 方向上每單位寬度的流量，m3/s；

N——y 方向上每單位寬度的流量，m3/s；

u——x 方向上的速度分量，m/s；

v——y 方向上的速度分量，m/s；

h——流動深度，m；

zb——原始床高程測得的侵蝕或沉積厚度，m；

θbx0——原始床面坡度x 分量，°；

θby0——原始床面坡度y 分量，°；

ib——侵蝕/沉積速度，m/s；

c——水流中泥沙濃度，kg/m3；

c*——河床中的最大泥沙濃度，kg/m3；

β——動量修正系數(shù)；

g——由于重力引起的加速度，m/s2；

ρT——泥沙密度，kg/m3；

τbx——x 方向的底部剪切應力，N/m2；

τby——y 方向的底部剪切應力，N/m2；

sb——河床中的飽和度；

q——滲透率，m/s。

通過考慮吸力的影響，侵蝕速度ib描述為：

式中：K——常數(shù)；

θb——河床斜率；

σ——泥沙質(zhì)量體積比；

Δτsuc——由于吸力引起的剪切應力增量，N/m2；

ρm——平衡沉積物密度；

dm——平均粒徑，mm；

c∞——平衡沉積物濃度；

ρw——孔隙流體密度。

壩面侵蝕的模擬結(jié)果如表3 所示：

表3 壩面侵蝕模擬結(jié)果

3 水槽實驗

水槽實驗是用模型模擬試驗，收集試驗參數(shù)，對工程實驗結(jié)果進行參照對比。試驗采用長500 cm、寬30 cm、深50 cm 的水槽（見圖1），壩體采用石英砂（泥沙混合物1 和2）制作，兩種類型的泥沙混合物用于不同粒度的結(jié)果比較。

圖1 水槽實驗示意圖（單位：cm）

圖2 和圖3 分別顯示了粒度分布曲線和壩體的詳細信息。泥沙混合物1 和2 的平均粒徑分別為1.4 mm 和1.04 mm，最大粒徑為4.75 mm，泥沙質(zhì)量體積比為2.65 g/cm3。

圖2 泥沙混合物粒度分布曲線

圖3 不同壩體的詳細信息

通過從湖泊上游端供給恒定的排水量來填充湖水。當水位上升漫頂而導致冰磧壩潰決時，不斷從湖的上游端向壩頂供水并侵蝕。然而，在冰磧壩滲漏潰決的情況下，通過從湖泊上游端供給恒定的排水量，湖水被填滿的深度約為16 cm。實驗條件見表4。

表4 實驗條件匯總表

4 成果分析

含水量反射計（WCR）測量的壩體中的水分運動如圖4 所示。圖5 顯示了湖泊中測得的水深隨時間的變化。van Genuchten 本構(gòu)關系的土壤參數(shù)值（泥沙混合物1 和2：α=6.3，η=3.1 和α=4.0，η=3.2）。Brooks 和Corey 本構(gòu)關系的土壤參數(shù)值（泥沙混合物1 和2：λ=2.5，ψ=-0.19 和λ=1.8，ψ=-0.19）。這均是使用兩種混合物的土壤-持水曲線擬合確定的。Ks=0.000 5 m/s 和θs=0.312 為泥沙混合物1 的測量值，Ks=0.000 25 m/s和θs=0.296 為泥沙混合物2 的測量值。

圖4 壩體中含水量反射計（WCR）1~9 的位置（單位：cm）

圖5 湖泊水深隨時間的變化

圖6 和圖7 分別顯示了案例IV 和案例V 壩體中的模擬水分分布和實驗水分分布。圖中，t=0 是湖泊上游端供水的開始時間，壩體中的水分運動是由于上游湖水的深度所致。通過使用Brooks 和Corey 關系計算的壩體水分運動速率比實驗值快。非飽和土壤中的導水率在Brooks 和Corey 關系式中高于van Genuchten函數(shù)。采用van Genuchten 函數(shù)本構(gòu)關系的模擬結(jié)果比采用Brooks 和Corey 本構(gòu)關系的模擬結(jié)果與試驗結(jié)果更吻合。儲水系數(shù)與滲透系數(shù)的關系對非飽和區(qū)域水分運動的計算非常重要。壩內(nèi)水分運動還取決于飽和導水率。泥沙混合物2 的細顆粒泥沙含量高于泥沙混合物1 的細顆粒泥沙含量，泥沙混合物2 的水力傳導率較小，水分運動較泥沙混合物1 的慢。

圖6 案例IV 泥沙混合物1 下模擬和實驗水分含量曲線的比較

圖7 案例V 泥沙混合物2 下模擬和實驗水分含量曲線的比較

圖8 給出了考慮吸力和不考慮吸力（Janbu 的簡化方法）的土體抗剪強度中滲流引起的冰磧壩潰決模擬滑動面與試驗滑動面的對比。利用van Genuchten函數(shù)關系式計算土壤中的水分運動，采用動態(tài)規(guī)劃法計算臨界滑動面。邊坡穩(wěn)定性分析中考慮吸力的模擬破壞面比Janbu 簡化方法更符合試驗結(jié)果。在非飽和土中，破壞一般表現(xiàn)為淺層破壞面。在邊坡穩(wěn)定性分析中不考慮基質(zhì)吸力的影響，可能會出現(xiàn)深部破壞面。因此，在非飽和壩體邊坡穩(wěn)定性分析中，需要考慮基質(zhì)吸力作用下的非飽和抗剪強度。泥沙混合物1 的破壞面比泥沙混合物2 更深，這可能是由于泥沙混合物2 中吸力的作用更大。

圖8 冰磧壩的模擬和實驗滑動面對比

為了計算冰磧壩破壞的潰決流量，以下采用van Genuchten 函數(shù)關系式進行滲流計算，并考慮非飽和土中的吸力進行邊坡穩(wěn)定性分析。圖9a 為溢流（案例I）導致冰壩潰決的水槽下游端潰決流量的模擬和試驗結(jié)果。此時，壩體內(nèi)初始含水率約為0.21%。圖9b 為壩內(nèi)初始含水率為5.54%的溢流冰磧壩潰決結(jié)果，壩內(nèi)初始含水率為5.54%時的洪峰流量較高。梯形壩形情況下潰決流量的相似結(jié)果如圖9c 所示。發(fā)現(xiàn)三角形壩形情況下的洪峰流量高于梯形壩形。三種情況下泥沙排出量的模擬結(jié)果與實驗結(jié)果一致。壩面侵蝕還考慮了土體中吸力的影響，溢流對冰磧壩的侵蝕也通過數(shù)值模擬模型得到了很好的解釋。

圖9 溢流潰決導致冰磧壩潰決的潰決流量

圖10 為泥沙混合物1～6 和1～7 在冰磧壩滲流情況下潰決流量的模擬和試驗結(jié)果。冰磧壩潰決后，湖水立即漫過壩面并侵蝕壩面，導致湖水快速下降。泥沙混合物1～6 的峰值流量高于泥沙混合物1～7。在冰磧壩潰壩滲流破壞的情況下，潰決流量的模擬結(jié)果與試驗結(jié)果一致。模擬結(jié)果與實驗結(jié)果的對比表明，本文采用的二維流動與侵蝕模型計算結(jié)果穩(wěn)定可靠，所提出的模型可用于計算冰川潰決導致冰湖潰決的滲流和溢流特征。

圖10 冰磧壩滲流破壞引起的模擬和實驗流量

5 結(jié)論

通過數(shù)值分析計算由滲流和溢流引起的冰磧壩潰決的特征。使用van Genuchten 函數(shù)計算的水分廓線與實驗水分廓線更加吻合。邊坡穩(wěn)定性分析中考慮吸力的模擬破壞面比不考慮吸力的Janbu 簡化方法給出的結(jié)果更精確。潰決流量和壩面侵蝕的結(jié)果也與試驗結(jié)果一致。湖泊潰決流量的突然釋放可能會造成下游洪澇等災害的發(fā)生。