王亞妮，劉德高，董千里，申國(guó)強(qiáng)

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

1 冰磧壩潰決成因

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

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

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

2 數(shù)值模型

2.1 滲流模型

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

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

Kx和Kz——x 和z 方向的水力傳導(dǎo)率，m/s；

C——比濕容量；

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

x——水平空間坐標(biāo)，m；

z——向上為正的垂直空間坐標(biāo)，m；

t——時(shí)間，s。

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

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

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

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

ψ——壓力水頭；

Ks——飽和導(dǎo)水率，m/s；

m——1-1/η。

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

式中：ψb——空氣進(jìn)入值；

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

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

通過(guò)滲流模型計(jì)算的不同案例含水率模擬值見(jiàn)表1。

表1 含水率模擬值匯總

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

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

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

li——每個(gè)切片底部的長(zhǎng)度，m；

α——每個(gè)切片底部的斜率；

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

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

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

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

式中：σn——總法向應(yīng)力，N/m2；

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

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

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

式中：Pi——第i 個(gè)切片的大氣壓強(qiáng)；

uαi——第i 個(gè)切片的孔隙氣壓；

mαi——第i 個(gè)切片底部的重量。

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

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

2.3 流動(dòng)和侵蝕模型

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

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

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

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

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

h——流動(dòng)深度，m；

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

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

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

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

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

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

β——?jiǎng)恿啃拚禂?shù)；

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

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

τbx——x 方向的底部剪切應(yīng)力，N/m2；

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

sb——河床中的飽和度；

q——滲透率，m/s。

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

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

θb——河床斜率；

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

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

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

dm——平均粒徑，mm；

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

ρw——孔隙流體密度。

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

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

3 水槽實(shí)驗(yàn)

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

圖1 水槽實(shí)驗(yàn)示意圖（單位：cm）

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

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

圖3 不同壩體的詳細(xì)信息

通過(guò)從湖泊上游端供給恒定的排水量來(lái)填充湖水。當(dāng)水位上升漫頂而導(dǎo)致冰磧壩潰決時(shí)，不斷從湖的上游端向壩頂供水并侵蝕。然而，在冰磧壩滲漏潰決的情況下，通過(guò)從湖泊上游端供給恒定的排水量，湖水被填滿的深度約為16 cm。實(shí)驗(yàn)條件見(jiàn)表4。

表4 實(shí)驗(yàn)條件匯總表

4 成果分析

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

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

圖5 湖泊水深隨時(shí)間的變化

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

圖6 案例IV 泥沙混合物1 下模擬和實(shí)驗(yàn)水分含量曲線的比較

圖7 案例V 泥沙混合物2 下模擬和實(shí)驗(yàn)水分含量曲線的比較

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

圖8 冰磧壩的模擬和實(shí)驗(yàn)滑動(dòng)面對(duì)比

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

圖9 溢流潰決導(dǎo)致冰磧壩潰決的潰決流量

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

圖10 冰磧壩滲流破壞引起的模擬和實(shí)驗(yàn)流量

5 結(jié)論

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