河流水冰沙耦合模型研究Ⅰ：原理和方法

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

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

1 研究背景

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

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

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

2 河冰水力學理論框架

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

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

圖1 河冰水力學理論框架

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

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

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

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

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

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

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

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

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

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

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

受溫升融水和降雨的影響，春汛引起的開河冰壩在北方河流最為常見，相應的凌汛災害也更嚴重?；诤颖鶆恿W理論的數(shù)學模型被成功用于日本渚滑川和加拿大圣約翰河等河流的冰壩過程模擬［54-55］，但開河相關的力學機制尚待完善。開河期碎冰塊間歇性的運動和停滯易誘發(fā)鏈式冰壩的形成和釋放，伴隨的凌汛洪水風險更高。開河冰壩最大的難題是冰壩坍塌釋放機理和碎冰塊與冰蓋間的相互作用機制，包括冰塊與冰蓋間的相互作用力及水對冰體的作用力等。此外，上游下泄的碎冰塊可能撞擊擠壓下游穩(wěn)定冰蓋，進而導致更多冰蓋破裂，相關力學機理還需進一步研究。

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

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

3 水冰沙耦合模型

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

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

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

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

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

3.3 二維水沙數(shù)值模塊二維水沙數(shù)值模塊包括二維淺水運動、泥沙輸運和河床變形三個部分?；谏蠈痈”拖聦铀w的雙層流體系，考慮河冰阻力和質量引起的水位流量變化，北方河流河冰影響下的二維淺水方程組為［69］：

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

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

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

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

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

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

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

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

式（9）—式（12）假設河冰顆粒為連續(xù)介質，采用基于拉格朗日場的光滑粒子法求解方程組。光滑粒子法不受網格的限制，能準確計算冰蓋前沿河冰下潛、冰塞體中冰塊的坍塌變形及冰塞的沖蝕和釋放等復雜過程，避免了復雜水冰作用下離散格式的不穩(wěn)定。

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

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

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

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

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

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

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

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

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

4 結論

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

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

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