积雪密度演变对北极积雪深度模拟的影响

尹豪，苏洁*，Bin Cheng

(1.中国海洋大学 物理海洋教育部重点实验室，山东 青岛 266100；2.中国高校极地联合研究中心，北京 100875；3.芬兰气象研究所，芬兰 赫尔辛基 00101)

1 引言

积雪的热力学变化在极区大气-冰雪-海洋耦合系统中起到至关重要的作用。近年来北极海冰呈现出冰厚变薄[1]、冰龄减小的趋势[2]。冰上积雪的年际变化趋势体现为融化提前[3-4]、积雪深度减小[5]。与20 世纪苏联漂流冰站数据相比，2009-2013 年的平均积雪深度呈现出减小的趋势：西北冰洋平均积雪深度从（35±9）cm 减小到（22±2）cm；而波弗特海、楚科奇海域平均积雪深度也由（33±9）cm 减小到（15±12）cm[5]。气候态意义下，整个北极冰上积雪的季节变化大体相似：秋季受气旋事件影响频繁，冰面积雪深度快速增加；冬季的大部分时间里降雪速率较低[6]。在空间分布上，以北冰洋大西洋扇区年降水、降雪量最大[7]。春季积雪深度与秋季海冰开始冻结时间呈负相关[5]。夏季海冰的衰退会导致上层海洋吸收短波辐射增加，增暖，进而使得秋季海冰生成时间推迟，更多的降雪直接落入开阔海洋中，导致冰上积雪总量减小。因此，翌年春季海冰将更早地暴露在阳光直射下，表面反照率进一步降低，强化反照率正反馈机制[8]。

积雪所具有的高反照率，可以有效减少表面对短波辐射的吸收，改变表面热力学平衡。此外，雪的热传导系数较低，在寒冷的天气条件下阻碍海冰与大气之间的热交换[9]，减慢海冰的增长。而在某些情况下，积雪也有助于海冰的质量增长。当冰上积雪重量超过海冰浮力时，会出现部分积雪被淹没生成雪泥，特定条件下形成雪冰[10]；春季积雪的融水可以重新冻结为层化冰[11-12]，这种积雪向冰的单向转化近年来在北极大西洋扇区也被发现[13-15]。在气候变暖的条件下，逐渐变薄的海冰对冰上积雪的变化更加敏感[16]，准确的冰上积雪量计算对海冰热力过程的模拟和预测至关重要。

与积雪质量密切相关的3 个参数为积雪深度、雪水当量和雪密度。积雪深度指积雪的厚度，由浮标或现场观测测量得到，同时也是海冰热力学模式的重要变量；雪水当量指积雪全部融化成水以后水层的厚度。积雪来源于大气固体降水，其沉积后质量变化的途径多种多样，如：融化、升华、风吹雪重分布以及雪向冰的单向转化等。此外，积雪随时间发生形态变化，从而导致雪密度发生改变，积雪应力变形的致密化率随雪的温度、密度和液态水含量呈指数变化[17]。在质量守恒的条件下，积雪的密度变化会导致积雪深度的变化（积雪压实增密，深度减小）。基于雪的形态特征，Anderson[18]建立了过程较为复杂的雪积累和消融一维模型（SNOW17）。模型针对积雪密度的处理中，除将3 种积雪变质（等温变质、等温度梯度变质和融化变质）的物理过程进行参数化外，还一并考虑了上层雪重量对下层雪的压实。雪密度随时间推移发生变化，这种效应在降雪发生后几周内尤为显著。Brun 等[19]在模型中引入了垂向水通量来描述液态水在雪层中的渗流，还考虑了与上层雪重量以及温度有关的积雪密度演变机制。SNOWPACK 模式中将积雪描述为三相共存的多孔介质，着重描述各相的体积含量以及如颗粒半径等的微观结构参数，且这些微观结构参数在温度梯度变质、等温变质作用下随时间变化[20]。Boone[21]对SNOW17 模型进行了简化，令积雪密度随时间的变化仅为应力变形压实与等温变质这两项之和。基于COSIMA（Coupled Snowpack and Ice Surface Energy and Mass Balance Model）模型[22]，Sauter等[23]将雪密度描述为冰密度、液态水密度以及空气密度基于每层中各自占比的加权平均，并使用了Boone[21]的密度变化参数方案。Liston 等[24]的Snow-Model-LG（Lagrangian snow-evolution model）模型中，考虑了压实、风、雪层中液态水的重新冻结以及层间水汽的转移对雪层密度的影响。

区别于这些复杂的积雪模式，海冰热力学模式中对积雪的考虑以其绝热效果及雪向冰的单向物质转化为主[25-27]。这些模式或将积雪密度假设为常量，或参照SNOW17 对积雪密度进行参数化。对短期过程研究，定常的积雪密度合理可行；对季节变化的研究，考虑随时间变化的雪密度方案的确能完善积雪热力学参数，提高模拟冰厚的准确性[12]。但在这些模式的计算中，随时间变化的雪密度只用来计算雪的热传导系数[12]，仅通过改变雪内、冰雪间热传导来影响冰雪热力学过程，而没有考虑在质量守恒的条件下积雪深度对雪密度变化的响应。在模式计算中，通常使用观测或大气再分析数据的降水/降雪数据作为外部强迫。在以往的一维高分辨率冰雪热力学模式HIGHTSI 中，基于Huwald 等[28]提出的常数块体雪密度均值(320～340 kg/m3)，将外强迫降雪雪水当量转化为积雪深度，这种假设可以减少降雪的不确定性以及雪密度变化对积雪深度的影响，但也忽略了不同次降雪对积雪产生的层化效应[12,27,29]。而观测数据显示，降雪密度在不同地区、不同天气条件下在50 kg/m3至250 kg/m3范围内变化。在Anderson[18]提出的降雪密度参数化方案中，降雪密度随气温变化，其数值限定在50～169 kg/m3。实际降雪密度总体来说要比模式的平均块体雪密度取值（320～340 kg/m3）小得多，后者实际上对应的是降雪以较小的密度沉降到现有雪层之后，经风吹雪、变质、应力变形等作用压实后的集成效果。

在以上研究的基础上，本文对HIGHTSI 中的两个物理过程进行改进：（1）降雪取更接近实际的密度，降落沉积后压实增密；（2）在质量守恒条件下考虑积雪深度对密度变化的响应。本文参照先进的拉格朗日冰上积雪模型（SnowModel-LG）[24]中雪密度随降雪事件、时间变化及其在质量守恒下对积雪深度产生影响的算法，将降雪以较实际的密度沉降到现有雪层，而后随时间推移压实致密的过程进行简化，引入到HIGHTSI 模式中。

2 数据来源

2.1 冰质量平衡浮标数据

冰质量平衡浮标（IMB）来自美国寒带研究与工程实验室的质量平衡浮标项目（http://imb-crrel-dartmou.org/），自2000年起，大量放置于北极冰区，如波弗特海、楚科奇海、拉普捷夫海等海域，用于测量、记录海冰质量平衡中的热力学变化[30-32]。IMB 可以提供近表面温度及气压、海冰厚度、冰上积雪深度以及垂直方向上从空气到积雪到海冰以及冰下海水的温度剖面。IMB 由3 部分组成：(1)配有气压计的数据控制与传输模块；(2)热敏电阻温度链，每10 cm 安装一个热敏电阻温度传感器，精度为0.1°C；(3)一根冻结在冰中的长杆，顶部安装有气温探针，上下两端分别安装两个不同的声学测距仪。长杆上端，安装一个专门用来测量积雪深度的声波测距传感器（sonic ranging sensor，型号为SR50A），俯视雪面并测量仪器与雪面之间的距离，以得到积雪深度变化以及融化季节雪、冰表面消融信息。长杆低端，安装有水下声呐测高计（underwater sonar altimeter，型号为PSA-916），仰视海冰底部，测量仪器与冰底之间的距离，以得到海冰底部生长、消融的信息。

为针对积雪密度变化对积雪深度的影响进行定量分析，本文只选取降雪积累时间段进行研究：对于9 月以前就已经放置了的浮标，以9 月初作为模拟的初始状态；而对于其他9 月后才放置、开始工作的浮标，模拟浮标记录始末所有时段，但仅对其中积累期的积雪深度结果进行统计计算。我们选取了降雪积累期冰面有积雪时段的15 个浮标，浮标漂移的轨迹如图1所示。浮标漂移大多位于北极中心区、波弗特海、楚科奇海以及格陵兰海海域。2010 年以前浮标时间间隔为2 h，2010 年及之后的浮标记录数据时间间隔为1 h。我们将GPS 经纬度位置数据统一插值为1 h 时间间隔，并据此来提取大气强迫数据。对浮标实测积雪深度时间序列进行质量控制，剔除异常点：对每个浮标的积雪深度序列进行24 h 滑动平均，与滑动平均值差大于0.03 m 的积雪深度数据点被剔除[33]。15 个浮标的模拟起止时间、位置、积雪深度信息见表1。

表1 15 个浮标模拟起止日期、位置、积雪深度Table 1 The starting and ending data,position,and snow depth of 15 buoys during the simulation period

图1 本文模拟时段内15 个IMB 的漂移轨迹Fig.1 Drift trajectories of 15 IMB during the simulation period

2.2 ERA-Interim 再分析数据

使用欧洲中期天气预报中心提供的再分析数据（ERA-Interim）作为外部大气强迫，包括风速、2 m 的气温和露点温度、总降水量、降雪量以及云量数据。该数据覆盖全球，本文使用的空间分辨率为0.5°×0.5°，时间分辨率为3 h。

除总降水、降雪量外的大气变量均为瞬时值，使用线性插值处理为1 h 的时间间隔。ERA-Interim 提供的总降水、降雪量为累计值。本文下载使用了3 h总降水、降雪累计值，线性内插平均分配到每小时。依据浮标时间、经纬度位置数据提取对应时间点位置临近的再分析数据格点作为该时刻大气强迫。

2.3 气候态积雪深度数据

基于1954-1991 年苏联冰漂移站长时间序列海冰表面积雪深度观测资料，Warren 等[34]通过积雪深度与空间位置的关系构建了含有6 个系数的二维二次多项式函数（以下简称W99）来描述北极多年冰的气候态冰上积雪深度分布。W99 文中给出了12 个月各自用来计算积雪深度与积雪密度的二次多项式的6 个系数、均方根差、线性趋势、年际变化与不确定度。给出某月北极地区的经纬度，就可根据其给出的二次多项式函数计算给定地点的积雪深度与密度。通过与W99 气候态积雪深度进行对比，简单地对模拟试验结果做出评价，这部分内容于4.2 节呈现。

3 模型介绍与试验设置

3.1 HIGHTSI 模式简单介绍

HIGHTSI 是一个一维冰雪热力学模式，通过数值求解冰雪内部温度变化的偏微分热传导方程，以得到冰雪随时间变化的垂直温度剖面，进一步计算冰雪质量平衡。在冰雪表面，考虑包括大气层结影响在内的热量收支平衡与融化过程；在冰雪交界面，温度连续，考虑雪向冰的单向转化；海冰底部考虑冰水间的热交换以及冻结、融化过程[25,27]。本小节对HIGHTSI 中的积雪过程做简单介绍。原HIGHTSI 模式中，积雪深度为

式中，hs,ini为初始积雪深度；hsf,tot为来自大气强迫的降雪总深度（总为正）；hsm,tot为总融雪量（考虑雪冰时包括雪泥的生成），总为负。

式中，SFtot为总大气强迫降雪的雪水当量；ρw为水的密度；ρSF为降雪密度，取常数。模式中积雪融化包括两部分：表面融化与内部融化。当表面温度大于融点时，使用表面吸收短波辐射、净长波辐射、感热潜热以及传导热提供相变潜热计算表面融化；当内部雪层达到融化温度时，使用上下雪层间传导热以及穿透雪内部的短波辐射提供相变潜热计算内部融化。

3.2 算法改进：积雪深度响应分层压实密度变化

拉格朗日冰上积雪模型（SnowModel-LG）考虑了风吹雪、雪升华等物理过程以及多层雪的密度、粒径以及雪热导率的季节演变[24]。由于同一次风暴事件中降雪的性质是相似的，该模式将雪的垂向结构依据风暴事件分层，每一层携带诸多物理参数，包括雪层的雪水当量、雪密度、积雪深度、粒径大小、导热率等，这些参数都随时间变化。

本文将其质量守恒下积雪深度响应分层压实密度变化的算法应用到HIGHTSI 中，考虑到模型的匹配与计算效率，对层设置进行简化。将集成的总积雪分为：（1）A 层，来自模型模拟初始即存在的积雪，其深度由IMB 实测或模式初值假设给定，初始密度设为330 kg/m3；（2）B 层,该层全部来自模拟时段内的大气强迫降雪，其深度由与该层雪水当量相对应的总雪层质量除以密度计算得到。分别追踪这两层雪水当量、雪密度以及作为计算结果的积雪深度随时间的变化。两层积雪深度的累计值为总积雪深度，具体计算过程如下：

对于雪层A 来说，只需要考虑融化造成的质量减少，而其他造成质量减少的效应，如升华等暂不予以讨论。t时刻雪层A 的积雪深度为

式中，下角标s 代表积雪深度；swe代表雪水当量；hswe,A(t-1)为t-1时刻结束时雪层A 的雪水当量；MA(t)为t时刻雪层A 的融化量的水当量（总为负值）；ρA(t) 为t时刻雪层A 的密度；ρw为水的密度。而在密度随时间变化的试验设置中，ρA初值为 330 kg/m3，随Anderson[18]提出的重力变形压实公式变化

式中，负荷雪重量Ws,A=hswe,A(t-1)+MA(t)+(t)+MB(t)=hswe,A(t)+hswe,B(t)；Ta为气温；Tf为冰点温度；常数C1=7 m-1·h-1，C2=21 g/cm3。

而对于雪层B 来说，除考虑融化外，还要考虑降雪质量增加。t时刻内，在计算雪层B 融化之前，将t时刻来自大气强迫的降雪量S F(t)（雪水当量）添加到雪层B 中。则考虑降雪后（角标′表示）雪层B 的雪水当量为

雪层B 积雪深度为

式中，ρsn为 降雪密度；雪层B 密度为

之后考虑融化并更新雪层积雪深度

总积雪深度为雪层A、B 在t时刻积雪深度之和，即式（3）与式（8）之和。

与原始HIGHTSI 相比，不同之处在于分别追踪了两层雪（A 与B）密度随时间的变化，且在质量守恒条件下考虑了积雪深度对随时间变化的雪密度的响应。

3.3 试验设置

表2给出了针对表1中15 个IMB 在降雪积累期内进行模拟试验的基本设置。试验T1 为标准试验，为原HIGHTSI 设置[12,27,29]，采用常数雪密度（330 kg/m3）将降雪的雪水当量转化成积雪深度；同样基于原模式设置，试验T2 采用更接近现实的降雪密度（200 kg/m3）将降雪的雪水当量转换为积雪深度。试验T1 与T2中B 层新雪的密度是固定值，初始积雪（雪层A）的密度是随时间变化的(公式（4）)，但如前文所述，密度变化只作用于积雪导热系数的改变，并没有考虑质量守恒下积雪深度对密度变化的响应。试验T3 使用本文3.2 节所述算法，降雪密度同样取200 kg/m3，积累后随时间致密的过程作用在A、B 两层，即A 和B 层均的雪密度均发生变化，且考虑了质量守恒下积雪深度对密度变化的响应。简单起见，试验暂未考虑雪冰和层状冰的计算。

表2 HIGHTSI 雪密度数值试验的设置Table 2 Configurations of HIGHTSI numerical tests

本文采用ERA-Interim 再分析数据提供的累计降雪作为模式降雪输入，处理方法见2.2 节所述，初始积雪深度和冰厚由浮标数据提供。

4 试验结果及分析

4.1 积雪深度模拟结果与IMB 的总体比较

考虑到IMB 观测取样的时间间隔，我们仅取与观测值相同时刻的模拟值进行比较，综合统计15 个浮标降雪积累期模拟结果与实测的偏差。对于模拟时段经历暖季融化的浮标，本文的统计计算只针对融化开始前积累期的积雪深度的时段，即图2中初始至灰色竖线之间的时间段。研究时段内总计14 948 个有效观测积雪深度数据。与观测相比，3 个试验模拟结果总体高估：分别有62.4%（T1）、85.6%（T2）以及67.5%（T3）的模拟积雪深度数据点相对实测偏高。试验T1 模拟积雪深度平均偏差最小，小于3 cm，均方根误差为8 cm，相关系数为0.74；Wang 等[12]在与T1 相同的密度设置下，模拟了斯瓦尔巴特峡湾固定冰生长情况，结果显示相对实测有正偏差2 cm，均方根误差6 cm、相关系数为0.83，本文试验T1 统计结果与之相比，偏差在可接受的范围之内。试验T3 积雪深度模拟结果的平均偏差与试验T1 接近，仅比试验T1 高出1 cm，但标准差在3 个试验中最小；平均绝对差与均方根差也与试验T1 的结果十分相近；而试验T2 相对实测的平均偏差与标准差都是最大的（表3）。这说明虽然新雪降落时取较小密度值的处理与实际观测更为符合，但若不考虑降雪后的雪的压实效应，一直以较小的降雪密度积累会造成积雪深度过大。

图2 15 个浮标轨迹积累期HIGHTSI 各试验模拟积雪深度与IMB 实测积雪深度（黑色散点）对比（灰色竖线表示浮标实测融化开始时间）Fig.2 Comparison of modelled snow depth against measured (black scattered points) along 15 IMB trajectories during the accumulation period (the measured melt onset is represented by the gray normal line)

各个浮标的模拟试验与观测的积雪深度随时间变化对比如图2所示。浮标积累期内ERA-Interim 大气强迫降雪量由图3给出，图2中模拟积雪深度快速增长与图3强降雪时期相对应，不同降雪密度设置影响单次降雪事件中模拟积雪深度的增幅。试验T2 和T3 使用较低降雪密度，单次降雪事件中模拟积雪深度的增幅较试验T1（较高新雪密度）更大。总体上，模拟能反映出降雪积累的季节变化过程。但多数情况下模拟的值偏高，这与表3统计结果一致。在3 个试验中，试验T1 与T3 模拟的积雪深度变化相近，而T2 模拟积雪厚度则明显高于两者。2009F、2012G、2013G 以及2013I 这4 个浮标试验T1 和T3 模拟结果与实测十分相似；2010E、2010F、2011I、2012L、2013A、2013H 以及2014F 共7 个浮标的试验T1 和T3 在模拟中后期有或多或少的高估，但明显优于试验T2 的结果；只有2010A 与2014E 这2 个浮标的模拟试验中T2 表现较好，T1 和T3 的模拟结果没有追赶上实测9 月中旬（前者）及10 月初（后者）实测积雪深度的快速增长；2012I 实测积雪深度约处在试验T2 与T1 及T3 模拟中间的平均位置，11 月初以前T1/T3 试验模拟结果更好，11 月初实测积雪深度快速增长，这之后T2 试验模拟结果与实测更为相近；2013B 的模拟3 个试验都与实测不一致，前期低估、后期高估：11 月初模拟试验中由ERA-Interim 再分析降雪数据（图3j）提供了大量降雪，但在浮标轨迹上的观测积雪深度不增反降，除降雪强迫的不准确性以外，也可能是由风吹雪等其他动力原因造成的。

图3 积雪积累期15 个浮标轨迹上ERA-Interim 再分析降雪量强迫数据（灰色竖线表示浮标实测融化开始时间）Fig.3 Snowfall as model input from ERA-Interim reanalysis data along 15 IMB trajectories during the accumulation period (the measured melt onset is represented by the gray normal line)

表3 15 个浮标实测积雪深度平均值及各试验的偏差统计结果比较Table 3 Comparison of statistical snow depth results of 3 tests against the measurement

风对冰上积雪的影响多种多样，包括风吹雪重分布、风吹雪升华以及风压实等。HIGHTSI 中没有风吹雪相关的参数化方案，以往使用该模式的研究中采用减少一定比例固态降雪外强迫的方式考虑风对降雪积累的影响[12,35]，本文未做相似的处理，也就是说未将风对雪的作用加以考虑。Snowmodel 中认为，当平均风速大于5 m/s 时，风就足以对积雪产生影响[24]。我们分别选取3 个浮标中ERA-Interim 平均风速相近且大于该阈值的一段时间进行进一步分析：（1）2012I浮标10 月中下旬至11 月初，这段时间内平均风速为（5.8±1.9）m/s，实测积雪深度从约0.27 m 快速上升至约0.4 m，同时段内ERA-Interim 提供的降雪量完全不足以重现这段时间内实测积雪深度的快速增长，因此推测很可能是由于风将附近的积雪吹过来造成的；（2）2012L 浮标10 月至11 月中旬这一段时间内，平均风速为（5.9±2.4）m/s，模拟积雪深度随强迫降雪量增加而不断增加，而实测在10 月份不增反减，这也可能是由风引起的；（3）2013I 浮标模拟初始至11 月这段时间内，强迫风速为（5.4±2.4）m/s，模拟与实测具有良好的一致性。理论上风对浮标中这3 段时间内积雪的影响相似，然而实测积雪深度相对不考虑风作用的模拟却呈现出了3 种截然不同的情况（模拟低估、模拟高估以及模拟与实测相近）。基于此，我们可以推测：在积雪模拟中，风作用的不确定性非常强。

Li 等[36]通过浮标实测冰雪厚度计算了103 个IMB 的干舷，结果显示其中仅有7 个浮标记录有负干舷存在，本文所使用的浮标2012I 就是其中之一。实际上，干舷为负为冰面被淹没的必要不充分条件。很多其他过程都会导致冰雪界面淹没，如卤水通过多孔海冰中相互连接的卤水通道垂向渗透，或海水通过变形冰脊断裂处到达积雪底部等[37]。HIGHTSI 模式考虑了雪冰以及层状冰的生成，认为当干舷为负时，淹没雪泥形成，冻结条件下转化成雪冰[12,38]，这里针对雪冰补充对浮标2012I 的模拟计算。结果显示，试验T1 中约5 cm 的积雪转化成了7 cm 的雪冰；而T3 中约4 cm 的积雪转化成了6 cm 的雪冰。其余浮标积雪较薄、海冰较厚，没有雪冰生成。

试验T3 模拟积雪深度相较T2 减小了，这是雪层A、B 积雪深度响应密度变化的结果。图4给出了T3 试验中A、B 层密度随时间演变。根据我们的试验设置，雪层A 密度只受到压实增密（公式（4））的影响，于积累期时间段内，A 层雪水当量不变，密度不断增大，因此，根据质量守恒，该层积雪深度不断减小。对于B 层，密度除随公式（4）变化之外，在发生降雪的时间点还会受到降雪密度加权平均效应（公式（7））的影响。公式（7）本质是分别以B 层原有雪与降雪的体积为权重对密度加权平均。由于降雪的密度更小，因此强降雪时段该层的密度会发生一定的减小，图4中B 层密度在总体增加的趋势之上叠加了锯齿形变化即是该效应的体现。其中，总体积雪密度增加的趋势是公式（4）压实增密的效应；局部强降雪时期呈现的短期下降过程是公式（7）的效应，曲线谷区宽度及下降幅度取决于降雪的量（图3）与持续时间。在积累期，由于外强迫降雪（雪水当量）不断增加的同时其雪层密度也不断增加，B 层的积雪深度总体表现为缓慢增长，其增幅明显小于T2。2010E、2012G、2012L、2014F 共4 个浮标的A 层积雪在模拟结束前就已融尽，因此A 层雪密度的时间序列较B 层短。对于这些模拟过程中经历暖季的浮标的B 雪层发生融化时密度快速下降，是快速融化与新雪降落加权平均的共同结果。

在试验T3 中，冰上积雪块体密度（图4黑线）为A、B 层总积雪质量与积雪深度的比值。初值均为330 kg/m3，即为预设的初始即存在的雪层A 的密度，这时尚没有新的降雪。初始积雪深度很小的浮标（如2012L、2013A、2013B、2013G、2014E 与2014F），块体密度受B 层密度变化的影响显著—因为块体密度这个量具有加权平均的意义，相较不断增长的B 层，A 层（初始即存在的积雪）的质量很小，故对雪层整体密度影响很小。无融化积累期时段内，每个浮标的总积雪质量（初始加大气强迫降雪量）在3 个数值试验中都是相同的，因此3 个试验中积雪深度的差异主要来自于块体密度的差异。2009F、2010A、2012I等浮标试验T1、T3 模拟积雪深度十分相近，几乎没有差异，这是因为试验T1 的块体密度始终为330 kg/m3，而这几个浮标在试验T3 中的块体密度十分接近330 kg/m3，于该值上下小幅度浮动（块体密度标准差小）。而对于那些T1、T3 模拟积雪深度差异较大的浮标，其初始积雪相对后续降雪较少，块体密度受B 层大气强迫降雪层密度变化影响较大（块体密度标准差也较大）。粗略统计15 个浮标试验T3 块体密度平均值约为305 kg/m3，这也是试验T1 与T3 15 个浮标平均意义下表现相近的原因。

在积累期，由于外强迫降雪（雪水当量）不断增加的同时其雪层密度也不断增加，试验T3 的B 层的积雪深度总体表现为缓慢增长，其增幅明显小于T2 的结果。可见，T3 显著地改善了T2 积雪深度过度积累的问题，这是积雪深度变化响应密度变化的结果。

我们进一步采用泰勒图（图5）再现积累期试验T2 和T3 积雪深度时间序列的标准差、相对于各浮标观测数据的均方根误差与相关系数。为将各浮标的试验误差整体呈现于一幅泰勒图中，我们首先以各观测点（图5中OBS 点，标准差=1，均方根误差=0，自相关系数=1）为基准，将各浮标模拟积雪深度序列均除以观测序列标准差作归一化处理。通过泰勒图我们可以直观地比较各试验模拟积雪深度与实测的符合程度。由于试验T3 与T1 的模拟结果相近，在图5中仅给出了T2 与T3 相对于观测数据的误差统计结果，以体现质量守恒条件下考虑压实致密效应的积雪深度算法（T3）对过度积累问题的明显改善。图中浮标序号与表4对应。由图5可见，除序号2（2010A）、7（2012I）、14（2014E）3 个浮标外，T3 的模拟积雪深度仅就标准差与均方根误差这两项指标而言，较T2 结果均有明显改善，虽然相关系数改变不大。如T2 均方根误差大于2 cm 的几个浮标试验中，以第3 个浮标（2010E）最为明显，T3 相对T2 均方根误差减小了5 cm，为T2 均方根误差的47.2%。而对于这3 个浮标来说，本身均方根误差和标准差都是这15 个浮标中较小的，试验T3 模拟结果相较T2 标准差减小了，且两个试验均方根误差相差在2 cm 以内。

雪积累期各浮标实测积雪深度平均值及标准差、模拟积雪深度以及误差统计结果详见表4。试验T1 大多数时间模拟结果高于实测，有6 个浮标（2010A、2012G、2012I、2013B、2013G、2014E）例外，有负偏差，其中2012G、2012I 与2013B 试验T1 与T3 模拟积雪深度与实测误差都在3 cm 以下。综合平均差、标准差、均方根差及绝对误差，15 个浮标中只有2010A与2014E 的试验T2 表现更好（这与图2呈现的一致），原因可能是强迫降雪量不足或实际降雪密度更小；对于其他浮标试验T2 持续以较小的降雪密度累计会造成积雪深度过大。应用考虑致密效应的积雪深度算法（试验T3）可以较大程度地避免这个问题，平均积雪深度减小了7 cm，误差由实测平均值的55.6%（T2）减小至16.4%（T3）。T3 相对T1 的相差仅为1 cm，模拟结果也较为相近。实际上，T1 所使用的平均雪块体密度本身也可以看作降雪沉积后经压实增密的时间平均结果，但从物理过程的处理上，显然T3 更为合理。

表4 雪积累期各浮标实测积雪深度平均值及标准差、模拟积雪深度以及误差统计结果Table 4 Statistical results including mean,standard deviation,and error of 3 tests against observation during the accumulation period

总的来说，就这15 个浮标的雪积累期而言，试验T3 所采用的改进算法不仅通过可变积雪密度的处理，避免了T2 由于持续以较小降雪密度积累造成的积雪深度高估的问题，又在标准试验T1 的基础上加入了雪压实致密对积雪深度影响的过程。

4.2 模拟结果与W99 气候态积雪深度数据的比较

3 个数值试验模拟积雪深度峰值出现的时间相近，挑选15 个浮标模拟积雪深度最大值出现的月份与其对应的经纬度，将3 个试验模拟结果、浮标实测以及W99 气候态积雪深度模拟结果对比。之所以选择模拟最大值，并非为了选择一个常规认识下的最大积雪深度数据去与气候态模拟结果对比，而是优先保证浮标轨迹上融化开始之前，积累期内的全部降雪都被纳入考虑。对于不经历融化的浮标，积累期内的全部积雪于模拟结束时全部纳入模型；而对于经历融化的浮标来说，融化开始前的积雪深度为积累期内全部降雪沉积得到的——对于两者来说，模拟积累期内全部降雪等同于模拟时间序列最大值。由表5可见，综合15 个浮标的平均值，试验T2 的积雪深度平均值明显较W99 偏高，试验T1 和T3 平均值与W99 气候态积雪深度平均值十分相近，但其中有3 个浮标（2010F、2012I 与2013A）试验T1 和T3 相对W99 的偏差都超过10 cm。前两个浮标T3 试验相对气候态误差在10 cm左右，而2013A 气候态积雪深度结果远高于模拟与实测结果。该浮标5 月所处位置于格陵兰岛与埃尔斯米尔岛之间，不在冰漂移站轨迹覆盖范围以内。近年来北极冰上积雪深度有负趋势[7]，这里计算的W99 气候态积雪深度结果并未考虑积雪深度趋势变化，可能较实际偏高。

表5 15 个浮标模拟峰值所在月份、模拟峰值时段内 3 个数值试验模拟积雪深度结果、IMB 实测积雪深度结果与 W99 气候态积雪深度对比Table 5 Months when modelled peak occurs,snow depth from 3 numerical tests,and observation results during the peak period against W99 climatology snow depth results

5 结论与讨论

本文以ERA-Interim 再分析数据作为外部大气强迫，基于SnowModel-LG[24]中的考虑压实致密的积雪深度的多雪层计算方法，将其简化为两层雪层，将新降雪以接近实际的较低密度降落，继而随时间压实增密，并在质量守恒条件下影响积雪深度的物理过程加入一维高分辨率热力学冰雪模式HIGHTSI。模拟研究了IMB 轨迹上积雪深度的演变，并将模拟结果与IMB 实测积雪深度进行对比验证。

结果显示，对积雪积累期而言，本文改进算法相较以往[12,27,29]使用了更接近观测的降雪密度，更加贴近现实，且得到了与以往算法相近的结果；避免了采用新雪密度转化降雪外强迫所导致的过度积累，以浮标观测为标准，分层积雪密度压实（试验T3）计算得到的积雪深度相对T2 减小了24.5%，平均绝对误差减小了5 cm。

北冰洋的积雪深度模拟依然是一项挑战。总体上，本文引入的降雪以更符合实际、较小的密度沉降，后经压实增密，雪层变薄的变化过程，对描述非融化期积雪积累过程有一定能力，但目前的模拟还存在若干问题。降雪发生时，强迫数据的不准确是误差的主要来源之一。本文使用的ERA-Interim 提供0.5°×0.5°的网格分辨率，而IMB 测量声波测距传感器与积雪表面之间的距离，这样的测量结果具有强烈的局地特征，且易被雪丘影响，即观测与强迫场的不同尺度带来不确定性。其次，对于经历融化期的浮标，模拟出现了融化过早、融化过快的问题。HIGHTSI 模式中雪物理过程还是较为简单，雪导热系数的参数化是否合理、对雪层中液态水的渗流过程[19]的量化能否对热力学过程有所改进都还有待检验。其他一些影响积雪密度和积雪质量变化的因素还未加以考虑[24]。以往的模式中也对风压实，风吹雪重分布以及风吹雪升华过程进行了参数化[39-40]。此外，风对降雪的积累也有极大的影响。气候态数据分析表明，即使在陆地上，早冬观测降水量与实际降雪积累量之间的相关性也很差，只有不到一半降水对雪的积累做出贡献[41]。风吹雪侵蚀、沉积过程只在1～500 m 水平距离尺度上起作用，因此在非高分辨率的模式中可以忽略，但风吹雪的升华在北极积雪质量收支模拟中有显著作用[24]。风对雪的作用是复杂的，基于4.1 节对风作用的分析，本文暂未考虑风的作用，尚未将风吹雪参数化纳入模式，风的效应还有待进一步研究。

另外，本文主要专注于对积雪深度算法的改进以及对随时间变化的积雪密度对模拟积雪深度的影响进行研究，没有对海冰模拟进行细致的讨论。冰上积雪不仅影响海冰热力学行为，还可以通过冰雪转化形成雪冰、层状冰对海冰生长做出贡献。HIGHTSI 模型中考虑了雪向冰的单向转换，模拟能较好地反映北极雪冰、层状冰形成的基本特征[11-14]。雪冰在南极海冰热力平衡中的重要性已是共识，近年来的实测与模式结果表明，北极雪冰生成潜力增大[14-15,42]。后续工作将会进一步研究雪向冰的单向转化对海冰的影响。

致谢：感谢林龙博士在浮标数据质量控制方面的讨论和指导。