基于離散元方法的河冰動力過程短期預測模型及應用

翟必垚 張寶森 季順迎 田治宗

摘 要：秋末冬初時黃河河道平穩(wěn)封河對于預防凌汛意義重大，而什四份子作為黃河內(nèi)蒙古河段主要首封河段備受重視。根據(jù)寒區(qū)河道具體冰情，基于河冰離散元方法，考慮河流表面冰的動力輸移和河流水動力學的相互耦合作用，建立了針對特定重點河段的短期預測模型，可根據(jù)實際水力要素、氣溫、現(xiàn)場冰情等條件進行河冰輸移、封河開河等情況的模擬和預測。模型考慮了河道的岸冰邊界以及上游來冰情況，包括來冰的密集度、冰塊體的尺寸分布、冰塊體的凍結強度等因素對河冰輸移和卡冰封河過程的影響。基于該預測模型，針對什四份子彎道在不同工況下的封河情況進行研究，結合現(xiàn)場情況探究了流速、來冰量、冰塊尺寸、岸冰尺寸等因素對卡冰封河的影響，從冰流量和河道輸冰量的角度解釋了各要素的影響機理，為特定河段的封河預測工作提供參考。

關鍵詞：河冰;封河;離散元方法;數(shù)值模型;什四份子彎道

中圖分類號：TV133;TV882.1

文獻標志碼：A

doi：10.3969/j.issn.1000-1379.2021.10.015

引用格式：翟必垚，張寶森，季順迎，等.基于離散元方法的河冰動力過程短期預測模型及應用[J].人民黃河，2021，43（10）：75-80.

Abstract： Based on the discrete element method， a short-term prediction model of the nature river ice dynamic process was established， considering the coupling effect between surface ice and the river hydrodynamics. It could simulate and predict the several river ice dynamic processes， such as， river ice transport， ice bridge formation， ice jam， break up and jam release. Based on the prediction model， the formation of ice bridge in Shensifenzi reach was simulated and the effects of flow velocity， ice supplement， ice size and border ice size on the jamming were investigated. The mechanism of each factor was explained by ice discharge and ice transport capacity of the channel. This paper provided a reference for the forecast work of ice bridge formation and river freeze up in natural rivers.

Key words： river ice; freeze up; discrete element method; numerical model; Shensifenzi reach

冰凌在我國北方寒冷地區(qū)的江河非常常見，是氣象水文、河道地形地勢和水力條件綜合影響的結果。嚴重的冰情會直接影響江河的航運交通、灌溉供水、水力發(fā)電，帶來巨大的經(jīng)濟損失，還會誘發(fā)凌汛災害，對人民生命財產(chǎn)造成巨大的損失[1-3]。秋末冬初時節(jié)，氣溫降低，江河出現(xiàn)流凌，隨著低溫的持續(xù)，會發(fā)生封河現(xiàn)象。一旦封河，應及時減小流量爭取平封，防止立封和冰塞的產(chǎn)生[4]。什四份子彎道是黃河內(nèi)蒙古河段高頻的首封河段，其地勢形態(tài)特征典型，上寬下窄，彎變急且角度變化大，是防凌重點關注的位置[5-6]。該河段已安裝配套了多種監(jiān)測設備進行冰情、氣候、水力要素等的實時監(jiān)測，并配合無人機、衛(wèi)星雷達定期觀察冰情發(fā)展，為冰情預測提供了充足的參考依據(jù)[7-9]。

國內(nèi)外的學者基于冰力學理論基礎開發(fā)了多種河冰熱力學、動力學過程模擬的數(shù)學模型和數(shù)值模型，較為系統(tǒng)的有RICE[10]、DynRICE[11-13]、RIVJAM[14]、ICEJAM[15] 、ICESIM[16]、RIVICE[17]、River1D[18]、HEC-RAS[19]等模型。已有的河冰模型從河冰動力學、熱力學機理出發(fā)，能夠很好地模擬長時段、長距離的河冰熱力生消、動力輸移等演化過程，對冰塞冰壩的出現(xiàn)位置及引起的水位、流量變化有較好的模擬結果[20-22]。近年來隨著人工智能技術的發(fā)展，統(tǒng)計模型、人工神經(jīng)網(wǎng)絡、模糊數(shù)學等方法也在冰情預測預報中得到了應用，并取得了很好的效果[23-25]。然而，以上模型和方法在實際應用中可能會受到長距離、多時段監(jiān)測數(shù)據(jù)不完備的限制，且關注點多集中在大中尺度的水位、流量和冰塞位置上，針對特定局部河段的冰情精細化模擬和預測相對較少。

黃河防凌監(jiān)測技術發(fā)展迅速，已由早年的監(jiān)測手段少、監(jiān)測指標單一發(fā)展形成了初步的冰情觀測技術體系，測量指標包含流速、水深、冰厚、水面冰密集度、水溫、冰溫等[21]，監(jiān)測手段擴展到了衛(wèi)星、雷達、無人機、遠程4G視頻拍攝等[26]。這些都為精細化預測模型的參數(shù)設置提供了依據(jù)。

本研究基于河冰離散元模型，建立面向真實河道的數(shù)值模型，根據(jù)特定天然河道具體冰情，對封河前期的水力要素、地勢地形、初始冰情等重要條件建立模擬參數(shù)，形成針對小范圍特定重點河段的短期預測模型;結合什四份子險工處已建立的完備監(jiān)測手段，對該河段河冰輸移及封河情況進行模擬研究，以期為防凌預報提供有效的參考。

1 基于河冰離散元方法的河冰動力過程短期

預測模型

1.1 河冰動力學短期預測模型的計算框架

結合特定河段的監(jiān)測數(shù)據(jù)，基于河冰離散元方法，建立特定河道的短期預測模型，可對特定局部河段的河冰動力學過程進行精細化模擬，具體模型框架如圖1所示。

由于熱力計算的輸入條件較為復雜，涉及日照、氣溫、水溫等條件，且熱力計算多針對長時間尺度的模擬、誤差不易控制，因此模型回避了熱力學計算，將河冰的熱力學生效結果轉化為動力學模擬的初始條件，如岸冰邊界、冰塊凍結強度等，再通過短期的動力學模擬以提高模型的適用性和準確性。完備的河段冰情資料包含氣象、水文、河道地形、冰情等條件，其可對模型的正確性進行驗證，也可對模型計算的初始參數(shù)進行校核，進而對特定的河冰現(xiàn)象如封河、開河、冰壩、潰壩等進行模擬計算，為河冰現(xiàn)象的預報提供參考。

1.2 河冰動力過程的離散元方法

河冰動力過程計算模型考慮了河冰動力學與水動力學的耦合作用[27]。河冰動力學采用了二維河冰輸移模型，可計算河冰受到的水流和風的拖曳力、水面梯度力、冰塊內(nèi)力、冰與河岸及建筑物的接觸力;河流水力學采用Shen根據(jù)表面冰層作用改進的二維淺水模型[11];河冰內(nèi)力主要采用離散元方法顯式逐步求解[28]。此外，該模型考慮了河冰離散單元的浮力、拖曳力以及水流與冰層的耦合[27]。

根據(jù)河冰在低溫條件下的凍結現(xiàn)象，建立了河冰離散元黏聚力模型，即計算單元之間建立黏結對，使得單元和單元之間可以相互黏結在一起，黏結對可以在單元之間傳遞拉力、壓力、剪切力、彎矩和扭矩，當傳遞的力超過破壞準則時，黏結對就會失效，宏觀表現(xiàn)為凍結冰的破碎。擴展多面體單元的黏聚力模型有兩種類型（見圖2）。一種是建立在共面接觸的黏結單元間，即必須是面面接觸模式才可建立黏結對[28-29]，黏結面積即為接觸面積。另一種是針對自然狀態(tài)下冰堆積體中的冰塊接觸隨機不規(guī)則的情況，建立的基于任意接觸方式離散元黏聚力模型，這種黏聚力模型在單元之間建立虛擬黏結面，虛擬黏結面的面積直接影響?zhàn)そY對傳遞的力的大小，其主要受到氣溫T和累計接觸時間ta控制，其表征了冰塊由接觸到凍結以及凍結面生長的過程。具體關系需要根據(jù)特定河道的現(xiàn)場試驗結果或者相關規(guī)范進行選取。

根據(jù)河冰在動力過程中表現(xiàn)出的離散分布特性，采用Voronoi切割算法結合擴展多面體離散單元形態(tài)對冰塊體進行構造。二維Voronoi切割算法是在一個二維平面內(nèi)隨機撒點，再通過隨機點將該二維平面劃分為隨機的任意多面形，對多面形賦予冰厚便可得到隨機分布的碎冰模型。碎冰的尺寸、形狀以及碎冰分布的密集度可以根據(jù)隨機點的分布參數(shù)進行調(diào)節(jié)[30]。

1.3 河冰動力過程離散元模擬的邊界條件

初始岸冰的構造主要依據(jù)現(xiàn)場觀測的岸冰形態(tài)，通過繪制岸冰邊界，再采用隨機三角單元對岸冰區(qū)域進行剖分，組成岸冰的三角冰塊初始狀態(tài)是相互黏結的，由此得到岸冰模型。在岸冰與河岸接觸的位置采用固接的方式將二者黏結。由于岸冰的熱力生成過程長、凍結時間久，通常有較大的強度，因此岸冰區(qū)域的冰單元黏結強度通常設置為較大值。

根據(jù)不同算例工況對冰塊補給方式的需求，模型中給出了兩種冰塊補給方法，分別為周期邊界法和上游供給法[31]。周期邊界法能夠大大節(jié)省程序計算內(nèi)存以及消除邊界對顆粒流動的影響，被廣泛應用于離散元模擬中。周期邊界法可應用在規(guī)則河道的兩岸以消除邊壁效應，也可應用于河道上下游邊界處，以節(jié)省程序預構造冰塊數(shù)量。這里主要針對后者的應用進行說明，其實質(zhì)為當冰單元從河道下游邊界流出后，會從上游重新進入河道區(qū)域。圖3為周期邊界法中的河冰輸移。為了保證邊界處的冰單元可以正確地搜索并接觸到另一端邊界處的冰單元，根據(jù)邊界附近的冰單元（Z1，Z2，Z3，Z4）在另一邊界處建立相應的虛擬冰單元（V1，V2，V3，V4）。如果河道為規(guī)則直河道，則可通過平移實現(xiàn);當河道上下游邊界方向不一致時，可采用四元數(shù)控制旋轉角度的方法通過平移和旋轉來實現(xiàn)。

上游供給法是另一種用來補給河道冰量的方法，即在河道的上游邊界或者特定位置進行輸冰。程序初始化階段預先構造足夠多的冰塊單元，開辟相應內(nèi)存但不參與計算。在河道上游邊界或者特定位置設置來冰邊界（見圖4），確定來冰邊界輸冰個數(shù)Nb和輸冰間隔dl0。邊界輸冰個數(shù)Nb是指輸冰邊界處一個時間間隔最多能夠輸送Nb塊碎冰，冰塊的輸冰間距dB0 =B/Nb;輸冰間隔dl0可以控制當輸送的冰塊運動多遠距離后開始進行下一輪輸冰。由于在輸入冰單元時必須保證輸入冰塊之間以及輸入冰塊與河道已有冰塊之間沒有接觸重疊，這里采用了最大顆粒球形包圍盒尺寸Dmax控制冰塊間隔，即dB0

周期邊界法多用于規(guī)則河道的算例研究，而上游供給法冰塊由上游邊界輸送，更符合實際情況，在自然河道算例中推薦采用。不論采用哪種冰塊補給方式，計算中都會根據(jù)初始時刻實際河道的冰分布情況，預先在河道內(nèi)布置一定量的冰塊，這種方式可減少冰塊由上游邊界流動到占據(jù)整個河道的時間，大大提高了計算效率。

2 黃河內(nèi)蒙古什四份子河段封河過程的離散元模擬

2.1 黃河內(nèi)蒙古什四份子河段概況

什四份子河段地處黃河內(nèi)蒙古河段下游（北緯40°17′39″，東經(jīng)111°2′53″），見圖5。該河段地形特征明顯，由上游的寬河道在入彎處逐漸束窄，河道走勢方向也發(fā)生了劇烈變化，因此該位置較易產(chǎn)生卡冰現(xiàn)象，歷史上多次成為秋冬封河期首封位置。

根據(jù)2019—2020年該河段封河過程的觀測資料，2019年什四份子封河時間為12月6日，封河流量約為787 m3/s，平均流速約為0.71 m/s，封河當日氣溫為1.4～-15.3 ℃，平均氣溫為-10.1 ℃;而2020年封河時間為12月8日，封河流量約為362 m3/s，當日氣溫為-7.2～-18.2 ℃，平均氣溫為-13.6 ℃，平均流速約為0.38 m/s。這兩年封河當日冰況差別較大，圖6為什四份子河段2019年和2020年封河過程的現(xiàn)場航拍圖。2019年什四份子河段封河之前上游來冰量較大、冰塊尺寸小且密集，河岸邊已有一定數(shù)量的岸冰生成;2020年封河之前來冰呈大塊片狀，冰塊稀疏但是尺寸較大，在持續(xù)低溫的作用下岸冰生長也更加明顯，與2019年情況相比更加束窄了河道。

2.2 短期封河的離散元數(shù)值預測

根據(jù)2019年和2020年封河前一天的水力情況和冰情，建立模型模擬封河現(xiàn)象，模擬參數(shù)見表1。模型中水力邊界條件上游為流量，下游為水位，平均糙率0.03。圖7為計算網(wǎng)格及河床高程。采用恒定流研究河冰的輸移特性，2019年工況設定初始流量787 m3/s，下游水位988.2 m，穩(wěn)定流態(tài)如圖8（a）所示，主槽的平均流速為0.7 m/s，上游輸入冰塊的平均面積為500 m2，厚度為0.2 m，密集度為0.8左右;2020年工況將初始流量調(diào)整為362 m3/s，得到的穩(wěn)定流態(tài)如圖8（b）所示，主槽平均流速約為0.4 m/s，輸入冰塊平均面積為5 000 m2，厚度為0.2 m，密集度為0.3左右。

圖9為這兩年封河模擬的結果。可見在相應水力條件和來冰量的情況下，冰塊在河道束窄位置卡冰，模擬結果與實際觀測現(xiàn)象（見圖6）較為吻合。從結果中可以看出岸冰的束窄作用決定了卡冰的位置并直接導致卡冰的發(fā)生，冰塊大小和流凌密度也影響了卡冰難易程度。2019年什四份子河段是由小尺寸高密度來冰在大流量下發(fā)生卡冰封河的，而2020年該河段是在小流量低流速的情況下因大尺寸的來冰以及岸冰對河道的束窄作用最終導致卡冰封河的。

圖10為2019年封河模擬結果中統(tǒng)計的彎道出口處（圖9中紅虛線所示）冰流量變化情況。在流凌初期河道將來冰順暢地輸往下游，隨著來冰量的持續(xù)增加，冰流量隨之增加并到達穩(wěn)定值附近，該穩(wěn)定值即反映了河道的輸冰能力。在1.0～2.3 h內(nèi)，主流繼續(xù)保持最大輸冰能力輸冰，出彎處冰流量基本不變，而超過河道輸冰能力的冰塊會在主流兩側停留、聚集，最終導致河冰在束窄位置出現(xiàn)擁堵。一旦出現(xiàn)卡堵的情況，河道輸冰能力會迅速降低，冰流量隨之減小，進而卡堵封河。由圖10可知不足0.5 h時間內(nèi)出彎處的冰流量由穩(wěn)定值降為0。2020年的冰流量變化情況與2019年較為類似，這里不做贅述。

2.3 河冰封河過程的主要影響因素分析

基于對實際工況的分析得知，影響封河的因素較為復雜，包括水力要素（主要是水流流速）、冰量、冰尺寸以及岸冰的邊界約束等。為了理解卡冰封河潛在的物理機制，通過數(shù)值模型對更多工況進行模擬，進而對封河的條件建立一定認識。

首先考慮無岸冰的工況。圖11是密集度為0.4和0.8的來冰輸移情況，分別代表低密度來冰量和高密度來冰量。當密集度較小時，冰塊主要在河道主槽輸移，只有極少量的碎冰停留在岸邊低速區(qū)。當密集度增大時，主槽的輸冰能力無法輸送大量來冰，會有部分冰塊被排擠至主槽兩側的低流速區(qū)域，在彎道凹岸處逗留，并且隨著來冰的持續(xù)增加，凹岸的冰塊不斷聚集，無法及時流出。二者都沒有出現(xiàn)卡冰封河現(xiàn)象。

當主槽的平均流速降到0.4 m/s左右時，輸移情況沒有太大變化，河道仍沒有發(fā)生卡冰封河現(xiàn)象。由于流速降低，主槽的輸冰能力隨之降低，因此有更多的冰塊在河道兩側停留。

在無岸冰的幾種工況下都未發(fā)生卡冰封河，可見岸冰的生長對卡冰封河有必要性。流凌初期在持續(xù)輸冰下，會有部分冰塊在凹岸側停留，來冰量越大，進入凹岸低速區(qū)的冰塊越多，這給什四份子彎道凹岸處沖擊岸冰和平滑岸冰的生長提供了基礎，進而為封河的發(fā)生提供了有利條件。由此得出：在流凌初期該處河道卡冰封河的可能性較小，這時應主要關注岸冰的生長發(fā)展趨勢。

其次，根據(jù)2019年、2020年兩種岸冰形態(tài)，探究在岸冰存在的條件下，不同流速和來冰量的封河情況。結果顯示，在同樣的條件下，密集度越大、流速越慢時，越容易卡冰。這里統(tǒng)計了兩種岸冰狀態(tài)下0.1、0.4、0.8、1.2 m/s幾種流速封河時的臨界來冰密集度，見圖12。較高的流速需要較高的臨界來冰密集度才能滿足卡冰封河的條件，并且岸冰二較岸冰一的尺寸更大，其對河道的束窄作用更明顯，因而對應的封河臨界來冰密集度較小，更易形成卡冰封河。

3 結 論

針對天然河道及工程條件，建立了基于河冰離散元方法的短期河冰預測模型，可準確表征實際河道的地形地勢、水力要素、冰況信息。將該方法應用于黃河內(nèi)蒙古段什四份子彎道處河冰輸移的模擬和研究中，結合2019年、2020年兩年封河前的工況條件進行封河預報模擬，得到的封河結果和卡冰位置與實際情況吻合較好。進一步對該河段的封河條件進行研究，建立了兩種特定岸冰形態(tài)和特定冰尺寸下不同水流速度對應封河臨界來冰密集度的關系，得到了一些定性的結論：河道的邊界直接影響卡冰封河的可能性，什四份子河段為由上游寬大到下游狹窄的形態(tài)，流冰容易在彎道凹岸的低速區(qū)停留，形成岸冰。岸冰的形成是卡冰的必要條件之一，岸冰對河道束窄得越嚴重，越容易造成卡冰;大尺寸的冰塊加速了卡冰封河的過程;流速決定河道的輸冰能力并影響卡冰的穩(wěn)定性，高流速下河道輸冰能力強，冰塊會沿主槽向下游輸移，或直接被帶入水底、流往下游，不易穩(wěn)定封河，低流速下河道輸冰能力弱，冰塊容易在束窄處聚集形成堆積之勢，進而穩(wěn)定，攔截來冰形成封河。本研究為河冰現(xiàn)象短期預測提供了新的方法，對封河條件的研究成果也可以為封河預警工作提供一定的參考。

參考文獻：

[1] 楊開林.河渠冰水力學、冰情觀測與預報研究進展[J].水利學報，2018，49（1）：81-91.

[2] 王軍，陳胖胖，江濤，等.冰蓋下冰塞堆積的數(shù)值模擬[J].水利學報，2009，40（3）：348-354.

[3] BELTAOS S. Hydrodynamics of Storage Release During River Ice Breakup[J]. Cold Regions Science & Technology， 2017，139： 36-50.

[4] 蔡琳.中國江河冰凌[M].鄭州：黃河水利出版社，2008：5-7.

[5] 羅紅春，冀鴻蘭，郜國明，等.黃河什四份子彎道冰期水流及冰塞特征研究[J].水利學報，2020，51（9）：1089-1100.

[6] 趙水霞，李暢游，李超，等.黃河什四份子彎道河冰生消及冰塞形成過程分析[J].水利學報，2017，48（3）：351-358.

[7] 武桂芝，張寶森， 李春江，等.陣列地質(zhì)雷達在黃河堤防隱患探測中的應用[J].人民黃河，2020，42（8）：113-116.

[8] 劉輝，冀鴻蘭，牟獻友，等.基于無人機載雷達技術的黃河冰厚監(jiān)測試驗[J].南水北調(diào)與水利科技（中英文），2020，18（3）：217-224.

[9] 張磊，牟獻友，冀鴻蘭，等.基于多波段遙感數(shù)據(jù)的庫區(qū)水深反演研究[J].水利學報，2018，49（5）：639-647.

[10] LAL W，SHEN H T. A Mathematical Model for River Ice Processes[J].Journal of Hydraulic Engineering， 1991，117（7）： 851-867.

[11] SHEN H T， SU J， LIU L. SPH Simulation of River Ice Dynamics[J].Journal of Computational Physics， 2000，165（2）： 752-770.

[12] SHEN H T ， WANG D S， LAL W. Numerical Simulation of River Ice Processes[J].Journal of Cold Regions Engineering， 1995， 9（3）：107-118.

[13] WANG D S， SHEN H T， CRISSMAN R D. Simulation and Analysis of Upper Niagara River Ice Jam Conditions[J].Journal of Cold Regions Engineering， 1995， 9（3）：119-134.

[14] BELTAOS S. Numerical Computation of River Ice Jams[J]. Canadian Journal of Civil Engineering， 1993，20（1）： 88-99.

[15] FLATO G M， GERARD R. Calculation of Ice Jam Thickness Profiles[C]//Proceedings of 4th Workshop on Hydraulics of River Ice. Montreal， Quebec， Canada， 1986（C3）：1-25.

[16] CARSON R， BELTAOS S， GROENEVELD J， et al. Comparative Testing of Numerical Models of River Ice Jams[J].Canadian Journal of Civil Engineering， 2011，38（6）：669-678.

[17] LINDENSCHMIDT K， XU J Y. Rivice-A Non-Proprietary， Open-Source， One-Dimensional River-Ice Model[J].Water， 2017， 9（5）：314.

[18] BLACKBURN J， SHE Y. A Comprehensive Public-Domain River Ice Process Model and Its Application to a Complex Natural River[J].Cold Regions Science and Technology， 2019， 163：44-58.

[19] BRUNNER G W. HEC-RAS（River Analysis System）[C]//North American Water and Environment Congress & Destructive Water. Washington D C： ASCE， 2002：3782-3787.

[20] SHEN H T. Mathematical Modeling of River Ice Processes[J].Cold Regions Science and Technology， 2010. 62（1）： 3-13.

[21] 郭新蕾，王濤，付輝.河渠冰水力學研究進展和動向[J].力學學報，2021，53（3）：655-671.

[22] BELTAOS S. Progress in the Study and Management of River Ice Jams[J].Cold Regions Science and Technology， 2008， 51（1）：2-19.

[23] GUO X， WANG T， FU H， et al. Ice-Jam Forecasting During River Breakup Based on Neural Network Theory[J].Journal of Cold Regions Engineering， 2018， 32（3）：1-8.

[24] WANG T，YANG K，GUO Y.Application of Artificial Neural Networks to Forecasting Ice Conditions of the Yellow River in the Inner Mongolia Reach[J].Journal of Hydrologic Engineering， 2008， 13（9）：811-816.

[25] MORSE B， HESSAMI M，BOUREL C. Mapping Environmental Conditions in the St. Lawrence River onto Ice Parameters Using Artificial Neural Networks to Predict Ice Jams[J].Canadian Journal of Civil Engineering， 2003， 30（4）：758-765.

[26] 張寶森，張興紅.冰凌監(jiān)測新技術研究進展[J].中國防汛抗旱，2017，27（6）：14-19，33.

[27] 翟必垚，劉璐，張寶森，等.基于離散元方法與水動力學耦合的河冰動力學模型[J].水利學報，2020，51（5）：617-630.

[28] 劉璐，尹振宇，季順迎.船舶與海洋平臺結構冰載荷的高性能擴展多面體離散元方法[J].力學學報，2019，51（6）：1720-1739.

[29] LIU L， JI S. Bond and Fracture Model in Dilated Polyhedral DEM and Its Application to Simulate Breakage of Brittle Materials[J].Granular Matter， 2019，21（3）：1-16.

[30] SOTOMAYOR O E，TIPPUR H V. Role of Cell Regularity and Relative Density on Elasto-Plastic Compression Response of Random Honeycombs Generated Using Voronoi Diagrams [J].International Journal of Solids and Structures， 2014，51（21-22）：3776-3786.

[31] BABIC M. Discrete Element Simulations of River Ice Transport[C]//Proceedings of IAHR the 10th International Symposium on Ice. Espoo， Finland： IAHR，1990：564-574.

【責任編輯 許立新】