基于实测冰芯数据的南极海冰厚度的数值模拟

白乙拉，张韩晗

（渤海大学 数理学院，辽宁 锦州 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.