过去几十年里，全球气候变化对冰川的影响引起了学术界广泛关注。根据IPCC（Intergovernmental Panel on Climate Change）AR6最新报告，2011—2020年全球地表温度比1850—1900年升高了1.1 ℃^［1］。这导致了全球山地冰川在1971—2019年间的物质亏损速率为（170±80） Gt·a^-1。进入21世纪后，全球山地冰川物质亏损速率更为显著，从2000—2009年间的物质亏损（240±9） Gt·a^-1增加到2010—2019年间的（290±10） Gt·a^-1［1-2］。冰川是地球上重要的淡水资源，冰川对全球水循环、海平面上升及生态系统稳定发挥着至关重要的作用。在全球气候变化背景下，青藏高原的气温快速升高导致冰川迅速退缩、冰湖明显扩张和冰川径流增加等现象频繁发生^［3-6］。帕隆藏布流域位于青藏高原东南部，气候变化显著影响流域内冰川变化。

藏东南是青藏高原冰川物质亏损最为严重的区域^［7-9］，冰川物质平衡的研究方法主要有花杆观测法、模型模拟法、大地测量法等。在区域尺度上，利用大地测量的方法观测冰川物质平衡较为常见，研究发现2000—2020年藏东南海洋性冰川物质平衡总体处于亏损的状态，平均物质平衡介于-0.66~-0.61 m w.e.·a^-1之间^［8-10］，且自2010年以来，冰川的物质亏损呈现加速态势，2010—2020年藏东南地区冰川平均物质平衡达到了-0.75m w.e.·a^-1［8］。藏东南的不同山区之间冰川物质平衡差异非常明显，根据多源遥感数据监测显示，念青唐古拉地区是藏东南冰川亏损最强烈的区域，在2010—2019年最高达-1.11 m w.e.·a^{-1［11-13］}，横断山区在2003—2009年平均冰川物质亏损为（-0.40±0.41） m w.e.·a^-1［14］。在模型模拟方面，Zhang等^［15］使用考虑表碛覆盖的能量平衡模型重建藏东南海螺沟流域7条海洋性冰川物质平衡，1952—2009年海螺沟流域平均物质平衡为-0.42m w.e.·a^-1；Yang等^［16］在藏东南24K冰川，在考虑有表碛覆盖和无表碛情况下使用度日模型对物质平衡计算结果的精度进行了评估，并使用能量平衡模型重建帕隆94号冰川1980—2010的物质平衡，结果表明，冰川的物质平衡呈现年际波动，且与大尺度的大气环流有关^［17］。综上分析，大地测量数据集只能计算近些年冰川物质平衡变化，不能反映冰川物质平衡的年际变化特征，且受限于影像的时期、自然条件，如藏东南云雨天气较多，会对不同时间段的影像质量造成一定的影响^［18］。能量平衡模型考虑精细化的输入参数，以提高重建和预测的精度，但仅限于单条冰川，且受限于输入数据的稀缺，不能在区域尺度上进行应用^［19-20］。相比而言，冰川温度指数模型可以重建区域尺度冰川的历史物质平衡，也可以分析区域冰川的年际变化特征、时空变化特征，具有区域可推广性^［21-22］。

帕隆藏布流域位于雅鲁藏布江流域下游，随着青藏高原气温快速升高，该区域冰川物质平衡已发生显著变化^［23］。目前，该流域冰川实测数据匮乏，大地测量数据对长时间历史序列重建存在局限性^［9］，较为精细的能量平衡模型不具有区域推广性，流域冰川年际变化特征以及对冰川的敏感性研究相对较少^［24］。因此，本文将对1980—2019年帕隆藏布流域冰川物质平衡进行模拟研究，明晰帕隆藏布流域冰川物质平衡年际变化特征以及对气候变暖的响应机制，对于了解该区域冰川物质平衡变化，评估区域水资源、防范水资源危机以及制定合理水资源管理策略具有重要的科学价值和现实意义。

帕隆藏布流域地处青藏高原东南部（图1），位于喜马拉雅山脉东段，念青唐古拉山北部、伯舒拉岭东侧，与南部的祁灵公山相邻。帕隆藏布发源于八宿县的然乌湖，向西南流经，与易贡藏布汇合后汇入雅鲁藏布江。流域地势东南高、西北低，地理范围为（95°07′~97°10′ E， 29°12′~30°67′ N）。流域属于高原温带季风半湿润气候，位于印度洋夏季季风进入青藏高原的重要水汽通道上，平均气温适中，年均气温在10~12 ℃之间；降水丰富，年均降水量为884.5 mm，主要集中在6—10月^［25］。湿润的水汽和较高的海拔使该区域成为海洋性冰川理想的发育场所，是我国除横断山以外海洋性冰川最重要和最集中的发育区^［26］。帕隆藏布流域拥有丰富的冰川资源（图2），根据RGI 6.0（Randolph Glacier Inventory），该流域分布有1 554条冰川，冰川总面积达1 799.4 km²，流域内最大的冰川是雅弄冰川（G096657E29334N），面积为179.58 km²，流域内冰川面积的平均值为1.15 km²。

本文使用的气象数据集为GSWP3_W5E5，由ISIMIP3b（Inter-Sectoral Impact Model Intercomparison Project 3b）生产，是OGGM v1.6系统数据集。时间覆盖的范围是1979－2019年，空间分辨率为0.5°，ISIMIP项目的固定模型输入分辨率为0.5°。W5E5是由ERA5数据经过偏差矫正之后得到的数据集，W5E5数据集通过应用WATCH（WFD）强迫数据方法进行偏差矫正，它通过应用一系列高程和月度偏差校正方法来改善再分析数据的准确，W5E5是由更高分辨率的ERA5聚合而成W5E5具有更高的时间分辨率（每小时）和更高的空间变异，对全球分布的13个FLUXNET2015站点评估表明，与使用原始的ERA5数据进行模型驱动相比，W5E5的偏差矫正后的月降水量在未校准的水文模型（WaterGAP）中分析时模拟出更合理的全球水文水量平衡成分^［27］。W5E5数据集可以用来计算冰川特定海拔高度的气温和降水，对于单条冰川，月气温和降水的时间序列从距离冰川最近的W5E5气象格点提取，然后根据特定的气温梯度［0.65 ℃·（100m）^-1］转换成相应海拔梯度的气温和降水，由降水矫正参数（pf）用来校准再分析降水数据，从而更精确获取高海拔地区降水特征。

冰川的边界信息来自RGI（Randolph Glacier Inventory） 6.0^［28］，由全球陆地冰测量空间（GLIMS）开发，2017年发布，用于模型的初始地形处理。RGI V6.0冰川边界的生产年份集中在2000年左右，在本文中，对模拟时间段内的冰川动力学过程进行简化，假设1980—1999年的冰川边界及冰川面积不发生变化，2000—2019年冰川边界及冰川面积发生变化^［29］。本文使用具有1弧秒分辨率的NASADEM，将根据冰川位置和轮廓下载并插值到本地网络，NASADEM是现代化和航天飞机雷达地形任务（SRTM）数据的相关产品，它依靠冰、云和陆地高程卫星（ICE Sat）地球科学激光高度计系统（GLAS）的激光雷达地面控制点来改进地表高程测量，从而提高了地理定位的准确性^［30］。OGGM中冰川厚度的计算方法采用冰通量的方法^［31］，冰川的厚度数据使用Farinotti等^［32］发布的全球冰川厚度数据集，当在level3~level5运行OGGM模型时，模拟区域的冰川厚度数据将进行校准，以匹配全球冰川厚度数据集在区域尺度上的值。

在本文研究中，OGGM（Open Global Glacier Model）是全球开放冰川模型，用于模拟冰川的物质平衡^［33］，OGGM是可以模拟冰川物质平衡、冰厚估算和动力学模拟的开源数值模型框架，是可以模拟冰川动力学的全球模型，该模型依靠浅冰近似来计算沿多条流线上冰川的深度积分通量^［34］。OGGM的计算成本合理，计算性能优越^［35］，可以在个人计算机上模拟区域冰川尺度。利用OGGM模型在全球冰川的模拟上，Zhao等^［36］利用模型重建高亚洲冰川物质平衡及预测到21世纪末高亚洲冰川径流的变化趋势，Afzal等^［34］在喀喇昆仑的Hunza地区重建了该流域1960—2020年冰川历史物质平衡，并利用CMIP6数据驱动模型计算到本世纪末该流域冰川径流的变化。Tang等^［37］利用模型量化青藏高原表面变暗对当地冰川的影响。

在本研究中，使用更新的OGGM v1.6。OGGM中的冰川物质平衡模块是Marzeion等^［38］开发的温度指数模型的拓展版本，能模拟流线冰川网格点历史物质平衡，计算公式如下：

式中：Mi(z)为z海拔高度带处的月度物质平衡的值；pf为降水偏差矫正因子；μ*为温度敏感性因子；Pisolidz为在z海拔高度处月固态降水的值；Ti为月空气温度；Tmelt为冰川冰消融的温度，一般将温度设置成0 ℃，ε为残差因子，降水的状态是根据雨雪分离的阈值来确定的，融点是冰融化时的温度，在本研究中，Tmelt设置为-1 ℃，我们假设固态降水量等于总降水的百分比，若Ti≤Tsolid（0 ℃），固态降水量为100%，若Ti≥Tliquid（2 ℃），固态降水量为0%，若Tsolid≥Ti≥Tliquid，线性插值得到固态降水的百分比^［39］，高海拔地区的降水由于实测数据的匮乏难以量化，因此采用模型的默认设置。

大尺度冰川模型模拟的最大挑战之一是对于模型的校准，主要原因在于降尺度后的气象要素不能准确描述冰川区的具体情况^［40］，且冰川模型不能精确地解决大尺度范围冰川的复杂过程，因此需要校准。大多数冰川模型都是使用实测物质平衡数据用于模型的校准^［41］，在一些研究中，模型参数被进一步微调以匹配由冰川观测数据外推得到的区域冰川物质平衡估计值^［42］。

帕隆藏布流域内，有长时间实测数据序列的冰川较少，考虑到实测数据的局限性，本文使用全球冰川物质平衡变化数据集（Accelerated global glacier mass loss in the early twenty-first century）^［10］用于OGGM模型的校准， 数据集包含2000—2020年单条冰川的物质平衡数据，本文选取2000—2020年冰川平均物质平衡数据集作为校准数据集，对OGGM模型模拟的流域进行校准。校准的参数包括：温度敏感性参数、降水校正参数、温度偏差参数。校准采用的方法Huss等^［43］采用的三步校准法，使用OGGM v1.6的‘mb_calibration_from_scalar_mb’函数进行校准，每一个参数都会根据历史模拟情况被设置在合理的范围，如温度敏感性因子被设置在的区间范围在（1.5，17），降水矫正参数被设置在的范围为（0.1，10），温度偏差参数被设置在（-15，15）。通过设置模型的校准顺序，根据大地测量数据首先校准温度敏感性因子，若校准后的温度敏感性参数不足以匹配大地测量数据，则依次校准降水参数，若校准后的降水参数依旧不能匹配大地测量数据，最后校准温度偏差参数，使得校准后的参数模拟的物质平衡以匹配大地测量数据，若在三步校准中出现提前匹配大地测量数据，校准程序提前结束。为了最大程度上减少超参数化对模型校准的影响，选择优先校准温度敏感性参数，从而使得校准后的参数可以更大程度上匹配冰川区实际的气象情况。

用校准后的参数去分别模拟2000—2010年和2010—2020年该流域1 554条冰川的物质平衡，再使用全球冰川物质平衡变化数据集2000—2010年及2010—2020年单条冰川的大地测量物质平衡数据用于对校准后的模拟结果进行验证（图3）。考虑到实测数据的局限性，验证数据包括：帕隆94号冰川2006—2018年实测数据、全球冰川物质平衡数据集，验证的范围包括：流域尺度、海拔带尺度、单条冰川尺度。

使用2000—2020年大地测量数据集校准后的OGGM物质平衡模型，模拟帕隆藏布流域2000—2010年与2010—2020年帕隆藏布流域冰川物质平衡，将模拟的数据分别与大地测量物质平衡（2000—2010年每一条冰川10年平均物质平衡与2010—2020年每一条冰川10年平均物质平衡）对比分析，在前10年，模拟结果与大地测量数据集的相关性达到了0.76，后10年达到了0.79。证实了模型在流域尺度上模拟的可靠性。

为了验证校准后的冰川模型在海拔梯度上模拟的可靠性，选取帕隆藏布流域内物质平衡实测冰川——帕隆94号冰川，帕隆94号冰川实测物质平衡的时间序列为2006—2018年，且帕隆94号冰川花杆的位置主要分布在海拔5 000~5 300 m，利用模型模拟2006—2018年帕隆94号冰川在5 000~5 300 m海拔带内物质平衡的平均值，与花杆数据对比，考虑到实测数据在每一年花杆的海拔高度和数量的差异性，将测杆数据统计成4个海拔带：5 030~5 100 m、5 100~5 200 m、5 200~5 300 m、5 300 m及以上，除5 300 m仅统计一个数据，每一个海拔带统计两个数据，利用校准后的模型模拟帕隆94号冰川2006—2018年在5 000~5 300 m海拔的物质平衡。模拟结果与实测结果显示出较高的一致性，证明校准后的模型在海拔尺度上模拟结果的可靠性（图4）。

为了验证校准后的物质平衡模型在单条冰川上模拟的准确性，选取帕隆94号冰川2008—2018年的实测物质用于模拟结果的验证。帕隆94号冰川的中流线分布在5 000~5 500 m，但是花杆实测数据主要集中在海拔5 000~5 300 m，这里使用外插的方法将冰川的物质平衡数据源外插到整个冰川流线，模拟结果显示，校准后的模型在单条冰川上的模拟表现较好，在趋势上，模拟数据与实测数据有较好的相关性（图5）。

考虑到不同冰川的规模和大小不同。在计算流域物质平衡时，为了提高流域物质平衡计算的准确性，在单条冰川基础上，计算流域冰川加权面积的物质平衡^［44］。

式中：M1代表第一条冰川的年均物质平衡；S1代表冰川的面积；Mn代表的是第n条冰川的年均物质平衡；Sn代表第n条冰川的面积。

通过上述公式，计算得到帕隆藏布流域在1980—2019年物质平衡为-0.41 m w.e.·a^-1， 在1980—1999年，1 554条冰川加权面积后，物质平衡为-0.25 m w.e.·a^-1，在2000—2019年物质平衡亏损更为严重，达到-0.56 m w.e.·a^-1，在2000—2010年物质亏损为-0.49 m w.e.·a^-1，2010—2019年物质亏损为-0.63 m w.e.·a^-1（图6）。将模型模拟结果与相关研究进行对比，表现出较好的一致性。

在空间上，分析了1980—1999年帕隆藏布流域冰川物质平衡，圆圈的大小代表了不同的冰川规模，不同的颜色代表了冰川物质平衡的不同区间，分析1980—1999年流域冰川物质平衡分布，流域在1980—1999大体上呈现负平衡，其中流域中部和西北部呈现轻微的正平衡，流域的东南部是该时间段内物质亏损最为强烈的区域，最高达到-1.7 m w.e.·a^-1。整体的物质亏损特征是：中部、西北部亏损较低，出现轻微正平衡，东南部亏损最为严重，在1980—1999年，流域冰川整体亏损为-0.25 m w.e.·a^-1［图7（a）］。

2000—2019年帕隆藏布流域全域冰川物质平衡亏损较为强烈，其中帕隆藏布流域东南部和亏损较为强烈，最高可达-2.7 m w.e.·a^-1，流域中西部冰川物质亏损相较于流域的其他区域，稍显平缓。相较于1980—1999年，流域西北部冰川的亏损速率有所加快，中部冰川依旧是全域亏损最为缓慢的区域。整体的物质亏损特征是：东南部亏损最为强烈，西北部冰川亏损速率较快，中部区域亏损最为缓慢。在2000—2019年，流域冰川整体物质平衡达到-0.56 m w.e.·a^-1［图7（b）］。

分析流域1980—2019年全域冰川物质平衡变化（图8），发现物质平衡亏损较为强烈的区域位于帕隆藏布流域东南部和西北部，流域中部和流域西部的冰川相较于其他区域而言，冰川物质亏损较小。分析不同规模冰川的变化情况，1980—1999年面积小于1 km²的冰川物质平衡在-1 000~0 mm w.e.·a^-1的占比为52%，在2000—2019面积小于1 km²的冰川物质平衡在-1 000~0 mm w.e.·a^-1的占比为80%，而面积在1~10 km²的冰川1980—1999年物质平衡在-1 000~0 mm w.e.·a^-1的占比为62%，在2000—2019的占比则为86%。分析表明，小冰川在21世纪初的20年内，退缩速率在持续加快。总体而言，1980—2019年全域冰川呈现亏损的态势，随着全球气候变暖和局部极端天气事件，该流域冰川将会处于持续亏损的状态中。

受印度洋夏季季风的影响，藏东南地区的海洋性冰川是对气候变化最敏感的区域。为了厘清流域冰川对气温和降水的敏感性，在其他参数不变的情况下，选取帕隆藏布流域冰川物质亏损最为强烈的时间段（2000—2019年），通过调整温度或者降水的单一参数实验，分析流域冰川物质平衡的变化趋势（图9）。实验中，将温度以0.5 ℃（K）为参数调整间隔，在其他参数不变的情况下，将温度在-1.5~1.5 ℃区间进行实验。同样，将降水以10%为参数调整间隔，将降水在-30%~30%区间进行实验，分析冰川的物质平衡对气温和降水变化的响应。

通过分析发现，在提高1 ℃的温度，流域全域冰川处于亏损之中，其中71.75%的冰川物质平衡变化在-1 000~-500 mm w.e.·a^-1，通过下降20%的降水，流域全域物质发生亏损，其中62.81%的冰川物质平衡变化在-450~-300 mm w.e.·a^-1。整体而言，帕隆藏布流域冰川，对气温的敏感性会大于对降水的敏感性。近些年随着局部气温上升，降水补给对冰川的影响小于气温变化带来的影响，且冰川对气温的敏感性远大于对降水的敏感性，气候变化下气温快速升高，降水对冰川的补给不足，造成冰川处于持续亏损之中。

帕隆藏布流域附近，分布有两个距离最近的国家气象站，分别是八宿和波密国家气象站（图10）。分析这两座气象站的历史气温和降水数据，发现在最近40年内，八宿站和波密站气温呈现不断上升的趋势，降水呈现轻微下降的趋势，其中八宿气象站从1980—2014年的年平均气温上升了超过1.5 ℃，波密气象站从1980—2019年的40年间，气温上升也超过了1.5 ℃。分析流域冰川物质平衡时间序列，1997年与2009年分别是全域冰川物质平衡亏损最小与最大的年份，通过对比八宿站和波密站的气温数据，发现1997年与2009年均为两站气温最低与最高的年份，说明气温的加速上升是造成帕隆藏布流域冰川快速亏损的主要原因。

分析两个国家气象站的降水数据，波密站1980—1999年这20年间的总降水量达18 554 mm，2000—2019年总降水量为16 669 mm，相比于前20年，波密气象站后20年的降水减少了1 885 mm，总量下降了10%。八宿气象站在1985—1999年总降水为4 115 mm，在2000—2014年为3 789 mm，后15年的总降水相较于前15年的总降水有轻微的减少。总体而言，在帕隆藏布流域的两个国家气象站的气温和降水长时间序列的数据表明，在1980—2019年，帕隆藏布流域气温呈现不断上升的趋势，相比于1980年，两个国家气象站的气象上升幅度都超过了1.5 ℃，两个国家气象站的降水都呈现出减少的趋势，其中波密站降水减少更为明显。

提取1980—2019年模型驱动数据集在研究区域的变化情况，分析气温和降水的空间变化率，全域气温呈现上升趋势，其中气温上升最为明显的区域在流域的西北部，最高可达0.19℃·（10a）^-1，流域中部是气温上升相对较小的区域，这与我们物质平衡的模拟结果变化情况一致。受印度洋夏季季风的影响，流域南部气温上升速率最小，北部气温上升速率最大。降水在研究区域的总体变化趋势为不断减少，但在帕隆藏布流域的减少幅度较小，平均为-40~-20 mm·（10a）^-1，这与在气象站观测到的结果接近，说明W5E5数据在帕隆藏布流域的数据表现与在两个气象站点观测到的结果具有较好的相关性（图11）。

通过分析再分析数据气温和降水变化率的空间分布情况，发现气温变化的空间分布情况相比降水更为复杂，存在明显的空间异质性。现提取流域格点数据，分析流域气温和降水的季节性变化特征。分别提取1980—2019年流域夏季（7—10月）与冬季（11月—次年2月）气温和降水数据。分析发现夏季气温存在明显的波动，且自21世纪以来，气温呈现明显的加速上升趋势，夏季降水的总量较大，但总量总体呈现下降趋势。冬季气温上升平缓，但自21世纪以来，呈现加速上升的趋势，冬季降水总量较夏季降水而言偏小，趋势变化不是特别明显，但总体而言，呈现略微下降的趋势（图12）。

利用OGGM模型模拟冰川的物质平衡，存在诸多的不确定性因素。首先，冰川对气候变化的响应剧烈，冰面裂隙、冰内冰下过程、冰崩等过程，都在冰川物质平衡变化中起到了一定的作用，但由于缺乏详细的野外监测数据和气象、地形等驱动数据，这些物理过程OGGM模型考虑较少，且对冰川边界的理想化考虑，使得模型模拟结果存在一定的不确定性。其次，模型在固液态降水的划分上采取了相对简单的双阈值法，冰川位于高海拔地区，实测资料匮乏，固液态降水阈值的选择是一个复杂的问题，模型采取双阈值法对固液态降水进行分离，可能会造成最高海拔地区降水量的误判，从而对物质平衡模拟的准确性造成一定的影响。最后，冰川表碛对冰川物质平衡有显著影响，藏东南地区冰川存在表碛覆盖的情况，表碛覆盖会对冰川的消融产生影响，然而目前大部分的冰川模型都没有考虑表碛覆盖对冰川物质平衡模拟的影响，可以将表碛参数的经验模型与冰川模型相结合，以量化表碛对冰川物质平衡的影响。

本文使用OGGM模型对藏东南地区帕隆藏布流域近些年冰川物质平衡进行模拟及对其物质亏损进行了归因分析，主要结论如下：

（1）1980—2019年帕隆藏布流域全域冰川物质平衡呈现不断亏损的状态，为-0.41m w.e.·a^-1，在2000—2019年物质平衡亏损更为严重，达到-0.56 m w.e.·a^-1，2000—2010年物质亏损为-0.49 m w.e.·a^-1，2010—2019年物质亏损为-0.63 m w.e.·a^-1。

（2）流域东南部和流域西北部是冰川亏损最为严重的区域，流域中部和西部冰川亏损相对更小。流域内小冰川的亏损速率在不断加快，面积小于1 km²的小冰川在1980—1999年冰川物质平衡在-1 000~0 mm w.e.·a^-1的比重为62%，在2000—2019年比重上升到80%。而面积在1~10 km²的冰川1980—1999年物质平衡在-1 000~0 mm w.e.·a^-1的占比为62%，在2000—2019的占比则为86%。

（3）通过对帕隆藏布流域冰川进行敏感性分析，冰川对气温变化的敏感性会大于对降水变化的敏感性。温度上升1 ℃，71.75%的冰川物质平衡变化在-1 000~-500 mm w.e.·a^-1，降水减少20%，62.81%的冰川物质平衡变化在-450~-300 mm w.e.·a^-1。

（4）通过分析流域内国家气象站点数据及再分析数据，发现1980—2019年气象站气温上升超过1.5 ℃，波密站2000—2019年总降水相较于前20年，减少了10%，流域降水整体呈现减少的趋势，再分析数据在冬季与夏季呈现与站点一致的趋势，气温的持续升高和降水的轻微减少导致帕隆藏布流域冰川处于持续亏损之中。