1. 引言
格陵兰岛位于大西洋与北极之间，经纬度59.46˚N~83.39˚N；73.08˚W~11.39˚W，岛屿面积216万平方公里，是世界最大的岛屿；格陵兰冰盖(GrIS)作为大小仅次于南极冰盖的世界第二大冰盖，是气候变化的重要调控因素之一。
GrIS巨大的蓄积率和显著的径流以及融水生成，使其成为一个高度动态变化的地方。GrIS的快速变化可能会影响全球海平面上升，并可能改变AMOC和全球气候 [1] [2] 。GrIS是正在发生的气候变化的一个指标，整个冰盖对气候的影响已经被许多研究所证明。1990s使用机载高度计测量的GrIS显示其边缘地区冰盖变薄，在中部区域有轻微的增厚 [3] 。随后的测量结果显示，沿海地区冰盖变薄加剧，并集中在被冰川出口占据的狭窄山谷中 [4] 。高度计数据研究显示 [5] [6] ，从1992年到2003/4年，2000米以上的GrIS内部出现了显著的(2~5 cm/yr)增长，这可能是由大气湿度和降水量的增加和(或)风暴路径的移动所致 [7] 。Velicogna等人 [8] 利用22个月的GRACE重力场首次从太空测量估算整个GrIS的质量变化。Ramillien等人 [9] 于2006年使用GRACE数据(利用新的大地水准面)计算了(2002年7月~2005年3月)这一时间段的GrIS质量变化及其对海平面上升的贡献。由于GRACE卫星数据的采样时间跨度非常短，尚无法区分年际振荡与气候变化相关的长期趋势。Velicogna等人 [10] 使用了一段较长的连续时间(2002~2009年)的GRACE数据分析格陵兰的冰盖质量变化。但是，这些开创性的分析估计之间存在许多差异和不确定性，对于影响格陵兰的冰盖质量变化的物理机制的探讨也比较少。
Hanna等人 [11] 将1990年代以来冰盖变薄融化及径流增加与夏季温度显著增加归因于全球气候变暖。如果GrIS全部融化，全球海平面将上升约7.4 m，另外研究也表明GrIS容易受到人类行为导致的气候变暖的影响 [12] 。自上世纪90年代初以来，随着大气变暖，GrIS地表径流增加，格陵兰表面融化过程约占整个GIS总质量损失的一半 [13] ，在表面温度长期升高的情况下，这一比例甚至更大 [14] 。然而，由于针对GrIS整体研究的数据跨度较短，尚未有研究提供一个更有说服力的多年代际视角，来研究冰盖如何应对长期气候变化，尤其是自20世纪70年代以来明显的全球变暖现象。
本文基于更新的GRACE数据重新对GrIS的时空分布特征进行分析，并使用区域模式数据结合多种研究方法探究控制冰盖质量变化的主要原因。
2. 数据资料与研究方法
2.1. 数据资料
格陵兰质量变化研究的主要数据是来自GRACE卫星的观测数据 [15] 。版本是GFZ RL05，水平分辨率1度，数据变量为液态水当量厚度，单位是厘米。从2003年1月到2016年12月，GRACE数据有13个月的缺测。格陵兰岛上的天气站(丹麦气象研究所，DMI)提供了环绕格陵兰岛的观测温度(GrT)站点数据，72个气象站中有68个站点的温度观测数据是从1958年开始的，其中的四个站点(Upernavik, Ilulissat, Nuuk, Tasiilaq)有长达百年尺度的持续观测数据。研究格陵兰表面能量变化时使用的向下短波辐射(SWD)、向下长波辐射(LWD)、向上长波辐射(LWU)、感热通量(SHF)、潜热通量(LHF)、反照率(AL)和云量(CC)等变量的日平均数据由区域气候模式Modèle Atmosphérique Régional (MAR)获得 [16] ，本文使用的MAR模式数据是由NCEP再分析数据驱动的，该模式在格陵兰岛运行的分辨率高达20 km，数据长度为1948~2017年。格陵兰周边的海表面温度(SST)数据来源于ERSST V3b，使用的数据时间长度为1979~2016。格陵兰的平均海平面气压数据是欧洲中期天气预报中心(ECWMF)的ERA-Interim的1979~2017年的1/8度分辨率月平均数据。
2.2. 研究方法
本文采用数理统计、线性拟合等方法对GrIS质量的时空特征进行分析，用质量收支平衡方法和表面能量平衡方法对控制GrIS质量变化的具体原因进行研究。最后使用经验模态分解方法(EMD) [17] 分析GrIS表面温度的不同时间尺度信号，并通过回归分析方法初步研究控制GrIS质量变化的动力机制。
3. 格陵兰冰盖质量变化特征
3.1. 格陵兰冰盖质量的时空变化
如图1所示，基于GRACE卫星数据显示的GrIS质量在21世纪初期呈现加速减少的趋势，其中2012年夏季还出现了极端融化事件。然而从2013年开始，GrIS质量损失的速度开始减缓。
Figure 1. The time series of GrIS mass change observed by GRACE satellite since 2003. The grey points represent the data of missing months filled with linear interpolation method
图1. GRACE卫星自2003年以来观测到的GrIS质量变化时间序列，灰点代表使用线性插值法填充的缺测月份的数据
从空间分布来看，自有GRACE卫星观测以来，2013年以前GrIS东南区域呈现质量加速损失的趋势，2013年以后GrIS西北区域保持与前十年几乎一致的变化速度，而其南部冰盖加速减损的趋势明显减缓，在格陵兰东部区域冰盖质量变化趋势甚至转变为轻微增加(图2)。尽管发生质量损失减速变化的时间段尚短(仅五年长度)，线性趋势的估计结果对于时间段的挑选十分敏感，GrIS质量加速损失的趋势减缓现象仍是十分明显的。
Figure 2. (a) The linear trend of GrIS melting observed by GRACE satellite in 2003-2012 and (b) the spatial pattern in 2013-2016, respectively. Except for the light grey dotted area in (b), all other spatial lattice points have passed 95% confidence test for linear trend estimation
图2. (a) GRACE卫星观测的GrIS融化的线性趋势分别在2003~2012年和 (b) 2013~2016年的空间型。除了(b)中的浅灰色打点区域外，其他所有空间格点对线性趋势的估计都通过95%置信检验
3.2. 格陵兰冰盖质量收支
质量收支平衡方法表明冰盖质量变化由表面质量平衡(SMB)和冰流量(D)决定，其中SMB由降水(雪)的压缩结冰量和融化或蒸发(升华)的冰雪量决定。总体来看，冰盖质量变化可以分成增加和减少两部分，增加的过程包括降雨、降雪、结冰，减少量则由表面融化、蒸发(升华)以及D控制，其中D与表面融化量密切相关 [18] 。观测结果与模型预测一致认为SMB而非D主导冰盖对21世纪海平面上升的贡献 [19] 。由于数据及方法的限制，本文暂不对D进行详细计算，主要对SMB进行分析研究。MAR模式数据中的SMB变量由降雪(SF)、降雨(RF)、径流(RU)、蒸发(SU)和地表水(SW)计算构成。所有冰盖质量变化控制因素中，SU和SW相比于其它量的变化可忽略不计，因此本文主要对SF、RF、RU进行分析研究。
模式数据中GrIS的SMB也呈现出与GRACE卫星观测数据一致的变化，2003~2012年SMB迅速减少，2013年开始有所回升。SMB的变化主要受RU控制，而RU主要由融化量决定(图3)，模式数据中降雪(雨)的影响并不明显。
图4显示了MAR模式数据中的SMB及其变量的夏季(JJA)季节变化。结果显示自2013年开始的SMB回升在夏季表现明显，并且SMB主要被RU(融化)所控制，RU与SMB的相关系数为0.99。夏季降雨降雪的量级相比融化径流小，并且没有明显变化趋势。
MAR模式数据中的SMB及其变量的冬季(DJF)季节变化如图5所示。冬季SMB主要受降雪控制，除了振幅在1997年以后减小，并无其他变化趋势。融化降雨和径流都接近于零。
Figure 3. The anomaly time series of melt, runoff, snowfall, rainfall and SMB in MAR model data from 1979 to 2017
图3. MAR模式数据中的融化、径流、降雪、降雨及表面质量平衡在1979~2017年的异常变化
Figure 4. The anomaly time series of summer melt, runoff, snowfall, rainfall and SMB in MAR model data from 1979 to 2017
图4. MAR模式数据中的夏季融化、径流、降雪、降雨及表面质量平衡的1979~2017年际变化
Figure 5. The anomaly time series of winter melt, runoff, snowfall, rainfall and SMB in MAR model data from 1979 to 2017
图5. MAR模式数据中的冬季融化、径流、降雪、降雨及表面质量平衡的1979~2017年际变化
通过对控制GrIS质量变化的SMB及其分量的变化趋势进行分析发现，降雨降雪并不是影响冰盖质量损失的主要因素，夏季融化量的改变才是控制SMB改变，进而使得冰盖质量变化的主要原因。
3.3. 格陵兰冰盖表面能量平衡
通过MAR模式数据的初步分析基本可以确定冰盖质量加速损失自2013年开始的减缓现象是被表面融化速度减缓控制。根据表面能量平衡方法，冰盖表面融化是由净热通量控制。表面能量平衡对应的热通量融化贡献计算公式为：
。其中
为净热通量(wm−2)；δt为单位时间间隔(s)；
是融解潜热(3.335 × 105 Jkg−1)；另外
是水的密度(1000 kgm−3)。对比MAR模式数据给出的融水量和利用表面能量平衡计算的净热通量对应的融化量，二者基本一致(图6)。
Figure 6. Variation of melting water (blue line) and melting water volume (red line) calculated by net heat flux from 1979 to 2017 based on MAR model data
图6. MAR模式数据的1979~2017年融水量变化(蓝线)及由净热通量计算的对应融化量(红线)
净热通量由辐射通量(短波辐射、长波辐射)和湍流热通量(感热通量、潜热通量)组成(图7)。规定热
Figure 7. Time series of surface net heat flux, short wave radiation, long wave radiation, sensible heat flux and latent heat flux in Greenland from 1979 to 2016
图7. 格陵兰岛1979~2016年地表净热通量及短波辐射、长波辐射、感热通量、潜热通量变化的时间序列
通量向下为正，一般净短波辐射向下，净长波辐射方向向上。1979~2016年格陵兰岛上的净热通量的变化与短波辐射的变化相似，感热通量与净长波辐射变化是负相关关系，潜热通量相比其他三个分量，量值和变化可以忽略不计。
GrIS表现出的质量损失是由于冬季净增加质量(结冰量)小于夏季净减少质量(融化量)形成的，而导致质量减少的速度减缓的可能原因有两个，分别是夏季融化量减少和冬季结冰量增加。
进一步对夏季融化时期的净热通量的各个分量进行分析，结果表明由MAR模式数据直接给出的融化质量(图8蓝线)与通过表面能量平衡方法用净热通量计算的能量对应融化质量(图8红线)的时间序列变化一致。图8(a)为净热通量及其各分量在夏季对融化的能量贡献的年际变化，夏季净热通量主要被净短波辐射控制，净短波辐射的变化趋势的改变控制了夏季融化质量在2012年以后的减少。而GrIS夏季的湍流热通量(潜热、感热)并没有明显的变化趋势。图8(b)显示的冬季净热通量的分量中主要变化部分是长波辐射和感热通量。净热通量及融化质量都没有明显的变化，冬季短波辐射接近于零。无论在冬季还是夏季，潜热通量对冰盖质量改变的能量贡献都很小。
Figure 8. Time series of ice sheet melting in Greenland in summer (a) and winter (b) based on net heat flux and its components (net short wave radiation, net long wave radiation, sensible heat flux and latent heat flux) and MAR model melt production
图8. (a)夏季和(b)冬季的格陵兰岛上根据净热通量及各分量(净短波辐射、净长波辐射、感热通量、潜热通量)计算的对应冰盖融化量及MAR模式给出的融化量的变化趋势
3.4. 格陵兰冰盖表面能量平衡的控制因素
结合上述分析可知，GrIS自2013年开始的质量损失速度减缓的现象主要是由于夏季融化质量的减少。GrIS的融化取决于表面热通量(能量平衡)，能量平衡又由表面特征属性以及气象条件控制，例如地表反照率、云量、近地面气温等，冰盖能量平衡过程中，各项表现是密切联系的。根据上一节的分析，夏季融化主要受到净短波辐射的控制。而净短波辐射由向下短波辐射和向上短波辐射组成，向下短波辐射主要受云量控制，向上短波辐射主要受表面反照率控制。图9展示的夏季冰盖上空云量及反照率时间序列都具有自2013年前后增加的趋势。当云量开始回升后，夏季向下的辐射通量减少，融化减弱，进而反照率也降低。
Figure 9. Trends of GrIS cloud cover (blue line) and albedo (red line) in summer
图9. GrIS夏季上空云量覆盖率(蓝线)及反照率(红线)变化趋势(数据来源于MAR)
4. 格陵兰冰盖质量变化动力机制的初步研究
由GRACE卫星数据展现出的冰盖表面融化的记录与站点观测的近地面空气温度变化一致，并常常作为冰盖融化研究的替代系数(Tedesco等，2011)。在格陵兰上有72个天气站点，其中的68个站点有从1958年开始的温度观测数据记录，剩下四个站点(Upernavik、Ilulissat、Nuuk、Tasiilaq)有百年尺度的长时间记录。尽管大部分站点位于海岸附近，站点的温度观测序列在十年和多年代际尺度上都有一致的变化，冰盖近地面温度升高时，质量减少的速度增大，反之对应(图10)。本文使用四个百年尺度夏季(JJA)近地面温度的平均作为夏季格陵兰近地面温度时间序列(GrT)。
Figure 10. Summer time series and mean values of near-surface temperature in Greenland; changes in summer melting of GrIS (black line)
图10. 夏季格陵兰站点温度及冰盖融化的变化趋势
Figure 11. Empirical mode decomposition (EMD) components (IMF) of GrT in summer of 1910-2016
图11. 1910~2016年夏季GrT的经验模态分解(EMD)各分量(IMF)
将GrT进行经验模态分解得到6个近乎正交的IMF分量(图11)，代表着GrT不同时间尺度的变化及周期，其中IMF6为长期的变暖趋势，IMF1~2为年代际时间尺度信号，IMF3~5为十年际和多年代际时间尺度信号。将GrT分解后的各分量投影到SLP和SST上的空间型呈对应关系(图12)，这是由于SST与SLP本身也有密切的联系。其中SLP的空间型表现的特征更为明显，IMF2~6的SLP投影都大体呈现出格陵兰岛与欧洲大陆西北部的两极结构，该结构与夏季北大西洋涛动类似。因此，夏季北大西洋涛动控制的气压环境(大气环流)改变可能对冰盖融化状况产生影响。
Figure 12. The spatial pattern of IMF components of GrT in summer that projected to sea level pressure field (SLP) and sea surface temperature field (SST) from 1982 to 2016, respectively
图12. 夏季GrT的6个IMF分量分别投影到1982~2016年的海平面气压场(SLP)和海表面温度场(SST)
5. 结语
本文基于GRACE重力卫星数据、MAR区域模式数据、DMI天气站点数据和ERA全球再分析数据对GrIS冰盖质量变化的时空变化特征进行分析，研究了控制GrIS质量变化的主要因素并对主要的动力过程进行初步探讨。得到的结论为GrIS在2003~2012年质量加速损失，从2013年开始，质量损失趋于平缓，主要减缓的区域在GrIS东南部。2013年开始的GrIS减缓事件的原因是云量增多、反照率增强，净短波辐射减少进而夏季融化减少，这可能与夏季北大西洋涛动改变的气压场有关。
研究过程中由于数据的限制，忽略了D对冰盖的影响，尽管D与表面融化量有关，但二者关系并非一成不变。未来的研究中除了进一步研究夏季北大西洋涛动的具体动力过程，还可以深入研究D对冰盖的作用。