托木尔峰青冰滩72号冰川表碛区夏季消融模拟研究

何 捷， 王璞玉,， 李宏亮， 李忠勤,， 周 平，牟建新， 余凤臣， 戴玉萍

（1.中国科学院西北生态环境资源研究院,冰冻圈科学国家重点实验室,甘肃 兰州 730000；2.中国科学院大学,北京 100049; 3.石河子大学理学院,新疆 石河子 832000）

20 世纪以来，全球变暖导致的冰川加速退缩已引起了广泛关注，山地冰川作为中亚干旱内陆区域的重要淡水来源，其消融变化与融水径流显著影响下游的农业和生产生活用水［1-3］。中国境内表碛覆盖型冰川共有1723 条，总面积为12974 km2，表碛所占面积为1494 km2（11.5%）；尤其是在我国天山，表碛覆盖型冰川数量最多，面积也最大，所占全部冰川面积比例达54%，表碛覆盖区域约481 km2，表碛区域占表碛覆盖型冰川面积的12.2%［4］。表碛的存在使得冰/雪表面与大气的能量交换变为表碛表面与大气的能量交换，其消融特征显著区别于裸冰或积雪表面。大量研究表明［5-7］，表碛厚度小于某一阈值时，由于表碛表面的反照率小于雪冰表面致使表碛吸收更多的太阳辐射，从而加速表碛下覆冰的消融；当表碛厚度大于这一阈值时，表碛表面盈余的能量传导至冰面时消耗过多，进而抑制下覆冰的消融。因此，研究表碛的物理特征（如厚度、导热系数、表面粗糙度）和能量平衡过程对揭示表碛覆盖型冰川的消融状况尤为重要。

表碛覆盖型冰川能量平衡模拟一直被学者们所关注。Nakawo 等［6］基于野外表碛观测数据计算了表碛下覆冰的消融量，得到了与实测消融量较为一致的模拟结果，并给出了通过表碛热阻估算消融量的方法；Nicholson 等［7］使用基于表碛内热传导模块的能量平衡模型，结合日尺度气象数据作为模型驱动数据，应用在阿尔卑斯山和斯瓦尔巴特群岛两条冰川，结果表明模拟的下覆冰消融量较好地反映了实测消融量。此后，表碛覆盖型冰川的模拟研究基于野外实测数据，开展了表碛层垂直结构内的水汽和热传导等更为细致的研究［8-9］。其中，Rounce等［10］通过对表碛导热系数、反照率和表面粗糙度进行了野外测量，在此基础之上模拟了表碛表面温度和下覆冰消融量，分析了模型参数的不确定性对模拟结果的影响程度。国内有关表碛覆盖型冰川的模拟研究起步相对较晚，Han 等［11］基于表碛热传导理论构建了适用于天山科其喀尔冰川表碛状况的能量平衡模型；Yang 等［12］利用表碛能量平衡模型分析了藏东南嘎隆拉冰川的气象和辐射分量特征，与非表碛覆盖型的帕隆4 号冰川相比，嘎隆拉冰川吸收了更多的太阳辐射和下行长波辐射，从而增加了上行长波辐射和湍流热通量的亏损进而导致消融量增多。此外，Zhang等［13-14］构建了基于表碛热阻系数的能量平衡模型，进一步发展了表碛厚度估算模型，并在藏东南表碛覆盖型冰川的分布区域得到了应用。基于上述研究发现，表碛覆盖型冰川的模拟研究主要集中在阿尔卑斯山脉和青藏高原地区，天山以及托木尔峰地区表碛覆盖型冰川的消融研究较少，限制了我们对该区域冰川消融速率和冰川融水资源变化的认识。因此，本文从冰川能量平衡过程出发，对青冰滩72 号冰川进行消融模拟研究，验证能量平衡模型和表碛属性参数在实际模拟中的适应性，完善对该区域表碛覆盖型冰川消融机制的认识。

天山托木尔峰地区作为亚洲内陆最大的冰川发育区，该区蕴藏着丰富的冰雪水资源，托木尔峰山区平均每年产流约63.4×108m3，其中约50%为冰雪融水［15］，对区域水资源和人民生产生活至关重要，同时对该区域冰川的现状和消融研究始终是国内外冰川学的研究热点［16］。鉴于此，本文选择天山托木尔峰青冰滩72号冰川作为研究对象，基于实测气象数据，结合野外表碛温度和消融量观测，使用表碛覆盖型冰川能量平衡模型开展表碛区冰川消融模拟研究，对表碛区的能量平衡特征和热传导状况进行分析，以期揭示表碛覆盖对冰川能量收支和消融过程的影响。

1 研究区概况

托木尔峰地区位于中国境内天山山脉的最西端，是亚洲内陆最大的冰川发育区之一。据中国冰川编目资料［17］，托木尔峰地区在我国境内共有冰川1858 条，总面积达4195 km2，是我国新疆阿克苏地区和伊犁地区重要的水资源，也是塔里木盆地北侧塔里木河主要支流的源区。托木尔峰地区为大陆性干旱气候，冰川区雪线附近年平均气温为-11～-7 ℃，降水主要来自大西洋和北冰洋的暖湿气流，集中在5—9月［18-19］。

青冰滩72 号冰川（41°45.51′N，79°54.43′E）发源于阿克苏河的支流库玛拉克河上游区域，地处托木尔峰南端，是一条朝南的山谷-冰斗冰川，由降雪和冰/雪崩补给（图1）。青冰滩72 号冰川面积为5.61 km2，冰舌末端至冰川顶部的海拔介于3560～5986 m之间，冰川长度约为7.4 km，冰川存在冰崖和冰面湖［20-21］。表碛主要分布于消融区冰舌末端和两侧，表碛覆盖面积为0.87 km2，占消融区面积的52%［22-23］。表碛厚度整体随海拔升高而减少，厚度最大处超过1 m，冰舌东侧和西侧表碛的平均厚度分别为7 cm和16 cm，且西侧的表碛覆盖面积大于东侧。72 号冰川的表碛临界厚度约为4 cm，超过该厚度的表碛区消融量随厚度增加而减少，其面积约为0.66 km2，占表碛总面积的76%；表碛厚度小于4 cm的区域消融量随表碛厚度的增加而变大，面积为0.21 km2［24-25］。与天山山脉其他冰川相比，该冰川运动速度较快，年均水平运动速度为47 m∙a-1［26］。

图1 研究区概况示意图Fig.1 Map of study area

2 数据与方法

2.1 气象数据

本文的气象数据来源于在青冰滩72 号冰川表碛区架设的自动气象站（AWS），型号为Davis Vantage Pro2 Plus。该气象站位于冰川冰舌东侧海拔3950 m 处，该处的表碛厚度为12 cm，主要由灰色和白色花岗岩及碎屑组成（图1）。观测的气象要素包括气温、相对湿度、风速、降水、入射/反射短波辐射等（表1），采集时间间隔为1 h。本研究使用的气象数据时段为2008 年7 月30 日14：00—8 月29 日22：00。由于复杂冰面气象环境的影响，观测期间自动气象站的数据记录在2008 年8 月6 日17：00—7 日12：00和8月8日14：00—16：00期间存在一定的缺失。需要说明的是，由于缺失的时间段内未能进行消融模拟，所以数据缺失时间段的消融量由平均消融量代替。同时，由于辐射传感器在2008 年8 月28 日12：00—8 月29 日22：00 之间无记录，因此，模拟的结束时刻提前了34个小时。

表1 自动气象站气象传感器和指标Tab.1 Meteorological sensors and their indicators of the automatic weather station

2.2 ERA5气象再分析数据

ERA5 是ECMWF 对过去全球气候和天气的第五代再分析产品，其水平分辨率为0.25°×0.25°，时间分辨率为1 h。本文选取了研究期内ERA5 云量因子数据分析不同天气状况下气象和能量条件的变化。不同的天气状况划分参考Van Den Broeke等［27］基于云量因子的划分方法，ERA5 的云量因子＜0.3时划分为晴天，＞0.7时划分为阴天，不考虑降水的影响，划分后通过逐小时的能量分量时间序列变化进行分析。数据来源于哥白尼气候变化服务(Copernicus Climate Change Service) 数据平台(https://cds.climate.copernicus.eu)。

2.3 表碛温度数据

表碛温度数据是表碛覆盖型冰川能量平衡模拟的关键参数和模拟结果的重要验证资料。研究人员于消融期间在冰川海拔3950 m 处的表碛区钻了一个2 m 的温度孔。该温度孔处表碛层厚度为13 cm，在表碛层内共布设了2 个温度探头。其中一个探头置于表碛表面以下约l cm 处，用于测量表碛表面的温度。由于表面温度的观测容易受到太阳辐射的影响，因此使用表碛表面以下1 cm 处的温度来近似代替表面温度，这与Rounce 等［10］测量喜马拉雅山脉Imja-Lhotse Shar 冰川表碛表面处温度的方法相一致。另一个温度探头布设在表碛表面以下10 cm 处，用于测量表碛层内部的温度。表碛以下冰体内，每隔20 cm 布设一个温度探头测量冰温。测温探头均为热敏电阻，由中国科学院西北生态环境资源研究院（原中国科学院寒区旱区环境与工程研究所）研制，精度为±0.05 ℃。野外试验中，人工使用万用表测量探头电阻率，进而计算表碛层和冰体内部的温度。共开展了6次观测，分别为7月31日、8月1日、8月3日、8月4日、8月7日和8月8日。

2.4 冰川实测消融量数据

研究人员在自动气象站附近表碛区布设了1根消融花杆（图1c），表碛厚度为12 cm。在2008 年7月29 日—9 月1 日期间，开展了26 次详细观测。具体的观测内容为花杆至冰川表面的垂直高度。消融量的计算公式为：

式中：A为消融量(m w.e.)；ρi为冰川冰的平均密度(900 kg·m-3)；Δm为观测期间花杆至冰面的垂直高度变化(m)。最终得到研究期内花杆的消融量值，本研究使用该数据作为冰川消融量模拟的验证资料。

2.5 表碛区能量平衡方程

本文使用表碛覆盖型冰川能量平衡模型（Debris Energy-Balance model，简称DEB 模型）进行单点尺度冰川能量和消融模拟，该模型开发较为成熟，物理意义明确，完整描述了表碛区冰川的能量交换过程，与度日模型相比能获取能量分量对消融过程的影响，且作为一款开源模型在阿尔卑斯山脉Miage 冰川和藏东南嘎隆拉冰川取得了很好的模拟效果［12,28］。其中，表碛覆盖型冰川能量平衡模型主要基于如下能量平衡等式：

式中：S为净短波辐射；L为净长波辐射；H为感热通量；LE为潜热通量；P为降水带来的热通量；G为传导热通量，即能量从表碛表面传输至冰面所需的损耗；M为冰川消融能量；单位均为W∙m-2。在本文的模型中，表碛表面温度和上行长波辐射整体带入能量平衡方程求解，表碛表面温度Ts不仅是计算各辐射通量的输入项，同时其数值也取决于辐射通量的大小。除太阳短波辐射和下行长波辐射外，其他能量通量均可以表示为表碛表面温度Ts的函数。该模型使用Crank-Nicolson数值迭代方案，计算时间步长为1 h，Ts既是能量平衡模型中的未知数，又是计算能量项的所需参数，通过将表碛分为数层，迭代计算每层表碛的温度，从而得出冰面的热通量［7,27］。由于表碛之间的缝隙中存在大量的空气和水汽，太阳辐射和湍流热通量在表碛内发生复杂的能量交换，能量平衡模型根据表碛内的温度变化表示能量传递过程。受太阳辐射和长波辐射的影响，表碛表面温度时刻发生变化，昼夜温差大。但根据表碛温度观测和自动气象站气温数据，表碛层温度在消融期内基本都为正温，因而可假定表碛底部与冰体交界处的温度一直处于融点（0 ℃）。与此同时，通过对表碛不同厚度处的温度进行观测，表碛内的温度曲线呈现出非线性变化，昼夜的能量传递方向也可能发生改变。因此，使用非定态热传导方程描述表碛内的温度变化：

式中：α为热扩散系数，表示非定态热传导过程中物体内部温度趋于均匀的能力；T为某一表碛层的温度(K)；z为每一表碛层的深度(m)；t为时间(h)；k为表碛导热系数(W∙m-1∙K-1)，ρd和cd分别为表碛的密度(1496 kg∙m-3)和热容量(948 J∙kg-1∙K-1)。

2.5.1 短波辐射 地表的短波辐射来源于太阳辐射，同时也是冰川表面的主要能量来源。云量对入射短波辐射有重要的影响，表碛反照率很大程度上也决定了出射短波辐射的大小。本文使用自动气象站的入射/出射短波辐射数据作为模型输入。

2.5.2 长波辐射 长波辐射包括地表发射的上行长波辐射和大气发射的下行长波辐射，由于自动气象站处没有长波辐射的实测数据，因此使用斯特藩-玻尔兹曼定律计算下行和上行长波辐射［29-30］。下行长波辐射Lin的计算公式为：

式中：σ为玻尔兹曼常数［5.67×10-8W∙(m2∙K4)-1］；Ta为大气温度(K)；εeff为有效辐射率，即考虑云层覆盖条件下的大气辐射率，采用Idso［31］和Unsworth等［32］的公式进行计算：

式中：εcs为晴空辐射率；ea为海平面水汽压(Pa)；n为云量，来自于ERA5再分析数据。

上行长波辐射Lout的计算同样采用黑体辐射定律，其不确定性主要取决于模拟的表碛表面温度：

式中：εd为表碛的比辐射率；Ts为表碛表面温度(K)。

2.5.3 湍流热通量 低层大气与下垫面之间因温度或湿度梯度产生的能量交换称为湍流热，包括感热通量和潜热通量。湍流热通量是冰川能量平衡的重要组成部分。通过计算整体理查森数和判断冰面稳定状态可以较好的模拟湍流热通量，因此，本文使用整体空气动力学方法计算感热通量，计算所需的表碛表面温度数据同样是模拟值，其基本思路为假设湍流热通量与温度梯度、湿度梯度和风速梯度存在相关关系［33-34］。感热通量是近地面气温与下垫面温度之间的差异引起的热量交换。潜热通量则是近地面大气与下垫面之间水汽相变产生的能量交换，包括蒸发和凝结。感热通量H和潜热通量LE的计算公式为：

式中：ρa为当地空气密度(kg∙m-3)；ca为空气热容量(J∙kg-1∙K-1)；Lv为蒸发潜热(J∙kg-1)，近地面气温为283 K 时其数值为2.476×106J∙kg-1；kvk为卡门常数(Von Karman constant)，是假设混合长和速度廓线的关系所引入的经验系数，实验室和近地面大气测定值在0.35～0.43 之间，本文取值为0.4；u为近地面风速(m∙s-1)；Ta为近地面气温(K)；Ts为表碛表面温度(K)；ea和es分别为近地面大气和表碛表面的水汽压(Pa)；z为自动气象站气温传感器距离地面的高度(2 m)；z0m为表碛表面（动量）的粗糙度(m)，定义为水平风速为0处距地物表面的高度，取值为0.016 m［28］；z0t和z0q分别为表碛表面的热量和水汽通量的粗糙度，二者取值与z0m相同。Φm、Φh和Φv分别为动量、热量和水汽的无量纲稳定性函数。在计算潜热通量的过程中，需定义表碛表面的相对湿度，本文将自动气象站处测得的降水量＞0.2 mm 时的表面相对湿度定义为100%，否则为0%［10］。

Brutsaert 等使用理查森数Ri来描述动量、热量和水汽通量的稳定性状态［33,35］，本文假设自动气象站处风速为0时，三者处于稳定状态，否则将处于非稳定状态。当三者在近地面处于稳定状态时，理查森数为正，可表示为：

当近地面处于非稳定状态时，理查森数为负，此时可表示为：

式中：理查森数Ri使用Brutsaert 描述的方法计算［33,36］，公式如下：

式中：g为重力加速度(m∙s-2)；Tm为气象站测得气温Ta与表碛表面气温Ts的平均值。

2.5.4 降水热量 降水与表碛之间的热量交换P采用Reid等［28］提供的方法计算：

式中：ρw为水的密度(kg∙m-3)；cw为水的比热（4.18×103J∙kg-1∙℃-1）；r为降雨速率（m∙s-1），Tr为降水温度(K)，由于缺乏实地观测，使用气温Ta代替降水温度。

2.5.5 传导热通量 一般将物体或系统内的温度差产生的热量传递现象称为热传导。导热系数k是计算传导热通量的关键参数，它用来描述物体传导热量的固有能力，热量传递的速率取决于温度梯度的大小和表碛的热特性，表碛的导热系数与其颗粒大小、密度、表碛内孔隙湿度和温度密切相关，其数值大小受到周围环境的影响，表碛内的温度和湿度越大，导热系数越大，传递到冰面的能量越多。模型通过对表碛进行等距离分层，每层表碛的厚度为0.01 m，使用傅里叶定律计算表碛表面和表碛/冰面交界处的传导热通量，使得每层表碛的温度变化可以近似表示为线性，具体计算公式如下：

式中：Gice为表碛底部传导至冰面的热通量(W∙m-2)；Td(N-1)和Tf分别为第N层表碛温度和冰体熔点温度（即冰面温度）；h为每层表碛的实际厚度(m)。

2.5.6 冰川消融 能量平衡模型模拟的表碛以下冰融化量近似为冰川消融量。冰川的消融量a由可供冰川消融的能量M求得［37］：

式中：ρi为冰密度(kg∙m-3)；Lf为融化潜热(J∙kg-1)。

3 结果与分析

3.1 模型验证

表碛表面温度反映了表碛与近地面大气之间的能量交换程度，表碛表面的能量主要以热传导的形式向下传递，基于表碛表面温度和表碛内热传导过程可估算出下覆冰的消融量，野外观测的表碛温度数据是模型的直接验证参数。本文选取了能量平衡模型模拟的表面温度和10 cm 深度处温度同实测值进行验证（由于8 月6—8 日有两段时间气象数据存在缺失，数据缺失的时间段内未显示能量平衡的模拟值），结果表明，表碛表面温度（图2a）的模拟值和实测值具有很好的一致性（R2=0.91, RMSE=1.78 ℃），同时表碛下10 cm 温度（图2b）的模拟值也与实测温度值拟合较好（R2=0.60, RMSE=0.48 ℃）。尽管表碛下10 cm 深度处的温度模拟精度有所降低，但从表碛覆盖型冰川的能量交换过程中可知，表碛表面温度通过近地面气象条件计算得到，进而基于表面温度计算表碛内不同深度处的温度，因此，表碛下10 cm 深的温度模拟结果包含了表面温度的计算误差。从表碛温度的时间序列变化中可以看出，表面温度的日波动较大，白天温度一般超过15 ℃，夜晚降温迅速，温度较低；表碛下10 cm 深度靠近冰面，日温差小，白天温度不超过5 ℃，夜晚温度接近0 ℃。基于野外观测的温度数据中，有4 d的正午表面温度超过15 ℃，10 cm 深处的实测值也相应较高；剩余2 d 温度低于10 ℃，10 cm 深处的温度同样较低。通过对气温和表碛表面温度进行相关性检验，发现二者存在较好的相关性（r=0.55，P<0.01），因此，表碛温度变化在一定程度上可以反映近地面气温的时间变化。

图2 模拟与实测的表碛表面温度和表碛下10 cm处温度Fig.2 Simulated and measured debris temperature at surface and depth of 10 cm

利用夏季模拟的消融与实测的消融进行对比，如图3 所示，结果显示，7 月30 日—8 月28 日期间，模拟的冰川累积消融量为0.39 m w.e.，实测的累积消融量为0.40 m w.e.。虽然研究期内，不同时间段的累积消融量模拟值普遍较实测值要小，但整体来看，模拟值和实测值两者之间具有较好的一致性，决定系数（R2）为0.92，均方根误差（RMSE）为±0.03 m w.e.。两者之间存在差异的原因可能主要有三个方面：其一，本文中表碛属性参数并未开展实地观测，是通过率定得到的；其二，模型没有涉及表碛内水汽的蒸发、再冻结过程，模拟结果与真实冰面能量收支状况存在一定差距；其三，该条冰川运动速度较快，存在较强的运动补给。模拟与实测的累积消融量在8 月15—19 日间均呈现出较为平缓的趋势，通过对比气温数据发现这一时段气温整体偏低，最低气温为-3.9 ℃，表明，此时气象条件不利于消融，以致消融亏损减少。尽管消融的观测时间只有一个月左右，但取得的观测结果精度较高。消融观测采用的是花杆法，研究人员每隔1～2 d 观测，充分记录了短时间尺度内消融的详细变化情况，为开展冰川能量消融模拟提供了良好的实测验证数据。

图3 模拟与实测的累积消融量Fig.3 Simulated and measured glacier melt

3.2 表碛区能量平衡特征

青冰滩72 号冰川夏季表碛区能量特征如图4和图5 所示。在冰面的能量交换过程中，净短波辐射是唯一的能量收入项。感热通量是最大的能量支出项（49.7%），其次分别为传导热通量（消融耗热）（25.8%），净长波辐射（19.8%）和潜热通量（4.6%），降水热量不足1%，可忽略其对冰面能量交换的影响。

图4 表碛表面辐射分量的日平均变化Fig.4 Mean diurnal variation of debris surface heat fluxes components

图5 表碛表面能量分量的逐小时变化Fig.5 Hourly variation of debris surface heat fluxes components

观测期内净短波辐射表现出明显的日变化，单日峰值可超过1000 W∙m-2，也可小于400 W∙m-2，这可能受到云量的影响，具体表现为云层的遮挡大幅减小了入射短波辐射透过大气到达冰川表面，使得净短波辐射曲线出现扰动或峰值削减。昼夜变化方面，从日出开始短波辐射持续增大，午后可超过600 W∙m-2，随之短波辐射持续减小直至日落。感热通量是冰面的最大能量支出项，这是由于太阳辐射能量加热表碛使表碛温度显著高于近地面气温而产生热量对流交换。午后感热通量的亏损高达-400 W∙m-2，值得注意的是，由于夜间表碛温度较低，使得感热通量处于轻微的盈余状态。传导热通量、净长波辐射和潜热通量均呈现出热量支出状态，其数值在白昼期间波动不大，曲线变化较为平稳，且亏损量不超过-200 W∙m-2。观测期间降水带来的热量几乎为0，所以，降水热量对冰面能量交换过程中的影响可忽略不计。表碛表面各辐射通量的日均变化特征较为一致，即在夜间能量波动小且趋近于0，日出前后收入和支出量持续增加，在午后绝对值达到最大，随后逐渐减小。

3.3 表碛区消融特征

基于青冰滩72号冰川表碛区的气象数据，结合能量平衡方程计算出模拟期内的日平均消融量为0.55 mm w.e.，研究期内冰川始终处于消融状态，日均最大消融发生在8月3日，消融量为0.88 mm w.e.，日均最小消融发生在8 月17 日，消融量为0.23 mm w.e.，从7 月30 日—8 月28 日，累积消融量呈现出持续减缓的趋势，为进一步研究冰川消融的变化特征，对模拟期间的日平均气温进行计算，发现每日的消融变化与气温变化具有高度的一致性（r=0.73，P＜0.01），最高的日均温（7.8 ℃）同样出现在8 月3日，日均温最低（-0.1 ℃）的一天为8月19日，这说明了表碛覆盖型冰川的消融变化特征同样受到环境气象条件尤其是气温的控制。

3.4 云量对表碛区气象和能量条件的影响

为进一步研究不同天气条件下表碛区的气象和能量特征，本文分别探讨了晴天和阴天天气中气象要素的日均变化，云量对各种气象要素都有不同程度的影响（图6），云量对入射短波辐射的影响主要体现在削弱其峰值，晴天正午左右入射短波辐射的数值高达850 W∙m-2，阴天时不足600 W∙m-2，云层对下行长波辐射的作用则完全相反，阴天的下行长波辐射始终大于晴天，但二者在一天之中的变化趋势保持一致。晴天气温略高于阴天，变化幅度在1.5～6 ℃之间。相对湿度并无明显的变化规律，阴天的相对湿度水平处于较高水平，平均相对湿度为70%左右，晴天上午时间段内相对湿度较小，仅为50%，随后相对湿度逐渐增大。一天之中风速的变化有较为明显的2个峰值，正午左右风速最小，快到日落时风速最大，阴天的风速变化与晴天相比较为平缓。降水主要发生在下午，阴天的降水量明显较大，最高可达0.7 mm·h-1。由此可以看出，冰面的气象条件在云量的影响下时间分布特征和数值大小都发生了相应改变，也会对表碛区的热量和能量交换造成一定的影响。

图6 不同天气条件下气象要素的变化Fig.6 Effects of different weather conditions on meteorological elements

通过云量区分表碛区的天气状况后，可以看出，不同天气条件对能量平衡方程的各辐射通量也产生了相应的影响（图7和表2）。云量对净短波辐射通量的影响较大，晴天正午时刻的净短波辐射超过750 W∙m-2，阴天的净短波辐射则削减为530 W∙m-2左右，与研究期内所有天气条件相比，晴天的平均净短波辐射增加了11.9%，阴天则减少了20.5%。同样的，晴天条件下感热通量的最大亏损量超过-500 W∙m-2，阴天则削减为-330 W∙m-2左右，晴天和阴天的平均感热通量分别为109.04 W∙m-2、72.30 W∙m-2。云量的增加虽然减弱了净短波辐射的收入和感热通量的支出，但二者的时间序列变化特征却几乎没有发生变化。阴天条件下净长波辐射、潜热通量和传导热通量与晴天时相比，其数值均减小了50%左右，晴天时潜热通量的平均亏损量较小，仅为-2.63 W∙m-2，净长波辐射和传导热通量的平均亏损量则在阴天条件下明显减少。由于观测期间降水稀少，降水热量可忽略不计，在不同天气条件下也无明显变化。晴天和阴天的平均消融量分别为0.58 mm w.e.和0.51 mm w.e.，阴天的消融量比晴天减少了12%。总体而言，云量的多少对能量通量的数值大小有显著影响，但未改变能量通量随时间的变化特征，阴天时净短波辐射的收入和其他能量项的亏损都有所减少，晴天时除潜热通量亏损减少外，净短波辐射的收入和其他能量项的亏损都在增加。

表2 不同天气条件下的平均能量平衡参数分量及其变化Tab.2 Mean fluxes and variation under different weather conditions

图7 不同天气条件对辐射通量的影响Fig.7 Effects of different weather conditions on radiation flux

4 讨论

4.1 敏感性分析

表碛属性作为表碛覆盖型冰川能量-消融模拟的重要不确定性因素，本文通过对模型中各表碛属性参数进行敏感性分析，可以确定影响表碛覆盖型冰川能量-消融模拟的关键参数，为后续进一步研究提供可靠的基础。在敏感性测试中，除表碛反照率（α=0.086）通过气象站的反射/入射短波辐射计算得出，其他表碛参数的取值均通过文献［28］获取。本研究分别选择表碛导热系数（k）、表碛表面粗糙度（z0）、表碛比辐射率（ε）和表碛反照率（α）4种表碛属性参数用于评估表碛属性参数对消融模拟的影响（表3）。在对青冰滩72 号冰川的实地观测和模拟中，表碛导热系数设定为0.94 W∙m-1∙K-1，当导热系数增加/减少0.2 W∙m-1∙K-1时，模拟消融量平均增加/减少17.8%；当导热系数增加/减少0.4 W∙m-1∙K-1时，消融量平均增加/减少35.7%。模型中表碛表面粗糙度设定为0.016 m，当粗糙度增加/减少0.01 m时，模拟的消融量减少4.6%或增加9.3%。表碛比辐射率的取值为0.94，当辐射率增加/减少0.03 时，模拟的消融量平均减少/增加4.3%；当比辐射率增加/减少0.06 时，消融量平均减少/增加8.7%。表碛对可见光的吸收能量较强，反照率接近于0，气象站处测量的平均表碛反照率为0.086，其数值增加0.1或0.2，模拟的消融量将减少5.6%或11.5%。除导热系数外，表碛的粗糙度、比辐射率和反照率与冰川的消融均呈现出负相关关系，表碛导热系数在能量平衡模型中表现为最敏感的参数，表碛反照率和表面粗糙度具有一定的敏感性，但二者在表碛区的变化幅度并不明显，此外，敏感性分析中表碛比辐射率的变化对消融结果的影响不大。因此，在对表碛覆盖型冰川进行消融模拟时，应高度关注表碛导热系数的实地观测值，以及表碛表面粗糙度和表碛反照率的取值，这些表碛参数的微小变动会对消融结果产生较大变化。由于本文未对所有表碛参数开展实地测量，其率定值可能与青冰滩72号冰川的实际数值存在一定的差距，未来应加强对表碛属性的实地观测，以确定其时空变化和不确定性程度。

表3 表碛参数的敏感性分析Tab.3 Sensitivity analysis of debris parameters

4.2 不同区域表碛覆盖型冰川对比

为进一步理解青冰滩72 号冰川消融区域的能量平衡特征，本文对比分析了不同区域山地冰川的能量和消融状况（表4）。西琼台兰冰川和科其喀尔冰川均位于托木尔峰地区，它们与青冰滩72号冰川有着较为相似的特点，冰川面积巨大，冰崖和冰面湖广泛发育。西琼台兰冰川面积超过100 km2［38］，冰面差异化消融现象显著，1978年在海拔4000 m处表碛区观测的能量收入来源中，净辐射收入占比为73%，感热通量占比23%，凝结潜热仅占4%左右，平均消融能量为77 W∙m-2［39］。位于科其喀尔冰川4200 m 处的观测站显示5—9 月的平均气温在1 ℃左右，而极端最高气温达到了12 ℃左右，太阳辐射强烈，有利于冰川消融，与此同时5—9 月集中了全年75%的降水，降低了冰川能量交换的水平［40］，2004年在该冰川表碛区模拟的消融能量为10.5 W∙m-2［41］。藏东南嘎隆拉冰川［12］的能量通量分布状况与地处天山的青冰滩72号冰川较为相似，仅有净短波辐射和潜热通量存在一定差异，可能原因是嘎隆拉冰川纬度较低接收了更多的短波辐射，同时受到东亚季风环流的影响冰面水汽含量较高从而产生大量的潜热交换。与青冰滩72号冰川相比，阿尔卑斯山脉的Miage 冰川接收了更多的太阳辐射，而净长波辐射的亏损比72号冰川高出90%以上，传导热通量的能量消耗也十分巨大，最终导致该冰川可供消融的能量不足5 W∙m-2［28］。喜马拉雅山脉的Lirung 冰川纬度低于30°N，冰川末端海拔接近4000 m，2013 年5 月该冰川可供消融的能量为123 W∙m-2［42］，消融程度远高于其他区域的冰川。喀喇昆仑山所处纬度较高，因此该地冰川接收的短波辐射显著低于其他区域，除此之外，Collier等［43］在喀喇昆仑山模拟的感热通量和传导热通量均为正值，而其他单点尺度的模拟研究则显示除净短波辐射外，其余能量一般为能量支出状态，这表明模型尺度和空间分辨率的差异对辐射通量和消融模拟结果具有重要的影响。总体来看，不同区域的冰川消融状况主要取决于当地的环流特征与地表气象条件，同样，不同的能量平衡模型以及气象要素参数化方案的差异也会影响模拟结果的准确性。

表4 不同区域表碛覆盖型冰川的能量分布状况Tab.4 Surface heat fluxes components of debris-covered glaciers in different regions

5 结论与展望

本文以天山托木尔峰青冰滩72 号冰川为例进行消融模拟研究，该冰川所在区域蕴藏着丰富的冰雪水资源，平均每年产流约63.4×108m3，其中约50%为冰雪融水，且分布着众多的表碛覆盖型冰川，这类冰川由于表碛覆盖形成差异化消融，对冰川物质平衡和径流评估带来一定的挑战，因此，本文基于表碛区自动气象站小时尺度数据作为输入，采用表碛覆盖型冰川能量平衡模型开展下覆冰消融模拟研究，结合野外实测消融和表碛温度数据进行验证，得出以下结论：

（1）研究期内能量平衡模型模拟的消融量为0.39 m w.e.，基于花杆实测的消融量为0.40 m w.e.，R2=0.92，RMSE=±0.03 m w.e.，能量平衡模型取得了较好的模拟效果，模型计算的表碛表面温度和表碛层内温度也与野外实测温度数据具有较好的一致性。

（2）在冰川表碛区的能量收支过程中，净短波辐射是唯一的能量收入项，是冰面的主要能量来源，其他能量项均处于亏损状态，感热通量是最大的能量支出项（49.7%），其次分别为传导热通量（消融耗热）（25.8%），净长波辐射（19.8%）和潜热通量（4.6%），降水热量不足1%，可忽略其对冰面能量交换的影响。

（3）不同天气条件对表碛区气象和能量变化有显著影响，但不改变它们的时间分布特征，入射短波辐射、下行长波辐射和相对湿度受云量的影响程度较大，阴天的云量显著削弱净短波辐射通量，其他辐射通量也较晴天减少50%左右，阴天的平均消融量比晴天减少了12%，需注意云量对冰川能量收支状况的影响。

（4）通过对表碛属性进行敏感性试验，模拟的消融量对导热系数的变化最为敏感，表碛反照率和表面粗糙度的变化量也不容忽视，而表碛比辐射率对消融模拟影响不大。

本文仅涉及冰面单点的消融模拟，由于冰川不同区域的气象条件及地形特征具有很大的差异性，为了厘清整条冰川乃至流域尺度的表碛覆盖型冰川的消融状况，未来应从以下三个方面加强大区域的能量消融和动力模拟：（1）完善对冰面气象要素的时空变化特征研究，以确定气象要素的参数化方案在大范围区域内取得最佳的全局效果；（2）加强对表碛厚度以及表碛导热系数的实地观测和模拟，以精确驱动消融模型；（3）构建冰川动力学参数化方案，并耦合能量平衡模型，更为细致地研究冰川动力和运动特征对冰川消融的影响。

致谢：感谢中国科学院青藏高原研究所杨威老师在模型方面给予的指导帮助！