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年?，F代冰川數據來源于國家冰川凍土沙漠科學數據中心（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］?；贏AR與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）?？傮w來看，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。