渤海海冰演化的三維數(shù)值模型

張慶河，張 娜,

(1. 天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300072；

2. 天津城市建設(shè)學(xué)院土木工程學(xué)院，天津 300384)

渤海海冰演化的三維數(shù)值模型

張慶河1，張 娜1,2

(1. 天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300072；

2. 天津城市建設(shè)學(xué)院土木工程學(xué)院，天津 300384)

基于非結(jié)構(gòu)化網(wǎng)格、有限體積三維海洋模型FVCOM建立了海洋-海冰耦合的渤海海冰數(shù)值模擬系統(tǒng)．模型中考慮了海冰動(dòng)力過(guò)程和海冰熱力過(guò)程，討論了適合渤海的冰厚度分類(lèi)規(guī)則、臨界冰厚、短波輻射在冰內(nèi)的透射率以及冰的熱傳輸系數(shù)的合理取值．利用美國(guó)國(guó)家環(huán)境預(yù)報(bào)中心(NCEP)一小時(shí)一次的再分析氣象數(shù)據(jù)作為驅(qū)動(dòng)對(duì)2003—2004年冬季渤海海冰過(guò)程進(jìn)行了模擬．模擬的渤海溫鹽計(jì)算值與實(shí)測(cè)值吻合較好，模擬的渤海冰蓋分布與同期衛(wèi)星遙感結(jié)果接近，模擬的遼東灣初冰日、終冰日的海冰面積、海冰厚度及其最大值和最大值對(duì)應(yīng)的時(shí)刻與實(shí)測(cè)值基本一致．模擬與實(shí)測(cè)結(jié)果對(duì)比表明，建立的三維海冰數(shù)值模型能比較合理地反映渤海動(dòng)力過(guò)程、氣象條件影響下的溫鹽場(chǎng)過(guò)程以及渤海海冰的生消演化過(guò)程，具有較高的計(jì)算精度和計(jì)算效率，適合用于渤海海冰中長(zhǎng)期數(shù)值模擬.

渤海海冰；有限體積海洋模型(FVCOM)；非結(jié)構(gòu)化網(wǎng)格；數(shù)值模擬

海冰是較高緯度地區(qū)海域出現(xiàn)的自然現(xiàn)象，我國(guó)的渤海和黃海北部海域，冬季常受西伯利亞寒潮的影響，氣溫降低，每年都會(huì)形成海冰，且容易造成災(zāi)害[1]．1969年2—3月間渤海曾出現(xiàn)歷史上罕見(jiàn)的冰情；2010年冬，渤海及黃海北部又發(fā)生了近30年來(lái)最嚴(yán)重冰情，海水結(jié)冰導(dǎo)致航道阻塞，船只、海上設(shè)施和海岸工程遭到破壞，港口碼頭封凍，水產(chǎn)養(yǎng)殖受損[2]．鑒于海冰災(zāi)害給人類(lèi)在海上的生命財(cái)產(chǎn)安全以及生產(chǎn)活動(dòng)造成的嚴(yán)重影響和危害，深入研究海冰生消演變規(guī)律、建立和完善海冰預(yù)警預(yù)報(bào)機(jī)制，對(duì)于海洋防災(zāi)減災(zāi)具有十分重要的意義．

海冰數(shù)值模擬是研究海冰生消變化規(guī)律的一種重要方法，在與氣象預(yù)報(bào)結(jié)合的基礎(chǔ)上可以進(jìn)一步進(jìn)行海冰預(yù)報(bào)．我國(guó)早期的海冰數(shù)學(xué)模型工作主要集中于海冰動(dòng)力學(xué)的研究[3-4]，只適宜做非常短期的海冰預(yù)報(bào)，而海冰的生長(zhǎng)與消融是同時(shí)受動(dòng)力條件和熱力條件影響的，特別是在進(jìn)行較長(zhǎng)期預(yù)報(bào)時(shí)必須綜合考慮熱力學(xué)影響．近年來(lái)隨著冰-海耦合技術(shù)的不斷提升，海冰數(shù)值模擬的時(shí)間尺度也從短期發(fā)展到中長(zhǎng)期．唐茂寧等[5]以MM5數(shù)值產(chǎn)品作為氣象輸入條件，利用結(jié)構(gòu)化網(wǎng)格的POM模式建立海冰-海洋耦合數(shù)值模式同時(shí)模擬動(dòng)力學(xué)和熱力學(xué)過(guò)程，對(duì)2003—2004年冬季渤海海冰過(guò)程進(jìn)行了模擬；劉煜等[6]基于冰-海耦合模型，模擬了從1997—1998年冬季至2008—2009年冬季渤海海冰的季節(jié)性變化．在國(guó)外，Herman等[7]基于結(jié)構(gòu)化網(wǎng)格的POM模型，同時(shí)考慮海冰熱力學(xué)和動(dòng)力學(xué)，模擬了波羅的海海冰的季節(jié)性變化；Gao等[8]基于有限體積法的FVCOM三維海洋模型開(kāi)發(fā)了非結(jié)構(gòu)化網(wǎng)格海冰-海洋耦合數(shù)值模型，該模型同時(shí)考慮海冰動(dòng)力學(xué)和熱力學(xué)過(guò)程，成功應(yīng)用于北極海域海冰的數(shù)值模擬．為了深入研究具有復(fù)雜幾何形狀的渤海的海冰演化規(guī)律，發(fā)展渤海海冰預(yù)報(bào)的數(shù)值模型，筆者將基于FVCOM模型建立渤海海冰演化的海冰-海洋耦合三維數(shù)值模式，并收集資料驗(yàn)證了模型的合理性．

1 數(shù)值模型

海冰的數(shù)值模擬涉及到三維水動(dòng)力學(xué)、海溫、鹽度以及海冰熱力學(xué)、海冰動(dòng)力學(xué)控制方程的求解及其耦合，圖1給出了海洋-海冰耦合數(shù)值模型的組成情況和各子模型間的變量傳遞．下面對(duì)涉及的模型進(jìn)行簡(jiǎn)要介紹．

圖1 海冰數(shù)值模型的組成Fig.1 Composition of the sea ice numerical model

1.1 FVCOM三維海洋模型

FVCOM是描述海洋運(yùn)動(dòng)的三維模型，其基于直角坐標(biāo)系的原始控制方程如下．

動(dòng)量方程

連續(xù)性方程

海溫方程

鹽度方程

密度方程

式中：x、y和z分別代表坐標(biāo)系的東向、北向和垂向坐標(biāo)軸；u、v和w分別為x、y和z軸的速度分量；T為海溫；S為鹽度；ρ為密度；p為壓強(qiáng)；f為科氏力系數(shù)；g為重力加速度；Km為垂向渦黏系數(shù)；Kh為垂向渦黏熱擴(kuò)散系數(shù)；Fu(Fv)、FT和FS分別代表水平動(dòng)量項(xiàng)、水平熱擴(kuò)散項(xiàng)和水平鹽擴(kuò)散項(xiàng)．控制方程(1)～(7)可實(shí)現(xiàn)垂向上的σ坐標(biāo)變換，即

式中σ在[-1,0]之間變化．在海表面z=b，σ=0，在海底zH=-，1σ=-；D為總水深，DHb=+，H和b分別為相對(duì)于0z=的底部水深和自由表面高度.

1.2 海冰模型的基本方程

海冰模型的基本方程主要包括海冰熱力學(xué)方程、海冰動(dòng)力學(xué)方程、海冰傳輸方程和海冰本構(gòu)方程．

1.2.1 海冰熱力學(xué)的控制方程

1) 來(lái)自大氣到達(dá)海冰上表面的凈熱通量

式中：NETF為凈熱通量；SHF為感熱通量；LHF為潛熱通量；LWF↑為冰面向大氣的長(zhǎng)波輻射；LWF↓為大氣逆輻射；SWF↓為來(lái)自大氣到達(dá)海冰表面的總太陽(yáng)輻射；α為冰面反照率；0I為短波輻射在冰內(nèi)的透射率．

2) 海水結(jié)冰的判斷條件

當(dāng)海水溫度低至臨界凍冰溫度f(wàn)TSμ=-時(shí)開(kāi)始形成海冰，μ為海水的凍結(jié)溫度與鹽度比例，μ= 0.054．用Ffrzmlt判斷凍結(jié)或融化，其表達(dá)式為

式中：wT為海表面溫度；wρ為海水的密度；pc為海水的比定壓熱容；mixh為海洋混合層深度；TD為海冰熱力學(xué)的計(jì)算時(shí)間步長(zhǎng)．當(dāng)frzmlt0F＞時(shí)海水開(kāi)始結(jié)冰，此時(shí)，冰底面(指海冰下表面)溫度為臨界凍冰溫度f(wàn)T.

3)海冰融化的判斷條件

冰底面的增長(zhǎng)和消融是由冰底面的傳導(dǎo)熱通量與海洋熱通量之間的凈熱通量控制的，用BOTF表示，其表達(dá)式為

式中：hc為熱傳輸系數(shù)；*u為海冰和海水之間的相對(duì)速度．

1.2.2 海冰動(dòng)力學(xué)的控制方程

直角坐標(biāo)系下，描述海冰運(yùn)動(dòng)的動(dòng)量方程為

式中：m為單位面積的冰和雪的混合質(zhì)量；us和vs為冰速在x和y方向上的分量；σij為內(nèi)部應(yīng)力張量，下標(biāo)1、2和i、j分別代表x和y方向上的分量；(τax,τay)和(τwx,τwy)分別代表x和y方向上的海表面風(fēng)應(yīng)力和海表面水應(yīng)力．表面風(fēng)應(yīng)力和表面水應(yīng)力的矢量形式可寫(xiě)為

式中：c為值從0～1的海冰密集度；su為表面速度矢量；u、C、ρ分別代表速度矢量、拖曳力系數(shù)和密度；下標(biāo)a和w分別代表空氣和水．

海冰本構(gòu)關(guān)系采用Hunke的彈黏塑性(EVP)模型[9]，其方程形式為

式中：E為由彈性模量定義的彈性參數(shù)；下標(biāo)kk表示Einstein求和；η為剪切黏度；ζ為體積黏度；p為冰強(qiáng)度；δij為Kronecker函數(shù)，i=j，δij=1，i≠j，δij=0；ε˙ij為冰應(yīng)變率．

1.2.3 海冰傳輸方程

海冰熱力學(xué)和動(dòng)力學(xué)模型求解的基本問(wèn)題就是計(jì)算冰厚度隨空間和時(shí)間的演化分布，這2個(gè)過(guò)程可通過(guò)海冰傳輸方程耦合起來(lái)，其標(biāo)準(zhǔn)對(duì)流方程形式為

1.3 模型的數(shù)值離散

FVCOM水動(dòng)力、溫鹽方程以及海冰模型控制方程均采用有限體積法進(jìn)行數(shù)值離散，潮流與溫鹽模型通過(guò)內(nèi)外模分離的數(shù)值離散方法，二維外模數(shù)值格式為基于非結(jié)構(gòu)化三角形網(wǎng)格的有限體積法，將連續(xù)性方程、動(dòng)量方程在控制體區(qū)域積分后，通過(guò)改進(jìn)的四階龍格庫(kù)塔方法求解，計(jì)算垂向平均的潮位和潮流．三維內(nèi)模的動(dòng)量方程和溫鹽方程采用簡(jiǎn)單的顯式與隱式相結(jié)合的差分格式，對(duì)流項(xiàng)采用二階精度的龍格庫(kù)塔時(shí)間推進(jìn)格式，垂向擴(kuò)散使用隱式格式求解，計(jì)算三維潮流和三維溫鹽，海冰模型的數(shù)值離散采用控制方程的積分形式．關(guān)于詳細(xì)的離散格式參見(jiàn)文獻(xiàn)[10]．

2 渤海模型參數(shù)討論與選取

2.1 水動(dòng)力、溫鹽模型的參數(shù)

渤海水動(dòng)力、溫鹽模型的主要計(jì)算參數(shù)及取值見(jiàn)表1．

表1 模型主要計(jì)算參數(shù)取值Tab.1 Values of main parameters in the model

2.2 海冰模型的參數(shù)

2.2.1 冰厚分類(lèi)規(guī)則

海冰傳輸方程通過(guò)離散每個(gè)網(wǎng)格內(nèi)的冰厚種類(lèi)求解，不同冰厚種類(lèi)的冰增長(zhǎng)率是不一樣的，因此冰厚的分類(lèi)規(guī)則對(duì)海冰模擬有重要影響．20世紀(jì)70年代，Thorndike等[11]將極地海冰視為包含無(wú)冰海水和厚冰脊的混合物，并首次提出冰厚度的分布函數(shù)．21世紀(jì)初，Lipscomb[12]給出了適合于北冰洋等極寒海域的冰厚分類(lèi)規(guī)則，通常包括無(wú)冰海水、一年冰、多年厚冰以及厚冰脊．渤海海冰的特點(diǎn)是冰期短(僅冬季有)、冰層薄(常冰年冰厚一般為0.20～0.40,m，最大冰厚0.60,m左右)，災(zāi)情最重的年份其最大堆積冰高度不超過(guò)3,m[13]．如果將北極海冰冰厚的分類(lèi)規(guī)則應(yīng)用于渤海，則會(huì)使薄冰與厚冰合并，不能準(zhǔn)確模擬．敏感性測(cè)試結(jié)果表明，冰厚度的分類(lèi)比較細(xì)時(shí)，模擬結(jié)果較合理，因此本文采用適合渤海的冰厚分類(lèi)規(guī)則為：0、0.10,m、0.30,m、0.60,m、1.00,m、1.50,m、2.10,m和2.80,m．敏感性計(jì)算表明，對(duì)于比該規(guī)則更細(xì)的劃分，模擬結(jié)果并沒(méi)有改善．

2.2.2 臨界冰厚

臨界冰厚的設(shè)置是為了保證熱力學(xué)計(jì)算的穩(wěn)定性，是影響冰模型計(jì)算結(jié)果的一個(gè)重要參數(shù)．在熱力學(xué)模式中，Hibler[14]提出在一個(gè)計(jì)算網(wǎng)格內(nèi)存在2類(lèi)冰厚分布，將薄冰和無(wú)冰海水歸為一類(lèi)，厚冰歸為一類(lèi)，針對(duì)北極海冰模擬給出了區(qū)分2類(lèi)冰的臨界冰厚為0.5,m．本文根據(jù)渤海海冰總體偏薄的特點(diǎn)，選取臨界冰厚為0.01,m．

2.2.3 熱傳輸系數(shù)

熱傳輸系數(shù)hc和垂向渦黏熱擴(kuò)散系數(shù)hK對(duì)于模式的計(jì)算結(jié)果影響較大，尤其是熱傳輸系數(shù)直接影響著海冰融化的速度．McPhee[15]在北極地區(qū)測(cè)得的適合該海區(qū)的海冰熱傳輸系數(shù)為0.006，Josberger[16]在白令海和弗拉姆海峽測(cè)量得到的為2×10-4～8×10-4；Haapala等[17]在對(duì)波羅的海海冰的研究中取ch為2.8×10-4；Shirasawa等[18]在加拿大Resolute Passage海的實(shí)驗(yàn)結(jié)果ch為2.3×10-3，本文選取ch為2.2×10-5[19]．敏感性測(cè)試研究結(jié)果表明，越小則海冰融化速度越慢．垂向渦黏熱擴(kuò)散系數(shù)hK的選取關(guān)系到模式的計(jì)算穩(wěn)定性，多次數(shù)值試驗(yàn)結(jié)果表明，選取1.0×10-4較合理．

2.2.4 短波輻射透射率

進(jìn)入冰內(nèi)的太陽(yáng)短波輻射能縮短融冰期、阻止海冰厚度增加，是影響冰模型計(jì)算結(jié)果的一個(gè)重要參數(shù)．海冰內(nèi)部可能含有大量氣泡、鹵水泡和泥沙等雜質(zhì)，因此太陽(yáng)輻射在海冰內(nèi)的透射率并不是一個(gè)恒定值．根據(jù)Beer定律[20]，本文的透射系數(shù)為

式中：I(z)為冰表面以下深度為z處的太陽(yáng)短波輻射；Fswabs為到達(dá)冰面除去反射部分的太陽(yáng)輻射，F(xiàn)swabs=(1-α)FSW↓；i0對(duì)于近紅外輻射取值為0，對(duì)于可見(jiàn)光輻射取值為0.7；κi表示冰內(nèi)太陽(yáng)輻射的衰減系數(shù)，對(duì)于可見(jiàn)光波長(zhǎng)設(shè)定為1.4,m-1．

3 渤海海冰模擬結(jié)果驗(yàn)證與分析

3.1 模型設(shè)置

3.1.1 空間網(wǎng)格與時(shí)間步長(zhǎng)

采用上面建立的模型對(duì)渤海海冰進(jìn)行了從2003年11月3日0：00至2004年2月29日23：00長(zhǎng)達(dá)4個(gè)月的數(shù)值模擬．海冰熱力學(xué)模擬的時(shí)間步長(zhǎng)與內(nèi)模式的時(shí)間步長(zhǎng)一致，為20,s．為了保證EVP本構(gòu)方程求解彈性波過(guò)程的穩(wěn)定性，海冰動(dòng)力學(xué)模擬的時(shí)間步長(zhǎng)設(shè)為熱力學(xué)時(shí)間步長(zhǎng)的1/2，為10,s；外模式的時(shí)間步長(zhǎng)為2,s．水平方向使用非結(jié)構(gòu)化三角形網(wǎng)格進(jìn)行劃分，并在渤海灣、遼東灣和萊州灣近岸將網(wǎng)格局部加密，水平方向的最小網(wǎng)格步長(zhǎng)為0.01°，最大網(wǎng)格步長(zhǎng)為0.055°．渤海局部(遼東灣)及渤海測(cè)站布置的網(wǎng)格劃分如圖2所示．垂向采用σ坐標(biāo)共劃分了10層，并在水面附近加密．

3.1.2 近岸動(dòng)邊界處理

根據(jù)實(shí)際的沿岸冰分布形態(tài)，在FVCOM三維干濕動(dòng)邊界處理技術(shù)的基礎(chǔ)上，設(shè)定潮間帶海冰不隨漲落潮移動(dòng)，只要海溫低于臨界凍冰溫度，則網(wǎng)格點(diǎn)干出時(shí)，依然有冰覆蓋，冰密集度與冰厚為網(wǎng)格點(diǎn)干出前一刻的數(shù)值，干出點(diǎn)網(wǎng)格的冰速為0．上述做法可以保證海冰求解的穩(wěn)定性，也接近實(shí)際情況．

3.1.3 模型初、邊值及外部驅(qū)動(dòng)條件

模型初值、邊界條件及外部驅(qū)動(dòng)條件如下．

(1) 模型計(jì)算的初始三維海溫場(chǎng)由美國(guó)國(guó)家環(huán)境預(yù)報(bào)中心(NCEP)發(fā)布的氣候預(yù)測(cè)系統(tǒng)(CFSR)一小時(shí)再分析海溫?cái)?shù)據(jù)插值而成，水平方向空間分辨率為0.3°[21]．初始鹽度場(chǎng)參考1998年10月實(shí)測(cè)渤海近岸秋季表層、底層鹽度分布[22]，結(jié)合NCEP一小時(shí)鹽度數(shù)據(jù)插值而成．

(2) 開(kāi)邊界潮位由潮汐預(yù)報(bào)軟件Chinatide提供[23]．

(3) 模型輸入條件是指風(fēng)、氣溫、相對(duì)濕度、云量、短波輻射、長(zhǎng)波輻射和海表面氣壓等全場(chǎng)逐時(shí)驅(qū)動(dòng)參數(shù)，其中氣象參數(shù)用來(lái)計(jì)算感熱通量、潛熱通量和凈熱通量等各種熱通量．外部驅(qū)動(dòng)場(chǎng)來(lái)自NCEP一小時(shí)再分析數(shù)據(jù)，該數(shù)據(jù)的水平方向空間分辨率為0.3°，時(shí)間間隔為1,h．

圖2 渤海局部(遼東灣)網(wǎng)格劃分及渤海測(cè)站布置Fig.2 Computational meshes for part of the Bohai Sea(the Liaodong Bay)and the observation stations in the Bohai Sea

3.2 模擬與實(shí)測(cè)結(jié)果的比較

為了驗(yàn)證模型的合理性，在渤海海域收集了計(jì)算時(shí)間內(nèi)的潮汐預(yù)報(bào)表潮位資料和各種實(shí)測(cè)資料，其中包括鲅魚(yú)圈和塘沽潮汐預(yù)報(bào)表潮位資料，葫蘆島和龍口實(shí)測(cè)海溫資料，大連鹽度實(shí)測(cè)資料，MODIS衛(wèi)星遙感資料以及JZ20-2平臺(tái)冰厚資料，圖2中標(biāo)出了資料收集點(diǎn)的位置．以下給出模擬結(jié)果與所收集資料的比較情況．

3.2.1 潮位模擬與潮汐表結(jié)果比較

水動(dòng)力模型計(jì)算了從2003年11月1日0：00至2004年2月29日23：00共4個(gè)月的渤海潮波，鲅魚(yú)圈和塘沽2個(gè)站點(diǎn)的潮位計(jì)算值與天文潮潮汐表的對(duì)比結(jié)果見(jiàn)圖3和圖4，圖中時(shí)間從2003年11月20日0：00至2003年11月30日0：00．由圖3和圖4可知，2個(gè)驗(yàn)潮站的計(jì)算潮位與潮汐表潮位吻合較好，水動(dòng)力結(jié)果的合理性為模擬渤海溫鹽及海冰奠定了基礎(chǔ)．

圖3 鲅魚(yú)圈計(jì)算水位與潮汐表對(duì)比Fig.3 Comparison between the calculated values and the water elevation from the tide table in Bayuquan

3.2.2 溫鹽模擬與實(shí)測(cè)結(jié)果比較和分析

圖5所示為葫蘆島、龍口海溫模擬值與實(shí)測(cè)值的對(duì)比情況，圖6所示為大連鹽度模擬值與實(shí)測(cè)值的對(duì)比情況．圖5和圖6中開(kāi)始時(shí)間為2003年11月9日，結(jié)束時(shí)間為2004年2月28日；受海面結(jié)冰影響，實(shí)測(cè)數(shù)據(jù)只到2003年11月30日．由圖5和圖6可知，驗(yàn)證點(diǎn)溫鹽模擬結(jié)果與實(shí)測(cè)結(jié)果有很好的一致性.

圖4 塘沽計(jì)算水位與潮汐表對(duì)比Fig.4 Comparison between the calculated values and the water elevation from the tide table in Tanggu

圖5 葫蘆島和龍口海表面溫度計(jì)算值與實(shí)測(cè)值Fig.5 Comparison between the calculated and measured sea surface temperatures for Huludao and Longkou

圖6 大連海表面鹽度計(jì)算值與實(shí)測(cè)值Fig.6 Comparison between the calculated and measured sea surface salinities for Dalian

為了直觀觀察渤海在海冰發(fā)生期的海溫平面分布以及沿水深的垂向分布，圖7給出了2004年1月25日模擬海表面溫度場(chǎng)．圖8以JZ20-2平臺(tái)所在位置為例給出不同日期垂向海溫分布模擬結(jié)果．從渤海溫、鹽的計(jì)算結(jié)果可知，由于冬季降水很少，徑流較弱，風(fēng)速較大引起海表面蒸發(fā)加速，致使鹽度略有上升，但總體變化不大．冬季海溫主要受冷空氣的影響呈現(xiàn)溫度不斷降低的趨勢(shì)，但海溫的日變化較氣溫平穩(wěn)，計(jì)算結(jié)果為冬季表層海溫南高北低(如圖5所示)，沿岸低外海高(如圖7所示)．由于冬季大風(fēng)使得表層與底層混合均勻，底層海溫分布與表層一致(如圖8所示)，這與以往的研究成果相吻合[24]．可以認(rèn)為，溫、鹽的計(jì)算值與實(shí)測(cè)值吻合較好，為海冰模擬奠定了基礎(chǔ)．

圖7 2004年1月25日渤海海表面溫度場(chǎng)分布Fig.7 Sea surface temperature distribution of the Bohai Sea on Jan. 25th，2004

圖8 JZ20-2平臺(tái)不同日期垂向海溫分布Fig.8 Vertical distribution of sea water temperature on different dates at JZ20-2 platform position

圖9 2003年12月7日渤海海冰數(shù)值模擬冰蓋分布與MODIS遙感影像冰蓋分布對(duì)比Fig.9 Comparison between the calculated ice cover and remote sensing image of MODIS in the Bohai Sea on Dec. 7th，2003

3.2.3 海冰模擬與衛(wèi)星遙感結(jié)果比較

渤海海冰冰蓋模擬的結(jié)果與MODIS衛(wèi)星遙感的同期冰蓋分布對(duì)比如圖9～圖12所示．根據(jù)MODIS遙感影像監(jiān)測(cè)，2003—2004年冬季為常冰偏輕年，初冰日為2003年12月6日，終冰日2004年2月24日，冰期76天．模式計(jì)算的初冰日與遙感監(jiān)測(cè)的日期吻合，模式計(jì)算的終冰日為2004年2月23日，只比遙感監(jiān)測(cè)的日期提前1,d．

圖13以遼東灣為例給出了海冰面積變化過(guò)程的模擬結(jié)果與衛(wèi)星遙感結(jié)果對(duì)比．從MODIS衛(wèi)星遙感結(jié)果看，遼東灣最大海冰面積發(fā)生于2004年1月25日，約為13,814,km2，整個(gè)渤海最大海冰面積發(fā)生在2004年1月25日，約為20,402,km2．模式計(jì)算的海冰面積最大時(shí)刻與遙感監(jiān)測(cè)的時(shí)間一致，但最大面積較監(jiān)測(cè)結(jié)果略偏小，遼東灣最大海冰面積約為12,432,km2，渤海最大海冰面積約為18,361,km2．

圖10 2004年1月25日渤海海冰數(shù)值模擬冰蓋分布與MODIS遙感影像冰蓋分布對(duì)比Fig.10 Comparison between the calculated ice cover and remote sensing image of MODIS in the Bohai Sea on Jan. 25th，2004

圖11 2004年2月2日渤海海冰數(shù)值模擬冰蓋分布與MODIS遙感影像冰蓋分布對(duì)比Fig.11 Comparison between the calculated ice cover and remote sensing image of MODIS in the Bohai Sea on Feb. 2nd，2004

圖12 2004年2月14日渤海海冰數(shù)值模擬冰蓋分布與MODIS遙感影像冰蓋分布對(duì)比Fig.12 Comparison between the calculated ice cover and remote sensing image of MODIS in the Bohai Sea on Feb. 14th，2004

圖13 遼東灣海冰數(shù)值模擬面積與MODIS 遙感影像面積對(duì)比Fig.13 Comparison between the calculated sea ice area and the remote sensing data of MODIS in the Liaodong Bay

就渤海灣和萊州灣而言，遙感監(jiān)測(cè)的渤海灣和萊州灣于2004年1月22日開(kāi)始有冰蓋分布，2004年2月2日渤海灣和萊州灣內(nèi)海冰全部融化；模式計(jì)算的渤海灣于2004年1月20日開(kāi)始結(jié)冰，1月30日冰期結(jié)束；萊州灣于2004年1月21日開(kāi)始結(jié)冰，1月31日冰期結(jié)束，模式計(jì)算的渤海灣與萊州灣冰期僅有短短的10,d，與遙感監(jiān)測(cè)的日期吻合較好．

圖14顯示了JZ20-2平臺(tái)附近海域的冰厚模擬值與實(shí)測(cè)值[3]對(duì)比情況，圖中初始時(shí)刻為2003年12月31日0：00．由圖14可知，JZ20-2平臺(tái)的模擬值與實(shí)測(cè)值的初冰日、終冰日吻合較好，模擬的冰厚最大值以及最大值對(duì)應(yīng)的時(shí)刻、冰厚變化趨勢(shì)與實(shí)測(cè)較為一致．

圖14 JZ20-2平臺(tái)冰厚模擬值與實(shí)測(cè)值對(duì)比Fig.14 Comparison between the calculated and the measured sea ice thicknesses at JZ20-2 platform

4 結(jié) 論

基于FVCOM三維海洋模型建立了海洋-海冰耦合的渤海海冰數(shù)值模擬系統(tǒng)，討論了適合渤海的冰厚分類(lèi)規(guī)則、臨界冰厚和短波輻射在冰內(nèi)的透射率以及冰的熱傳輸系數(shù)的合理取值，利用NCEP一小時(shí)再分析氣象數(shù)據(jù)作為驅(qū)動(dòng)對(duì)2003—2004年冬季渤海海冰過(guò)程進(jìn)行了模擬，得出如下結(jié)論：

(1) 葫蘆島與龍口2個(gè)近岸實(shí)測(cè)站點(diǎn)的海表面溫度、大連的海表面鹽度的驗(yàn)證結(jié)果表明計(jì)算值與實(shí)測(cè)值吻合較好．由2004年1月25日渤海海表面溫度分布及JZ20-2平臺(tái)位置不同時(shí)刻垂向海溫分布結(jié)果可知，冬季表層海溫呈南高北低、沿岸低外海高特征．由于冬季大風(fēng)使得表層與底層混合均勻，海溫沿垂向分布非常均勻，底層海溫與表層幾乎相同．

(2) 海冰模擬結(jié)果與MODIS衛(wèi)星遙感的同期冰蓋分布、海冰面積及JZ20-2平臺(tái)附近海域的冰厚實(shí)測(cè)值對(duì)比，初冰日、終冰日吻合較好，模擬的冰面積和冰厚最大值及最大值對(duì)應(yīng)的時(shí)刻，以及冰面積、冰厚變化趨勢(shì)均與實(shí)測(cè)較為一致．

(3) 模擬和實(shí)測(cè)結(jié)果比較表明，本文建立的渤海三維海冰數(shù)值模型能比較合理地反映渤海動(dòng)力過(guò)程、氣象條件影響下的溫鹽場(chǎng)過(guò)程，再現(xiàn)了渤海海冰的發(fā)生、發(fā)展和消融過(guò)程．基于FVCOM三維海洋模型，綜合考慮海冰熱力學(xué)和動(dòng)力學(xué)過(guò)程的海洋-海冰耦合模式，具有與氣象預(yù)報(bào)相結(jié)合進(jìn)行海冰預(yù)測(cè)的潛力，非結(jié)構(gòu)化網(wǎng)格可實(shí)現(xiàn)對(duì)復(fù)雜岸線模擬及對(duì)關(guān)心區(qū)域局部加密，未來(lái)為精細(xì)化模擬特別是港口海冰發(fā)展模擬打下了基礎(chǔ)，值得深入研究和推廣．

［1］ 王相玉，袁本坤，商 杰，等. 渤黃海海冰災(zāi)害與防御對(duì)策[J]. 海岸工程，2011，30(4)：46-55.

Wang Xiangyu，Yuan Benkun，Shang Jie，et al. Disaster of sea ice and defense countermeasures in the Bohai Sea and the Yellow Sea[J]. Coastal Engineering，2011，30(4)：46-55(in Chinese).

［2］ 孫 劭，蘇 潔，史培軍. 2010年渤海海冰災(zāi)害特征分析[J]. 自然災(zāi)害學(xué)報(bào)，2011，20(6)：87-93.

Sun Shao，Su Jie，Shi Peijun. Features of sea ice disaster in the Bohai Sea in 2010[J]. Journal of Natural Disasters，2011，20(6)：87-93(in Chinese).

［3］ 季順迎，岳前進(jìn)，趙 凱. 渤海海冰動(dòng)力學(xué)的質(zhì)點(diǎn)網(wǎng)格法數(shù)值模擬[J]. 水動(dòng)力學(xué)研究與進(jìn)展，2003，18(6)：748-760.

Ji Shunying，Yue Qianjin，Zhao Kai. Numerical simulation of sea ice dynamics with particle-in-cell approach in the Bohai Sea [J]. Journal of Hydrodynamics，2003，18(6)：748-760(in Chinese).

［4］ 王瑞學(xué)，季順迎，岳前進(jìn). 海冰動(dòng)力學(xué)數(shù)值模擬中改進(jìn)的PIC方法[J]. 計(jì)算力學(xué)學(xué)報(bào)，2006，23(4)：440-446.

Wang Ruixue，Ji Shunying，Yue Qianjin. Modified PIC method for numerical simulation of sea ice dynamics[J]. Chinese Journal of Computational Mechanics，2006，23(4)：440-446(in Chinese).

［5］ 唐茂寧，劉欽政，劉 煜，等. 渤海海冰季節(jié)演變的數(shù)值模擬[J]. 海洋預(yù)報(bào)，2010，27(2)：48-52.

Tang Maoning，Liu Qinzheng，Liu Yu，et al. The numerical simulation of sea ice in the Bohai Sea[J]. Marine Forecasts，2010，27(2)：48-52(in Chinese).

［6］ Liu Yu，Liu Qinzheng，Su Jie，et al. Seasonal simulations of a coupled ice-ocean model in the Bohai Sea and North Yellow Sea since the winter of 1997/1998[C]//Proceedings of the Twenty-First(2011) International Offshore and Polar Engineering Conference. Maui，Hawaii，USA，2011：19-24.

［7］ Herman A，Jedrasik J，Kowalewski M. Numerical modelling of thermodynamics and dynamics of sea ice in the Baltic Sea[J]. Ocean Science，2011，7：257-276.

［8］ Gao Guoping，Chen Changsheng，Qi Jianhua，et al. An unstructured-grid，finite-volume sea ice model：Development，validation，and application[J]. J Geophys Res，2011，116(9)：2010JC006688.

［9］ Hunke E C. Viscous-plastic sea ice dynamics with the EVP model：Linearization issues[J]. J Comput Phys，2001，170：18-38.

［10］ Chen C S，Cowles G，Beardsley R C. An Unstructured Grid，F(xiàn)inite-Volume Coastal Ocean Model：FVCOM User Manual[R]. USA：School of Marine Science and Technology，University of Massachusetts-Dartmouth，2006.

［11］ Thorndike A S，Rothrock D A，Maykut G A，et al. The thickness distribution of sea ice[J]. J Geophys Res，1975，80(33)：4501-4513.

［12］ Lipscomb W H. Remapping the thickness distribution in sea ice models[J]. J Geophys Res，2001，106(C7)：13989-14000.

［13］ 白 珊，劉欽政，李 海，等. 渤海的海冰[J]. 海洋預(yù)報(bào)，1999，16(3)：1-9.

Bai Shan，Liu Qinzheng，Li Hai，et al. Sea ice in the Bohai Sea of China[J]. Marine Forcasts，1999，16(3)：1-9(in Chinese).

［14］ Hibler W D. A dynamic thermodynamic sea icemodel[J]. J Phys Oceanogr，1979，9：815-846.

［15］ McPhee M G. Turbulent heat flux in the upper ocean under sea ice [J]. J Geophys Res，1992，97(C4)：5365-5379.

［16］ Josberger E G. Bottom ablation and heat transfer coefficients from the 1983 marginal ice zone experiment [J]. J Geophys Res，1987，92(C7)：7012-7016.

［17］ Haapala J，Lepparanta M. Simulating the Baltic Sea ice season with a coupled ice：Ocean model[J]. Tellus A，1996，48：622-643.

［18］ Shirasawa K，Ingram G R. Currents and turbulent fluxes under the first-year sea ice in resolute passage，Northwest Territories，Canada[J]. J of Marine Syst，1997，11：21-32.

［19］ 蘇 潔，吳輝碇，劉欽政，等. 渤海冰-海洋耦合模式[J]. 海洋學(xué)報(bào)，2005，27(1)：19-26.

Su Jie，Wu Huiding，Liu Qinzheng，et al. A coupled ice-ocean model for the Bohai Sea[J]. Acta Oceanologica Sinica，2005，27(1)：19-26(in Chinese).

［20］ Hunke E C，Lipscomb W H. CICE：the Los Alamos Sea Ice Model Documentation and Software User’s Manual Version 4. 0 LA-CC-06-012[R]. USA：T-3 Fluid Dynamics Group，Los Alamos National Laboratory，2008.

［21］ 張 寅，羅亞麗，管兆勇. NCEP全球預(yù)報(bào)系統(tǒng)在ARM SGP站點(diǎn)預(yù)報(bào)大氣溫度、濕度和云量的檢驗(yàn)[J]. 大氣科學(xué)，2012，36(1)：170-184.

Zhang Yin，Luo Yali，Guan Zhaoyong. Temperature，relative humidity，and cloud fraction predicted by the NCEP global forecast system at the ARM SGP site during 2001—2008：Comparison with ARM observations[J]. Chinese Journal of Atmospheric Sciences，2012，36(1)：170-184(in Chinese).

［22］ 程濟(jì)生. 黃渤海近岸水域生態(tài)環(huán)境與生物群落[M].青島：中國(guó)海洋大學(xué)出版社，2004.

Cheng Jisheng. The Ecological Environment and Biological Communities in the Coastal Waters of Bohai Sea and Yellow Sea[J]. Qingdao：Ocean University of China Press，2004(in Chinese).

［23］ 李孟國(guó)，鄭敬云. 中國(guó)海域潮汐預(yù)報(bào)軟件Chinatide的應(yīng)用[J]. 水道港口，2007，28(1)：65-68.

Li Mengguo，Zheng Jingyun. Introduction to Chinatide software for tide prediction in China seas[J]. Journal of Waterway and Harbor，2007，28(1)：65-68(in Chinese).

［24］ 宋文鵬. 渤海冬、夏季溫鹽場(chǎng)結(jié)構(gòu)及其海流特征分析[D]. 青島：中國(guó)海洋大學(xué)海洋環(huán)境學(xué)院，2009.

Song Wenpeng. The Analysis of the Structure of T-S and the Current Characteristics in Baohai Sea During Winter and Summer[D]. Qingdao：College of Physical and Environmental Oceanography，Ocean University of China，2009(in Chinese).

A Three-Dimensional Numerical Model of Sea Ice Evolution in the Bohai Sea

Zhang Qinghe1，Zhang Na1,2
(1. State Key Laboratory of Hydraulic Engineering Simulation and Safety，Tianjin University，Tianjin 300072，China；2. School of Civil Engineering，Tianjin Institute of Urban Construction，Tianjin 300384，China)

An ice-ocean coupled sea ice numerical simulation system of the Bohai Sea，including the dynamic processes and thermodynamic processes of sea ice，was developed based on the unstructured grid three-dimensional finitevolume coastal ocean model(FVCOM). Reasonable values of ice thickness，critical thickness for thinnest ice，ratio of shortwave radiation penetrating into the ice and heat transfer coefficient were discussed. The sea ice evolution of the Bohai Sea in 2003—2004 winter was simulated with the hourly time-series atmospheric data of NCEP(National Centers for Environmental Prediction)as a driving force. The simulated results of sea water temperature and salinity of the Bohai Sea are in good agreement with the measured data. The simulated ice coverage of the Bohai Sea is close to the image obtained by satellite remote sensing in the same period. The simulated sea ice area and sea ice thickness of Liaodong Bay are basically the same as the measured data for the start date，the end date，the maximum value and the corresponding time. It is shown that the presented three-dimensional numerical model of sea ice can reasonably simulate the process of hydrodynamic，sea water temperature and salinity，and the evolution of sea ice forced by meteorological conditions. The model can be applied to medium- and long-term sea ice numerical simulation with high accuracy and computational efficiency.

sea ice in the Bohai Sea；finite-volume coastal ocean model；unstructured grid；numerical simulation

O242；P731.15

A

0493-2137(2013)04-0333-09

DOI 10.11784/tdxb20130408

2012-12-04；

2012-12-13.

國(guó)家自然科學(xué)基金創(chuàng)新研究群體科學(xué)基金資助項(xiàng)目(51021004)；天津市自然科學(xué)基金重點(diǎn)資助項(xiàng)目(12JCZDJC30200).

張慶河（1966— ），男，博士，教授，qhzhang@tju.edu.cn．

張 娜，coastlab@163.com.