1984—2019年念青唐古拉山中段冰川ELA变化估算及特征分析

李亚鹏，张 威，柴 乐，唐倩玉，葛润泽，孙 波

（1.辽宁师范大学地理科学学院，辽宁大连 116029；2.东华理工大学地球科学学院，江西南昌 330013）

0 引言

政府间气候变化专门委员会（IPCC）第六次评估报告对观测到的气候系统变化显示，人类活动估计造成了全球升温高于工业化前水平约1.0℃，可能区间为0.8～1.2℃。如果继续以目前的速率升温，全球升温可能会在2030年至2052年达到1.5℃（高信度）。同时，报告显示了过去二十年以来，几乎全球范围内的冰川都呈现退缩状态（高信度）［1-2］。冰川过度融化将会在区域和全球范围内造成严重影响［3-5］。

ELA是气候对冰川系统影响的关键参数，与当地气候条件密切相关。特别是固态降水和气温，即当气温升高，固态降水量减少时，ELA呈上升趋势，反之下降［6］。在全球变暖的背景下，ELA变化的研究应该引起重视［7-8］。ELA的计算主要是通过现代冰川卫星监测或数学模型模拟或者利用经验和统计方法等进行估算［9-13］，随着人们对ELA对气候响应的逐步充分认识，更多的方法也被应用在ELA研究中，如利用ELA变化量与冰川物质和能量平衡的关系式进行计算等［14-16］。尽管如此，由于青藏高原及其周边地区长时间尺度连续观测数据仅有几条冰川（主要有天山乌鲁木齐河源1号冰川和唐古拉山的小冬克玛底冰川），其余大多数监测冰川的平衡线高度只有近几年的观测资料，时间序列相对较短。同时，高原气象站点的分布主要集中于东缘，且部分冰川分布区没有站点检测数据，这在很大程度上限制了不同区域冰川平衡线高度变化状况及其气候敏感性的对比分析与研究［17-18］。

念青唐古拉山作为西南季风暖湿气流进入高原的通道和藏东南海洋型冰川区向唐古拉山大陆型冰川区的过渡区域，是青藏高原东南缘重要的山谷冰川分布区［19-20］。对于念青唐古拉山区域内各冰川相关研究表明不同流域冰川变化及其特征差异显著［21-30］。因此，在前人研究基础之上，本研究利用自然分域法，选取念青唐古拉山中段冰川分布区作为研究点，通过遥感影像、冰川编目数据以及气象数据提取与分析，利用气温与固态降水数据和ELA作为输入参数，尝试建立多元线性回归方程，重建研究区1984—2019年间冰川ELA，讨论并分析其年际尺度下的ELA变化趋势、区域异同性以及气候响应。

1 研究区概况

念青唐古拉山中段冰川作用区位于怒江流域以南，雅鲁藏布江大拐弯西北部的麦曲流域、霞曲流域及易贡湖以西的易贡藏布流域（图1）。地势南高北低，海拔5 500 m以上，相对高差大于1 000 m，平均海拔3 600 m。气候属高原温带半湿润气候类型。年极端最低气温为-40℃，年平均气温为-1℃，年均降水量为400～700 mm［31-32］。冰川编目资料显示［33］，区 域 内 共 有500条 现 代 冰 川，面 积 为2 716.608 km2，以中小型冰川（面积<5 km2）为主，约占区域整体冰川的78%，规模相对较大（面积≥10 km2）的冰川有47条，平均海拔范围为3 000～6 000 m，主峰孔嘎峰海拔6 490 m。

图1 研究区位置（a）与冰川分布（b）Fig.1 Location of the study area（a）and distribution of glaciers（b）［The base map is drawn with 30-m ASTER（https：//www.gscloud.cn），the surveying and mapping area is bounded by Nyainqêntanglha Mountains［19］，the modern glacier data is from RGI 6.0［34］，and basin data is from BasinATLAS_v10_lev07（https：//hydrosheds.org/downloads）］

2 数据与方法

2.1 数据来源

气象数据来源于美国国家气候数据中心（ftp：//ftp.ncdc.noaa.gov/pub/data/noaa/isd-lite/）和中国气象数据网（http：//data.cma.cn）。由于研究区没有气象站点，通过选取周边邻近地区8个气象站点，分别为丁青、类乌齐、昌都、洛隆、波密、八宿、察隅和林芝，利用标准克里金插值获取研究区数据，时间范围为1984—2019年。现代冰川数据来源于国家冰川冻土沙漠科学数据中心（www.ncdc.ac.cn）第二次冰川编目资料与Randolph Glacier Inventory（http：//www.glims.org/RGI/）数据相互补充。

2.2 气温与固态降水量提取

临界气温法区分固液态降水是冰川物质平衡模型中最常用的方法。主要包括双临界气温分离法与单临界气温分离法。双临界气温分离法只适用于我国的一些干旱地区，单临界气温法具有普遍适用性且具有较高精度［35-37］。本研究中采用单临界气温法［38］。

式中：Ps为月固态降水量（mm）；P为月降水量（mm）；T为月平均气温（℃）；Tl为液态降水临界气温值（℃）；Ts为固态降水临界气温值（℃）。为获取研究区固态降水数据，首先参考康尔泗等［38］与张太刚等［39］临界气温方法划分固液态降水的方案，分别取0℃和2℃作为固液降水的划分阈值。研究区气象资料显示，冰川作用现代冰川区固态降水临界气温为-4.2～-0.6℃（图2，折线），平均固态降水量为140～614 mm（图2，条带）。

图2 1984—2019年临界气温与固态降水量变化Fig.2 Variations of critical air temperature and solid precipitation from 1984 to 2019

2.3 ELA计算

积累区面积比率法（accumulation area ratio，AAR）与面积-高程平衡率法（area-altitude balance ratio，AABR）是目前最广泛使用的冰川ELA计算方法［40-41］。基于AAR与AABR法基本原理，本研究使 用Pellitero等［41］开 发 的 自 动 计 算ELA工 具（https：//github.com/cageo/Pellitero-2015）。该工具计算ELA过程主要包括两个部分：输入需要重建冰川表面的数字高程模型（digital elevation model，DEM）和定义AABR与AAR值，即可快速计算出单条冰川或多条冰川数据集。

具体操作步骤如下：①输入需要重建冰川表面的DEM。设置用于冰川面积计算的等高线区间，本研究设置为50 m。②定义AABR与AAR值。为获取适合本研究区的AABR与AAR值，需要对其进行确定。AAR值的确定利用Kern等［42］系统总结分析了全球不同地区的46条山谷和冰斗冰川的面积S与其稳定状态下的AAR值之间的关系式，为AAR=0.0648lnS+0.483。由此，计算得出本研究区域AAR值 为0.74±0.02。AABR值 的 确 定 是 基 于Rea［43］系统评估了全球冰川的物质平衡提出的适用于全球范围内的BR平均值为1.75±0.71。鉴于同一地区有多个冰川，Benn等［44］建议对一个BR范围进行计算，然后选择ELA标准差最低的一个。利用该方法，计算得出本研究区的BR均值为1.78±0.68。③使用Osmaston［45］提出的最低标准差方法来确定冰川最佳ELA。返回结果为“ELA_val⁃ues_AAR_and_AABR.txt”。

3 结果与讨论

3.1 ELA重建模型与验证

Singh等［46］利用多元线性回归统计分析法对1979—2013年喜马拉雅山脉的纳拉杜冰川平衡线高度的统计估算及其趋势分析，Zekollari等［47］使用该方法描述了瑞士的Morteratsch冰川的表面质量平衡，Sakai等［48］利用多元线性回归统计分析了亚洲高山冰川对近期气候变化的响应。基于前人研究，利用编目资料数据与研究区已有相关遥感影像解译数据［49］。利用ELA calculation toolbox计算工具，对研究区内1990年、1999年、2011年及2019年的ELA进行计算（表1）。以气象站点数据（气温与固态降水量）和ELA数据作为输入参数，建立1984—2019年间念青唐古拉山中段现代冰川ELA变化多元线性回归方程。

表1 念青唐古拉山中段ELA计算结果Table 1 Calculation results of ELA in middle section of the Nyainqêntanglha Mountains

式中：X1为月平均固态降水量（mm）；X2为月平均气温（℃）。

利用有冰川纪录的年份对1984—2019年间研究区冰川ELA模拟结果进行验证，冰川编目资料［32］记 录 了1989—1990年、1999—2000年 及2010—2011年研究区冰川变化特征，对应ELA分别为5 140.1 m、5 403.5 m和5 408.8 m，利用上述方程［式（2）］计算ELA结果分别为5 182.5 m、5 375.0 m和5 390.0 m，ELA模拟值与测量值具有高度一致性，确定系数为0.9977。

同时，为进一步对比验证模拟结果的有效性，通过对气象站1984—2019年的年平均气温、6—8月平均气温、年降水量和6—8月平均降水量进行分析，参考我国西部17条现代冰川的物质平衡以及气象观测资料［50］，建立冰川物质平衡线处年降水量（PELA）与6—8月平均气温（TSO）的拟合方程［式（3）］，并结合综合因子法［51］［式（4）］，来计算念青唐古拉山中段现代冰川ELA。

选取距离研究区域较近的丁青与洛隆气象站数据，站点海拔分别为3 873 m和3 640 m，6—8月平均气温T1=13.1℃，年降水量P1=534 mm，气温梯度∂T/∂Z=0.68℃·（100m）-1，降水梯度∂P/∂Z=29.6 mm·（100m）-1。带入式（3），计算得出海拔增加量ΔH=1 495 m，念青唐古拉山中段ELA为5 368 m，模拟结果为5 282 m，相对误差仅为0.16%，两者基本一致。

3.2 ELA年际变化及特征

通过重建念青唐古拉山中段1984—2019年间研究区冰川ELA。结果显示，研究区多年平均ELA为5 389 m a.s.l.，总体呈上升趋势（图3，直线），90%置信区间（虚线）的确定系数为0.7830。ELA年变化研究显示出波动变化特征（图3，折线），ELA波动范围为5 360～5 420 m a.s.l.，波动幅度约为60 m，平均上升速率为1.57 m·a-1。ELA年变化量正值与负值分别为18年和17年，基本呈现出交替状态（图3，条带）。其中，ELA年际变化量正值总和为319.9 m，负值变化总和为291.2 m。

图3 1984—2019年念青唐古拉山中段ELA及其年变化量Fig.3 ELA and its annual variation in middle section of the Nyainqêntanglha Mountains from 1984 to 2019

3.3 不同流域的ELA变化

利用流域自然分区法，选取HydroSHEDS and WaterGAP v2.2数 据 库 中 的BasinATLAS_Cata⁃log_v10数据，将研究区分为三个流域［图1（b）］，分别为霞曲流域［图4（a）］、易贡藏布流域［图4（c）］和麦曲流域［图4（e）］。霞曲流域、易贡藏布流域和麦曲流域现有现代冰川分别为140、256和104条。冰川表面高程范围分别为3 925～6 713 m a.s.l.、2 689～6 639 m a.s.l.和4 185～6 443 m a.s.l.。1984—2019年霞曲流域、易贡藏布流域和麦曲流域ELA变化范围分别为5 178～5 492 m a.s.l.、4 855～5 120 m a.s.l.和5 150～5 485 m a.s.l.［图4（b），（d），（f）］。

图4 三个流域冰川表面高程与ELA变化Fig.4 Variations of glacier surface elevation and ELA in three river basins（The black line segment in（b），（d），（f）is the average ELA elevation）

在气候变暖的背景下，ELA变化总体都呈上升状态，但不同流域ELA变化具有差异性。其中，麦曲流域ELA变化幅度最大，平均上升335 m，霞曲流域和易贡藏布流域分别上升265 m和314 m，上升速率分别9.57 m·a-1、7.57 mm·a-1和8.97 m·a-1。同时，统计分析显示霞曲流域，易贡藏布流域和麦曲流域年平均ELA高程分别为5 335 m a.s.l.、4 987 m a.s.l.和5 317 m a.s.l.，呈现出西北—东南方向降低，北坡高于南坡。

受印度季风，流域位置及冰川朝向等因素影响，各流域现代冰川区水热组合的差异化特征，进而影响着冰川ELA变化。从流域分布来看，霞曲流域冰川区位于研究区西北方向，冰川朝向以NW为主，易贡藏布流域位于研究区东南方向，冰川朝向以SE为主，麦曲流域位于研究区东北方向，冰川朝向以NE为主。季风路径（图5）显示出易贡藏布流域处于印度季风迎风坡，获得的水汽相较于处于背风坡雨影区的麦曲流域和霞曲流域更多。

从各流域冰川区气温与降水状况看，冰川区整体气温介于-5.6～-3.9℃，平均值为-4.8℃［图5（a）］；降水介于502～564 mm，平均值为547 mm［图5（b）］。易贡藏布流域冰川区气温为-5.5～-3.9℃，均值为-4.9℃；降水为502～554 mm，均值为547 mm；麦曲流域冰川区气温为-4.6～-4.1℃，均值为-4.3℃；降水为510～555 mm，均值为539 mm；霞曲流域现代冰川区气温为-5.5～-4.4℃，均值为-4.8℃；降水为549～557 mm，均值为553 mm。对比各流域冰川ELA和冰川区气温降水组合，发现相较于处于背风坡的麦曲流域，处于迎风坡的易贡藏布流域的水热组合状况更优，更利于冰川发育。

从冰川规模来看，易贡藏布流域冰川规模最大，最大冰川长度为33 km［图5（d）］，总面积为1 238.809 km2，麦曲流域次之，霞曲流域最小，最大冰川长度分别为24 km［图5（b）］和19 km［图5（f）］，总面积分别为760.817 km2和716.982 km2［33］。相关研究表明，冰川对气候变化的反馈作用［52］以及冰川变化对气候的响应具有滞后效应［53］，使得尽管相较于易贡藏布流域冰川区，霞曲流域水热组合状况更优，但冰川规模相对较小，ELA变化对气候响应更敏感。这也是造成短时间尺度范围同一气候区域冰川ELA变化差异性的重要原因。

图5 标准克里金插值处理后的冰川区气温与降水Fig.5 Air temperature and precipitation in glacier area after standard Kriging interpolation

3.4 ELA变化对气候变化的响应

姚檀栋等［54］对高亚洲冰川总体退缩研究表明，气温升高引起的冰川消融量增加和粒雪线上升，是近百年来全球山地冰川呈现总的退缩趋势的主要原因。同时，Caidong等［55］研究发现，当温度变化为1℃时，ELA变化为140 m。段克勤等［56］对小冬克玛底冰川ELA的气候敏感性进行了试验，发现在气候状态基础上，假设其他气候因子不变，当气温增加1℃将引起冰川ELA大幅度上升160 m，而增加20%的降水量引起的ELA仅降低10 m，这也说明冰川ELA对气温变化较敏感。

本研究通过1984—2019年间平均气温、固态降水量与冰川ELA统计，对念青唐古拉山中段冰川ELA变化对气候的响应进行分析（图6）。总体来看，ELA随气温和固态降水量变化都呈显著相关性，确定系数R2分别为0.5873和0.9540。同时，ELA响应气温与固态降水量变化结果显示，研究区冰川ELA主要受气温变化控制，随气温变化1℃，ELA波动幅度为126.02 m，而受固态降水量变化影响，ELA波动仅为0.2492 m。并且，随着气温的持续升高，念青唐古拉山中段冰川ELA上升，冰川将继续处于退缩状态。

图6 ELA随气温与固态降水量的变化Fig.6 Variation of ELA with air temperature and solid precipitation

4 结论

通过建立念青唐古拉山中段冰川作用区多元线性回归方程，对1984—2019年间ELA变化进行了重建。主要结论如下：

（1）研究时段内平均ELA为5 360 m a.s.l.。

（2）ELA总体呈上升趋势，上升速率为1.57 m·a-1。ELA年变化显示出波动变化特征，波动范围为5 360～5 420 m a.s.l.，上升幅度为60 m。

（3）受流域位置、冰川朝向及冰川规模影响，研究区ELA高程呈差异性，霞曲流域、易贡藏布流域和麦曲流域多年平均ELA高程分别为5 335 m a.s.l.、4 987 m a.s.l.和5 317 m a.s.l.，平均上升幅度分别为265 m、314 m和335 m，上升速率分别7.57 m·a-1、8.97 m·a-1和9.57 m·a-1。各流域平均ELA高程对比显示出西北—东南方向降低，北坡高于南坡的变化特征。

（4）冰川区多年ELA变化的气候响应分析显示，ELA变化主要受气温控制，随气温变化1℃，冰川ELA总体波动126.02 m。