白乙拉，張韓晗

（渤海大學(xué) 數(shù)理學(xué)院，遼寧 錦州 121013）

0 引 言

南極海冰是引起全球氣候變化的冷源的重要組成部分，極冰特別是海冰的變化對(duì)氣候的影響是當(dāng)今世界氣候研究計(jì)劃的前沿課題。近年來(lái)，許多學(xué)者針對(duì)極地海冰對(duì)氣候的影響及其時(shí)空變化特征和長(zhǎng)期變化趨勢(shì)做了大量的研究工作［1－5］。此外，依據(jù)海冰熱力學(xué)模式、動(dòng)力學(xué)模式、熱力學(xué)—— 動(dòng)力學(xué)耦合模式等進(jìn)行海冰數(shù)值模擬以及海冰熱力學(xué)系數(shù)的參數(shù)辨識(shí)研究也引起學(xué)者們的關(guān)注［6－8］。筆者根據(jù)中國(guó)第22次南極科學(xué)考察現(xiàn)場(chǎng)觀測(cè)數(shù)據(jù)，并重點(diǎn)考慮了太陽(yáng)輻射、云量、鹽度等因素對(duì)海冰厚度變化的影響，選取2006年4月1日0點(diǎn)至6月30日24點(diǎn)南極中山站內(nèi)拉灣海冰觀測(cè)點(diǎn)的實(shí)測(cè)氣溫、太陽(yáng)輻射、云量等現(xiàn)場(chǎng)觀測(cè)氣象數(shù)據(jù)，以及現(xiàn)場(chǎng)采集的冰芯分析數(shù)據(jù)作為海冰厚度數(shù)值模擬的原始數(shù)據(jù)，對(duì)南極中山站內(nèi)拉灣附近海域海冰的厚度變化進(jìn)行了數(shù)值模擬，并與現(xiàn)場(chǎng)觀測(cè)的海冰厚度變化數(shù)據(jù)進(jìn)行了對(duì)比。

1 模型建立

海冰熱力學(xué)模型一般由一維熱傳導(dǎo)方程描述。設(shè)南極海冰的比熱、密度、導(dǎo)熱系數(shù)分別為C、ρ、λ，取冰層表面一點(diǎn)為坐標(biāo)原點(diǎn)O，過(guò)原點(diǎn)O垂直向下為x軸的正向，T（x，t）表示t時(shí)刻冰層x點(diǎn)處的溫度，冰層和海水交界面為x＝S（t）。則南極海冰厚度數(shù)值模擬的熱力學(xué)方程可表示為［7］

海冰冰域方程：

冰面溫度：

冰水交界面條件：

初始冰厚：

其中：

上述各公式中，t為時(shí)間，S（t）為t時(shí)刻冰厚，C和ρ分別為海冰的比熱（J／kg·℃）和密度（kg／m3），λ（T）為熱傳導(dǎo)系數(shù)（W／m·℃），L為單位質(zhì)量海冰的熔解潛熱（J／kg），F(xiàn)w為海冰的熱通量。式（1）中的q（x，t）為冰內(nèi)熱源項(xiàng)，在熱傳導(dǎo)過(guò)程中起著較為重要的作用，是太陽(yáng)輻射、介質(zhì)物理特性等的函數(shù)。式（7）中S為海冰鹽度（‰）。式（8）中I0為冰層的透射率，這里N是云量參數(shù)，其取值從0到1。式（9）中ri為消光系數(shù)。其中導(dǎo)熱系數(shù)公式（7）是由第22次南極科考現(xiàn)場(chǎng)觀測(cè)數(shù)據(jù)辨識(shí)得到［10］，云量取值出自第22次南極科考現(xiàn)場(chǎng)實(shí)測(cè)氣象數(shù)據(jù)，海冰鹽度、密度取值引自第22次南極科考現(xiàn)場(chǎng)采集的冰芯分析數(shù)據(jù)，其余各參數(shù)表達(dá)式及其取值引自文獻(xiàn)［7－9］。

2 邊界條件

在上述問(wèn)題中，海冰的上邊界條件，即海冰表面溫度由大氣和冰面之間的熱交換所決定。筆者考慮到冰表面溫度的變化主要與大氣溫度有關(guān)，忽略冰表面溫度與冰內(nèi)熱傳導(dǎo)二者的耦合過(guò)程。根據(jù)前人研究成果可知，冰表面溫度與大氣溫度二者近似成線性關(guān)系［7］。

依據(jù)前人的經(jīng)驗(yàn)公式，嘗試采用多種不同的函數(shù)描述南極中山站內(nèi)拉灣附近冰表面溫度與大氣溫度的關(guān)系，通過(guò)最小二乘法，建立相應(yīng)的辨識(shí)冰表面溫度與大氣溫度關(guān)系的最優(yōu)化模型辨識(shí)，應(yīng)用MATLAB軟件編程求解該優(yōu)化模型，得到南極中山站內(nèi)拉灣附近冰表面溫度與大氣溫度的關(guān)系為

其中T1為冰表面溫度，Ta為當(dāng)?shù)卮髿鉁囟取?/p>

根據(jù)辨識(shí)得到的南極海冰表面溫度與大氣溫度的線性關(guān)系（10），將各觀測(cè)時(shí)刻的大氣溫度換算為相應(yīng)時(shí)刻的海冰表面溫度，從而得到海冰生消熱力學(xué)問(wèn)題的上邊界條件。

海冰的下邊界條件，即冰水交界面的溫度，由于海冰地理位置及鹽度情況不同，其結(jié)冰點(diǎn)也會(huì)略有不同，筆者根據(jù)現(xiàn)場(chǎng)觀測(cè)情況結(jié)冰點(diǎn)取為－1.81℃。

3 數(shù)值解法

對(duì)方程（1）的定解區(qū)域做網(wǎng)格剖分，取空間步長(zhǎng)為h，時(shí)間步長(zhǎng)為τ，Tj，n表示xj＝j(luò)h處時(shí)刻tn＝nτ的溫度。

圖1 計(jì)算網(wǎng)格示意圖

由于冰的厚度是隨時(shí)變化的，所以冰水交界面是不確定的。因此在每一步時(shí)間計(jì)算中都要確定冰厚S（t）。而自由邊界不一定恰在網(wǎng)格節(jié)點(diǎn)上，所以在自由邊界附近采用非等距步長(zhǎng)的方法來(lái)處理。設(shè)時(shí)刻t的自由邊界位置為S（t），這時(shí)S（t）已經(jīng)定下來(lái)，jh＜S（t）＜（j＋1）h，除了jh，S（t），（j＋1）h這3個(gè)點(diǎn)外，采用向前差分與Crank－Nicolson格式對(duì)方程（1）進(jìn)行離散，整理便可得到

設(shè)（S（t），nτ）點(diǎn)為B點(diǎn)，設(shè)B點(diǎn)至（jh，nτ）點(diǎn)的距離為Pn，在（j－1）至（j＋1）h點(diǎn)上采用Lagrange三點(diǎn)插值，選取（（j－1）h，Tj－1，n），（jh，Tj，n），（（jh＋Pnh），TB）三點(diǎn)進(jìn)行插值，經(jīng)計(jì)算可導(dǎo)出

TB為冰水交界面的溫度。利用（11）與向前差分即可得到自由邊界附近式（1）的離散方程

利用式（12）再結(jié)合拉格朗日中值定理，可以導(dǎo)出式（3）的離散方程

上述方程的數(shù)值求解方法參見(jiàn)文獻(xiàn)［11－12］。

4 數(shù)值模擬

在2006年3月到12月期間，中國(guó)第22次南極科學(xué)考察對(duì)南極中山站內(nèi)拉灣附近海域固定冰冬季物理性質(zhì)進(jìn)行了現(xiàn)場(chǎng)觀測(cè)，獲得了大量的有關(guān)氣溫、云量、太陽(yáng)輻射、冰溫、冰蓋厚度等現(xiàn)場(chǎng)觀測(cè)數(shù)據(jù)。從2006年3月末開(kāi)始，對(duì)冰面下0cm、6cm、12cm、18cm、24cm直至冰下2 750cm每間隔一定距離的海冰和海水溫度進(jìn)行了連續(xù)觀測(cè)，直到12月下旬結(jié)束，觀測(cè)期間采樣間隔為30min。針對(duì)不同經(jīng)緯度、不同厚度的海冰，并每間隔一定時(shí)間對(duì)不同位置的冰芯進(jìn)行采樣分析，得到海冰剖面的各種冰芯數(shù)據(jù)。此外，利用16套熱電阻絲手動(dòng)測(cè)量裝置，對(duì)冰雪厚度變化也進(jìn)行實(shí)測(cè)，每套間隔15m，4月初開(kāi)始測(cè)量，每天測(cè)量一次，測(cè)量持續(xù)到10月初，獲得了南極中山站內(nèi)拉灣附近海域海冰冬季生長(zhǎng)的現(xiàn)場(chǎng)觀測(cè)數(shù)據(jù)。

筆者利用2006年4月初至6月末海冰觀測(cè)點(diǎn)的實(shí)測(cè)氣溫、冰芯分析數(shù)據(jù)以及相應(yīng)的太陽(yáng)短波輻射、云量實(shí)測(cè)數(shù)據(jù)，作為海冰厚度數(shù)值模擬的原始數(shù)據(jù)。對(duì)于給定時(shí)間段起始時(shí)刻海冰層各觀測(cè)點(diǎn)的實(shí)測(cè)溫度數(shù)據(jù)，經(jīng)線性插值得到網(wǎng)格各節(jié)點(diǎn)的初始溫度，以各時(shí)刻的實(shí)測(cè)大氣溫度由公式（10）換算得到的海冰表面溫度，經(jīng)插值得到的各時(shí)間節(jié)點(diǎn)的溫度數(shù)據(jù)作為上邊界條件，由于冰層厚度是隨時(shí)變化的，下邊界是冰水交界面，下邊界處各時(shí)刻的溫度為南極海冰冰點(diǎn)溫度TB，根據(jù)南極海冰實(shí)測(cè)溫度數(shù)據(jù)，冰點(diǎn)溫度TB取值為－1.81℃。由冰層不同點(diǎn)采集的冰芯樣本分析測(cè)得的海冰密度、鹽度數(shù)據(jù)也略有不同。筆者根據(jù)現(xiàn)場(chǎng)采集的冰芯分析數(shù)據(jù)插值得到網(wǎng)格各節(jié)點(diǎn)的密度、鹽度數(shù)據(jù)。數(shù)值模擬結(jié)果與實(shí)測(cè)數(shù)據(jù)對(duì)比情況如圖2所示。

圖2 2006年南極中山站內(nèi)拉灣實(shí)測(cè)冰厚與模擬冰厚對(duì)比

5 結(jié) 論

筆者根據(jù)第22次南極科考實(shí)測(cè)數(shù)據(jù)以及現(xiàn)場(chǎng)采集的冰芯分析數(shù)據(jù)，考慮了太陽(yáng)輻射、云量、鹽度等多種因素對(duì)海冰厚度變化的影響，選取2006－04－01T00∶00—06－30T24∶00南極中山站內(nèi)拉灣的實(shí)測(cè)氣溫、冰芯分析數(shù)據(jù)以及相應(yīng)的太陽(yáng)輻射、云量實(shí)測(cè)數(shù)據(jù)，作為海冰厚度數(shù)值模擬的原始數(shù)據(jù)，通過(guò)最小二乘法辨識(shí)氣溫與冰表面溫度的關(guān)系，用該關(guān)系式（10）換算得到的冰表面溫度作為上邊界條件，依據(jù)現(xiàn)場(chǎng)實(shí)測(cè)數(shù)據(jù)取結(jié)冰點(diǎn)－1.81℃作為下邊界條件，對(duì)南極中山站附近海冰厚度變化進(jìn)行了數(shù)值模擬。由圖2可以看出數(shù)值模擬結(jié)果與南極海冰厚度實(shí)測(cè)數(shù)據(jù)吻合良好，表明采用文章方法進(jìn)行極地海冰厚度的數(shù)值模擬是可靠而有效的。

致謝 大連理工大學(xué)海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室李志軍教授，提供了中國(guó)第22次南極科學(xué)考察獲得的海冰現(xiàn)場(chǎng)觀測(cè)數(shù)據(jù)和冰芯分析數(shù)據(jù)，在此表示衷心的感謝！

［1］解思梅，鄒武，王毅.南極海冰的長(zhǎng)期變化趨勢(shì)［J］.海洋預(yù)報(bào)，1996，13（3）：21－29.

［2］唐述林，秦大河，任賈文，等.極地海冰的研究及其在氣候變化中的應(yīng)用［J］.冰川凍土，2006，28（1）：91－100.

［3］馬麗娟，陸龍驊，卞林根.南極海冰的時(shí)空變化特征［J］.極地研究，2004，16（1）：29－37.

［4］雷瑞波，李志軍，竇銀科，等.南極中山站附近固定冰生消過(guò)程觀測(cè)［J］.水科學(xué)進(jìn)展，2010，21（5）：708－712.

［5］卞林根，林學(xué)春.近30年南極海冰的變化特征［J］.極地研究，2005，17（4）：233－244.

［6］劉欽政，白珊，黃嘉佑，等.一種冰－海洋模式的熱力耦合方案［J］.海洋學(xué)報(bào)，2004，26（6）：13－21.

［7］CHENG B，VIHMA T，PIRAZZINI R，et al.Modelling of superimposed ice formation during the spring snowmelt period in the Baltic Sea［J］.Ann Glacio，2006，44（1）：139－146.

［8］SHI Liqiong，BAI Yila，LI Zhijiang.Preliminary results on relationship brtween thermal diffusivity and porosity of sea ice in the Antarctic［J］.Chin J Polar Sci，2009，20（1）：72－80.

［9］YEN Y C.Review of thermal properties of snow，ice and sea ice［R］.CRREL Report 81－10，Hanover，N.H.：U.S.Army，Corps of Engineers，Cold Regions Research and Engineering Laboratory，1981.

［10］白乙拉，關(guān)博，劉丹.南極海冰導(dǎo)熱系數(shù)優(yōu)化辨識(shí)［J］.內(nèi)蒙古民族大學(xué)學(xué)報(bào)：自然科學(xué)版，2010，25（4）：361－364.

［11］肖建民，金龍海，謝永剛，等.寒區(qū)水庫(kù)冰蓋形成與消融機(jī)理分析［J］.水利學(xué)報(bào)，2004（6）：80－85.

［12］馬振寧，高明.金屬材料裂紋尖端溫度場(chǎng)和熱應(yīng)力場(chǎng)的數(shù)值模擬［J］.沈陽(yáng)師范大學(xué)學(xué)報(bào)：自然科學(xué)版，2006，24（1）：34－37.