中國極地冰蓋數值模擬研究進展

趙勵耘，唐學遠，張通，冷偉，艾松濤，劉成彥，王玉哲，劉巖，岳超

（1.北京師范大學全球變化與地球系統科學研究院，北京 100875；2.自然資源部極地科學重點實驗室中國極地研究中心，上海 200136；3.北京師范大學地表過程與資源生態國家重點實驗室，北京 100875；4.中國科學院數學與系統科學研究院，北京 100190；5.武漢大學中國南極測繪研究中心，湖北 武漢 430079；6.南方海洋科學與工程廣東省實驗室（珠海），廣東 珠海 519080；7.山東師范大學地理與環境學院，山東 濟南 250014）

0 引言

極地冰蓋（圖1）對全球熱量平衡和海平面上升具有重要的影響。IPCC第六次評估報告指出，氣候變化評估中最大的不確定性之一來自對極地冰蓋對未來海平面上升貢獻量的估計，格陵蘭冰蓋物質損失對20世紀全球海平面變化貢獻了29%，南極冰蓋物質損失區域差異性大，其中西南極冰蓋底部由于受到海洋環流影響融化速率相對較快［1-2］。此外，東南極尤其是威爾克斯盆地（Wilkes Basin）也發現潛在的不穩定性［3］。

圖1 南極冰蓋（a）；北極的格陵蘭冰蓋、斯瓦爾巴群島和冰島（b）［在圖中標注了我國建立的臺站（長城站、中山站、昆侖站、泰山站、羅斯海新站、黃河站）以及本文提到的一些研究區域或冰川位置（威爾克斯盆地、伊麗莎白公主地、埃默里冰架、普利茲灣、Totten冰川、雅各布港冰川、瓦特納冰帽）。Dome A與昆侖站非常近，Austre Lovénbreen冰川和Pedersenbreen冰川與黃河站距離非常近，不再標注］Fig.1 Antarctic Ice Sheet（a）；Greenland Ice Sheet，Svalbard Islands and Iceland in the Arctic（b）［The stations established by China（Great Wall Station，Zhongshan Station，Kunlun Station，Taishan Station，Ross Sea New Station and Yellow River Station）and some study regions or glacier locations mentioned in this paper（Wilkes Basin，Princess Elizabeth Land，Amery ice shelf，Prydz Bay，Totten Glacier，Jakobshavn Glacier，Vatnaj?kull ice cap）are marked in the figure.Dome A，very close to Kunlun Station，Austre Lovénbreen Glacier and Pedersenbreen Glacier almost next to Yellow River Station，are not marked］

由于極地冰蓋氣候惡劣、難以接近，評估當前極地冰蓋變化的主要手段是基于遙感監測，在部分區域結合野外實地測量。遙感技術可以對冰蓋實現全覆蓋、連續性監測。地面實地測量難以大范圍實施，但是可以更精確地探測冰蓋內部構造和底部地形等關鍵物理信息。遙感技術和地面實地測量等多源觀測為研究極地冰蓋變化提供了必備的數據資料（表1）。

表1 冰蓋關鍵動力過程、模型中的相關變量和主要觀測手段Table 1 Key dynamic processes of ice sheet，relevant variables in the model and main observational approaches

基于多源觀測數據作為模式邊界條件和外部驅動，可以通過數值求解冰流動力學模型，詳細描述冰流的運動或預測冰蓋的未來演化。冰蓋數值模擬通常包括兩種：一種是診斷模擬（diagnostic modeling），用以診斷冰蓋在穩態假設和特定邊界條件及參數下的動態平衡狀態；也可以用來測試冰蓋對某些參數或邊界條件的敏感性以及微小變化對模型行為的影響，從而確定可能反饋機制的性質。診斷模型可用于提高對控制特定冰流行為的參數的理解，或研究一般冰蓋中一個或多個物理參數與邊界條件的重要性。另一種是預估模擬（prognostic modeling），用以描述冰蓋運動的過程，通過調整模型并使其包括不同的物理過程，以便更好地將模擬預估與實際情況相匹配，從而更好地預測冰蓋的未來演變。預估模型既可直接用于預估（或預測）冰蓋的未來變化；也可用于重建冰蓋的歷史演化過程，從而有助于準確理解冰蓋未來可能發生的變化。

如果要考慮到作用在冰上的所有應力分量，full-Stokes模型可能是最合適的工具。在某些情形下忽略部分應力分量，可以得到簡化近似方程。基于full-Stokes模型或者簡化近似的冰流模型，各國研究團隊陸續研發了一些用于模擬冰蓋狀態或者預估其未來變化的數值求解方案，即冰流模式［4］。雖然我國在冰蓋數值模擬方向的研究起步較晚，但是目前也形成了三個自主研發的冰流模式——多溫型陸地冰模型PoLIM［5］、三維并行冰蓋演化模型［6］和高精度并行有限元冰蓋模型［7-9］。PoLIM為二維Blatter-Pattyn高階近似熱-力耦合冰流模型，以焓方法計算冰川熱力場，包含冰下水文模型，采用有限差分方法求解計算。張懷等開發了三維并行冰蓋演化模型，進行了格陵蘭冰蓋的模擬演化［6］。冷偉等研發的三維full-Stokes和高階近似冰蓋模式求解器，使用了適用于大規模數值模擬的快速求解方法［7-8］，并采用了保持局部質量守恒的數值離散方法以提高計算精度［9］，具有熱-力耦合和處理觸地線進退過程的功能［10-11］。

冰蓋與大氣、基巖、海洋接觸，在邊界處存在復雜的過程，包括冰蓋表面物質平衡、冰蓋水文過程、冰下水熱環境、冰架崩解、冰架底部消融和再凍結等［12-13］（表1、圖2）。這些過程都會對冰流系統的物質平衡和穩定性產生重要的影響，是最近冰蓋動力學研究領域的熱點和難點。目前在理論和觀測上對上述關鍵動力過程和機制的認識都存在很大不足，嚴重制約了模式的發展和模擬結果的置信度。盡管國際冰蓋模式發展時間較長，但是其對上述過程多采取簡單的參數化處理或者尚未考慮，其模擬結果仍存在很大的不確定性［4］。特別的，近年來開展的冰蓋-海洋耦合模式比較計劃，不同耦合模式對海洋強迫的模擬結果差異很大。因此，模式的改進需要以冰蓋典型區域長期持續的高精度精細測量作為基礎支撐以提升對冰蓋關鍵過程和機制的理解。

圖2 冰蓋關鍵動力過程示意圖（本示意圖的地形更適合西南極冰蓋，而不適用于東南極和格陵蘭冰蓋大部分地區）Fig.2 Schematic diagram of key dynamic process of ice sheet（The topography in this diagram is more suitable for West Antarctic Ice Sheet than most region in East Antarctic and Greenland Ice Sheet）

前些年已有數篇針對極地冰蓋數值模擬的評述文章［14-15］。近十余年來，我國冰川學家在極地冰蓋數值模擬方面產出了不少新的研究成果。本文簡要介紹冰蓋數值模擬方法，進而梳理了我國在極地冰蓋數值模擬方面的新進展，并對未來的研究前景進行了展望。

1 冰蓋數值模擬方法簡述

開發冰蓋動力學數值或解析模型的主要用途是為了更好地了解冰蓋的流動行為以及它們如何對外部環境的強迫做出響應。這些模型通常基于一些描述冰流的基本定律或假設，經過簡化后可以得到解析或數值解。冰蓋動力學模型包括動力學方程和初始條件以及邊界條件。動力學方程主要包括冰蓋的質量、動量、能量守恒定律和本構方程Glen流動定律。目前最為完整的動力學方程是熱-力耦合的三維full-Stokes方程，它將描述質量、動量守恒的Stokes方程和描述能量守恒的熱傳導方程進行耦合，適用于各種情形的冰蓋模擬。在Stokes方程的基礎上進行不同程度的簡化近似，可以得到靜水壓力近似（hydrostatic approximation）、一階近似（first order approximation）、淺冰蓋近似（shallow ice approximation）、淺冰架近似（shallow shelf approximation）等［16-17］。在full-Stokes方程中忽略垂直切應力τzx和τzy沿著運動方向的導數，可以得到靜水壓力近似。在靜水壓力近似基礎上，忽略速度豎直分量的水平導數，即假設速度豎直分量vz僅依賴于豎直坐標z和時間t，vz=vz（z，t），得到一階近似。在一階近似基礎上進一步簡化，把冰蓋看作是平行底面的剪切流，水平切應力τxz和τyz對水平方向的動量平衡起到重要作用，忽略法向應力偏量和豎直切應力τxy，可以得到淺冰蓋近似模型；而如果應用在冰架上，假設冰架底面切應力為零，流速的水平分量不依賴于深度，即vx=vx（x，y，t），vy=vy（x，y，t），則得到淺冰架近似模型。淺冰蓋近似和淺冰架近似可以結合成混合模型（hybrid model），應用到整個冰蓋的模擬［18-20］。

黏性系數是動力學方程的關鍵參數，直接影響著冰流的速度大小，并且敏感地依賴于冰溫和冰晶組構等因素。因此，熱-力耦合的動力學方程的模擬結果相對更接近實際情形。full-Stokes方程及其近似方程是建立在冰各向同性的假設之上，由于實際冰蓋冰晶組構各向異性等原因，冰流動力學方程中的重要參數（例如冰的黏性增強因子、冰底摩擦系數）可以通過求解反問題來反演，使得表面冰流速的模擬值與觀測值誤差最小［21-22］。原則上，對于full-Stokes冰蓋模型而言，隨著計算機技術的發展使得完整應力分布的充分計算已成為可能并得到逐步實現。但對于百年以上的積分，某些大陸尺度的冰蓋演化模擬在數值計算上仍不可行。

技術上，冰蓋動力學模型求解的空間域是基于地形數據建立的研究對象的數字高程模型。基于冰厚、海平面高程和冰底高程數據，可以判斷出冰蓋觸地部分和冰架部分。在研究區域的邊界（冰面、冰底、冰崖、側面等）要給出合理的邊界條件，主要包括：（1）冰面溫度、表面物質平衡；（2）冰底觸地部分的熱力平衡方程（包括地熱通量等）、運動狀態等；（3）冰蓋底部的壓融點溫度和物質平衡等；（4）冰崖處的海水壓強和崩解過程等。其中，表面物質平衡和冰架底面物質平衡體現的是大氣和海洋對冰蓋作用的強迫場。冰蓋表面物質平衡的獲取途徑較多，包括實地測量、再分析資料和區域氣候模式的模擬結果等。冰架底面物質平衡的直接監測條件受限（主要是回聲測深儀［23］和鉆孔），缺乏基礎數據。以往和目前的研究多采用參數化方案或者海洋模式的模擬結果進行估算。如何給出較為精確的冰架底面物質平衡是目前研究的主要難點之一。另外，冰蓋底部的運動狀態很難監測，而且對冰下水熱環境變化非常敏感等，具有較大的不確定性。以往研究中的常用方法是在假定的滑動定律形式下通過數據同化反演底部滑動系數［24-25］。冰崖處的崩解過程對南極冰蓋的穩定性非常重要，也是南極冰蓋物質損失很大的不確定性來源。但是崩解過程如何在冰蓋模型中精確地刻畫是迄今尚未較好解決的難題。在以往的冰蓋模擬研究中常見的處理方式有：冰崖的位置保持固定，或者利用冰裂隙深度、冰面徑流、應變率、應力、冰厚等對崩解速率進行參數化［4，26-28］。

冰蓋模型的初始條件是確定起始時間冰蓋的狀態。如何得到符合冰蓋實際情況的初始狀態，即模式的初始化，是一個研究難題。前兩年國際專門針對模式初始化開展了模式比較計劃initMIP［29-30］。常見的方法有長期的預熱模擬（spin-up）和基于數據同化方法的穩態模擬，即融合觀測數據和數值模擬得到同化的模式參數和較為接近實際情形的初始狀態。同化的參數在預估冰蓋未來演化狀態中是否仍然適用，是有待商榷的問題。

根據基本物理定律而構造的用來模擬冰蓋狀態或者預估其未來變化的一套流體力學和熱力學偏微分方程組及其求解方案，稱之為冰流模式。針對冰蓋研發的冰流模式也稱為冰蓋模式。冰蓋模式通過數值求解帶有初始條件和邊界條件的物理模型，模擬冰蓋的速度場、溫度場、應力場等主要物理量的時空分布以及冰蓋形態和觸地線位置的變化，從而有效地理解冰蓋動力演化過程。冰蓋模式計算的大致流程見圖3。

圖3 冰蓋動力學數值模式計算流程圖Fig.3 Flow chart of numerical simulation of ice sheet model

2 中國在南極冰蓋數值模擬領域的主要進展

我國早期開展的南極冰蓋數值模擬研究主要是一些理想化的冰蓋診斷模擬試驗，探索冰蓋運動變化特征。唐學遠等［31］利用基于淺冰近似模型的三維冰蓋模式GLIMMER，計算了方形冰蓋積分區域的流場特征量，研究了冰蓋在大尺度長時間序列條件下對氣候變化的反饋，考察了穩定態下冰蓋演化的周期性行為。唐學遠等［32］采用三維冰蓋模式GLIMMER，對南極冰蓋開展了理想化冰蓋數值試驗，給出了冰蓋處于穩定狀態下的各種演化特征曲線，并將該模式從三維情形簡化為二維情形，在耦合溫度場條件下模擬了二維冰流的演化，探討了簡化后的二維模型在模擬Dome A地區年代-深度-古積累率關系的應用前景。Zhang等［11］比較了三維full-Stokes冰蓋模式求解器FELIX-S［7-9］和冰流模式Elmer/Ice兩個不同Stokes冰蓋模式對簡化的海洋性冰蓋的模擬結果，研究表明：對于診斷模擬，使用相同的幾何形狀和網格時，兩個模型給出了相似的結果（差異＜2%），原因是兩個模型的結構是相似的；對于預測模擬，發現FELIX-S（Elmer/Ice）接地線相對更后退（前進），結果與在診斷試驗結果中觀察到的微小差異一致，這是由于兩個模型中底部邊界條件的不同選擇造成的［11］。

在過去四十年，我國南極科學考察以中山站、昆侖站、泰山站等為支撐平臺，在東南極伊麗莎白公主地和Dome A地區逐步實施了冰芯鉆探和強化觀測［33-34］。沿中山站—Dome A斷面，已逐步建設了星-空-地一體化的業務監測空間網絡。基于這些極地考察區域的觀測優勢，我國在南極內陸冰蓋典型區域開展數值模擬研究，取得了一些國際水平的研究成果。

Dome A在昆侖站附近，位于東南極冰蓋分冰嶺的中心，被認為是獲取古老冰芯記錄的理想地點［35］。Sun等［36］利用冰流模式Elmer/Ice，基于各向同性和簡單的各向異性的冰晶組構的假設條件，對Dome A深冰芯所在區域冰溫和冰齡的空間分布進行診斷模擬，模擬結果表明昆侖站可能提供過去60～70萬年垂直鉆探的高分辨率記錄。Zhao等［37］采用冰流模式Elmer/Ice，將Wang等［38］發現的各向異性冰晶組構隨深度交替現象進行參數化，并利用雷達等時層的年代信息作為約束條件，重新模擬了Dome A地區的冰溫和冰齡空間分布，估算了昆侖站冰芯位置的底部融化速率，指出超過一百萬年的古老冰有望在距離昆侖站5～6 km以內、冰底以上200 m深度處找到。Wang等［39］利用內部等時層和Vostok冰芯年代資料，得到Dome A六條等時層位上的深度-年代關系，并將其輸入到Dansgaard-Johnsen冰流模式，揭示了Dome A地區和昆侖站深冰鉆探點在16萬年以來6個不同歷史時期的古積累率的時空變化。張良甫等［40］采用基于full-Stokes方程的三維熱力-動力耦合冰流模式Elmer/Ice，模擬了甘布爾采夫山脈的速度場和溫度場。泰山站是中國于2014年新建立的科考站。Tang等［41］對泰山站的冰川與氣象條件開展了綜合的分析，并采用一維垂向模型模擬了冰齡和冰底溫度。

在過去十幾年，我國對東南極的埃默里冰架（Amery Ice Shelf）和普里茲灣（Prydz Bay）開展了長期的觀測，對以這兩個區域為代表的南極冰蓋邊緣開展了數值模擬研究。基于高分辨率的海洋-海冰-冰架耦合數值模式，Liu等［42］首次命名了普里茲灣東岸流（PBECC），繼而診斷了PBECC的海洋熱輸送能力及其對Amery冰架底部質量平衡的影響［43］。Cheng等［44］利用自主研發的冰架海洋邊界層模式對Amery冰架底部的邊界層冰晶羽流過程開展了深入的研究和數值模擬，提高了我國自主研發冰架海洋邊界層模式的技術水平。季青原等［45］利用PISM冰蓋模式模擬了Amery冰架的速度場，模擬結果總體上與觀測結果較為接近。Li等［46］利用úa冰流模式模擬和分析了Amery冰架在三次崩解事件前后的變化，發現崩解事件對Amery冰架的穩定性影響很小；也對Amery冰架開展了反演和預測試驗，模擬了冰架解體極端情景下的冰流加速狀況［47］。Sun等［48］基于冰蓋模式BISICLES，采取多種冰架底面融化率參數化方案，對南極包含Totten冰川的奧羅拉（Aurora）流域開展動力學模擬和預估研究。

3 中國在北極的冰蓋數值模擬領域主要進展

我國在北極地區的冰蓋模擬研究主要集中在格陵蘭冰蓋和斯瓦爾巴群島、冰島等地區的冰川。在針對整個格陵蘭冰蓋的數值模擬方面，Yan等［49］利用冰蓋模式SICOPOLIS，探討了不同現代氣候強迫和初始化方案對于現代格陵蘭冰蓋模擬結果的影響，預估了21世紀格陵蘭冰蓋物質平衡的變化以及冰蓋融化對全球海平面上升的貢獻。

表面融水和溢出冰川動力過程的物質輸送是造成格陵蘭冰蓋物質損失的兩大主要因素。格陵蘭冰蓋表面消融產生大量融水，形成冰面水系并輸送到冰裂隙、豎井、冰面湖，下滲到冰蓋內部和底部，通過促進底部滑動和提升冰溫，對冰蓋運動速率與穩定性產生影響。格陵蘭冰蓋溢出冰川末端大多延伸到狹長的峽灣，通過峽灣連通到開闊的海域。海洋環流的變化帶動峽灣水體的變化，對冰川的物質平衡和動力過程起到重要的影響。運用數值模擬預估格陵蘭冰蓋表面徑流和溢出冰川的演化是理解和評估格陵蘭冰蓋未來物質平衡變化的重要方面。

在針對格陵蘭冰蓋表面徑流的模擬研究中，我國在不同的時空尺度上開展了研究。Smith和Yang等［50］利用衛星遙感與無人機影像，提取了一個典型集水區的冰面水系，結合實地測量的冰面河出口位置徑流量，提出了一種測算冰面徑流量的新方法，研究結果表明，目前常用的區域氣候模式過高地估算了冰面實際徑流量，冰面融水輸送的方式和效率對于冰蓋物質平衡以及動力過程具有重要影響。Moore等［51］利用冰蓋表面能量-物質平衡模型，研究平流層注入硫酸鹽氣溶膠的地球工程試驗G4對格陵蘭冰蓋表面徑流的影響，模擬結果顯示在G4情景下冰蓋表面徑流比RCP4.5情景下減少了20%。Yue等［52］利用度日模型，開展了21世紀和22世紀不同氣候情景下格陵蘭冰蓋表面徑流的模擬和預估研究。

在溢出冰川動力過程的數值模擬方面，我國主要對格陵蘭冰蓋最快的溢出冰川——雅各布港冰川進行了數值模擬。Guo和Zhao等［53］采用冰蓋模式BISICLES和同化的參數方案，重現了該冰川在2004—2013年間由于海水增暖造成的冰流速度加速等行為過程，在季節尺度上模擬的冰流速度、冰架末端位置變化、觸地線位置變化與觀測值具有很好的一致性，并且在RCP4.5情景下預估了冰川的動力學演化。

除了格陵蘭冰蓋，我國科研工作者依托黃河站為支撐平臺，對北極斯瓦爾巴群島的冰川等開展了實地觀測和數值模擬有效結合的研究。Ai等［54］利用Elmer/Ice冰流模式對北極黃河站附近的Austre Lovénbreen冰川開展了數值模擬，發現該冰川除了以前實測發現的快速冰流區，還有另一個速度更快的區域。基于該模擬發現，2016年起開展了高精度GPS加密測量，證實了該發現，新發現的快速冰流區比以前測得的快速流速區流速要快8%左右。此外，以歷史觀測記錄和Elmer/Ice冰流模式為參照，艾松濤等［55］在忽略動力學因素的條件下，利用Arc-GIS軟件在冰川輪廓內的每個像素處計算冰面高程變化，進而模擬山地冰川面積、冰厚和體積的長期變化，并探討了該方法的可行性和適用條件。在2015/2016年的北極異常升溫背景下，黃河站Austre Lovénbreen冰川和Pedersenbreen冰川的流速都有較大增加。Wang等［56］通過模擬發現，底部滑動是冰流增速的主要影響因素，直接原因是冰下水量（subglacial water）的增加。隨著冰川融水滲透到冰川內部，還可能出現冬季緩慢釋放至冰川末端形成冰湖的現象［57］。通過冰川生命周期的模擬，斯瓦爾巴群島上類似Austre Lovénbreen冰川的小冰川很可能在一百余年后徹底消失［58］。另外，Yue等［59］應用冰蓋模式PISM，模擬和預估了冰島冰帽在21世紀物質平衡和動力過程的演化。

4 總結與展望

回顧過去十幾年，我國在極地冰蓋數值模擬領域的研究取得了明顯的進步。我國已經自主研發了三個冰蓋模式，并在持續改進。同時，在極地冰蓋及其典型區域上應用冰蓋模式模擬，得到了一些重要的科學認識，取得了具備國際水平的研究成果。經過數十年的極地科考活動，我國在極地冰蓋的監測方面形成了強化監測優勢區域，其中包括南極的Dome A、Amery冰架、伊麗莎白公主地和北極的黃河站等，積累了大量寶貴的數據資源。近年來在格陵蘭冰蓋典型區域也開展了空-天-地一體化監測［60］。這些都為研究冰蓋變化關鍵過程提供了良好的數據平臺。

與此同時，我國在極地冰蓋數值模擬領域也面臨著一些問題，制約著當前和未來的工作，其中包括：研究人員和人才很少，沒有形成實力較強的團隊；擁有自主知識產權的冰蓋模式在研發上處于起步階段，需要持續投入；我國在國際冰蓋模式比較計劃中的參與度不足；研究成果的數量和質量盡管有很大進步，但是涉及入海冰川或冰架崩解和消融的模擬研究較少；國際上已有研究團隊探索了大氣環流模式與冰蓋模式的耦合［61-64］以及海洋模式與冰蓋模式的耦合［65-67］，提高模式模擬的精度，而我國目前鮮有海-冰-氣耦合模式方面的研究。

針對這些問題，需要提高對冰蓋模式模擬研究領域的重視程度，持續加大投入力度，加強國內研究院所的交流合作，有效結合實地觀測、遙感監測和數值模擬，形成研究合力；加強與大氣、海洋等其他圈層領域互動，推進海-冰-氣多圈層耦合和相互作用研究，提升極地冰蓋模擬在氣候變化科學領域的重要性；在我國傳統觀測優勢區域，通過觀測、模擬等綜合研究手段開展聯合攻關，有望在冰蓋關鍵動力過程和機制的科學認識上有所突破，提升冰蓋模式對物理過程的刻畫能力，產出具有中國特色的世界先進水平的研究成果。展望未來，可能需要在以下幾個方面加強研究力量：

（1）充分利用多源數據和區域氣候模式，構建長時間高分辨率時空連續的冰蓋表面物質平衡，準確刻畫區域尺度空間變化特征；開展格陵蘭冰蓋表面融水和冰面水系的觀測和模擬，關聯冰面-冰下融水，探索冰面融水對冰蓋運動的影響。

（2）開展冰蓋典型流域的動力學模擬，研究海洋強迫對入海冰川末端或者冰架的作用機理，改進冰架消融和崩解在冰蓋模式中的實現，探索入海冰川或者冰架與海洋之間的相互作用和耦合模擬，加深對關鍵動力學機制的理解。

（3）研究冰蓋冰下水系分布、冰下水文過程及其對冰蓋不穩定性的影響，開展冰下水文過程模擬，在冰蓋模式中考慮冰下水文過程的影響。