潘佳佳，Hung Tao Shen，郭新蕾

（1.流域水循環(huán)模擬與調控國家重點實驗室，中國水利水電科學研究院，北京 100038；2.Department of Civil and Environmental Engineering，Clarkson University，Potsdam，NY 13699-5710）

1 研究背景

北方河流冬季流凌能引起嚴重岸灘刮擦侵蝕和崩塌破壞，伴隨著河冰運動的洪水威脅著河道堤防安全。例如，我國黃河和黑龍江等河流的凌汛時有發(fā)生，給人民生命和財產安全帶來巨大隱患［1-3］。黃河部分河段遇到不利氣溫和大水條件時，流凌或冰壩能造成堤防潰決，進而引起漫堤洪水，形成嚴重凌災［4-7］。此外，流凌刮擦不僅影響岸灘穩(wěn)定性，也能破壞水工建筑物，危及沿岸建筑或房屋安全。研究冬季河冰運動和岸灘刮擦侵蝕的作用機制對防凌減災具有重要意義。

國內外很少有研究關注河冰運動對岸灘的刮擦侵蝕機制，已有的部分研究僅定性分析了河冰運動規(guī)律或岸灘侵蝕現象［8-9］。Lawler 在南威爾士伊爾斯頓河的兩個彎道上開展了長達兩年的連續(xù)監(jiān)測，揭示該河道彎道的岸灘侵蝕后退主要發(fā)生在冬季［10］。該研究表明冬季河道水沙運動過程對岸灘演變具有重要影響，但沒有從機理上解釋引起岸灘冬季侵蝕的原因。Zabilansky等在1998—1999年冬季對美國中部密蘇里河的7個代表性斷面開展現場觀測，發(fā)現這些斷面均存在顯著的流凌刮擦侵蝕［11］。該研究揭示了河冰運動與冬季岸灘侵蝕后退的聯系，解釋了冬季河道岸灘侵蝕頻發(fā)的誘因。他們的研究主要側重于現場觀察分析，并沒有從機理上解釋河冰運動對岸灘侵蝕的促進機制。Beltaos 等在加拿大阿薩巴斯卡河下游的冰塞觀察中也發(fā)現冰塞體的釋放能引起巨大的泥沙輸運，并伴隨著嚴重的河岸刮擦侵蝕［12］。這些研究均強調河冰運動對岸灘侵蝕的影響，但缺少對河冰與岸灘作用力的分析，也沒有分析岸灘崩塌破壞機制。

關于河岸崩塌破壞的研究大部分集中在無河冰影響的工況。Simon等將土力學坡體穩(wěn)定性分析中的靜力平衡理論應用到河岸穩(wěn)定分析，建立堤腳沖刷和岸坡崩塌調整模型［13-14］。該模型能成功模擬河道水流沖刷引起的堤腳淘刷和岸坡失穩(wěn)變化過程，為堤岸穩(wěn)定性分析提供了有效的數值模擬技術。Duan 和Julien 采用岸坡平行后退假設，建立耦合岸灘演變和水沙運動的平面二維數學模型，并成功用于長周期彎道演變機理研究［15］。該研究假設初始河岸坡度處于臨界穩(wěn)定坡度，但這種假設并不符合實際河道復雜的岸坡條件，不能準確模擬河岸演變過程［16］。Abderrezzak 等采用統(tǒng)一泥沙休止角判別岸坡穩(wěn)定與否，在滿足泥沙質量守恒的基礎上調整岸坡坡面，并成功模擬實驗條件下的河岸演變過程［17］，但統(tǒng)一的泥沙休止角并不能反映不同土體性質和含水量的岸坡穩(wěn)定性［18］。Zech等采用兩個不同的泥沙休止角分別代表水上和水下岸坡的臨界失穩(wěn)坡度，并采用較小的休止角重新分布崩岸土體［19］。該方法適合自然環(huán)境下非均質岸坡的穩(wěn)定分析，成功模擬了實驗條件下潰壩引起的岸坡崩塌侵蝕過程。這些研究為岸灘侵蝕和崩塌破壞提供了有效的技術手段，但沒有涉及河冰運動對岸坡穩(wěn)定性的影響分析。

綜上所述，目前的研究缺少對河冰運動與岸灘侵蝕相互作用的分析。本文提出一種北方河流二維水冰沙耦合數學模型，旨在研究冬季北方河流復雜的水冰沙相互作用機制。研究采用耦合水動力過程、河冰運動過程、泥沙運動、河床演變和岸灘侵蝕的方法，模擬了實驗條件下潰壩引起的岸坡崩塌侵蝕過程、冰蓋下的岸灘侵蝕及冰塞冰壩形成與釋放引起的岸灘崩塌和河床變形，有利于揭示河冰運動與岸灘侵蝕間的耦合作用機制。

2 水冰沙耦合數學模型

以Hung Tao Shen河冰研究團隊開發(fā)的二維河冰水沙動力學模型為基礎［20-21］，潘佳佳等建立了一種北方河流二維水冰沙耦合數學模型［22-23］。該模型耦合了二維水沙數值模塊、河冰動力學模塊和岸灘侵蝕模塊，能模擬冰塞冰壩等極端工況下的水位流量波動、河冰聚集和釋放、流凌對岸灘和河床的刮擦侵蝕、泥沙輸移和岸灘的崩塌侵蝕過程［24］。模型的框架示意見圖1。在給定地形、水位、流量和氣溫等初始條件和邊界條件下，水冰沙耦合模型先采用具有迎風特性的Petrov-Galerkin（SUPG）型有限元法計算三角形非結構網格上的水位、流速、流量、水溫、泥沙運動和河床變形，然后將計算結果傳遞給河冰動力學模塊［25-30］；利用無網格的光滑粒子法（SPH）計算河冰分布、冰速、冰厚、面密度和冰濃度等信息，將計算結果傳遞給岸灘侵蝕模塊［23］；采用雙泥沙休止角法判斷不同含水層岸坡的穩(wěn)定與否，在質量守恒的基礎上調整河岸崩塌后的坡面及失穩(wěn)土體堤腳的再分布；最后將調整后的岸灘邊界反饋給河冰動力學模塊和二維水沙數值模塊，校正地形變化后的水位流量、河冰運動、泥沙輸運和河道變形?；谝陨?個模塊間的信息傳遞和數據反饋實現水冰沙的耦合模擬。

圖1 二維水冰沙耦合數學模型框架

二維水冰沙耦合數學模型旨在解決河流全季節(jié)的水沙運動、河床及河岸演變問題。其中，二維水沙數值模塊負責計算河冰影響下的水位流量波動、泥沙輸運及河床沖淤變化，河冰動力學模塊主要模擬水流驅動下的河冰運動、堆積、沖蝕及釋放過程，而岸灘侵蝕模塊提供流凌刮擦下的河道邊界變化。該模型適用于以下研究：（1）復雜河道邊界和地形突變下的急緩流交替過程；（2）降雨降雪、太陽輻射和不同風場氣溫影響下的水溫變化過程；（3）冬季低氣溫影響下的水內冰、錨冰、岸冰和浮冰的生長、消融及變化規(guī)律；（4）冰蓋下流冰花的侵蝕、輸移和沉降過程；（5）開河和封河過程中的冰塞冰壩產生、發(fā)展和釋放機制；（6）河流冰蓋的熱力增長消融和武開河過程；（7）河冰影響下的泥沙運動和河床侵蝕規(guī)律；（8）流凌刮擦和冰塞冰壩影響下的岸灘崩塌侵蝕及再平衡過程。具體的理論背景和基本方程參考北方河流水冰沙耦合研究［24］。

水冰沙耦合數學模型考慮了河冰的形成、發(fā)展、輸移、堆積、釋放及熱力生長消融過程，能計算冬季河流水和冰的相互作用，也能研究夏季或無冰河流水、沙、河床與河岸間的相互作用機理，還能用于北方河流全季節(jié)的水冰沙耦合過程模擬，將水沙理論和河冰理論相融合［22］。

3 模型驗證

二維水冰沙耦合數學模型［20-21］的河冰動力學模塊已被多個實驗算例驗證，并成功用于黃河、加拿大皮斯河和美國尼亞拉加河等多個河流的河冰過程模擬［20，25，28］。本文進一步驗證了二維水冰沙耦合數學模型的岸灘侵蝕模塊，重點采用Zech等試驗條件下潰壩引起的河岸侵蝕后退資料檢驗該模塊的準確性與可靠性［23］。試驗設置在一均勻順直水槽，斷面形態(tài)為半梯形斷面，具體地形與試驗布置見圖2。河床為平坡，水槽長12 m，寬0.495 m。梯形斷面底部寬0.15 m，初始邊坡為50°，坡高為0.16 m。初始時刻閘門設置在距離上游邊界5 m處，位置見圖3。閘門上游河道和河岸初始時刻是不可沖刷的，閘門開啟后可沖刷，水位高0.15 m；其下游是干河床，河道與河岸由1.8 mm的均勻沙組成，相對容重為2.615。該均勻沙在干燥狀態(tài)下泥沙休止角為82°，大于初始邊坡角度；其水下泥沙修止角為30°，小于初始邊坡角度。因此岸坡在無水條件下是穩(wěn)定的，在水位淹沒以下部分是不穩(wěn)定的。該試驗展示了潰壩水流淹沒引起的岸灘崩塌侵蝕過程，被用于河岸侵蝕模型的驗證，具有較好的代表性。

圖2 潰壩試驗布置（改自Zech等［19］）

在試驗中，閘門設置在x=0m處，閘門下游0.5 m 和1.5 m 分別設有斷面地形高程記錄儀，用于提供不同時刻斷面侵蝕后退的試驗數據。二維河冰動力學模型采用非結構的三角形網格，以適應自然河流復雜的地形邊界變化。沿流線方向網格精度為0.05 m，橫斷面方向網格精度為0.016～0.1 m。網格精度在關心的坡面區(qū)域為0.016 m，在主河道區(qū)域相對較粗。計算的時間步長為0.001 s。

圖3展現不同時刻模擬的潰壩水位過程。當閘門在t=0s突然開啟后，水跌引起的潰壩水波快速向下游傳播，在t=1s時經過第一個測站0.5 m處；在t=3s時，水流頭部到達第二個測站1.5 m處。5 s后水流傳播到閘后3 m，但并沒有遇到下游邊界。因此下游邊界的影響可以忽略不計。

圖3 不同時刻模擬的潰壩水位過程

圖4對比了本文模型、Zech 等［19］和Swartenbroekx 等［31］計算的斷面侵蝕崩塌過程與實測資料。在t=3s時，Zech等模擬的坡面侵蝕過大，而堤腳附近淤積偏少；Swartenbroekx等模擬的坡面侵蝕也偏大，而堤腳淤積偏小。本文所采用模型計算的0.5 m和1.5 m兩個橫斷面地形高程與實測資料吻合良好，模型計算精度優(yōu)于其他兩個模型。在t=5s時，Zech 等模型計算的兩斷面坡面崩塌面積大于實測值；Swartenbroekx 等模擬的坡面侵蝕面積較試驗值偏小，堤腳淤積量也偏小。水冰沙耦合數學模型模擬的斷面形態(tài)與試驗資料最為吻合。對比分析可知，本文提出的模型能準確模擬潰壩水流引起的岸坡侵蝕和崩塌過程，而Zech等和Swartenbroekx等模擬的結果存在較大誤差。這主要是因為Zech等采用一維水沙數值模型，而斷面形態(tài)僅通過5個點來捕捉坡面侵蝕和堤腳淤積，橫斷面的網格精度較低。Swartenbroekx 等采用平面二維數學模型，但其輸沙模型采用平衡輸沙模式，不能準確模擬潰壩水流引起的急劇坡面侵蝕和崩塌破壞過程。本文采用非平衡的水沙運動模塊、河冰動力學模塊及岸灘崩塌侵蝕模塊相耦合的數學模型，能更準確地揭示岸灘侵蝕、失穩(wěn)坡面變化和坍塌土體再分布規(guī)律。

圖4 不同時刻和位置下模擬的岸坡侵蝕過程與實驗資料的比較

4 算例應用

4.1 恒定均勻流模擬在潰壩水槽實驗尺度下驗證二維水冰沙耦合數學模型后，該模型被進一步應用于實際尺度冰蓋冰塞影響下的河床沖淤和岸灘侵蝕研究。概化的數值計算區(qū)域為順直、均勻梯形河道，代表性橫斷面形態(tài)及正常水位見圖5。河底寬12 m，計算區(qū)域寬34 m，河岸坡度為1∶4.33，坡頂高度為2.2 m。河道長為400 m，河床底坡為0.0005。河岸與河道由粒徑為0.5 mm 的均勻沙組成。計算區(qū)域的非結構三角形網格，共計節(jié)點5025個，三角形單元9600個。沿流線方向網格精度為2 m，橫斷面方向網格精度在0.9～2 m 間，其中坡面處網格精度較高。網格沿河道中心線y=17m對稱，以保證非結構網格計算的準確性。數學模型計算的時間步長為0.05 s，模型所采用的主要參數見表1。上游邊界設為恒定的流量30m3/s，下游邊界條件設為正常水深1.6 m。在忽略泥沙運動和河冰運動時，模擬的恒定均勻流水深和流速分布如圖6所示。模擬的水深和流速沿軸線對稱，且分布均勻，驗證了該模型在均勻對稱河道水流計算中的可靠性。斷面平均流速為0.98m/s，弗勞德數為0.29。該恒定均勻流模擬為以下河冰運動和岸灘侵蝕模擬提供了可靠的初始條件。

表1 河流水冰沙耦合模型關鍵參數

圖5 均勻梯形斷面形態(tài)及水位示意

圖6 恒定均勻流下模擬的水深和流速分布

4.2 冰蓋對河道沖淤的影響在圖6計算的恒定均勻流下，分別設置有無冰蓋的計算工況，上游泥沙邊界采用循環(huán)邊界，即下游出口流出多少泥沙，上游邊界就提供多少泥沙。循環(huán)泥沙邊界適用于長河道的水沙模擬。圖7為明渠和冰蓋覆蓋兩種工況下100 h 和200 h 模擬的水位和斷面地形。模擬200 h后梯形河道達到沖淤平衡狀態(tài)，冰蓋條件下穩(wěn)定的水位高于明流工況。這主要是因為冰蓋增加了河道邊界的濕周和總體阻力，河道斷面平均流速下降，因此考慮冰蓋時的河道水位更高。4個不同斷面均顯示冰蓋下水位更高，水面下岸坡的沖刷起點更高，但冰蓋下岸坡總體沖刷量小于明流工況，相應主河道的泥沙淤積也偏小。這主要是因為冰蓋下的流速比明流工況低，岸坡和床面的水流拖曳力低于明流工況，因此岸坡的沖刷和主河道的淤積小于明流條件。對比100 h和200 h的模擬結果，下游河道比上游更快達到沖淤平衡。計算結果顯示二維水冰沙耦合數學模型能模擬河道中冰蓋對岸坡沖刷和河道淤積的影響。

圖7 明流和冰蓋影響下不同斷面模擬的水位及河床沖淤對比

4.3 冰塞形成和釋放過程為了進一步驗證二維水冰沙耦合數學模型在冰塞冰壩等極端工況下的可靠性，研究采用第3節(jié)工況下的梯形斷面河道及河岸泥沙條件，模擬了沖積河流冰塞冰壩的形成和釋放過程。在河道中間設置攔冰柵并假設冰塞形成過程中河床為可沖刷而河岸固定，在冰塞釋放過程中河床與河岸均可沖刷。圖8顯示了穩(wěn)定冰塞形成后浮冰顆粒的分布、河道地形的沖淤變化、縱剖面的水位、上下冰面及河床高程分布。圖中，黑色方塊示意浮冰，地形沖淤變化中藍色示意沖刷，紅色示意淤積，兩條白線示意左右岸堤腳位置。模擬結果顯示卡冰處冰體堆積嚴重，最大冰厚超過0.12 m，冰塞引起上游水位壅高，上下水位差達0.04 m。此外，冰塞體前端冰厚最大，并導致斷面過流面積減小，局部流速增大，河道主槽出現顯著沖刷。模擬結果與Beltaos等在加拿大阿薩巴斯卡河觀測的沖積河流冰塞形成過程相符合［32］，證實了模型對冰塞形成中水冰沙耦合模擬的可行性。

以圖8模擬的冰塞體為初始條件，在t=0 時刻撤去攔冰柵，進一步采用二維水冰沙耦合數學模型模擬冰塞釋放過程中的河冰運動、河道沖淤變化、岸坡崩塌和再平衡、水位流量波動過程。模擬的不同時刻水位、地形和河冰分布結果見圖9。計算2 s后，儲存在卡冰點上游的浮冰和水體以類似潰壩水體的洪水波向下游傳播。下游岸坡在淹沒和水流沖刷下發(fā)生河岸侵蝕、崩塌和堤腳淤積，而上游河道流凌及水流也導致岸坡上部沖刷坍塌及堤腳淤積。隨著冰塞的釋放，冰塞體前端水位及冰厚減小，并向下游平滑過渡。在3至6 s，上游儲存的浮冰和水體繼續(xù)向下游釋放，流凌刮擦引起的岸灘崩塌更為顯著，坍塌的土體堆積于堤腳附近，部分泥沙被輸移到下游。計算的12 s，冰塞釋放的浮冰流出計算區(qū)域，上游壅水部分回落，水面附近的岸坡出現顯著侵蝕和崩塌，而兩側堤腳出現大量淤積。其中上游的岸坡沖刷和堤腳淤積大于下游河道。這是合理的，因為上游河道壅水高于下游，冰塞體刮擦影響區(qū)域也較下游河道大。

圖8 冰塞形成引起的水位壅高及河床沖刷模擬結果

圖9 不同時刻模擬的冰塞釋放過程及水位、地形和河冰分布

圖10顯示了上游至下游4個典型斷面0、3、6、9和12 s計算的斷面地形沖淤變化過程。4個斷面均顯示水面下存在岸坡沖刷侵蝕，水面以上的岸坡也出現崩塌侵蝕，大量失穩(wěn)土體落淤在堤腳附近。從上游到下游，沖刷和崩岸的土體面積逐漸減小，而堤腳淤積的土體面積相應也減少。這主要是因為水位從上游到下游逐漸減小，河冰大量堆積在上游河道，流凌刮擦和水流侵蝕范圍從上游向下游遞減。此外，水面以上岸坡再平衡的坡角大于水面以下岸坡的平衡坡角。計算結果顯示水冰沙耦合模型能準確模擬水面上下不同含水層的穩(wěn)定坡面。圖8—10 證實該模型能有效模擬極端冰塞下的岸灘刮擦侵蝕、岸坡崩塌、失穩(wěn)土體堤腳落淤和冰塞釋放過程，能進一步應用于其它水冰沙問題研究。

圖10 冰塞釋放引起的岸灘侵蝕和斷面地形變化

5 結論

本文初步建立二維水冰沙耦合數學模型，由水沙數值模塊、河冰動力學模塊和岸灘侵蝕模塊3部分組成。模型采用非結構的有限元法計算水沙運動和河床沖淤變化，能適應自然河流復雜的河道邊界；采用無網格的光滑粒子法計算河冰運動、堆積和釋放過程，能模擬極端冰塞冰壩的形成、發(fā)展和釋放過程；采用雙泥沙休止角法計算岸灘的侵蝕、崩塌和再平衡過程，能揭示岸坡的周期性侵蝕規(guī)律。研究采用潰壩水流引起的岸灘侵蝕和崩塌實驗驗證了模型的準確性和可靠性，并將模型應用于冰蓋和冰塞形成與釋放引起的河床沖淤與岸灘崩塌侵蝕分析。模擬結果顯示二維水冰沙耦合數學模型能揭示流凌對岸灘的刮擦侵蝕、河岸崩塌及失穩(wěn)土體再分布規(guī)律，為極端冰塞冰壩下的水冰沙相互作用研究提供了新工具。二維水冰沙耦合模型創(chuàng)新性地融合河冰理論和水沙理論，同時滿足無冰和有冰河流水沙和河冰模擬需求，為河流全季節(jié)的水冰沙耦合問題研究提供了有力支撐。

北方河流冬季水流運動、河冰輸運、泥沙輸移、河床變化與岸灘侵蝕存在復雜的水冰沙相互耦合作用，本文在數值模擬和理論計算方面開展了詳細研究，還需通過原型觀測和物理模型實驗等多種手段深入研究水冰沙耦合作用機理。關于北方河流河冰與河岸的相互作用及泥沙運動的實測資料和理論分析欠缺，亟需開展更多系統(tǒng)水冰沙耦合作用研究。