潘佳佳，Hung Tao Shen，郭新蕾

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

1 研究背景

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

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

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

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

2 水冰沙耦合數(shù)學(xué)模型

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

圖1 二維水冰沙耦合數(shù)學(xué)模型框架

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

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

3 模型驗(yàn)證

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

圖2 潰壩試驗(yàn)布置（改自Zech等［19］）

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

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

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

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

圖4 不同時刻和位置下模擬的岸坡侵蝕過程與實(shí)驗(yàn)資料的比較

4 算例應(yīng)用

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

表1 河流水冰沙耦合模型關(guān)鍵參數(shù)

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

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