1 引言

气候变化情景下，冰盖表面物质平衡可以直接影响海平面变化，是研究南北极冰盖物质输入和输出的重要参数^[1]。据研究，由于西南极的冰架快速融化和崩塌，南极冰盖对海平面的贡献已由1992~2001年间的0.08(-0.10至0.27) mm/a上升至2002~2011年间的0.40(0.20至0.61) mm/a^[2]。与此同时，有研究指出南极冰盖物质损失的增强，部分被东南极降雪的增加所抵消^[3~5]。通过冰芯记录和气候模拟发现南极的降雪量随着气温的上升而增加；通过卫星雷达和激光测高研究发现近几十年来东南极的降雪有所增加^[6]，特别是Droning Maud Land地区，2009年后发生了多次高降雪事件^{[7, 8]}。当然，也有研究指出虽然过去几十年中南极海岸地区降雪增强，但整个南极的平均降雪并没有发生长期变化^[9~12]。

冰芯是提供南北极降雪历史记录的重要手段之一，也是校正冰雷达等时层^[13~15]、驱动冰盖流动模式和定年模式^[16]、以及评估区域气候模式^[17]的重要基础资料。然而，相比于南极广袤的面积，目前的冰芯记录仍然十分稀少，特别是南极冰盖高原。由于局地气候背景的不同，南极某些地区的冰芯记录显示过去一百年中降雪持续增加^[18]，另外一些地区则处于减少趋势^[19]。从全南极尺度来看，虽然南极冰盖雪积累率有所增加，但趋势并不显著^[20]。

一般来讲，随着纬度升高和距海岸距离的加大，南极冰盖降雪逐步降低，所获取的冰芯记录分辨率同时下降，因此很难同时获取时间分辨率高、时间序列长的冰芯。在近岸地区，因为冰盖厚度较薄，年降雪量高，因此很少有超过200年的冰芯记录。为了更好地研究南极冰盖物质平衡和气候变化，冰盖物质平衡和海平面变化研究计划(Ice Sheet Mass Balance and Sea Level，简称ISMASS)多次呼吁加强东南极冰盖近岸地区冰芯记录的获取(http://www.scar.org/ssg/physical-sciences/ismass)。

受到冰盖/冰川运动的影响，冰芯定年的误差会随着深度的增加而增加，因为冰流的水平运动会造成年层向四周扩展，因此而变薄。一般来说，在东南极内陆冰盖地区，由于应力变形引起的水平延展率在0.01 % /a；若冰芯的长度不超过100m，相对于冰盖的厚度就可以忽略，相对误差也不会超过5 % ^[13]。但是在运动较快的冰川、近岸冰盖或冰架地区，即使浅冰芯也会受到很大影响，需要进行校正。

过去几十年中，大多数冰芯研究都集中在气温变化和环境变化上(同位素、可溶性离子和痕量金属记录)^[21]，降雪记录恢复则少有出现。即使冰芯降雪记录，很多研究也未考虑冰流作用^[20]，可能造成了不同程度的误差，对于评估南极冰盖表面物质平衡产生了无法估量的影响。

2002年，中国南极考察队在东南极冰盖近岸地区LGB69钻取了一支102m的冰芯。本文选用该冰芯，通过两种冰流模型消除了冰盖运动的影响，恢复其过去300年历史降雪记录。本研究所得结果对过去全球变化研究计划2000年工作组(PAGES 2k)具有重要帮助。

2 数据与方法 2.1 冰芯采集

2002年1月，中国南极科学考察第18次考察队，在东南极冰盖Lambert冰川流域LGB69(70°50′06.6″S，77°04′28.9″E；1850m a.s.l.，图 1)，采用中国科学院冰冻圈科学国家重点实验室自主研发的BXZ-100电动机械钻钻取了一支102m的冰芯。该点距南极大陆海岸约为160km，据气象站观测年均积累率约为70cm/a，风速约为8.5m/s，主风向为东北偏东风57度，表面分布有高度介于0~20cm的雪丘^{[22, 23]}。冰芯钻取期间，用游标卡尺测量冰芯直径，记录每段长度，并用电子秤称重，放入塑料桶内依顺序编号，最后装入保温箱保存。冰芯钻取和搬运严格按规范操作，避免污染和融化。

2.2 分析与定年

定年是对冰芯资料进行研究的第一步，也是其他研究工作的基础。以海洋飞沫来源为主的离子如Cl^-、Na⁺和K⁺受到海洋表面环境的影响，在冰芯中呈现出明显的季节循环，可以在高积累地区用于定年；氢氧稳定同位素比率与降水发生时的气温有关，也是定年的重要指标；火山爆发释放大量的SO₂，约一个月后被氧化成硫酸，并在空气中留存很长时间，所以能随着大气环流传输到南极地区，一定时间后(取决于喷发地点，一般不超过3年)沉降到南极雪冰中，是长时间序列冰芯定年的重要手段；放射性物质参考层、理论模型、古温度记录等也可以用于冰芯定年^[24]。

为了能获得准确的冰芯年代记录，考察队回国后将冰芯送至美国缅因大学进行分样处理，并利用离子色谱仪(Dionex IC)分析了主要可溶性离子(Na⁺、K⁺、Mg²⁺、Ca²⁺、Cl^-、NO₃^_、SO₄^2-和MSA)。根据火山事件法，Li等^[25]确定了冰芯中的主要年层深度；之后我们根据主要海盐离子的季节性变化，利用年层法确定了该冰芯的详细年代，该冰芯顶部年层为2001年，底部年层为1708年(图 2a)。季节性参数所确定的定年结果与冰芯中保存的历史重大火山记录年代^[25]一致，通过与全球火山喷发计划进行对比^[25]，可知定年误差约为± 2a。

2.3 降雪量校正方法

如前言所述，沉积后过程不仅仅会影响垂直方向上的雪层密度，还会造成雪层的侧向变形。如果不考虑冰流所导致的雪层侧向变形，会低估冰芯年层厚度，进而影响对降雪/积累率的重建，这种误差会随着深度的增加而增大，也与冰流速度、底部滑动速度等因素有关。目前有3种方法/模型可用于校正该误差：1) 幂律模型(power-law model)^[26]；2) 分段线性模型(a piece-wise linear model，后称D-J模型)^[27]；3) 完全线性模型(fully linear model，后称Nye模型)^[28]。其中幂律模型需要测量钻孔的水平位移，我们未在LGB69冰芯钻取的时候测量该参数，因此，本研究使用D-J模型和Nye模型进行计算。

2.3.1 Nye模型

Nye模型假设在冰盖的任一垂直面上垂直应变率全都相同，同时忽略冰盖底部的融化作用，并假定冰层运动为水平平行的。因此，冰层的初始厚度(λ_H)与减薄后的厚度的比值和冰层所处的高度之间存在直接联系。关系式如下：

(1)

公式(1) 中y和H分别是冰层现在和初始的距离基岩高度。由此，我们可以得到在该模型中单位时间内的变化为：

(2)

τ为1年。整合上式得到关于时间的函数如下：

(3)

由此，结合公式(1)，我们可以得到不同年代的冰层距离基岩的高度：

(4)

因此，冰盖表面下H-y处的冰层年龄为：

(5)

2.3.2 D-J模型

该模型假设水平速度廓线v与据基岩高度(y)、水平位移距离(x)之间存在一定的关系，表示为：

(6)

结合Glen法则(Glen's law)^[29]，可以得到：

(7)

公式(7) 中，σ为剪切应力；k′和n分别为取决于温度的常数；y为冰层离底部的高度。只有当任意高度(y)纵向的应变率相比于剪切应力都较少时，公式(7) 才需要进行修正。然而，利用Weertman公式^[30]对公式(7) 进行转换发现，我们可以用一个简单的近似来对公式(7) 进行求解。近似结果如下：

(8)

Dansgaaard和Johnsen^[27]研究发现，冰盖内部水平运动速度会随着深度的增加而增加，直到约400m，水平速度基本不再增加。因此，在计算中，h一般设定为400m，在400m以上认定f(y)为定值。

另一方面，根据冰的不可压缩性可以得到：

(9)

由此垂直速率的部分可以计算为：

(10)

并且，结合公式(8) 可以得到：

(11)

当y=h时，

(12)

当y=H时，

(13)

(14)

将式(12)、(13) 和(14) 代入公式(11) 可以得到：

(15)

将上式的结果与Nye模型的结果进行比较可以发现，当h≤y≤H时，可以用Nye模型的方程来表示，在该范围内Nye模型的直线与y轴的交点为h/2。分别用H-h/2和y-h/2来代替H和y，代入到式(1) 可以得到：

(16)

(17)

而在范围为0≤y≤h时，有式(15) 可以得到：

(18)

当间距不断减少，h和t_h无限趋近与y和t，即h→y和t_h→t。其中t_h为当y=h时，可由式(17) 计算得。由此：

(19)

或者，

(20)

同时结合式(18) 和式(19) 可得：

(21)

由上文介绍可知，Nye模型为经典模型，其假设冰层的初始厚度λ_H与减少后的厚度的比值与冰层所处的高度之间呈线性关系，即不同高度层上冰层扩散速率均为恒定值。D-J模型则假设水平扩散速率与高度并不是简单的线性关系，认为水平速度廓线与距基岩高度、水平位移距离x之间存在一定的函数关系。相比Nye模型，D-J模型更接近实际情况，并通过观测及Glen法则^[29]发现水平扩散速率与高度之间的呈曲线关系。

3 结果与讨论 3.1 年降雪量/雪积累率恢复

如上所述，在钻取LGB69冰芯时，我们记录了冰芯直径、每段长度以及重量，进而我们可以获取该冰芯的密度随深度的变化。通过图 3可以看出，LGB69冰芯的表层雪密度约为400kg/m³，受到沉积后作用特别是密实化作用的影响，在12.9m深度达到第一临界密度(550kg/ m³)，密实化过程由机械压密阶段转变为塑性变形和再结晶阶段；然后密度增大减缓，在80m深度达到第二临界密度(830~840kg/ m³)，冰内气泡封闭，雪转换为冰，密度变化趋势进一步减小。

利用拟合公式对密度进行插值，结合冰芯定年深度，我们可以计算得到LGB69的逐年降雪序列。之后，使用Nye模型和D-J模型，我们计算得到了LGB69站点1708~2001年间年降雪量的校正序列(表 1和图 4a)，该区域平均年降雪量为261.7±69.7kg/ m²a。根据经验，我们将过去293年划分为不同的时期(图 4b，表 1)，分别具有不同的降雪量(以D-J模型校正为例)：1708~1863年，年降雪量为264.8±63.4kg/ m²a；1863~1884年，年降雪量为236.4±50.3kg/ m²a；1884~1917年，年降雪量为319.6±84.6kg/ m²a；1917~1937年，年降雪量为293.5±66.4kg/ m²a；1937~1971年，年降雪量为207.5±48.1kg/ m²a；1971~2001年，年降雪量为240.2±58.8kg/ m²a。从图 4b和表 1可以看出，1708~1863年间，LGB69年降雪量非常稳定，之后则变化快速，特别是1884~1917年间，该处降雪突然升高，显著高于平均水平。主要原因可能有两个：一是全球气温快速变化使得中低纬度向南极输送水汽的能力发生快速变化^[32]；二是LGB69地处冰盖边缘，其降水与普里兹湾附近的海冰面积变化有关^[33]，海冰面积减小有利于更多水汽输送至南极冰盖，反之亦然。同时，南极海冰变化还受到南极环状模(SAM)强度变化的影响，与全球气温升降并不是严格正相关关系^[34]。

在南极冰盖海岸地区，年降雪量主要受到中高纬大气环流和局地气候环境的影响，由气旋风暴活动引起的降水事件影响尤为显著^[35]。杨佼等^[36]通过分析1979~2001年间该区域降雪极高年与极低年的大气环流场，发现位于南印度洋高纬地区的南印度洋低压带(SIOL)的强度和活动是控制Lambert流域降水的主要气候系统。由此可以推论，1863年之后，南印度洋高纬地区的气旋活动变化剧烈，在SIOL增强的年份，气压加深并向南扩张，影响范围扩大，海洋性气团的向极输送增强，LGB69地区降水增多；反之亦然。近年来，有研究指出海冰面积变化可能影响冰盖降水，主要机制是海冰面积减少有利于水汽蒸发和输送^[33]；同时，也有研究认为，向极输送大气活动的增强(同时海冰面积减小)起到了最重要的作用而不是海冰面积减小^[32]。有关研究可能解释为何LGB69地区降雪量多变，但仍未在Lambert流域和Prydz湾开展。

3.2 误差分析

为了进一步检验不同方法反演降雪量历史记录的误差，我们计算得到了Nye模型和D-J模型的校正序列与未校正序列的差及相对误差。如图 5所示，自冰盖表面至下部，由于冰盖运动造成的年层减薄越来越明显，最大误差可达20.0kg/ m² a(Nye模型)和22.5kg/ m² a(D-J模型)，最大相对误差可达5.6 % (Nye模型)和6.3 % (D-J模型)，模型修正差异明显。需要指出的是，由于算法的不同，两个模型之间的差异可达0.7 %。

4 结论

结合火山沉积记录和可溶性离子季节循环，本文重建了LGB69钻取的102m冰芯的年降雪量/雪积累率历史记录，发现该冰芯最深处年龄为1708年。之后，通过Nye模型和D-J模型分别对该序列进行了校正，消除了由于侧向应力所造成了误差。与常规重建记录相比，该序列精度非常高，对于南极冰盖长期物质平衡估算具有重要参考意义。

与东南极冰盖其他冰芯记录相比，Lambert冰川流域降雪历史变化具有明显不同。Thomas等^[37]重建了公元1000~2000年间南极降雪历史记录，发现南极地区降雪在1900~2000年间整体处于增加的趋势，特别是唐宁毛德地和南极半岛地区，其降雪量快速增加。在流域尺度上，同处Lambert冰川流域的LGB65冰芯记录显示1900年以来局地年降雪量缓慢增加^[23]，与LGB69地区的快速下降具有明显区别。与LGB69相似的只有维多利亚地的几支冰芯^[37]，究其原因，可能与LGB69及维多利亚地更多地受到附近海冰快速变化有关，详细机制还需要进一步分析。

东南极冰盖面积广袤，不同地区降雪量历史记录的恢复，对于评估南极与海平面的相互作用具有重要影响。本研究发现，在LGB69附近区域，冰盖运动造成的降雪量低估在102m深度可达5.6 % (Nye模型)和6.3 % (D-J模型)，随着深度的增加，该误差还会继续增大。因此，在使用南北极冰芯特别是快速冰流区的冰芯反演历史降雪记录时，应考虑该误差造成的影响。

致谢: 感谢同行评审专家和编辑部老师提出的宝贵修改意见!感谢评审专家和编辑部老师提出的宝贵修改意见！本研究得到国家自然科学基金杰出青年基金项目“极地冰冻圈”(批准号：41425003)、国家海洋局南北极环境综合考察与评估专项项目(批准号：CHINARE2016-2020) 和中国气象科学研究院基本科研业务费项目(批准号：2017Z002) 共同资助。