潘佳佳，Hung Tao Shen，郭新蕾，王 濤

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

1 研究背景

我國(guó)北方河流如黃河、黑龍江和松花江等每年都有超過(guò)100天的冰期［1］，而冬季河冰運(yùn)動(dòng)對(duì)泥沙輸運(yùn)和河道演變的影響常常被忽略。一方面受地區(qū)和時(shí)間的限制，河冰影響下的水沙問(wèn)題是季節(jié)性過(guò)程，不及明渠水沙研究更具代表性，常常被研究人員忽略［2-4］。另一方面冰期河流涉及冰體堆積釋放、水位壅高、流量波動(dòng)、河床沖淤變化和岸灘崩塌侵蝕等多種過(guò)程，存在復(fù)雜的水冰沙相互作用。北方河流水冰沙耦合作用機(jī)理問(wèn)題是水力學(xué)、河冰動(dòng)力學(xué)和河流動(dòng)力學(xué)的交叉方向，涉及的物理因素多［4-6］，問(wèn)題復(fù)雜，是河冰領(lǐng)域研究的前沿和難點(diǎn)。

冬季河冰過(guò)程對(duì)北方河流水沙運(yùn)動(dòng)的影響至關(guān)重要。全球氣候變化和人類(lèi)活動(dòng)影響下，極端冰塞冰壩發(fā)生的可能性更大。冰塞冰壩能引起上游河道水位迅速抬高，流凌刮擦割蝕岸灘能導(dǎo)致堤岸崩塌破壞，由此引發(fā)的凌汛洪災(zāi)嚴(yán)重威脅北方河流冬季輸水安全和河流管理。河冰不僅影響泥沙運(yùn)動(dòng)和河道變化，還顯著影響水體溫度和含氧量，例如錨冰和冰蓋的形成會(huì)壓縮水生物生長(zhǎng)繁殖空間，進(jìn)而影響水生態(tài)環(huán)境［7］。這些河冰過(guò)程吸引了眾多學(xué)者的關(guān)注，并在水內(nèi)冰、岸冰、錨冰、冰蓋、冰塞和冰壩等方面取得長(zhǎng)足進(jìn)步［8-11］，但缺少耦合水沙運(yùn)動(dòng)和河冰動(dòng)力過(guò)程的研究［3-4］。目前岸灘的崩塌侵蝕研究主要基于無(wú)冰條件［12-14］，不能滿足北方河流季節(jié)性岸灘侵蝕的研究需求。

傳統(tǒng)的河流動(dòng)力學(xué)主要研究水流、泥沙與河道邊界（包括河床和河岸）三者之間的復(fù)雜相互作用，而目前河冰動(dòng)力學(xué)研究主要分析河冰運(yùn)動(dòng)與水流間的相互作用。本文首先簡(jiǎn)要總結(jié)了河冰水力學(xué)的研究進(jìn)展，詳細(xì)綜述了有關(guān)河冰過(guò)程的理論成果和已有研究的不足。這些研究是水冰沙耦合分析的基礎(chǔ)，為系統(tǒng)的河冰模擬奠定了堅(jiān)實(shí)的理論框架。在綜合各河冰數(shù)學(xué)模型優(yōu)點(diǎn)的基礎(chǔ)上，本文提出河流水冰沙耦合數(shù)學(xué)模型［4，15］，既能模擬無(wú)冰河流水沙運(yùn)動(dòng)和岸灘演變過(guò)程，也能計(jì)算冬季水冰變化，還能模擬北方河流全季節(jié)及河冰生消全過(guò)程的水冰沙耦合過(guò)程。該研究有助于分析流凌和冰塞冰壩引起的洪水過(guò)程及堤岸潰決問(wèn)題，為北方河流防凌減災(zāi)提供技術(shù)支撐。

2 河冰水力學(xué)理論框架

Beltaos 從河冰熱力生長(zhǎng)和消融過(guò)程、水內(nèi)冰、錨冰及岸冰的發(fā)展、冰塞冰壩形成和釋放引起的水位流量波動(dòng)等方面詳細(xì)綜述了相關(guān)理論、數(shù)學(xué)模型以及原型觀測(cè)的顯著進(jìn)步［16］，但也強(qiáng)調(diào)關(guān)于氣候變化和人類(lèi)活動(dòng)影響下的冰塞冰壩形成機(jī)理和預(yù)測(cè)預(yù)報(bào)仍有很多不足。Shen［17］詳細(xì)綜述了河冰全過(guò)程所涉及的熱力、動(dòng)力和水力過(guò)程，指出河冰研究是水力學(xué)、冰體力學(xué)、熱力學(xué)和河流動(dòng)力學(xué)等多學(xué)科的交叉領(lǐng)域，所包含的物理過(guò)程復(fù)雜，關(guān)于河冰的理論和數(shù)學(xué)模型在過(guò)去30年有了長(zhǎng)足進(jìn)步，能協(xié)助解決天然河流和實(shí)際工程中涉冰的防洪、發(fā)電、航道、生態(tài)及環(huán)境問(wèn)題。河冰數(shù)學(xué)模型的建立與發(fā)展為北方河流冬季用水安全和管理提供了有力的技術(shù)工具。

受冬季低氣溫的影響，北方河流水體會(huì)釋放熱量，進(jìn)而生成過(guò)冷卻水和水內(nèi)冰［18］。新生成的水內(nèi)冰具有較強(qiáng)的吸附性，能黏附于河床、水工結(jié)構(gòu)物及河岸，形成錨冰和岸冰；而大量聚集的水內(nèi)冰上浮到水面后可以形成表面浮冰。隨著氣溫的進(jìn)一步降低，河流低流速區(qū)域的表面冰蓋會(huì)出現(xiàn)聚集，直至整個(gè)斷面被河冰覆蓋，形成冬季河流的首封位置［19］。在封河過(guò)程中，由于上游徑流的變化和氣溫的波動(dòng)，在彎道、收縮斷面以及阻水建筑物附近能形成冰體堆積。堆積體前大量上游來(lái)冰的聚集和堵塞會(huì)縮小斷面過(guò)流面積、顯著增加河道阻力，進(jìn)而引起冰體卡塞點(diǎn)上游水位的壅高，形成封河冰塞［20］。在冬末春初氣溫回升時(shí)，河道冰蓋發(fā)生熱力消融和動(dòng)力破壞，進(jìn)而出現(xiàn)開(kāi)河。在開(kāi)河過(guò)程中，如果水流平穩(wěn)，河面上的冰蓋最終會(huì)完全消融，這種平穩(wěn)的熱力開(kāi)河為“文開(kāi)河”［21］。如果上游融雪融冰產(chǎn)流大或降雨較多，在過(guò)水?dāng)嗝媸芟薜暮拥罆?huì)出現(xiàn)大量上游來(lái)冰堆積，持續(xù)的冰體聚集和上游水位壅高會(huì)產(chǎn)生開(kāi)河冰壩。嚴(yán)重的冰壩能引起洪水漫堤或堤岸潰決，進(jìn)而誘發(fā)凌汛災(zāi)害。冰壩引起的河面冰蓋破裂和釋放常常伴隨著急劇水位流量波動(dòng)，這種水冰動(dòng)力誘發(fā)的開(kāi)河為“武開(kāi)河”［11］。河冰的生長(zhǎng)、運(yùn)動(dòng)和釋放會(huì)影響水工建筑物結(jié)構(gòu)穩(wěn)定性，也能刮擦侵蝕河床岸灘，進(jìn)而影響沖積河流的輸沙過(guò)程和河道演變。詳細(xì)的河冰過(guò)程和河冰水力學(xué)理論框架見(jiàn)圖1，主要包括水體失熱、河流產(chǎn)冰、封河、開(kāi)河及河冰影響等五個(gè)方面。

圖1 河冰水力學(xué)理論框架

2.1 水體失熱北方河流水體的熱交換包括徑流和支流的能量匯入、下游出流的能量輸出、空氣與水體的熱交換、水體動(dòng)能和勢(shì)能的轉(zhuǎn)換及河床與水體的熱交換。冬季水體失熱主要是河流與空氣間的熱交換，占水體熱循環(huán)的90%以上［22］。持續(xù)的水體失熱是河冰形成的前提。河床一般在溫暖的季節(jié)儲(chǔ)存能量，在冬季向水體釋放熱量［23］。河床與水體間的熱交換在水體熱循環(huán)中占比較小，但顯著影響河床上錨冰的生長(zhǎng)和釋放過(guò)程［24-25］。水體運(yùn)動(dòng)中動(dòng)能和勢(shì)能轉(zhuǎn)換的熱能一般較小，可以忽略［15］。因此，河流水體熱循環(huán)主要考慮水面或冰面與空氣間的熱交換。

針對(duì)河流表面熱交換，Shen 考慮太陽(yáng)輻射、水體長(zhǎng)波輻射、蒸發(fā)熱損失、降雨降雪熱損失和水體與空氣間的熱傳導(dǎo)等因素，建立了冬季北方河流表面熱交換的準(zhǔn)確計(jì)算方法［26］。除了太陽(yáng)輻射是水體吸收熱量外，其它熱傳導(dǎo)均為河流向外界釋放熱量，進(jìn)而促進(jìn)水體產(chǎn)冰。該方法能提供詳細(xì)的水體熱交換過(guò)程及各影響因子的權(quán)重，但需要詳細(xì)的太陽(yáng)輻射角、降雪量、河道經(jīng)緯度、空氣透射程度、云層狀況、氣溫、水溫和風(fēng)速風(fēng)向等各項(xiàng)資料，不便于實(shí)際工程應(yīng)用。為了簡(jiǎn)化計(jì)算，Shen等提出了近似的線性公式計(jì)算河流表面的熱交換［27-28］

式中：φ為河流表面損失的凈熱通量，W/m2；φs為太陽(yáng)輻射的熱通量，W/m2；β為與氣候有關(guān)的經(jīng)驗(yàn)參數(shù)，W/m2；α為水面或冰面與空氣間的熱傳導(dǎo)系數(shù)，W/(m2℃)；Ts為河流表面水面或冰面的溫度，℃；Ta為空氣的溫度，℃。針對(duì)長(zhǎng)期大尺度的河冰數(shù)值模擬，式（1）可簡(jiǎn)化為［27-28］

式中α′為考慮太陽(yáng)輻射和其它熱損失的綜合熱傳導(dǎo)系數(shù)，W/(m2℃)。

2.2 產(chǎn)冰當(dāng)河流水面溫度降到0℃后，水體開(kāi)始產(chǎn)冰。隨著水面的持續(xù)失熱，水體溫度能出現(xiàn)一定的過(guò)冷卻到達(dá)-0.01℃的量級(jí)，大量的水內(nèi)冰產(chǎn)生并釋放潛熱以彌補(bǔ)水體熱量損失，接著水體溫度回升并維持在結(jié)冰溫度［18，29-30］。湍流下水體過(guò)冷卻和水內(nèi)冰產(chǎn)生計(jì)算的理論公式相對(duì)完善，具體見(jiàn)文獻(xiàn)［18，30-31］。然而，水內(nèi)冰上浮、懸浮和下沉的運(yùn)動(dòng)機(jī)理尚不明確，有待研究［32］。

在河流的低流速區(qū)如河岸附近易形成靜態(tài)的岸冰，進(jìn)而促進(jìn)斷面冰蓋的產(chǎn)生。Matousek［33］總結(jié)了不同類(lèi)型流冰的原型觀測(cè)資料，給出了靜態(tài)冰蓋和運(yùn)動(dòng)薄冰形成的水溫和流速條件，冰蓋的生長(zhǎng)過(guò)程與湍流強(qiáng)度有關(guān)。水流速度過(guò)大或者水溫不夠低時(shí)冰蓋均不能發(fā)展［34］。岸冰的橫向發(fā)展速率與表面流冰密度成正比，與水流拖曳力和流冰與岸冰摩擦力的大小密切相關(guān)［17，35］。

錨冰是一種黏附在河床上的冰，相對(duì)浮冰的觀測(cè)難度更大。錨冰的生長(zhǎng)一般在水體達(dá)到結(jié)冰溫度之后，在河面完全冰封之前。當(dāng)錨冰受到的浮力大于冰體與河床間的吸附力時(shí)會(huì)從河床釋放，進(jìn)而上浮到水面。錨冰底部因吸收太陽(yáng)輻射而融化時(shí)也會(huì)上浮釋放。冷凍水槽試驗(yàn)顯示錨冰的生長(zhǎng)和釋放能顯著影響河床的綜合糙率［36］。錨冰的釋放能搬運(yùn)所吸附的河床泥沙，甚至輸運(yùn)幾十公斤的大型卵石或石塊，并在融化后釋放泥沙［37］。最近Pan等［25］的研究顯示錨冰的形成與釋放能引起河床高程的變化、斷面流量的波動(dòng)和河道整體糙率的急劇變化，進(jìn)而造成水深和流量高達(dá)30%的增長(zhǎng)波動(dòng)。

2.3 封河當(dāng)河道出現(xiàn)卡塞點(diǎn)和初始冰蓋后，上游的來(lái)冰會(huì)在卡塞處堆積。如果水流速度較低，斷面弗勞德數(shù)低于一定臨界值時(shí)，水面浮冰會(huì)以平鋪上溯的方式發(fā)展（Juxtaposition mode）；如果水流速度較大，浮冰所受的水流拖曳力能拖動(dòng)冰塊下潛或翻轉(zhuǎn)到冰蓋下，冰蓋會(huì)以水力加厚的模式（Hydrau?lic thickening mode）向上游發(fā)展，又被稱(chēng)為窄河冰塞［16］；當(dāng)下潛的浮冰過(guò)多，冰塞體的內(nèi)摩擦力小于水流拖曳力和上游壓力時(shí)，冰塞體會(huì)出現(xiàn)坍塌和力學(xué)加厚，進(jìn)而以更厚的冰塞向上游發(fā)展（Mechani?cal thickening mode），這被稱(chēng)為寬河冰塞，大部分河流的冰塞屬于該類(lèi)型［38］；如果水流速度進(jìn)一步增大，斷面弗勞德數(shù)超過(guò)上限臨界值時(shí)，上游來(lái)冰在水流拖曳作用下下潛并沿冰塞體底部滑移，最終沖蝕到下游河道，且不能停留堆積在冰塞體下，因此冰蓋不會(huì)向上游發(fā)展［17］。這會(huì)在上游河段出現(xiàn)明流區(qū)域，并不斷產(chǎn)生水內(nèi)冰。水內(nèi)冰在下游冰蓋下的輸移和堆積可進(jìn)一步形成懸冰壩［29］。

經(jīng)典的河冰水力學(xué)理論采用靜力平衡原理分析了冰蓋的形成和臨界冰厚，將水力學(xué)理論應(yīng)用到冰塞問(wèn)題分析［39］。Michel 將能量守恒方程應(yīng)用于冰塞頭部計(jì)算，建立窄河冰塞下臨界流速與冰厚的關(guān)系，進(jìn)一步量化了冰塞的發(fā)展過(guò)程［40］。Uzuner and Kennedy 采用連續(xù)方程和動(dòng)量方程分析了寬河冰塞下的水流特征，提出非恒定流下寬河冰塞的冰厚計(jì)算方法［41］。通過(guò)自然河流中大量冰塞的原型觀測(cè)，Beltaos 與Fan 等［42-43］收集了平衡和非平衡冰塞資料，提出冰塞對(duì)水體的阻力與冰厚成正比，并指出冰塞釋放引起的岸灘刮擦侵蝕和河床沖刷遠(yuǎn)大于夏季汛期的河道沖刷［44-45］。沈洪道團(tuán)隊(duì)針對(duì)冰塞過(guò)程中的非恒定水流運(yùn)動(dòng)，采用冰水雙層流連續(xù)介質(zhì)假定，建立一二維河冰動(dòng)力學(xué)模型，能模擬河冰的產(chǎn)生、輸運(yùn)、堆積、消融及冰塞的發(fā)展釋放過(guò)程，為河冰研究提供了有效的技術(shù)手段［17，46-50］。

2.4 開(kāi)河相比于封河研究，冬末春初的開(kāi)河研究較少。關(guān)于冰蓋橫縫和縱縫的開(kāi)裂過(guò)程及碎冰蓋的卡塞位置和時(shí)間仍是河冰研究的技術(shù)瓶頸。隨著冬末氣溫的升高和太陽(yáng)輻射的增強(qiáng)，冰蓋消融，冰體強(qiáng)度下降。在上游洪水波的作用下，陡峭河段的冰蓋破裂成塊，并向下游輸運(yùn)。類(lèi)似封河時(shí)的冰塞形成，開(kāi)河堆積的冰塊也能形成冰塞或冰壩。國(guó)內(nèi)外關(guān)于冰塞和冰壩的定義存在較大差異。國(guó)外將開(kāi)河和封河過(guò)程中浮冰或流動(dòng)冰塊堆積形成的冰體雍塞均定義為冰塞（Ice jam），將穩(wěn)定冰蓋下流冰花或流冰塊的懸浮堆積稱(chēng)為懸冰壩（Hanging ice dam）［16-17，29］。在國(guó)內(nèi)，一般將封河過(guò)程中產(chǎn)生冰體堆積和擁堵定義為封河冰塞，將開(kāi)河時(shí)期冰塊堆積和擁堵定義為開(kāi)河冰壩。受流量和來(lái)冰量的影響，一般開(kāi)河冰壩引起的水位和流量波動(dòng)比封河冰塞大［21］。

Shulyakovskii［49］觀測(cè)了大量河流開(kāi)河過(guò)程，采用開(kāi)河和封河時(shí)期的水位差作為冰蓋破碎與開(kāi)河的標(biāo)準(zhǔn)，建立冰蓋強(qiáng)度與開(kāi)河水位和封河水位的經(jīng)驗(yàn)公式。Beltaos 通過(guò)大量加拿大河流的原型觀測(cè)，進(jìn)一步提出開(kāi)河與封河期的水位比值與冰厚、冰封寬度、冰體強(qiáng)度和水流拖曳力有關(guān)，建立了開(kāi)河冰蓋破裂和移動(dòng)的概化模型［50］。但這些經(jīng)驗(yàn)公式具有較大的主觀因素，不利于其他河流的應(yīng)用。國(guó)內(nèi)開(kāi)河預(yù)報(bào)主要基于人工智能算法［10］。陳守煜等［51］基于BP神經(jīng)網(wǎng)絡(luò)模型，預(yù)報(bào)了黃河寧蒙河段開(kāi)河日期。王濤等［52］基于GIS 地理信息系統(tǒng)和神經(jīng)網(wǎng)絡(luò)預(yù)報(bào)模型，建立黃河冰情預(yù)報(bào)專(zhuān)家系統(tǒng)，能有效預(yù)報(bào)寧蒙河段開(kāi)河日期，預(yù)報(bào)期精度在10 d左右。王軍等［53］采用BP神經(jīng)網(wǎng)絡(luò)模型模擬了實(shí)驗(yàn)條件下彎道冰塞的水位壅高，能較好地預(yù)測(cè)河冰影響下的水位變化。

受溫升融水和降雨的影響，春汛引起的開(kāi)河冰壩在北方河流最為常見(jiàn)，相應(yīng)的凌汛災(zāi)害也更嚴(yán)重?；诤颖鶆?dòng)力學(xué)理論的數(shù)學(xué)模型被成功用于日本渚滑川和加拿大圣約翰河等河流的冰壩過(guò)程模擬［54-55］，但開(kāi)河相關(guān)的力學(xué)機(jī)制尚待完善。開(kāi)河期碎冰塊間歇性的運(yùn)動(dòng)和停滯易誘發(fā)鏈?zhǔn)奖鶋蔚男纬珊歪尫牛殡S的凌汛洪水風(fēng)險(xiǎn)更高。開(kāi)河冰壩最大的難題是冰壩坍塌釋放機(jī)理和碎冰塊與冰蓋間的相互作用機(jī)制，包括冰塊與冰蓋間的相互作用力及水對(duì)冰體的作用力等。此外，上游下泄的碎冰塊可能撞擊擠壓下游穩(wěn)定冰蓋，進(jìn)而導(dǎo)致更多冰蓋破裂，相關(guān)力學(xué)機(jī)理還需進(jìn)一步研究。

2.5 河冰影響冬季河冰運(yùn)動(dòng)和堆積可能影響水工建筑物運(yùn)營(yíng)及安全［56-58］。封河期生成的大量水內(nèi)冰能吸附在攔冰柵和取水管口，造成取水口堵塞，進(jìn)而影響冬季供水和輸水安全［25］。流凌運(yùn)動(dòng)直接撞擊或刮擦橋墩，嚴(yán)重的能引起橋墩混凝土表面脫落，影響橋梁結(jié)構(gòu)的安全性。此外，橋墩結(jié)構(gòu)所在的斷面容易造成浮冰堆積，促進(jìn)冰塞冰壩的形成，引起凌汛洪水［59］。嚴(yán)重的冰壩能在橋梁或其它阻水建筑物迎冰面形成較大推力，造成水工結(jié)構(gòu)的破壞。

Ettema 和Kempema 詳細(xì)綜述了河冰對(duì)沙質(zhì)河床地形、泥沙輸運(yùn)、河岸侵蝕的影響，強(qiáng)調(diào)封河期和開(kāi)河期河冰運(yùn)動(dòng)對(duì)泥沙運(yùn)動(dòng)和河道演變的影響最為劇烈［60］。錨冰的釋放能搬運(yùn)大量泥沙，甚至輸運(yùn)非冰期無(wú)法啟動(dòng)的卵石［37］。冰塞冰壩的形成和釋放會(huì)導(dǎo)致河冰刮擦割蝕河岸，嚴(yán)重的能導(dǎo)致堤岸崩塌破壞，誘發(fā)凌汛洪水。另一方面，岸冰的凍結(jié)能避免流凌直接刮蝕河岸，進(jìn)而保護(hù)河岸［2，6］。受限于冬季惡劣的天氣條件、儀器裝備的不足和理論的缺失，目前仍缺少關(guān)于河冰對(duì)河床沖淤影響的研究，無(wú)法定量評(píng)估河冰對(duì)岸灘的刮擦侵蝕、冰蓋下的泥沙輸運(yùn)和冰塞冰壩作用下的河道演變規(guī)律。

3 水冰沙耦合模型

3.1 河冰數(shù)學(xué)模型的基礎(chǔ)自從1960年代Pariset和Hausser 將水力學(xué)基本方程引入到河冰冰塞分析，眾多研究者完善發(fā)展了冰塞理論，并促進(jìn)河冰數(shù)學(xué)模型的產(chǎn)生和發(fā)展［16-17，39，41，61］。基于靜力平衡理論，F(xiàn)lato和Gerard開(kāi)發(fā)了一維恒定流下的非平衡冰塞模型（ICEJAM）［62］，Beltaos 建立了考慮冰塞中水體滲流的寬河冰塞模型（RIVJAM）［63］。隨后，靜力平衡冰塞理論被包含在非恒定一維商業(yè)模型中，包括HECRAS和MIKE11-ICE等［64-65］。這些一維模型能預(yù)測(cè)冰塞冰壩等極限條件下的冰厚和水位上限，為冬季凌汛防治提供了依據(jù)，但它們忽略了河冰運(yùn)動(dòng)的動(dòng)力過(guò)程，并不能判斷冰塞和冰壩形成的條件和位置。

Shen 等將一維非恒定圣維南方程與河冰運(yùn)動(dòng)、質(zhì)量守恒和面密度平衡方程相結(jié)合，奠定了河冰動(dòng)力學(xué)模型的基礎(chǔ)，能有效模擬河冰的運(yùn)動(dòng)以及冰塞的動(dòng)態(tài)發(fā)展過(guò)程［46］。隨后，Lal 和Shen 建立了一維非恒定流河冰數(shù)學(xué)模型（RICE），考慮水溫變化、河冰熱力生長(zhǎng)過(guò)程及河冰運(yùn)動(dòng)與水流過(guò)程的相互影響，能模擬水內(nèi)冰、表面浮冰、岸冰和冰塞過(guò)程［66］。Shen 等在RICE 的基礎(chǔ)上進(jìn)一步開(kāi)發(fā)了RICEN 模型，擴(kuò)展了復(fù)雜河網(wǎng)模塊、錨冰模塊和基于輸冰率的水內(nèi)冰輸移模塊，并成功應(yīng)用于黃河下游和美國(guó)尼亞拉加河上游河道的河冰模擬［27］。隨后，該一維河冰模型被升級(jí)為CRISSP1D和RICEE，能完整考慮冬季河冰全生命周期的產(chǎn)生、發(fā)展和演變過(guò)程［48，67-68］。Pan 等基于CRISSP1D 模型，開(kāi)發(fā)了考慮錨冰生長(zhǎng)和釋放的錨冰洪水波模塊，揭示了錨冰引起的河床地形抬高、斷面流量變化和河床綜合糙率變化，其中錨冰影響下的河床糙率變化是水位和流量波動(dòng)的主要因素［25］。

為了考慮復(fù)雜地形和河岸的影響，Shen 等開(kāi)發(fā)了二維河冰動(dòng)力學(xué)模型（DynaRICE），采用基于歐拉場(chǎng)的有限元法計(jì)算二維淺水方程，采用基于拉格朗日式的無(wú)網(wǎng)格光滑粒子法計(jì)算河冰運(yùn)動(dòng)，并應(yīng)用于密西西比河中游河段冰塞過(guò)程模擬，為寒區(qū)河流的防凌減災(zāi)提供了有效技術(shù)支撐［69］。隨后，Liu等開(kāi)發(fā)了河冰全生命周期的二維數(shù)學(xué)模型（Two-dimensional Comprehensive River Ice Simulation System，CRISSP2D），能模擬（1）復(fù)雜地形下的急緩流過(guò)程；（2）太陽(yáng)輻射及風(fēng)場(chǎng)影響下的水溫升降規(guī)律；（3）水內(nèi)冰的生成及錨冰、岸冰和浮冰過(guò)程；（4）冰蓋下的浮冰輸移和沉降過(guò)程；（5）冰塞冰壩的產(chǎn)生、發(fā)展和釋放過(guò)程；（6）冰蓋的熱力增長(zhǎng)和消融及武開(kāi)河過(guò)程［70］。Knack和Shen耦合二維河冰模塊與水沙數(shù)值模塊，將CRISSP2D 發(fā)展為Hydro-Thermal-Ice-Sediment River Dynamics Model，能模擬北方河流河冰影響下的推移質(zhì)和懸移質(zhì)泥沙運(yùn)動(dòng)及河床變形［4］。此后，Pan 和Shen 將CRISSP2D 發(fā)展為RICES2D，考慮河冰影響下的泥沙輸運(yùn)、河床沖淤變化及岸灘侵蝕等，為北方河流水冰沙耦合機(jī)理研究奠定了基礎(chǔ)［71］，也是本系列文章的前提。

3.2 數(shù)學(xué)模型的構(gòu)建基于沈洪道團(tuán)隊(duì)長(zhǎng)期開(kāi)發(fā)的河冰動(dòng)力學(xué)模型和河冰水力學(xué)理論［17］，本文建立了二維水冰沙耦合數(shù)學(xué)模型，以研究冰蓋和流凌作用下的泥沙輸移、河床演變及岸灘崩塌侵蝕問(wèn)題。該模型耦合了水動(dòng)力過(guò)程、熱力學(xué)過(guò)程、河冰運(yùn)動(dòng)、泥沙輸運(yùn)、岸灘侵蝕和河床演變過(guò)程，模型框架及功能見(jiàn)圖2。該模型包括二維水沙數(shù)值模塊、河冰動(dòng)力學(xué)模塊和岸灘侵蝕模塊三部分。其中，二維水沙數(shù)值模塊采用基于非結(jié)構(gòu)網(wǎng)格的有限元法計(jì)算非恒定水流運(yùn)動(dòng)、推移質(zhì)輸沙率和河道沖淤變形，河冰動(dòng)力學(xué)模塊利用無(wú)網(wǎng)格的光滑粒子法計(jì)算河冰的運(yùn)動(dòng)、分布、水力和動(dòng)力加厚過(guò)程，而岸灘侵蝕模塊采用雙泥沙休止角法計(jì)算岸坡崩塌和土體再分布。通過(guò)給定耦合時(shí)間下的信息傳遞和反饋，調(diào)整變化河冰條件下的水流和泥沙過(guò)程，進(jìn)而實(shí)現(xiàn)河流冰期水冰沙耦合過(guò)程的模擬。以下分節(jié)給出三個(gè)模塊的控制方程和計(jì)算方法。

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

3.3 二維水沙數(shù)值模塊二維水沙數(shù)值模塊包括二維淺水運(yùn)動(dòng)、泥沙輸運(yùn)和河床變形三個(gè)部分。基于上層浮冰和下層水體的雙層流體系，考慮河冰阻力和質(zhì)量引起的水位流量變化，北方河流河冰影響下的二維淺水方程組為［69］：

式中：x、y 為平面直角坐標(biāo)，m；t 為時(shí)間變量，s；qtx、qty為兩個(gè)坐標(biāo)方向上的單寬流量，m2/s；t′i為淹沒(méi)在水下的冰厚，m；H 為總水深包括河冰的影響，m；Ht為單寬流量對(duì)應(yīng)的有效水深，m；C為河冰的面密度或河冰面積占整個(gè)河面的比例；Txx、Tyx、Txy、Tyy分別為不同方向和作用面上的切應(yīng)力沿水深的積分值，N/m；ρ為水體密度，kg/m3；τsx、τsy分別為兩個(gè)方向上的表面冰對(duì)水體的切應(yīng)力，Pa；τbx、τby分別為兩個(gè)方向上的床面切應(yīng)力，Pa；η為水面高程，m。

淺水方程組采用基于三角形非結(jié)構(gòu)網(wǎng)格的有限元方法求解，能準(zhǔn)確描述復(fù)雜地形下不規(guī)則河道邊界。該有限元法基于具有迎風(fēng)特性的皮德羅夫-蓋爾金（Petrov-Galerkin）算法，能有效捕捉干濕邊界并適應(yīng)急緩流變化的條件，甚至能模擬潰壩引起的洪水波傳播過(guò)程，在多個(gè)國(guó)家眾多河流開(kāi)展了應(yīng)用［4，17］。

Knack 和Shen 在總結(jié)已有實(shí)驗(yàn)和原型資料的基礎(chǔ)上，提出適應(yīng)冰蓋條件和流凌影響下的推移質(zhì)輸沙率公式［3］?；谠摴胶头瞧胶夥蔷鶆蛏撤匠?，考慮河冰影響的泥沙質(zhì)量守恒方程為［72］：

式中：qbk為k粒徑泥沙單寬推移質(zhì)輸沙率，m2/s；k為粒徑分組；ubk為k粒徑推移質(zhì)運(yùn)動(dòng)速度，m/s；zb為河床高程，m；p′m為河床泥沙的孔隙率；α′bx，α′by為推移質(zhì)在水流和重力沿岸坡分量綜合作用下的方向向量分量。根據(jù)不平衡輸沙的經(jīng)驗(yàn)公式，推移質(zhì)的運(yùn)動(dòng)方程為：

式中：Lbk為k粒徑泥沙不平衡輸移長(zhǎng)度，m；q*bk為k粒徑泥沙平衡輸沙下的單寬推移質(zhì)輸沙率，m2/s。岸坡上的泥沙運(yùn)動(dòng)受重力和水流拖曳力的綜合影響，推移質(zhì)泥沙運(yùn)動(dòng)方向向量的修正計(jì)算公式為［15，72］：

式中：αbx、αby為水流方向向量的分量；β為考慮岸坡和床面泥沙性質(zhì)的經(jīng)驗(yàn)性修正系數(shù)，可以由實(shí)測(cè)資料率定或經(jīng)驗(yàn)公式計(jì)算［73-76］；φx、φy分別為床面沿兩個(gè)方向的傾角。

3.4 河冰動(dòng)力學(xué)模塊河冰動(dòng)力學(xué)模塊的主要控制方程為河冰質(zhì)量守恒方程、河冰面密度連續(xù)性方程和動(dòng)量方程［64］。河冰厚度和面密度由河冰連續(xù)性方程計(jì)算，河冰運(yùn)動(dòng)速度由動(dòng)量方程計(jì)算。河冰的動(dòng)量方程、質(zhì)量守恒方程、河冰質(zhì)量及河冰面密度的連續(xù)性方程分別為：

式中：M為單位面積冰的質(zhì)量，kg；為冰速向量，m/s；為冰內(nèi)部的相互作用力，N；為風(fēng)作用在冰上的力，N；為水流作用在冰上的力，N；為重力在水平方向的分力，N；Em為單位面積冰的增長(zhǎng)速率，由熱力生長(zhǎng)消融或外部源匯控制，kg/s；ρi為冰的密度，kg/m3；ti為冰的厚度，m；Ra為機(jī)械作用引起的冰面積變化，s-1；Ea為冰面積的熱力生長(zhǎng)消融速率，s-1。為了計(jì)算河冰之間的內(nèi)應(yīng)力，河冰模塊采用黏彈塑性應(yīng)力應(yīng)變方程計(jì)算冰塊的內(nèi)應(yīng)力，采用摩爾-庫(kù)倫屈服準(zhǔn)則計(jì)算冰塞冰壩中冰塊的屈服及冰與河岸的臨界摩擦力［77］。

式（9）—式（12）假設(shè)河冰顆粒為連續(xù)介質(zhì)，采用基于拉格朗日?qǐng)龅墓饣Ｗ臃ㄇ蠼夥匠探M。光滑粒子法不受網(wǎng)格的限制，能準(zhǔn)確計(jì)算冰蓋前沿河冰下潛、冰塞體中冰塊的坍塌變形及冰塞的沖蝕和釋放等復(fù)雜過(guò)程，避免了復(fù)雜水冰作用下離散格式的不穩(wěn)定。

3.5 岸灘侵蝕模塊自然河流會(huì)出現(xiàn)周期性的岸坡沖刷變陡、河岸失穩(wěn)破壞、坍塌土體落淤再平衡的演變過(guò)程。大部分研究采用泥沙休止角判斷非黏性沙岸坡的穩(wěn)定與否，通過(guò)原型觀測(cè)或參數(shù)率定確定岸坡穩(wěn)定的臨界坡角，即泥沙休止角［78-80］。Zech 等針對(duì)淹沒(méi)和非淹沒(méi)岸坡不同含水量的岸坡穩(wěn)定性，提出雙泥沙休止角法：采用稍大的動(dòng)泥沙休止角判斷岸坡的失穩(wěn)與否，采用靜泥沙休止角決定崩岸坡面再平衡的穩(wěn)定坡面，采用兩組不同的動(dòng)、靜泥沙休止角分別計(jì)算水面上下的崩岸過(guò)程［81］。本文岸灘侵蝕模塊采用雙泥沙休止角法計(jì)算冰蓋及流凌刮擦影響下的岸灘侵蝕、崩塌和再平衡過(guò)程［71］。計(jì)算方法如下。

（1）臨界失穩(wěn)點(diǎn)判斷。如圖3計(jì)算岸坡各個(gè)網(wǎng)格點(diǎn)的坡角，判斷坡角是否大于相應(yīng)的動(dòng)泥沙休止角，并確定臨界崩塌破壞點(diǎn)的位置。采用式（13）計(jì)算坡面上各點(diǎn)的坡角：

式中αi為橫斷面i點(diǎn)的坡角。通過(guò)對(duì)比各點(diǎn)坡角與動(dòng)泥沙休止角的大小，確定臨界失穩(wěn)點(diǎn)位置，該點(diǎn)以上為不穩(wěn)定坡面，以下為穩(wěn)定坡面。圖3中NT為橫斷面網(wǎng)格點(diǎn)總數(shù)，Nc為臨界失穩(wěn)點(diǎn)。

圖3 橫斷面坡角及臨界失穩(wěn)點(diǎn)示意

（2）確定崩塌失穩(wěn)土體質(zhì)量。以臨界失穩(wěn)點(diǎn)為起點(diǎn)、靜泥沙休止角為坡角，確定崩塌后穩(wěn)定的坡面，計(jì)算崩岸土體的質(zhì)量，具體示意見(jiàn)圖4。圖中，Ns為失穩(wěn)坡面上端起點(diǎn)；Ne為結(jié)束河床調(diào)整的下端網(wǎng)格點(diǎn)；As為崩岸破壞土體面積，m2；Af為崩岸土體重新分布的面積，m2；Δzai為不同網(wǎng)格點(diǎn)崩塌岸灘調(diào)整的高程變化，m；Δzbi為不同網(wǎng)格點(diǎn)崩岸土體再分布后的岸灘調(diào)整高程變化，m。崩塌失穩(wěn)土體的面積采用下式計(jì)算：

（3）失穩(wěn)土體堤腳再分布：在保證土體質(zhì)量守恒的基礎(chǔ)上（As=Af），確定落淤土體分布位置，更新岸坡再平衡后各點(diǎn)的地形高程，見(jiàn)圖4。通過(guò)試算法確定再平衡后堤腳的落淤土體分布。先假定落淤土體以靜泥沙休止角為坡角平鋪在臨界失穩(wěn)點(diǎn)Nc以下，計(jì)算相應(yīng)的再分布土體面積A'f。如果試算的土體面積等于崩塌的土體面積As，則假定的穩(wěn)定坡面為最終落淤坡面。如果As≠A'f，則采用式（15）平行調(diào)整臨界失穩(wěn)點(diǎn)以下的落淤高程變化。

圖4 失穩(wěn)坡面崩塌及再平衡示意

岸灘侵蝕模塊通過(guò)坡面各點(diǎn)坡角與動(dòng)泥沙休止角的對(duì)比，判斷坡面的穩(wěn)定與否；利用靜泥沙休止角確定崩岸坡面，計(jì)算崩塌土體的質(zhì)量；再通過(guò)土體質(zhì)量守恒計(jì)算堤腳落淤土體分布，進(jìn)而為水沙計(jì)算提供更新的河床地形。

3.6 各模塊耦合步驟水冰沙耦合模型實(shí)施步驟和計(jì)算方法是：（1）在滿足顯格式穩(wěn)定性的時(shí)間步長(zhǎng)下運(yùn)行二維水沙數(shù)值模塊，直到給定的耦合時(shí)間，提供計(jì)算區(qū)域的水位、流速、水溫、水流拖曳力、挾沙力和床面高程等水沙要素，并傳遞給河冰動(dòng)力學(xué)模塊；（2）運(yùn)行河冰動(dòng)力學(xué)模塊獲得耦合時(shí)間的冰厚、冰速、冰濃度及冰摩擦力等河冰要素；（3）將前兩步計(jì)算的結(jié)果傳遞給岸灘侵蝕模塊，并計(jì)算水冰沙影響下的岸灘侵蝕速率、岸坡穩(wěn)定坡面及崩岸土體再分布位置，直至給定的耦合時(shí)間；（4）將第三步計(jì)算結(jié)果反饋給二維水沙數(shù)值模塊，校正更新岸坡和地形下的水位、流量和輸沙過(guò)程；（5）重復(fù)以上4 個(gè)步驟至計(jì)算結(jié)束，通過(guò)3 個(gè)模塊的信息傳遞和反饋實(shí)現(xiàn)水冰沙耦合模擬。

4 結(jié)論

受季節(jié)性的河冰影響，北方河流的輸水安全和河道演變需考慮水流、河冰、泥沙和邊界的復(fù)雜相互作用。本文首先總結(jié)了近幾十年來(lái)河冰研究的主要進(jìn)展和基本理論，包括冬季水體失熱、河流產(chǎn)冰、封河過(guò)程、開(kāi)河機(jī)理和河冰影響等五個(gè)方面。基于上述理論基礎(chǔ)，本文提出了二維河流水冰沙耦合數(shù)學(xué)模型，創(chuàng)新性地將復(fù)雜水流變化、河冰運(yùn)動(dòng)、泥沙輸運(yùn)、河床沖淤變形及岸灘崩塌侵蝕過(guò)程耦合起來(lái)。該模型分為二維水沙數(shù)值模塊、河冰動(dòng)力學(xué)模塊和岸灘侵蝕模塊。二維水沙數(shù)值模塊采用非結(jié)構(gòu)的有限元法計(jì)算非恒定水流運(yùn)動(dòng)、泥沙輸移和河道沖淤變化，能適應(yīng)復(fù)雜河道地形條件。河冰動(dòng)力學(xué)模塊采用無(wú)網(wǎng)格的光滑粒子法求解，能模擬冰塞冰壩等復(fù)雜河冰聚集和釋放過(guò)程。岸灘侵蝕模塊采用雙泥沙休止角法計(jì)算，能模擬不同含水層岸坡的崩塌和再平衡過(guò)程。新模型將河冰理論和水沙理論創(chuàng)新性地融合起來(lái)，既滿足河冰模擬需求，也適應(yīng)水沙問(wèn)題研究，能模擬北方河流全季節(jié)和河冰全過(guò)程的水冰沙相互作用。二維水冰沙耦合模型的初步建立為北方河流防凌減災(zāi)和河流管理提供了有力的技術(shù)支撐。

二維水冰沙耦合模型主要針對(duì)沙質(zhì)河道和岸灘而開(kāi)發(fā)，尚不包括黏性土質(zhì)岸坡和河床沖淤計(jì)算模塊，也不包括岸灘植被分析功能，不適于黏性岸坡穩(wěn)定性分析和密集植被影響的河道模擬。針對(duì)非黏性土和黏性土混合的復(fù)雜河道條件，還需進(jìn)一步拓展該模型的水沙數(shù)值模塊和岸灘侵蝕模塊。

關(guān)于河流水冰沙耦合模型研究的系列文章分為兩篇，本文主要介紹原理和方法，模型具體驗(yàn)證與應(yīng)用將在下篇給出。