基于實測冰芯數據的南極海冰厚度的數值模擬

白乙拉，張韓晗

（渤海大學 數理學院，遼寧 錦州 121013）

0 引 言

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

1 模型建立

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

海冰冰域方程：

冰面溫度：

冰水交界面條件：

初始冰厚：

其中：

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

2 邊界條件

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

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

其中T1為冰表面溫度，Ta為當地大氣溫度。

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

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

3 數值解法

對方程（1）的定解區域做網格剖分，取空間步長為h，時間步長為τ，Tj，n表示xj＝jh處時刻tn＝nτ的溫度。

圖1 計算網格示意圖

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

設（S（t），nτ）點為B點，設B點至（jh，nτ）點的距離為Pn，在（j－1）至（j＋1）h點上采用Lagrange三點插值，選取（（j－1）h，Tj－1，n），（jh，Tj，n），（（jh＋Pnh），TB）三點進行插值，經計算可導出

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

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

上述方程的數值求解方法參見文獻［11－12］。

4 數值模擬

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

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

圖2 2006年南極中山站內拉灣實測冰厚與模擬冰厚對比

5 結 論

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

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

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

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

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

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

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

［6］劉欽政，白珊，黃嘉佑，等.一種冰－海洋模式的熱力耦合方案［J］.海洋學報，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］白乙拉，關博，劉丹.南極海冰導熱系數優化辨識［J］.內蒙古民族大學學報：自然科學版，2010，25（4）：361－364.

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

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