山地冰川變化的數值模擬及其在亞洲高山區的應用

段克勤，石培宏，何錦屏

（陜西師范大學地理科學與旅游學院，陜西 西安 710119）

0 引言

在全球氣候日益變暖背景下，全球山地冰川普遍呈減薄后退狀態［1］，引發社會各界對冰川這一固體水庫的儲存量及其變化的關注。亞洲高山區（青藏高原及其周邊高海拔地區），因其巨大的海拔高度和特殊的地理位置，擁有中低緯度最多的冰川。這些冰川是亞洲大江大河的重要水源地，冰川退縮的直接后果就是淡水資源的減少［2-3］，而且冰川的快速退縮還引起了諸如冰湖潰決等自然災害［4-5］，以及河川徑流的季節性改變［6-8］。亞洲高山區正處于快速變暖過程中［9-10］，造成冰川大幅度退縮［11-15］。在當前全球持續變暖背景下，明確亞洲高山區冰川未來變化趨勢及幅度，對揭示亞洲高山區水塔失衡失穩的科學內涵、過程機理、減緩和應對這一演化趨勢，具有十分重要的科學意義和應用價值。

中國科學家自20世紀50年代開始研究亞洲高山區冰川以來，不論是定點觀測研究，還是利用衛星遙感資料進行空間上的研究，迄今已取得了巨大成績，特別是中國第一和第二次冰川編目的完成［15］，以及長期觀測冰川變化的數據積累［16］，在亞洲高山區冰川變化方面做了大量的工作，擴展和豐富了亞洲高山區的冰川研究內容，為深入研究亞洲高山區冰川變化奠定了基礎。然而，同國際同行相比，我國的冰川研究相對偏重于實地觀測和統計分析，而缺乏對冰川變化的機理與理論研究，特別是利用數值模擬方法對冰川變化的研究，與發達國家還存在較大差距。

隨著研究的深入，為明確冰川-氣候-水文之間的定量聯系，以及揭示冰川對氣候變化響應機理，并預測冰川的未來變化趨勢，應用物理模型對冰川變化進行模擬和預測是必然的發展趨勢。在老一輩冰川學家施雅風和謝自楚等，以及秦大河、姚檀棟和丁永建等冰川學家的大力推動下［17-19］，開展冰川變化的數值模擬研究已是我國冰凍圈研究的重要領域和趨勢［18，20］。雖然國內學者最近在冰川變化模擬方面取得了顯著研究進展［20-25］，然而如何利用完全的物理模型定量的對冰川變化進行全面深入的分析，并在此基礎上進行更科學的預測，與發達國家相比，還存在較大差距。因此，推動冰川變化的數值模擬研究，不僅可以深刻揭示山地冰川變化規律，而且可以推動中國的冰川研究向數值化方向發展，以縮小與發達國家的差距，相關研究急需加強。

本文介紹了冰川數值模擬的基本原理和方法，綜述了近年來亞洲高山區冰川數值模擬的研究現狀，基于目前的研究與認識，提出亞洲高山區冰川數值模擬研究的薄弱環節并進行了展望。

1 冰川變化數值模擬的基本原理

冰川是特定氣候在特殊地形下的產物。作為一個巨大的非牛頓流體，重力使其總處于持續運動狀態。由物質守恒定律，可以得到冰川表面厚度變化的控制方程：

式中：H為冰川厚度；V是冰川運動速度；t為時間；b為物質平衡。

式（1）表示冰川某處的厚度變化由兩部分決定：一是發生在此處冰川表面的物質平衡，二是此處冰流通量的散度變化。在重力作用下冰川總處于運動之中，造成冰川物質處于不斷調整和分配中。當冰川物質平衡為零時，消融區損失的物質會由冰川上部積累區向下輸送物質而得到補充，使冰川總處于動態平衡中。但當冰川消融加劇，消融區擴大而積累區縮小時，向下輸送的物質會減少，冰川消融區因得不到足夠物質補給，會造成冰川退縮減薄，反之則冰川變厚并引起末端前進（圖1）。因此，冰川對氣候的響應過程就是在地形約束下，冰川物質平衡過程和冰川動力過程相耦合的動態演變過程，只有從冰川表面的物質平衡和冰川流動兩方面考慮，才能完整刻畫冰川的動態變化。

圖1 冰川動態變化示意圖（A.縱剖面圖；B.俯視圖）Fig.1 Sketch map of glacier dynamic change（A.vertical section view；B.bird view）

基于此，模擬冰川變化的經典方案如圖2所示。要對一個冰川進行全面的研究，需從三個方面入手：一是冰川變化的氣候背景，即冰川及其周邊的氣象條件是什么；二是冰川對氣候響應過程，這不僅要充分考慮發生在冰川表面的物質平衡過程，還要考慮冰川流動引起的物質輻合輻散過程；三是冰川變化的影響問題，如冰川融水對河川徑流的影響。

圖2 冰川變化模擬的經典方案流程圖Fig.2 Flow chart of glacier variation simulation

2 冰川物質平衡模擬進展與方法

在式（1）中，冰川物質平衡b是指冰川表面各種相態水收支的代數和，是冰川發育水熱條件的綜合反映，是銜接冰川變化與氣候變化之間關系的重要參量。冰川物質平衡觀測主要基于冰川表面花桿/雪坑法，在過去幾十年對中國西部典型冰川進行了大量的物質平衡觀測，積累了寶貴的資料［11，16］。隨著多源遙感資料日益豐富，利用大地測量法獲取冰川變化信息，已成為研究冰川物質平衡變化的主要手段［16，26］。近年來發展的三維激光掃描和無人機等新技術，也為研究冰川物質平衡變化提供了新的重要手段［20］。

但要深入研究冰川物質平衡對氣候的響應及未來變化趨勢，必須依賴于基于物質平衡變化過程的數學模型。目前主要有兩類模型用來研究冰川物質平衡：一是基于參數化的溫度指數模型［27-28］（以下稱：溫度指數模型）；二是詳細描述冰川表面物質能量平衡過程的模型（以下稱：能量平衡模型）［21，29］。能量平衡模型所需要參數較多且不易觀測，目前只在觀測要素全面的個別冰川得到應用，國內一般都是依據氣象觀測計算單點的物質平衡［30］。相比較，溫度指數模型參數化程度較高，應用比較廣泛［31-32］。但利用物質能量平衡模型可深入揭示冰川物質平衡變化的原因及其對氣候響應方式。為了更加準確描述冰川的時空變化，如考慮地形因素及其對太陽輻射的遮蔽等，國際上已發展出分布式物質能量平衡模型［29，33］。模型也從冰川表面物質平衡的簡單模擬，向考慮冰川表面不同介質（例如積雪、粒雪、冰川冰和表磧物等）和不同層位（將表面能量平衡和冰面以下雪冰過程）相耦合的復合模型方向發展，使物理基礎更加健全和精細，如COSIMA［34］和COSIPY［35］模型。當前國際上另一個發展動態是利用再分析資料（例如EAR5、NCEP等）和區域氣候模式結果作為輸入驅動數據，通過求解物質能量平衡方程，模擬分析一個流域或更大范圍的冰川物質平衡變化，而不局限單個冰川［36］。并在此基礎上，基于氣候預估數據（如CMIP6），利用物質能量平衡方程實現評估區域尺度冰川變化［37］。這方面的研究非常有助于從整體上認識大范圍和流域尺度內冰川變化及其水文水資源變化。

國內能量物質平衡的研究源自20世紀50年代開始的山地冰川冰/雪能量平衡與轉化研究［38］。近年來隨著參數化方案的不斷完善，能量物質平衡模型在亞洲高山區快速被應用。如在珠穆朗瑪峰東絨布冰川［39］、祁連山七一冰川［40］和老虎溝12號冰川［23］、念青唐古拉山西布冰川和扎當冰川［41］、帕米爾［30］以及藏東南帕隆藏布4號冰川［42］等開展了基于能量物質平衡模型的冰川模擬研究。以上工作驗證了能量平衡模型能夠精準地模擬冰川消融過程。但是目前國內物質能量平衡的研究，主要是以單點或單條冰川為主。雖然基于度日模型的流域/區域研究已經在祁連山北大河流域和阿爾泰山地區已有應用［43-44］，但對整個流域或一個區域內冰川分布式物質能量平衡模擬研究工作還尚未系統開展。值得盡快開展這一方面的研究，以縮小與國際同行之間的差距，并為評價冰川變化對水資源的影響奠定堅實基礎。

能量物質平衡模型可用下面的公式描述［29］：

式中：b為物質平衡；QM為冰川表面消融耗熱；Lm為冰的融化潛熱；f為冰川融水再凍結量；QE為升華或蒸發耗熱；L為升華或蒸發潛熱；Ps為固態降水；G為太陽輻射；α為反射率；Lnet為凈長波輻射通量；QH為顯熱通量；QL為潛熱通量；QG為地面熱通量；QR為降水帶來的熱量。這些物理量通常都難以直接獲得，只能通過參數計算，目前這方面的理論研究比較充分，有完備的計算方法，關鍵是如何獲得研究區氣象要素的空間分布。對分布式物質能量平衡模型而言，除氣象資料外，還必須考慮諸如地形對輻射的遮蔽等因素［21，29］。

驅動分布式物質能量平衡模型，需要空間變化數據。氣象臺站資料比較稀疏，而同化后的再分析資料如NCEP和EAR5，分辨率又比較低，不能直接用于分布式物質能量平衡模型。因此，需要把氣象資料和再分析資料降尺度到研究區域。方法一，把氣象站資料和再分析資料利用統計方法（例如偏差校正）插值到空間格點［29，31，45］；方法二，利用氣候模式對再分析資料，以及全球氣候模式的模擬結果進行動力降尺度，把各氣象要素（溫度、降水、輻射、濕度、風速和氣壓）降尺度到研究區域［46-47］。受限于觀測資料和再分析資料自身不確定性的限制，降尺度（統計或動力）在面向單條冰川或流域/區域冰川模擬時存在一定的誤差［20，48］。但對于長時間尺度冰川物質平衡的模擬，使用降尺度后的數據能夠較好重現觀測結果［49-51］。

最后由數字化地形圖確定冰川范圍，劃分冰川區計算網格，確定每個網格的高程、坡向和坡度，結合網格化的各氣象數據，代入物質能量平衡方程，計算每一個格點的物質平衡變化，探討冰川物質平衡對氣候變化的響應。

3 冰川流動模擬進展

冰川流動會造成物質遷移，改變冰川的水熱、邊界條件，使冰川形狀不斷發生調整，并決定冰川末端的變化速率。因此，冰川的流動極大增加了冰川變化研究的復雜性。

冰川冰通常認為是不可壓縮、黏性的能導熱的非牛頓流體，則式（1）中的描述冰川的動量守恒方程為：

式中：P為壓力；η為動力學黏性系數；ρ為冰的密度。

在20世紀60—70年代，對冰川流動的研究主要采用淺冰近似模型，如采用一維淺冰近似流線模型模擬了挪威Nigardsbreen冰川末端變化［54］。20世紀80—90年代，二維冰流模型在山地冰川動力學過程中得到了廣泛使用［55-58］。進入21世紀后，隨著計算能力的提高和理論研究的深入，應用有限元方法求解三維的Stokes方程模擬冰流過程取得了進一步進展，如用非常復雜的三維冰流模型成功模擬了瑞士Rhonegletscher冰川和Grosser Aletschgletscher冰川的變化歷史，并對冰川變化進行了預測［59］。整體而言，盡管現有動力學模式都是建立在冰川流動定理與物質守恒原理之上，但使用了不同的關系式和簡化方法，區別在于是選擇全Stokes方程進行計算，還是采用全Stokes方程的簡化近似模型進行計算［60-64］。

在國內早期針對新疆天山烏魯木齊河源1號冰川進行過頻率響應［65］和主流線厚度簡單模擬［66］，以及對新疆伊犁河流域冰川的一些理論性探討研究［67］。近年來，段克勤等［22］利用有限元算法對二維的Stokes方程求解，計算了烏魯木齊河源1號冰川的流動速度，發現模擬結果與實測速度非常吻合。李慧林［68］利用“淺冰近似（SIA）”對天山烏魯木齊河源1號冰川進行了模擬和預測。

國內冰川流動模型研究近年最大的亮點是張通等［25］和王玉哲等［63］發展了PoLIM冰川流動模型，其核心控制方程為：

式中：u是水平速度；W是冰川寬度；s是冰川表面高程；η是黏性系數，定義如下：

即冰的有效黏性系數由冰川的流速場和流動系數A決定，而A由冰川溫度場（T）計算獲得：

而冰川溫度場（T）可由熱傳導擴散方程計算得到：

式中：K是熱傳導系數；ρ是冰川冰密度；c是冰的熱容量；w是垂向速度；E是冰體內部的形變熱。

該模型的優勢在于沿冰川流動方向考慮縱向應力梯度，可以模擬地形起伏。同時在橫剖面對側向拉力進行參數化，考慮了山體對冰川的阻力影響。此外，還考慮了熱力耦合模擬，可以獲得冰川內部的熱學特征。相比三維模型，所需輸入數據較少，適用于山地冰川的模擬。與其他動力學模型同樣面臨的困難在于，一些關鍵數據和參數化方案對于冰川動力過程十分重要，但很難觀測獲取，例如冰川底部運動速度、冰川底部溫度及含水量組成等。利用該模型研究了珠穆朗瑪峰的東絨布冰川的溫度和速度分布［25］，以及祁連山老虎溝12號冰川的熱力特征等工作［63］。此外，冷偉等［69］發展的三維Full-Stokes模式對格陵蘭冰蓋的模擬也取得了很好的結果。這些都極大促進了國內冰川流動模型的進一步應用和發展。

整體而言，冰川動力學研究需把流體力學和數值分析方法應用到冰川變化中，計算過程復雜，目前還不能達到普適的應用階段。與發達國家相比，我國冰川動力學的研究相對薄弱，尚未參與國際冰蓋和冰川模型比較計劃。從學科發展的趨勢和提高冰川學理論角度，國內對冰川的流動模型研究亟需加強，以滿足冰川變化預測研究的需求。

4 冰川變化預測研究進展

從冰川變化的理論角度，若由式（2）獲得冰川的物質平衡，由式（4）和（5）獲得冰川的運動速度，代入式（1），就可研究冰川的演化趨勢。但在實際應用中，由于氣候變化的不確定性和冰川變化過程的復雜性，對未來氣候情景下山地冰川的演化趨勢，仍存在較大的不確定性。例如IPCC第五次評估報告顯示在RCP2.6～RCP8.5情景下，至21世紀末全球不同規模山地冰川的消融量在15%～85%之間［70］，而最新IPCC第六次評估報告認為30%～80%的山地冰川將消融殆盡［1］。由此可見，對于冰川預測存在較大的不確定性。

在全球尺度上，2006年Raper等［71］首次利用度日模型模擬了全球1°×1°網格內的冰川物質平衡梯度，預估了21世紀全球冰川物質虧損及其對海平面的貢獻。Hock等［48］等利用溫度指數模型重建了1961—2004年全球冰川年時間尺度的物質平衡序列，以評估冰川在該時間段內對海平面上升的貢獻。Kraaijenbrink等［13］采用Huss模型的方法預估了全球氣溫較工業革命前上升1.5℃對亞洲高山區冰川變化的影響。Huss和Hock［8］預估了全球尺度下冰川變化對冰川區徑流的影響，并預估了對不同冰川區徑流峰值出現的時間。綜合現有研究看，目前公認有三個模型可以較好地模擬冰川在未來氣候情景下的變化，分別由Marzeion等［49］、Radi?等［50］和Huss等［51］開發。然而，由式（1）可見，因冰川流動會使物質從積累區向消融區遷移，使冰川對氣候的響應具有滯后性。從物質和能量守恒的角度，要模擬和預估冰川的真實變化，必須考慮冰川流動引起的冰川形態變化。雖然在Marzeion模型和Radi?模型中，采用“面積-體積-長度”比率法估算了冰川末端的進退，但并未嚴格考慮冰川厚度隨時間的變化，本質上缺乏冰川形態變化對冰川流動的反饋機制。相比較，Huss模型雖然由質量守恒角度對冰川每個高程帶的厚度和面積變化進行了參數化處理［72］，但由于設定的參數眾多，且部分參數物理機制不明確，造成大范圍的應用存在較大的不確定性。Clarke等［73］基于“淺冰近似（SIA）”理論，提出了一套考慮三維冰流過程的模擬框架，但地形等邊界輸入條件要求過高，模型初始化比較困難，普適性并不理想。Zekollari等［74］在Huss等［51］開發的GloGEM模型的基礎上，增加了冰流模型模擬冰川動力過程，來模擬歐洲阿爾卑斯山冰川對未來氣候變化的響應。Rounce等［75］開發了PyGEM模型來模擬預估21世紀冰川徑流的變化。雖然上述模型可以實現對全球尺度冰川變化的模擬，但由于模型之間的物理基礎和參數化有較大的差別，造成不同模型的預測結果不盡相同。為促進冰川模型的發展和減小冰川變化預測的不確定性，Hock等［76］于2015年發起了“冰川模型比較計劃（Glacier Model Intercomparison Project，簡稱GlacierMIP）”，截至2022年已進行到了第三階段（GlacierMIP3）。

冰川厚度作為冰川形態的關鍵參數，是冰川動力學數值模擬的重要輸入參數，但獲取大范圍的冰川厚度存在較大困難，通常采取參數化或反演方案估算冰川的厚度。為比較不同冰川厚度反演模型的精準度和減小不確定性，Farinotti［77］2017年發起了“冰厚反演模型比較計劃（Ice Thickness Models Intercomparison eXperiment，簡稱ITMIX）”，2020年業已推行到第二階段（ITMIX2）。

國內早在2001年施雅風先生使用經驗統計方法預估在氣候變暖情景下，祁連山北麓河地區小冰川將迅速衰退，流域內1 km2左右的小冰川將會在近二三十年大量消亡［78］。但國內對亞洲高山區冰川整體變化預測的研究較為缺乏，且現有缺乏物理基礎的模型在數值模擬與預測方面還存在較大的不確定性。除了模型自身的不確定性外，驅動數據的質量和可靠性對于模擬結果有著直接的影響。特別高海拔地區的觀測非常有限。對高分辨率、適應性好的ERA5等再分析資料和CMIP5等預估數據進行降尺度，驅動物質平衡和冰川流動模型是最行之有效的手段［51］。國際上，近年來已經在努力開展這方面的研究［75，79］，國內類似研究還比較欠缺。特別是對流域尺度冰川的數值模擬與預測是國內工作的難點［20］，使用再分析資料是一個很好的出口。此外，比CMIP5考慮物理機理更加完善的CMIP6數據的發布，對于山地冰川未來變化的準確預估也提供了新的契機。

5 研究展望

冰川變化是一個極其復雜的過程，冰川數值模擬則是對該復雜過程的高度概化?；诋斍皣鴥韧獾难芯楷F狀，以及未來的學科發展趨勢，可以看出發達國家對冰川的研究已經完全數值化，從冰川物理角度全方位進行數值模擬，使冰川學研究上升到一個新的理論高度，相關研究在冰川變化、冰川水資源和海平面上升方面有重大的應用價值。與其相比，我國在山地冰川模型的應用和建模方面開展了初步的探索工作，總體上還處于模型的應用或者使用大型計算軟件（如COMSOL、Elmer/Ice）對冰川進行數值模擬。這些方法的局限性很明顯：①僅能對單個冰川進行模擬；②與其他模式或者模塊耦合的可能性很小。因此，山地冰川的數值模擬研究必須加強。為此，提出以下展望：

（1）度日模型存在較大的空間異質性，且不能夠解釋冰川變化的物理機制和過程，引入能量物質平衡，特別是耦合多層過程的能量平衡過程并考慮全地形影響的模型，對于亞洲高山區冰川物質平衡過程至關重要。

（2）亞洲高山區山地冰川規模遠小于極地冰蓋，基于冰川變化的物理基礎，“淺冰近似原理”并不適用于山地冰川。因此，高階Stokes方程的求解是山地冰川模型建立的首選。

（3）明確對山地冰川未來變化趨勢做出準確預測，是我們長期的戰略任務和目標。通過耦合分布式物質能量平衡模型和冰川流動模型，對冰川進行全物理過程的數值模擬，才能深刻探討冰川變化對氣候的響應機理并預測冰川的未來變化，同時推動中國冰川學研究向數值化方向發展。

（4）在氣候情景下，流域尺度或區域尺度冰川數值模擬研究是目前工作的重點和急需突破的方向。在現有冰川編目、冰川觀測的基礎上，如何將分布式物質平衡模型和簡化的動力學模式結合，預測和評估未來冰川變化的可能狀態及其對冰川水資源的影響，為政府決策、國際氣候談判提供科學數據。在冰川變化數值模擬方面做出一流的研究成果，逐步確立我國在該研究領域的國際領先地位。

（5）冰川數值模擬的門檻相對較高，造成當前國內進行冰川數值模擬工作的科研人員極其不足。從業者不僅要具備冰川變化的基本知識，還需具備深厚的數學和物理知識，特別是流體力學、計算數學和編程方面的知識。對此，一方面要加強研究生的數學和物理知識，另一方要有組織地進行數值模式的研發。