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

李亞鵬，張 威，柴 樂(lè)，唐倩玉，葛潤(rùn)澤，孫 波

（1.遼寧師范大學(xué)地理科學(xué)學(xué)院，遼寧大連 116029；2.東華理工大學(xué)地球科學(xué)學(xué)院，江西南昌 330013）

0 引言

政府間氣候變化專門委員會(huì)（IPCC）第六次評(píng)估報(bào)告對(duì)觀測(cè)到的氣候系統(tǒng)變化顯示，人類活動(dòng)估計(jì)造成了全球升溫高于工業(yè)化前水平約1.0℃，可能區(qū)間為0.8～1.2℃。如果繼續(xù)以目前的速率升溫，全球升溫可能會(huì)在2030年至2052年達(dá)到1.5℃（高信度）。同時(shí)，報(bào)告顯示了過(guò)去二十年以來(lái)，幾乎全球范圍內(nèi)的冰川都呈現(xiàn)退縮狀態(tài)（高信度）［1-2］。冰川過(guò)度融化將會(huì)在區(qū)域和全球范圍內(nèi)造成嚴(yán)重影響［3-5］。

ELA是氣候?qū)Ρㄏ到y(tǒng)影響的關(guān)鍵參數(shù)，與當(dāng)?shù)貧夂驐l件密切相關(guān)。特別是固態(tài)降水和氣溫，即當(dāng)氣溫升高，固態(tài)降水量減少時(shí)，ELA呈上升趨勢(shì)，反之下降［6］。在全球變暖的背景下，ELA變化的研究應(yīng)該引起重視［7-8］。ELA的計(jì)算主要是通過(guò)現(xiàn)代冰川衛(wèi)星監(jiān)測(cè)或數(shù)學(xué)模型模擬或者利用經(jīng)驗(yàn)和統(tǒng)計(jì)方法等進(jìn)行估算［9-13］，隨著人們對(duì)ELA對(duì)氣候響應(yīng)的逐步充分認(rèn)識(shí)，更多的方法也被應(yīng)用在ELA研究中，如利用ELA變化量與冰川物質(zhì)和能量平衡的關(guān)系式進(jìn)行計(jì)算等［14-16］。盡管如此，由于青藏高原及其周邊地區(qū)長(zhǎng)時(shí)間尺度連續(xù)觀測(cè)數(shù)據(jù)僅有幾條冰川（主要有天山烏魯木齊河源1號(hào)冰川和唐古拉山的小冬克瑪?shù)妆ǎ?，其余大多?shù)監(jiān)測(cè)冰川的平衡線高度只有近幾年的觀測(cè)資料，時(shí)間序列相對(duì)較短。同時(shí)，高原氣象站點(diǎn)的分布主要集中于東緣，且部分冰川分布區(qū)沒(méi)有站點(diǎn)檢測(cè)數(shù)據(jù)，這在很大程度上限制了不同區(qū)域冰川平衡線高度變化狀況及其氣候敏感性的對(duì)比分析與研究［17-18］。

念青唐古拉山作為西南季風(fēng)暖濕氣流進(jìn)入高原的通道和藏東南海洋型冰川區(qū)向唐古拉山大陸型冰川區(qū)的過(guò)渡區(qū)域，是青藏高原東南緣重要的山谷冰川分布區(qū)［19-20］。對(duì)于念青唐古拉山區(qū)域內(nèi)各冰川相關(guān)研究表明不同流域冰川變化及其特征差異顯著［21-30］。因此，在前人研究基礎(chǔ)之上，本研究利用自然分域法，選取念青唐古拉山中段冰川分布區(qū)作為研究點(diǎn)，通過(guò)遙感影像、冰川編目數(shù)據(jù)以及氣象數(shù)據(jù)提取與分析，利用氣溫與固態(tài)降水?dāng)?shù)據(jù)和ELA作為輸入?yún)?shù)，嘗試建立多元線性回歸方程，重建研究區(qū)1984—2019年間冰川ELA，討論并分析其年際尺度下的ELA變化趨勢(shì)、區(qū)域異同性以及氣候響應(yīng)。

1 研究區(qū)概況

念青唐古拉山中段冰川作用區(qū)位于怒江流域以南，雅魯藏布江大拐彎西北部的麥曲流域、霞曲流域及易貢湖以西的易貢藏布流域（圖1）。地勢(shì)南高北低，海拔5 500 m以上，相對(duì)高差大于1 000 m，平均海拔3 600 m。氣候?qū)俑咴瓬貛О霛駶?rùn)氣候類型。年極端最低氣溫為-40℃，年平均氣溫為-1℃，年均降水量為400～700 mm［31-32］。冰川編目資料顯示［33］，區(qū) 域 內(nèi) 共 有500條 現(xiàn) 代 冰 川，面 積 為2 716.608 km2，以中小型冰川（面積<5 km2）為主，約占區(qū)域整體冰川的78%，規(guī)模相對(duì)較大（面積≥10 km2）的冰川有47條，平均海拔范圍為3 000～6 000 m，主峰孔嘎峰海拔6 490 m。

圖1 研究區(qū)位置（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 數(shù)據(jù)與方法

2.1 數(shù)據(jù)來(lái)源

氣象數(shù)據(jù)來(lái)源于美國(guó)國(guó)家氣候數(shù)據(jù)中心（ftp：//ftp.ncdc.noaa.gov/pub/data/noaa/isd-lite/）和中國(guó)氣象數(shù)據(jù)網(wǎng)（http：//data.cma.cn）。由于研究區(qū)沒(méi)有氣象站點(diǎn)，通過(guò)選取周邊鄰近地區(qū)8個(gè)氣象站點(diǎn)，分別為丁青、類烏齊、昌都、洛隆、波密、八宿、察隅和林芝，利用標(biāo)準(zhǔn)克里金插值獲取研究區(qū)數(shù)據(jù)，時(shí)間范圍為1984—2019年?，F(xiàn)代冰川數(shù)據(jù)來(lái)源于國(guó)家冰川凍土沙漠科學(xué)數(shù)據(jù)中心（www.ncdc.ac.cn）第二次冰川編目資料與Randolph Glacier Inventory（http：//www.glims.org/RGI/）數(shù)據(jù)相互補(bǔ)充。

2.2 氣溫與固態(tài)降水量提取

臨界氣溫法區(qū)分固液態(tài)降水是冰川物質(zhì)平衡模型中最常用的方法。主要包括雙臨界氣溫分離法與單臨界氣溫分離法。雙臨界氣溫分離法只適用于我國(guó)的一些干旱地區(qū)，單臨界氣溫法具有普遍適用性且具有較高精度［35-37］。本研究中采用單臨界氣溫法［38］。

式中：Ps為月固態(tài)降水量（mm）；P為月降水量（mm）；T為月平均氣溫（℃）；Tl為液態(tài)降水臨界氣溫值（℃）；Ts為固態(tài)降水臨界氣溫值（℃）。為獲取研究區(qū)固態(tài)降水?dāng)?shù)據(jù)，首先參考康爾泗等［38］與張?zhí)珓偟龋?9］臨界氣溫方法劃分固液態(tài)降水的方案，分別取0℃和2℃作為固液降水的劃分閾值。研究區(qū)氣象資料顯示，冰川作用現(xiàn)代冰川區(qū)固態(tài)降水臨界氣溫為-4.2～-0.6℃（圖2，折線），平均固態(tài)降水量為140～614 mm（圖2，條帶）。

圖2 1984—2019年臨界氣溫與固態(tài)降水量變化Fig.2 Variations of critical air temperature and solid precipitation from 1984 to 2019

2.3 ELA計(jì)算

積累區(qū)面積比率法（accumulation area ratio，AAR）與面積-高程平衡率法（area-altitude balance ratio，AABR）是目前最廣泛使用的冰川ELA計(jì)算方法［40-41］?；贏AR與AABR法基本原理，本研究使 用Pellitero等［41］開(kāi) 發(fā) 的 自 動(dòng) 計(jì) 算ELA工 具（https：//github.com/cageo/Pellitero-2015）。該工具計(jì)算ELA過(guò)程主要包括兩個(gè)部分：輸入需要重建冰川表面的數(shù)字高程模型（digital elevation model，DEM）和定義AABR與AAR值，即可快速計(jì)算出單條冰川或多條冰川數(shù)據(jù)集。

具體操作步驟如下：①輸入需要重建冰川表面的DEM。設(shè)置用于冰川面積計(jì)算的等高線區(qū)間，本研究設(shè)置為50 m。②定義AABR與AAR值。為獲取適合本研究區(qū)的AABR與AAR值，需要對(duì)其進(jìn)行確定。AAR值的確定利用Kern等［42］系統(tǒng)總結(jié)分析了全球不同地區(qū)的46條山谷和冰斗冰川的面積S與其穩(wěn)定狀態(tài)下的AAR值之間的關(guān)系式，為AAR=0.0648lnS+0.483。由此，計(jì)算得出本研究區(qū)域AAR值 為0.74±0.02。AABR值 的 確 定 是 基 于Rea［43］系統(tǒng)評(píng)估了全球冰川的物質(zhì)平衡提出的適用于全球范圍內(nèi)的BR平均值為1.75±0.71。鑒于同一地區(qū)有多個(gè)冰川，Benn等［44］建議對(duì)一個(gè)BR范圍進(jìn)行計(jì)算，然后選擇ELA標(biāo)準(zhǔn)差最低的一個(gè)。利用該方法，計(jì)算得出本研究區(qū)的BR均值為1.78±0.68。③使用Osmaston［45］提出的最低標(biāo)準(zhǔn)差方法來(lái)確定冰川最佳ELA。返回結(jié)果為“ELA_val?ues_AAR_and_AABR.txt”。

3 結(jié)果與討論

3.1 ELA重建模型與驗(yàn)證

Singh等［46］利用多元線性回歸統(tǒng)計(jì)分析法對(duì)1979—2013年喜馬拉雅山脈的納拉杜冰川平衡線高度的統(tǒng)計(jì)估算及其趨勢(shì)分析，Zekollari等［47］使用該方法描述了瑞士的Morteratsch冰川的表面質(zhì)量平衡，Sakai等［48］利用多元線性回歸統(tǒng)計(jì)分析了亞洲高山冰川對(duì)近期氣候變化的響應(yīng)?；谇叭搜芯?，利用編目資料數(shù)據(jù)與研究區(qū)已有相關(guān)遙感影像解譯數(shù)據(jù)［49］。利用ELA calculation toolbox計(jì)算工具，對(duì)研究區(qū)內(nèi)1990年、1999年、2011年及2019年的ELA進(jìn)行計(jì)算（表1）。以氣象站點(diǎn)數(shù)據(jù)（氣溫與固態(tài)降水量）和ELA數(shù)據(jù)作為輸入?yún)?shù)，建立1984—2019年間念青唐古拉山中段現(xiàn)代冰川ELA變化多元線性回歸方程。

表1 念青唐古拉山中段ELA計(jì)算結(jié)果Table 1 Calculation results of ELA in middle section of the Nyainqêntanglha Mountains

式中：X1為月平均固態(tài)降水量（mm）；X2為月平均氣溫（℃）。

利用有冰川紀(jì)錄的年份對(duì)1984—2019年間研究區(qū)冰川ELA模擬結(jié)果進(jìn)行驗(yàn)證，冰川編目資料［32］記 錄 了1989—1990年、1999—2000年 及2010—2011年研究區(qū)冰川變化特征，對(duì)應(yīng)ELA分別為5 140.1 m、5 403.5 m和5 408.8 m，利用上述方程［式（2）］計(jì)算ELA結(jié)果分別為5 182.5 m、5 375.0 m和5 390.0 m，ELA模擬值與測(cè)量值具有高度一致性，確定系數(shù)為0.9977。

同時(shí)，為進(jìn)一步對(duì)比驗(yàn)證模擬結(jié)果的有效性，通過(guò)對(duì)氣象站1984—2019年的年平均氣溫、6—8月平均氣溫、年降水量和6—8月平均降水量進(jìn)行分析，參考我國(guó)西部17條現(xiàn)代冰川的物質(zhì)平衡以及氣象觀測(cè)資料［50］，建立冰川物質(zhì)平衡線處年降水量（PELA）與6—8月平均氣溫（TSO）的擬合方程［式（3）］，并結(jié)合綜合因子法［51］［式（4）］，來(lái)計(jì)算念青唐古拉山中段現(xiàn)代冰川ELA。

選取距離研究區(qū)域較近的丁青與洛隆氣象站數(shù)據(jù)，站點(diǎn)海拔分別為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），計(jì)算得出海拔增加量ΔH=1 495 m，念青唐古拉山中段ELA為5 368 m，模擬結(jié)果為5 282 m，相對(duì)誤差僅為0.16%，兩者基本一致。

3.2 ELA年際變化及特征

通過(guò)重建念青唐古拉山中段1984—2019年間研究區(qū)冰川ELA。結(jié)果顯示，研究區(qū)多年平均ELA為5 389 m a.s.l.，總體呈上升趨勢(shì)（圖3，直線），90%置信區(qū)間（虛線）的確定系數(shù)為0.7830。ELA年變化研究顯示出波動(dòng)變化特征（圖3，折線），ELA波動(dòng)范圍為5 360～5 420 m a.s.l.，波動(dòng)幅度約為60 m，平均上升速率為1.57 m·a-1。ELA年變化量正值與負(fù)值分別為18年和17年，基本呈現(xiàn)出交替狀態(tài)（圖3，條帶）。其中，ELA年際變化量正值總和為319.9 m，負(fù)值變化總和為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變化

利用流域自然分區(qū)法，選取HydroSHEDS and WaterGAP v2.2數(shù) 據(jù) 庫(kù) 中 的BasinATLAS_Cata?log_v10數(shù)據(jù)，將研究區(qū)分為三個(gè)流域［圖1（b）］，分別為霞曲流域［圖4（a）］、易貢藏布流域［圖4（c）］和麥曲流域［圖4（e）］。霞曲流域、易貢藏布流域和麥曲流域現(xiàn)有現(xiàn)代冰川分別為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 三個(gè)流域冰川表面高程與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變化總體都呈上升狀態(tài)，但不同流域ELA變化具有差異性。其中，麥曲流域ELA變化幅度最大，平均上升335 m，霞曲流域和易貢藏布流域分別上升265 m和314 m，上升速率分別9.57 m·a-1、7.57 mm·a-1和8.97 m·a-1。同時(shí)，統(tǒng)計(jì)分析顯示霞曲流域，易貢藏布流域和麥曲流域年平均ELA高程分別為5 335 m a.s.l.、4 987 m a.s.l.和5 317 m a.s.l.，呈現(xiàn)出西北—東南方向降低，北坡高于南坡。

受印度季風(fēng)，流域位置及冰川朝向等因素影響，各流域現(xiàn)代冰川區(qū)水熱組合的差異化特征，進(jìn)而影響著冰川ELA變化。從流域分布來(lái)看，霞曲流域冰川區(qū)位于研究區(qū)西北方向，冰川朝向以NW為主，易貢藏布流域位于研究區(qū)東南方向，冰川朝向以SE為主，麥曲流域位于研究區(qū)東北方向，冰川朝向以NE為主。季風(fēng)路徑（圖5）顯示出易貢藏布流域處于印度季風(fēng)迎風(fēng)坡，獲得的水汽相較于處于背風(fēng)坡雨影區(qū)的麥曲流域和霞曲流域更多。

從各流域冰川區(qū)氣溫與降水狀況看，冰川區(qū)整體氣溫介于-5.6～-3.9℃，平均值為-4.8℃［圖5（a）］；降水介于502～564 mm，平均值為547 mm［圖5（b）］。易貢藏布流域冰川區(qū)氣溫為-5.5～-3.9℃，均值為-4.9℃；降水為502～554 mm，均值為547 mm；麥曲流域冰川區(qū)氣溫為-4.6～-4.1℃，均值為-4.3℃；降水為510～555 mm，均值為539 mm；霞曲流域現(xiàn)代冰川區(qū)氣溫為-5.5～-4.4℃，均值為-4.8℃；降水為549～557 mm，均值為553 mm。對(duì)比各流域冰川ELA和冰川區(qū)氣溫降水組合，發(fā)現(xiàn)相較于處于背風(fēng)坡的麥曲流域，處于迎風(fēng)坡的易貢藏布流域的水熱組合狀況更優(yōu)，更利于冰川發(fā)育。

從冰川規(guī)模來(lái)看，易貢藏布流域冰川規(guī)模最大，最大冰川長(zhǎng)度為33 km［圖5（d）］，總面積為1 238.809 km2，麥曲流域次之，霞曲流域最小，最大冰川長(zhǎng)度分別為24 km［圖5（b）］和19 km［圖5（f）］，總面積分別為760.817 km2和716.982 km2［33］。相關(guān)研究表明，冰川對(duì)氣候變化的反饋?zhàn)饔茫?2］以及冰川變化對(duì)氣候的響應(yīng)具有滯后效應(yīng)［53］，使得盡管相較于易貢藏布流域冰川區(qū)，霞曲流域水熱組合狀況更優(yōu)，但冰川規(guī)模相對(duì)較小，ELA變化對(duì)氣候響應(yīng)更敏感。這也是造成短時(shí)間尺度范圍同一氣候區(qū)域冰川ELA變化差異性的重要原因。

圖5 標(biāo)準(zhǔn)克里金插值處理后的冰川區(qū)氣溫與降水Fig.5 Air temperature and precipitation in glacier area after standard Kriging interpolation

3.4 ELA變化對(duì)氣候變化的響應(yīng)

姚檀棟等［54］對(duì)高亞洲冰川總體退縮研究表明，氣溫升高引起的冰川消融量增加和粒雪線上升，是近百年來(lái)全球山地冰川呈現(xiàn)總的退縮趨勢(shì)的主要原因。同時(shí)，Caidong等［55］研究發(fā)現(xiàn)，當(dāng)溫度變化為1℃時(shí)，ELA變化為140 m。段克勤等［56］對(duì)小冬克瑪?shù)妆‥LA的氣候敏感性進(jìn)行了試驗(yàn)，發(fā)現(xiàn)在氣候狀態(tài)基礎(chǔ)上，假設(shè)其他氣候因子不變，當(dāng)氣溫增加1℃將引起冰川ELA大幅度上升160 m，而增加20%的降水量引起的ELA僅降低10 m，這也說(shuō)明冰川ELA對(duì)氣溫變化較敏感。

本研究通過(guò)1984—2019年間平均氣溫、固態(tài)降水量與冰川ELA統(tǒng)計(jì)，對(duì)念青唐古拉山中段冰川ELA變化對(duì)氣候的響應(yīng)進(jìn)行分析（圖6）?？傮w來(lái)看，ELA隨氣溫和固態(tài)降水量變化都呈顯著相關(guān)性，確定系數(shù)R2分別為0.5873和0.9540。同時(shí)，ELA響應(yīng)氣溫與固態(tài)降水量變化結(jié)果顯示，研究區(qū)冰川ELA主要受氣溫變化控制，隨氣溫變化1℃，ELA波動(dòng)幅度為126.02 m，而受固態(tài)降水量變化影響，ELA波動(dòng)僅為0.2492 m。并且，隨著氣溫的持續(xù)升高，念青唐古拉山中段冰川ELA上升，冰川將繼續(xù)處于退縮狀態(tài)。

圖6 ELA隨氣溫與固態(tài)降水量的變化Fig.6 Variation of ELA with air temperature and solid precipitation

4 結(jié)論

通過(guò)建立念青唐古拉山中段冰川作用區(qū)多元線性回歸方程，對(duì)1984—2019年間ELA變化進(jìn)行了重建。主要結(jié)論如下：

（1）研究時(shí)段內(nèi)平均ELA為5 360 m a.s.l.。

（2）ELA總體呈上升趨勢(shì)，上升速率為1.57 m·a-1。ELA年變化顯示出波動(dòng)變化特征，波動(dòng)范圍為5 360～5 420 m a.s.l.，上升幅度為60 m。

（3）受流域位置、冰川朝向及冰川規(guī)模影響，研究區(qū)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高程對(duì)比顯示出西北—東南方向降低，北坡高于南坡的變化特征。

（4）冰川區(qū)多年ELA變化的氣候響應(yīng)分析顯示，ELA變化主要受氣溫控制，隨氣溫變化1℃，冰川ELA總體波動(dòng)126.02 m。