劉煜，吳輝碇

（國家海洋環(huán)境預(yù)報中心，海洋災(zāi)害預(yù)報技術(shù)重點實驗室，北京 100081）

第2講 海冰動力學(xué)

劉煜，吳輝碇

（國家海洋環(huán)境預(yù)報中心，海洋災(zāi)害預(yù)報技術(shù)重點實驗室，北京 100081）

簡單回顧了海冰動力學(xué)的發(fā)展，介紹了海冰漂移、形變和堆積等動力學(xué)過程，對比分析了不同海冰本構(gòu)關(guān)系在海冰流變學(xué)的應(yīng)用。

海冰；動力學(xué)；流變學(xué)；本構(gòu)關(guān)系

1 引言

影響海冰發(fā)生和發(fā)展的物理過程主要分為熱力學(xué)、動力學(xué)兩種過程。大氣和海洋間的熱量交換造成海水凍結(jié)和海冰形成。海冰厚度增長和消融的長期趨勢主要依賴于大氣、海洋和海冰間的能量交換，這些過程是構(gòu)成海冰熱力學(xué)模式的基礎(chǔ)。在極地和區(qū)域海域，海冰漂移和形變主要依賴于動力過程，受其影響，海冰冰脊和水道的形成可以在短時間內(nèi)明顯改變海冰厚度及分布。海冰動力過程主要取決于氣-冰-海間動量交換以及冰內(nèi)應(yīng)力，上述過程構(gòu)成了海冰動力學(xué)模式的基礎(chǔ)。

海冰運動主要受大尺度氣-冰-海相互作用影響，是決定海冰分布和形變的主要因素。海冰動力學(xué)主要包括海洋、大氣動力強迫下海冰的運動變化、動量傳輸、冰內(nèi)相互作用以及海冰堆積、脊化、破碎、斷裂等相關(guān)過程，主要分為以下3個方面：（1）海冰漂移的動量平衡；（2）決定海冰內(nèi)應(yīng)力和形變、冰強度關(guān)系的流變學(xué)；（3）由海冰分布所決定的海冰強度。海冰運動受到以下5種力作用：風(fēng)應(yīng)力、科氏力、冰內(nèi)應(yīng)力、海流應(yīng)力和海面高度梯度力。海冰內(nèi)應(yīng)力的處理歷來是計算海冰動力學(xué)模式的難點和關(guān)鍵。

早期海冰動力學(xué)研究主要探討海冰自由漂移的規(guī)律，并不注重海冰之間相互作用的研究[1-3]。其后海冰流變學(xué)的復(fù)雜性逐漸成為動力學(xué)研究的重點，主要包括：將海冰作為線性粘性流體[4-5]、牛頓粘性流體[6]或塑性物質(zhì)。20世紀(jì)70年代，北極冰動力學(xué)聯(lián)合試驗（Arctic Ice Dynamics Joint Experiment，ATDJEX）首次提出彈-塑性流變學(xué)[7]，于是非線性塑性流變學(xué)發(fā)展起來[8-9]。1979年，Hibler提出粘-塑性（V-P）流變學(xué)海冰模式，該模式已經(jīng)作為眾多海冰模擬和研究的基礎(chǔ)，同時，它也成功地應(yīng)用于北極海冰季節(jié)循環(huán)特征的數(shù)值模擬。

在V-P模式中，方程采用迭代方法求解[5,10]，有效地解決了因計算穩(wěn)定性限制而被迫在海冰密集區(qū)采用的時間步長強約束條件。當(dāng)計算網(wǎng)格加密時，收斂速度非常緩慢。

粘-塑性海冰本構(gòu)關(guān)系主要因模擬極區(qū)海冰變化對全球氣候影響而建立。但是因該模型較好地描述了大、中尺度的海冰由于受重疊、堆積、斷裂等作用導(dǎo)致的顯著流變學(xué)特征，同時省略了計算中海冰彈性變形作用產(chǎn)生的影響，將冰應(yīng)力僅作為應(yīng)變率函數(shù)，因此對大、中尺度下中長期海冰數(shù)值模擬和預(yù)報具有穩(wěn)定性好，計算量較小等優(yōu)點，能很好地再現(xiàn)大應(yīng)變率下海冰所具有的塑性流動行為。目前，這種基于粘-塑性本構(gòu)關(guān)系的動力學(xué)模式在包括極區(qū)及副極區(qū)在內(nèi)的大中尺度海冰數(shù)值模擬及預(yù)報中得到了非常廣泛的應(yīng)用。

自然界中，普遍存在冰隙或冰間水道，大尺度海冰一般承受不了拉應(yīng)力。為此，對粘-塑性模型中的橢圓屈服函數(shù)進(jìn)行平截，可以避免拉應(yīng)力產(chǎn)生[5]。此種方法可以簡化粘-塑性模型在輻散過程中對拉伸和剪切應(yīng)力的計算，實際情況與計算結(jié)果相當(dāng)吻合。Zhang[11]在Baltic海冰數(shù)值模擬中也采用了平截橢圓屈服函數(shù)的方法，取得很好的計算結(jié)果。

因為采用了半隱式求解方法，V-P模式在1d時間尺度內(nèi)，對于因明顯變化的強迫場導(dǎo)致的風(fēng)應(yīng)力反映不太準(zhǔn)確[12]。Hunke等[13]在粘-塑海冰模式中使用條件共扼梯度法改進(jìn)模式計算效率。另外在流變學(xué)中加入了彈性特征，即為彈-粘-塑（E-V-P）海冰模式，E-V-P海冰模式中彈性并不是物理意義的彈性，只是為了計算方便而引入的。這種處理方法避免了以往彈-塑模式[14-15]的復(fù)雜性。和V-P模式的對比試驗表明[16]：兩種模式在長時間尺度積分時，效果基本一致。模式結(jié)果與北極浮標(biāo)資料對比表明，E-V-P模式對短時間內(nèi)強天氣尺度強迫的反映更加準(zhǔn)確和迅速。在高密集度海冰覆蓋海域，E-V-P模式彈性波對海冰狀態(tài)的改變并不明顯。兩種模式冰內(nèi)應(yīng)力基本相同：在海冰高密集度海域表現(xiàn)為粘-塑特性，在低密集度海域則顯示出近似于自由漂流的特征。在氣候系統(tǒng)模式中，E-V-P模式比V-P模式具有更好應(yīng)用前景。

2 海冰漂移

（1）海冰漂移的動量平衡

位于海洋和大氣之間的海冰，受到因冰-氣界面湍流運動產(chǎn)生的大氣拖曳力、海冰-海水之間的拖曳力、冰區(qū)內(nèi)部不同冰塊間和冰塊內(nèi)部的張力以及海面因傾斜造成的應(yīng)力和海冰運動所產(chǎn)生的科氏力影響。海冰動量平衡表示如下：

式中：m為海冰在單位面積內(nèi)的質(zhì)量。τ→a和τ→w分別是大氣和海水對海冰的應(yīng)力。H是海洋表面動力高度，-mg?H為傾斜海面梯度力，主要因為海冰的重力位勢存在水平的差異而產(chǎn)生的應(yīng)力。F→是由于冰內(nèi)應(yīng)力的變化所引起的內(nèi)力。

通過觀測冬季風(fēng)應(yīng)力和水應(yīng)力，Hunkins[17]做出了海冰動量平衡圖（見圖1）。在海冰受力平衡中，風(fēng)應(yīng)力、水應(yīng)力、科氏力和冰內(nèi)力是主要的。風(fēng)應(yīng)力和水應(yīng)力的量級為0.1 N/m2，如果3 m厚冰的冰速為0.1 m/s，則科氏力是0.05 N/m2。3 m厚的冰15 min內(nèi)可以加速到0.1 m/s，產(chǎn)生慣性力0.1 N/m2。不過慣性振蕩引起的加速度項會很大。因為每km冰速變化0.1 m/s需0.1 N/m2力，因此動量平流項很小。海面動力高度梯度項很小，在局地問題研究中可以忽略不計。

研究者們對上述各項的量級大小持有不同的看法。Rothrock[18]對式（1）進(jìn)行尺度分析：海水應(yīng)力和風(fēng)應(yīng)力比時間變化項大3個數(shù)量級，后者可以忽略。在海冰動力學(xué)研究中發(fā)現(xiàn)，各種力平衡的觀測數(shù)據(jù)和研究結(jié)果較為一致。風(fēng)應(yīng)力和海水應(yīng)力是最主要作用力。海冰速度的局地變化、海面傾斜引起的梯度力以及非線性平流項的影響相對比較小，Coriolis力對海冰速度方向影響較大，而對海冰速度值影響較小。在不同條件下，海冰內(nèi)應(yīng)力所產(chǎn)生的影響差異較大[19]。近年來，對于海冰內(nèi)應(yīng)力（本質(zhì)是海冰流變學(xué)）的研究得到了廣泛關(guān)注。在動量平衡中，冰內(nèi)力的作用不容忽視。圖2表示在大氣邊界層的風(fēng)場和冰下海洋Ekman層海流的作用下，通過冰內(nèi)應(yīng)力響應(yīng)海冰的漂移。

在海冰受力平衡機制中穩(wěn)定平均流的作用較小，僅為風(fēng)驅(qū)動作用的百分之幾。然而在長時間的情況下，流和水位梯度的作用都不能忽視。圖3是3個北極漂流站兩年位置變化的觀測和模擬，清楚表明海流和水位梯度的作用。不考慮冰內(nèi)應(yīng)力作用，即所謂“自由漂移”?？疾齑蟪叨群１\動對地轉(zhuǎn)風(fēng)的響應(yīng)時[20]，發(fā)現(xiàn)北冰洋海冰運動對數(shù)月長時間平均海洋環(huán)流和地轉(zhuǎn)風(fēng)具有同等作用；對數(shù)天到1個月短時間尺度地轉(zhuǎn)風(fēng)決定了冰速方差的70%。Martinson等[21]在研究南大洋浮冰塊漂移時，發(fā)現(xiàn)冰漂移速度大約為風(fēng)速的3%，約向其左偏23°。

圖1 根據(jù)觀測所作的動量平衡（引自文獻(xiàn)[17]）

圖2 冰在大氣風(fēng)和海流驅(qū)動下通過冰內(nèi)應(yīng)力響應(yīng)而漂移

（2）海冰與大氣及海洋的動量交換

在海冰受力平衡中，風(fēng)應(yīng)力和海水應(yīng)力是最大項。冰氣界面的動量湍流垂直輸送造成風(fēng)應(yīng)力。風(fēng)應(yīng)力的量值既與冰表面特征有關(guān)，也與大氣層結(jié)狀況關(guān)系密切。風(fēng)應(yīng)力τa包括切向力τ→at和法向力τ→an。τ→an是風(fēng)作用在水面之上冰側(cè)面的力（見圖4），τ→at為冰與氣流間摩擦力。τ→at表示為平均風(fēng)速垂直梯度的函數(shù)：

鑒于表面層對數(shù)風(fēng)速分布，上式轉(zhuǎn)換為：

式中：Ca=(kln(z+z0)/z0)2，k=0.4是Karman常數(shù)，z0為冰面粗糙高度，一般經(jīng)驗確定，取值為5×10-2±2.5×10-2cm。就上述z0和10 m高的風(fēng)，阻尼系數(shù)Ca=1.64×10-3，其離散度約10%。如不直接給定z0，上述估計的切向應(yīng)力誤差不小于10%。海冰漂移速度比風(fēng)速小兩個量級，所以式（3）中的冰速可以忽略。

上述估計的冰面動力粗糙高度適合于平整冰面。通常覆蓋的冰面包括不同程度的堆積冰，從而導(dǎo)致較大的空氣阻尼系數(shù)。這時，法向應(yīng)力起重要作用。堆積冰迎風(fēng)一側(cè)的單位垂直面積受力可用類似于式（3）的公式表示。式中阻尼系數(shù)Ca須考慮風(fēng)向與冰脊軸線的交角φ的訂正（見圖4）。

圖3 北極3個測流站兩年位置變化的觀測和模擬（引自文獻(xiàn)[22]）

圖4 作用于冰面風(fēng)應(yīng)力示意圖

式中：法向應(yīng)力的系數(shù)Can尚缺野外現(xiàn)場觀測資料驗證。如果堆積冰脊的總長為l，高為h，則堆積冰脊對風(fēng)總的阻尼可表示為：

若冰脊側(cè)面和水平面的交角為ψ，那么所有堆積冰脊的底面積可表示為：

對于漂移冰塊面積A，其所受法向的風(fēng)應(yīng)力為：

式中：Ah/A是堆積冰與冰蓋的面積比，用Nh表示。系數(shù)0.5Cntanφ等于堆積冰的阻尼系數(shù)3.4×10-2，則Ca=3.4×10-2cosφ。由此可得，冰對風(fēng)總的阻尼為：

海洋的湍流渦動混合輸送造成了冰水間的應(yīng)力。假定Ekman螺旋層深度與渦動粘性有關(guān)并在邊界層中線性減小[23]，得到：

式中：Vd是解平衡動力高度和科氏力得到的海洋地轉(zhuǎn)速度，f為科氏參數(shù)。在考慮了邊界層的非線性特性同時可以更詳盡地計算水應(yīng)力。分析前蘇聯(lián)和美國北極浮冰漂移站的實驗結(jié)果表明，冰下某鄰近層摩擦切應(yīng)力的垂直變化很小，可表示為：

式中：Cw=(k/ln((z+z1)/z1))2；ρw為海水密度，V→W是z處流速，z1為冰下水的動力粗糙度，Karman常數(shù)取k=0.4。有人取得較小，如0.1，z1=2.65cm變化范圍為0.2～10 cm。一年冰的z1為0.8 cm，冬季平均冰厚為1.5～2.0 cm。當(dāng)k=0.4，z1=1cm，z=0.5～1m，Cw=0.01時，冰脊?jié)撊胨虏糠直人嫔喜糠执蟮枚?，產(chǎn)生的法向應(yīng)力相對重要很多。

大尺度冰場承受的風(fēng)應(yīng)力和水應(yīng)力：風(fēng)場和流場根據(jù)海面氣壓場和海水動力高度場，通過地轉(zhuǎn)關(guān)系得到地轉(zhuǎn)風(fēng)Vg和地轉(zhuǎn)流Vw，同時考慮風(fēng)和相對流速的扭轉(zhuǎn)偏角?和θ，則：

根據(jù)線性關(guān)系，式中的Ca′和Cw′取為常數(shù)；而非線性關(guān)系Ca′和Cw′為風(fēng)和流的函數(shù)，即：

式中：Ca和Cw均為無因次量，Mcphee[24]分別取0.001 2和0.005 5作為典型值。

（3）冰漂移數(shù)值模擬

式（1）指出決定海冰漂移的外界強迫主要由風(fēng)和流驅(qū)動，可將冰模式和海洋模式耦合并與大氣模式聯(lián)接，或提供大氣風(fēng)場，模擬大氣和海洋對冰漂移的強迫作用。

利用渤海1月海面平均風(fēng)應(yīng)力場，模擬不同冰覆蓋情況下的渤海冰漂移場，結(jié)果表示冰漂移方向偏向風(fēng)的右方10°～20°，遼東灣中部冰速最大，達(dá)到風(fēng)速的4%，其他海域約為風(fēng)速的2%～4%，當(dāng)然也包含風(fēng)生環(huán)流的影響。

將海冰動力學(xué)模式應(yīng)用于渤海并與潮流模式耦合，模擬風(fēng)與潮共同作用下海冰漂移。圖5取自遼東灣4個計算點1個潮周期的冰漂移軌跡，實線是只在潮流作用下的冰漂移軌跡，虛線是風(fēng)和潮共同作用的模擬結(jié)果，非常相似于在遼東灣破冰船跟蹤流冰塊路徑的特點，明顯地反映了渤海流冰隨橢圓主軸來回漂流，以及風(fēng)與潮共同作用的特點。Zhang等[25]利用渤海海冰模式和二維渤海潮流模式耦合研究了渤海潮流對冰漂移的影響。

3 海冰形變和冰脊、水道的形成

海冰動力學(xué)特性非常依賴于冰強度，隨冰厚特性而變化，而凍結(jié)、融化、平流和形變過程造成冰厚的變化，即冰厚改變是熱力過程和動力過程共同作用的結(jié)果。海冰動力過程主要由于形變通過脊化而形成厚冰，通過水道形成而產(chǎn)生開闊水和薄冰。

圖5 風(fēng)與潮流共同作用下渤海海冰漂移軌跡的模擬結(jié)果

為了表示地球物理尺度造脊過程，Thorndike等[27]引用冰厚分布函數(shù)，通過薄冰再分布形成不同類型厚冰。在冰厚分布方程中引用形變再分布函數(shù)ψ，表示冰非均勻運動引起水道和壓力脊形成的動力學(xué)過程，在模式中參數(shù)化地表示。形變再分布函數(shù)也是應(yīng)變率張量的函數(shù)，如果冰各向同性，則ψ僅依賴與冰速散度和切變有關(guān)的兩個應(yīng)變率張量的恒量。假如一塊冰被一些不規(guī)則裂縫分隔成若干準(zhǔn)剛體塊，平行于裂縫的位移不影響冰厚分布，但垂直于裂縫位移會引起類似于脊和水道形成的局地形變和冰厚再分布。

冰速場散度引起冰脊和水道形成是顯而易見的，圖6也明顯地表示純切變作用下仍能出現(xiàn)有些裂縫斷開，有些閉合。和純輻合、輻散形成脊和水道一樣，切變也會改變厚度分布。應(yīng)指出實際冰塊內(nèi)既有氣泡又有鹵水，以及其他雜質(zhì)。海冰復(fù)雜成分既影響它的熱力學(xué)特性，又決定它的力學(xué)強度，影響它的動力學(xué)特性。

圖6 純切變形變時水道和壓力脊形成示意圖

海冰模式中采用參數(shù)化方法表示非均勻運動引起水道和冰脊形成的動力學(xué)過程。將上述模式與大氣模式聯(lián)結(jié)，由大氣模式和邊界層模式提供冰模式海面風(fēng)場，圖7是1993年1月13日的平均堆積冰厚和冰漂移軌跡的3 d和5 d遼東灣預(yù)報結(jié)果，明顯地表示在強東北風(fēng)驅(qū)動下，沿遼東灣西北沿岸海域堆積冰的形成及5 d內(nèi)迅速發(fā)展的過程，預(yù)示了該海域堆積冰的形成。該模式還能成功地預(yù)報和模擬出水道和遼東灣底開闊水域的形成。

3 海冰流變學(xué)

海冰運動的動量平衡分析和數(shù)值試驗表明，海冰內(nèi)應(yīng)力對海冰漂移和形變非常重要。海冰流變學(xué)主要探討海冰切變和抗壓斷裂的線性和非線性性質(zhì)、切變和抗壓強度的相對量級、冰區(qū)內(nèi)冰塊隨機碰撞產(chǎn)生的無約束壓力項及冰強度與冰厚、冰塊大小特征的關(guān)系等。

（1）本構(gòu)關(guān)系

海冰的變形與海冰內(nèi)應(yīng)力關(guān)系密切，根據(jù)形變場可以確定內(nèi)應(yīng)力。在早期海冰季節(jié)變化數(shù)值模擬中采用速度相關(guān)方法近似這種作用：Parkinson等[3]采用了迭代訂正方法，即所謂“空化模型”[9]；更合理方法是根據(jù)海冰流變學(xué)，即采用海冰形變和厚度與冰內(nèi)應(yīng)力關(guān)系的本構(gòu)定律進(jìn)行科學(xué)計算。

海冰是非連續(xù)介質(zhì)，在足夠大的空間尺度上，可以被近似地描述為二維連續(xù)介質(zhì)，最普遍的本構(gòu)關(guān)系表示為[28]：

圖7 1993年1月13日渤海遼東灣堆積冰形成及冰漂移

式中：σij是二維應(yīng)力張量，εij是二維應(yīng)變率張量，δij是Kronecker算子，ζ和η是塊體和切變粘性系數(shù)，εkk=trεij=ε11+ε22，p是壓力項，p、ζ和η通常為應(yīng)變率張量的恒量的函數(shù)。其正應(yīng)力和偏應(yīng)力分別為：

這種本構(gòu)關(guān)系既包含粘性流變學(xué)又包含塑性流變學(xué)。在粘性流變學(xué)中，只有形變才能有應(yīng)力，流量的耗散是不可逆的。在塑性流變學(xué)中，非形變或能量耗散可逆形變（即彈性形變）也可維持應(yīng)力。在海冰研究中，通常采用流變學(xué)的特殊情形：

（a）線性粘性模型

對于線性粘性，應(yīng)力與應(yīng)變率成線性關(guān)系，p、ζ和η是常數(shù)。

理想固體塑性：應(yīng)力狀態(tài)不定或與應(yīng)變率無關(guān)，應(yīng)力狀態(tài)是固定的，ζ和η與應(yīng)變率恒量成非線性關(guān)系；對零應(yīng)變率狀態(tài)不能確定ζ和η，此時為剛體運動，由外力平衡決定內(nèi)應(yīng)力。

海冰的線性粘性流變學(xué)允許兩種應(yīng)力狀態(tài)，即切變粘性（Newton粘性，η=常數(shù)，ζ=p=0）和塊體粘性（ζ=常數(shù)，η=p=0）。

對于切變粘性

其應(yīng)變率張量為：

或：

應(yīng)變張量為：

式中：Si為位移分量。

對塊體粘性，介質(zhì)對膨脹和壓縮有阻尼，對切變無阻尼，因此

既包含切變粘性又包含塊體粘性的更為一般的線性粘性情況

或

式中：j是單位張量。Glen提出海冰對壓縮的阻尼大，而對拉伸的阻尼小。

忽略粘性梯度，則內(nèi)應(yīng)力F為：

即

該式已經(jīng)應(yīng)用于海冰漂移模式中，通過附加屈服判據(jù)考慮當(dāng)海冰密集度接近于1時的塑性性質(zhì)。

（b）彈-塑性模型

線性彈性本構(gòu)定律也被用于海冰動力學(xué)研究。理想的彈性過程是等溫過程，當(dāng)外力撤除時可回復(fù)到初始狀態(tài)。由于大尺度海冰變形幾乎總是不可逆，對海冰而言，線性彈性甚至非線性彈性假定通常都不合理，通常海冰表現(xiàn)為塑性性質(zhì)。Coon等[7]提出彈-塑性本構(gòu)關(guān)系，采用彈-塑性材料的屈服判據(jù)確定是否出現(xiàn)塑性性質(zhì)。

即：當(dāng)F(σ1,σ2)＜0時為彈性，當(dāng)F(σ1,σ2)=0時為塑性。

式中：F(σ1,σ2)為屈服條件。

AIDJEX研究組曾提出新屈服判據(jù)：

在彈-塑性模式中，假定彈性性質(zhì)是線性和各向同性的，則屈服函數(shù)為圓屈服函數(shù)：

于是得到：

這里

該關(guān)系常用于拉格朗日模式中，數(shù)學(xué)處理比較復(fù)雜。

圖8 海冰塑形流變學(xué)的屈服曲線

（c）粘-塑性模型

對于理想塑性本構(gòu)定律，應(yīng)力狀態(tài)位于某種固定屈服曲線之上或其內(nèi)，即圖8表示的橢圓屈服曲線。對塑性流變學(xué)來說，應(yīng)力狀態(tài)和應(yīng)變率成唯一關(guān)系，最一般的關(guān)系為給定應(yīng)力狀態(tài)，應(yīng)變率之比是垂直表面的矢量之比，即所謂的法向流定律。對于剛體運動，應(yīng)力狀態(tài)位于橢圓內(nèi)，上述兩種線性粘性定律允許大的抗壓強度（即正壓力）。然而，這種特殊塑性定律僅允許小的抗拉強度。Hibler[5]將海冰當(dāng)作非線性粘性可壓縮流體，提出粘-塑性本構(gòu)定律，這是最一般情況的流變學(xué)。冰內(nèi)壓力與冰厚、密集度相關(guān)，可表示為：

式中：η和ζ分別表示非線性切變和塊體粘性，可寫作：

因此冰內(nèi)應(yīng)力可表示為：

采用粘-塑性流變學(xué)對于強非線性海冰動力學(xué)數(shù)值模擬是很有效的方法，其原因在于它能有效避免彈-塑性流變學(xué)中彈性波所引起的問題[29]，在固定歐拉網(wǎng)格上可以更容易導(dǎo)出差分方程。當(dāng)模式采用隱式方案時，可以選取較大時間步長，這樣非常有利于模擬海冰長期季節(jié)性變化。Pritchard[14]指出，粘-塑性流變學(xué)比彈-塑性流變學(xué)包含著更普遍的本構(gòu)關(guān)系，同時粘-塑性流變學(xué)也適用于冰邊緣帶研究。

圖9表示了幾種流變學(xué)的比較情況。對于線性粘性，應(yīng)力與形變率成正比；塑性流變學(xué)無論對于大的形變率還是小的形變率，冰塊都允許相當(dāng)大的應(yīng)力，另外它還允許冰能強抗壓和切變形變，但只允許幾乎沒有形變的膨脹。圖9還表明：

圖9 幾種流變學(xué)的比較情況

（1）在彈-塑性情況下，冰抗壓應(yīng)力不變，而且對輻散應(yīng)變無阻尼；

（2）在粘-塑性情況下，具有類似性質(zhì)，但應(yīng)力狀態(tài)由形變率而不是由應(yīng)變率來決定；

（3）上述兩種塑性和線性粘性的不同在于，即使對于小的應(yīng)變率，應(yīng)力也可相當(dāng)大，并且除了很小的形變率外，應(yīng)力與形變率無關(guān)。另外，在彈-塑性情況下，無任何冰的相對運動，強應(yīng)力也能維持。

（d）海冰碰撞流變學(xué)

Shen等[30-31]和Lu[32]應(yīng)用粒狀流體理論于海冰邊緣帶的海冰運動，并假定海冰運動象“剛體微粒流”那樣，系統(tǒng)內(nèi)通過碰撞來傳遞動量消耗動能。Shen等[30]對格陵蘭海資料分析發(fā)現(xiàn)，觀測結(jié)果與根據(jù)碰撞流變學(xué)預(yù)報的冰速變化之間相關(guān)很大。

碰撞流變學(xué)假定流冰塊為均勻圓盤，由冰速隨機變化產(chǎn)生碰撞，碰撞頻率依賴于冰速大小和密集度，碰撞性質(zhì)用回復(fù)（或非彈性）系數(shù)表征，碰撞產(chǎn)生的應(yīng)力受動量變率控制，在系統(tǒng)內(nèi)轉(zhuǎn)換，因此應(yīng)力與流冰塊質(zhì)量成正比。為閉合系統(tǒng)，還需要確定隨機速度變化面，根據(jù)冰形變的能量耗散，通提供冰速變化面的機械能方程來進(jìn)行。因此假如形變?yōu)?，則變化面和應(yīng)力為0。

應(yīng)力與質(zhì)點碰撞轉(zhuǎn)換平均動量的關(guān)系為：

式中：σmn是垂直于Xm的面上，作用于Xn方向的應(yīng)力。pm是垂直于Xm的單位面積的質(zhì)點數(shù)，f為每個質(zhì)點的碰撞頻率是方向每次碰撞的平均動量輸送。為了利用粒子流理論得到冰應(yīng)力場，需要知道以下3個因子：所考慮面上的流冰塊數(shù)，每塊的碰撞頻率，每次碰撞的動量輸送。因為冰塊的大小和形狀變化多樣，所以又假定：所有冰塊具有相同形狀，即盤形；所有冰塊具有同樣的物理性質(zhì)（密度、鹽度、孔隙率和溫度等）；在同一計算網(wǎng)格上，冰塊具有相同厚度。

（e）空化流體模型

冰內(nèi)應(yīng)力依賴于表示冰速（形變）和外強迫之間關(guān)系所選取的特殊參數(shù)?？栈黧w流變學(xué)是一種特殊參數(shù)化方法，在數(shù)值模擬中考慮海冰最主要特性。它是在考慮海冰大尺度性質(zhì)基礎(chǔ)上經(jīng)過某些簡化假定而得到的。其基本假定是：將冰當(dāng)作理想的兩相介質(zhì)，一相為幾乎沒有強度的開闊水，另一相為相當(dāng)硬的冰；這種理想化介質(zhì)沒有切變和抗拉強度，當(dāng)輻散時允許開闊水的覆蓋面積自由地增加（海冰密集度減?。?，而當(dāng)壓縮時，由于“剛體”位相的原因，這種介質(zhì)具有某種阻尼。但冰位相不是完全剛體，它有可壓縮強度，當(dāng)超過此強度，海冰可以輻合和增厚，這種可壓縮強度是海冰厚度合密集度的函數(shù)。這樣的理想性質(zhì)非常便于數(shù)學(xué)處理。

（2）不同流變學(xué)的對比試驗

（a）粘-塑性流變學(xué)參數(shù)敏感性試驗

粘-塑性流變學(xué)已較廣泛地應(yīng)用于不同海冰模式，其橢圓屈服曲線各種參數(shù)的重要性已引起許多學(xué)者重視。各模式所采用參數(shù)差別也較大，且至今尚無定論[33]，應(yīng)用上述模式，采用海冰粘-塑性流變學(xué)模擬加拿大東海岸拉布拉多海域海冰漂移，并試驗橢圓屈服曲線主要控制參數(shù)冰強度參數(shù)P*對模擬冰速的影響。當(dāng)P*減小，模擬的冰速增大；反之，P*增大，冰速減小。

（b）不同流變學(xué)對冰邊緣帶（Margrinal Ice Zone，MIZ）海冰的動力學(xué)影響

Zhang[34]利用理想的MIZ海冰模式，采用浮冰塊碰撞，Hibler粘-塑性和Mohr-Coulomb 3種流變學(xué)本構(gòu)關(guān)系計算冰內(nèi)應(yīng)力。第1種海冰看作為一種微粒流體，而后兩者看作為塑性流體。給定風(fēng)場強迫，模擬結(jié)果顯示冰流變學(xué)對MIZ冰漂移有明顯影響（見圖10）。

（c）不同屈服曲線對浮冰線性運動學(xué)特性的影響

Wang[35]根據(jù)浮冰塊內(nèi)應(yīng)力場的特性分析發(fā)展了一種壓縮浮冰的屈服曲線觀測方法。分析指出屈服曲線斜率和線性運動學(xué)特性（Linear Kinematic Features，LKFs）的交叉角有關(guān)。因此根據(jù)這種交叉角測量可以反演屈服曲線。Wang等[36]進(jìn)一步利用2-level粘-塑性海冰模式研究不同模式網(wǎng)格分辨率和不同屈服條件對浮冰塊線性運動學(xué)特性（LKFs）模擬的影響。采用的網(wǎng)格分辨率包括10.5和2 km，采用3種屈服曲線（橢圓、修正的Coulomb及彎曲菱形屈服曲線，見圖11）。應(yīng)用NCEP/NCAR再分析風(fēng)場資料，并采用RADARSAT GPS海冰運動資料做比較。結(jié)果顯示應(yīng)用彎曲的菱形屈服曲線，LKFS能很好地再現(xiàn)，而應(yīng)用修正的Coulomb屈服曲線，增加網(wǎng)格分辨率僅能稍改進(jìn)LKFS模擬結(jié)果。圖12采用2 km格距，穩(wěn)定風(fēng)場3種不同屈服曲線模擬的最大切應(yīng)變率(MSSR)和用RGPS處理得到MSSR的比較。

圖10 冰邊緣帶定常風(fēng)場不同流變學(xué)本構(gòu)關(guān)系對MIZ冰運動的影響

圖11 模擬中所采用的3種屈服曲線橢圓（實線），修正的Coulomb（虛線）及彎曲菱形（點虛線）

圖12 1997年12月5—11日MSSR場的比較（網(wǎng)格距為2 km；穩(wěn)定風(fēng)場引自文獻(xiàn)[36]）

4 討論

在上述不同的海冰動力學(xué)本構(gòu)模型中，應(yīng)用最廣泛的是Hibler的粘-塑性本構(gòu)關(guān)系，它不僅應(yīng)用于極區(qū)及邊緣區(qū)域的大尺度海冰數(shù)值模擬，還應(yīng)用于像Baltic海、渤海等中尺度海冰的數(shù)值研究及數(shù)值預(yù)報。Hunke等在粘-塑海冰模式基礎(chǔ)上發(fā)展了彈-粘-塑海冰模式，可以采用較大時間步長的顯式差分計算，進(jìn)而可用于海冰動力學(xué)的并行計算，在氣候系統(tǒng)模式中得到了越來越廣泛應(yīng)用。另外，基于海冰碰撞流變學(xué)的顆粒流模型可以應(yīng)用于小尺度冰塊間的碰撞、堆積、斷裂和重疊等過程的數(shù)值模擬，其優(yōu)點在于物理意義明確，計算精度高，但因其計算量很大，從而限制了它在大、中尺度海冰動力學(xué)的數(shù)值模擬。

[1]Felzenbaum A I.The theory of steady drift of ice and the calculation of the long period mean drift in the central part of the Arctic Basin[J].Problems of the North,1961,2:13-44.

[2]Bryan K,Manabe S,Pacanowski R C.A global ocean-atmosphere climate model.Part II.The oceanic circulation[J].Journal of Physical Oceanography,1975,5(1):30-46.

[3]Parkinson C L,Washington W M.A large-scale numerical model of sea ice[J].Journal of Geophysical Research,1979,84(C1):311-337.

[4]Hibler III W D.Differential sea-ice drift.II.comparison of mesoscale strain measurements to linear drift theory predictions[J].Journal of Glaciology,1974,13(69):457-471.

[5]Hibler III W D.A dynamic thermodynamic sea ice model[J].Journal of Physical Oceanography,1979,9(4):815-846.

[6]Campbell W J.The wind driven circulation of ice and water in a polar ocean[J].Journal of Geophysical Research,1965,70(14):3279-3301.

[7]Coon M D,Maykut G A,Pritchard R S,et al.Modeling the pack ice as an elastic-plastic material[J].AIDJEX Bulletin,1974,24:1-105.

[8]Pritchard R S,Coon M D,McPhee M G.Simulation of sea ice dynamics during AIDJEX[J]. Journal of Pressure Vessel Technology,1977,99(3):491-497.

[9]Flato G M,Hibler III W D.Modeling pack ice as a cavitating fluid[J].Journal of Physical Oceanography,1992,22(6):626-651.

[10]Holland D M,Mysak L A,Manak D K,et al.Sensitivity study of a dynamic thermodynamic sea ice model[J].Journal of Geophysical Research,1993,98(C2):2561-2586.

[11]Zhang Z H.Comparisons between observed and simulated ice motion in the northern Baltic Sea[J].Geophysica,2000,36(1-2):111-126.

[12]H?kkinen S,Mellor G L.Modeling the seasonal variability of a coupled Arctic ice-ocean system[J].Journal of Geophysical Research,1992,97(C12):20285-20304.

[13]Hunke E C,Dukowicz J K.An elastic-viscous-plastic model for sea ice dynamics[J].Journal of Physical Oceanography,1997,27(9):1849-1867.

[14]Pritchard R S.An elastic-plastic constitutive law for sea ice[J].Journal ofApplied Mechanics,1975,42(2):379-384.

[15]Colony R,Pritchard R S.Integration of elastic-plastic constitutive laws[J].AIDJEX Bulletin,1975,30:55-80.

[16]Hunke E C,Zhang Y.A comparison of sea ice dynamics models at high resolution[J].Monthly Weather Review,1999,127(3):396-408.

[17]Hunkins K L.Substance eddies in the Arctic Ocean and baroclinic instability[M]//Weller G,BoWling S A.Climate of the Arctic.Fairbanks:UniversityAlaska,1975:398-406.

[18]Rothrock D A.The steady drift of an incompressible Arctic ice cover[J].AIDJEX Bulletin,1973,21:49-77.

[19]Steele M,Zhang J L,Rothrock D,et al.The force balance of sea ice in a numerical model of the Arctic Ocean[J].Journal of Geophysical Research,1997,102(C9):21061-21079.

[20]Thorndike A S,Colony R.Sea ice motion in response to geostrophic winds[J].Journal of Geophysical Research,1982,87(C8):5845-5852.

[21]Martinson D G,Wamser C.Ice drift and momentum exchange in winter Antarctic pack ice[J].Journal of Geophysical Research,1990,95(C2):1741-1755.

[22]Hibler III W D,Tucker W B III.Some results from a linearviscous model of the Arctic ice cover[J].Journal of Glaciology,1979,22(87):293-304.

[23]Reed R J,Campbell W J.Theory and observations of the drift of ice station alpha[R].Washington:Department of Meteorology,University of Washington,1960:255.

[24]McPhee M G.Ice-ocean momentum transfer for the AIDJEX ice model[J].AIDJEX Bulletin,1975,29:93-111.

[25]Zhang Z H,Wu H D.Numerical study on tides and tidal drift of sea ice in the ice-covered Bohai Sea[M]//Yu Z,Tang C L,Preller R H,et al.Sea Ice Observation and Modeling.Proceedings of Beijing’93 International Symposium on Sea Ice.Beijing:China Ocean Press,1994:34-46.

[26]Wu H D,Bai S,Zhang Z H,et al.Numerical simulation for dynamical processes of sea ice[J].Acta Oceanologica Sinica,1997,16(3):303-317.

[27]Thorndike A S,Rothrock D A,Maykut G A,et al.The thickness distribution of sea ice[J].Journal of Geophysical Research,1975,80(33):4501-4513.

[28]Glen J W.Thoughts on a viscous model for sea ice[J].AIDJEX Bulletin,1970,2:18-27.

[29]Colony R,Rothrock D A.A perspective of the time-dependent response of the AIDJEX model[M]//Pritchard R S.Sea Ice Processes and Models.Washington:University of Washington Press,1980:124-133.

[30]Shen H H,Hibler W D,Lepp?ranta M.On applying granular flow theory to a deforming broken ice field[J].Acta Mechanics,1986,63(1-4):143-160.

[31]Shen H H,Hibler III W D,Lepp?ranta M.The role of floe collisions in sea ice rheology[J].JournalofGeophysical Research,1987,92(C7):7085-7096.

[32]Lu Q M.On mesoscale modelling of the dynamics and thermodynamics of sea ice[R].Copenhagen:Technical University of Denmark,1988:50-73.

[33]吳輝碇,白珊,張占海.海冰流變學(xué)[M]//蔣伯城,張鎖春.高科技研究中的數(shù)值計算.長沙:國防科技大學(xué)出版社,1995:40-49.

[34]Zhang Z H.Effect of different rheologies on sea ice dynamics[J].Report Series in Geophysics,University of Helsinki,1998,40:39-72.

[35]Wang K G.Observing the yield curve of compacted pack ice[J].Journal of Geophysical Research,2007,112(C5):C05015,doi:10.1029/2006JC003610.

[36]Wang K G,Wang C X.Modeling linear kinematic features in pack ice[J].JournalofGeophysicalResearch,2009,114(C12):C12011,doi:10.1029/2008JC005217.

Sea ice dynamics

LIU Yu，WU Hui-ding
（National Marine Environment Forecasting Center,Key Laboratory of Research on Marine Hazards Forecasting,Beijing 100081 China）

In this paper,the development of sea ice dynamics is briefly reviewed.The dynamic processes of sea ice drift,deformation and accumulation are introduced.The applications of different sea ice constitutive equations in sea ice rheology are compared and analyzed.

sea ice；dynamics；rheology；constitutive equations

P731.15

A

1003-0239（2017）05-0099-12

10.11737/j.issn.1003-0239.2017.05.011

2017-09-22；

2017-10-09。

國家重點研發(fā)計劃（2016YFC1401502）；國家自然科學(xué)基金（41676189；41506109）。

劉煜（1973-），男，副研究員，碩士，主要從事海冰預(yù)報和研究工作。E-mail：liuyu@nmefc.gov.cn