基于實(shí)測冰芯數(shù)據(jù)的南極海冰厚度的數(shù)值模擬

白乙拉，張韓晗

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

0 引 言

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

1 模型建立

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

海冰冰域方程：

冰面溫度：

冰水交界面條件：

初始冰厚：

其中：

上述各公式中，t為時間，S（t）為t時刻冰厚，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)過程中起著較為重要的作用，是太陽輻射、介質(zhì)物理特性等的函數(shù)。式（7）中S為海冰鹽度（‰）。式（8）中I0為冰層的透射率，這里N是云量參數(shù)，其取值從0到1。式（9）中ri為消光系數(shù)。其中導(dǎo)熱系數(shù)公式（7）是由第22次南極科考現(xiàn)場觀測數(shù)據(jù)辨識得到［10］，云量取值出自第22次南極科考現(xiàn)場實(shí)測氣象數(shù)據(jù)，海冰鹽度、密度取值引自第22次南極科考現(xiàn)場采集的冰芯分析數(shù)據(jù)，其余各參數(shù)表達(dá)式及其取值引自文獻(xiàn)［7－9］。

2 邊界條件

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

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

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

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

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

3 數(shù)值解法

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

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

由于冰的厚度是隨時變化的，所以冰水交界面是不確定的。因此在每一步時間計算中都要確定冰厚S（t）。而自由邊界不一定恰在網(wǎng)格節(jié)點(diǎn)上，所以在自由邊界附近采用非等距步長的方法來處理。設(shè)時刻t的自由邊界位置為S（t），這時S（t）已經(jīng)定下來，jh＜S（t）＜（j＋1）h，除了jh，S（t），（j＋1）h這3個點(diǎn)外，采用向前差分與Crank－Nicolson格式對方程（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)插值，選?。ǎ╦－1）h，Tj－1，n），（jh，Tj，n），（（jh＋Pnh），TB）三點(diǎn)進(jìn)行插值，經(jīng)計算可導(dǎo)出

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

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

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

4 數(shù)值模擬

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

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

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

5 結(jié) 論

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

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

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

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

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

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

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

［6］劉欽政，白珊，黃嘉佑，等.一種冰－海洋模式的熱力耦合方案［J］.海洋學(xué)報，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)化辨識［J］.內(nèi)蒙古民族大學(xué)學(xué)報：自然科學(xué)版，2010，25（4）：361－364.

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

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