基于GRACE/GRACE-FO和Swarm衛星研究2002—2020年南極和格陵蘭島冰蓋質量時空變化

陳威,鐘敏,馮偉*，王長青，李文浩，梁磊

1 中國科學院精密測量科學與技術創新研究院，大地測量與地球動力學國家重點實驗室,武漢 430077 2 中國科學院大學地球與行星科學學院,北京 100049 3 中山大學測繪科學與技術學院，廣東珠海 519082 4 南京工業大學測繪科學與技術學院，南京 211800

0 引言

在全球氣候變暖的背景下，冰川消融不僅會導致全球海平面上升，還會使得極端天氣事件頻發(Church et al.,2013;Meredith et al.,2019).南極和格陵蘭島作為全球最大的兩個大陸冰蓋，約占地球總冰川面積97%、總冰量99%(Grinsted,2013).若兩極冰蓋完全消融，全球海平面大約上升68 m(Alley et al.,2005;Church et al.,2013).準確估計南極和格陵蘭島冰蓋的質量變化不僅有助于解釋全球海平面變化的成因以及對全球水循環的貢獻，還有助于理解海洋-大氣之間的相互作用和全球氣候變化.

自2002年以來，GRACE(Gravity Recovery and Climate Experiment)重力衛星在監測南極和格陵蘭島冰蓋質量變化中發揮了不可替代的作用.Ramillien等(2006)利用GRACE較早地給出了2002—2005年期間南極和格陵蘭島冰蓋質量變化，指出由于研究時間段較短很難從年際項中分離出趨勢信號.Barletta等(2008)顧及到冰川均衡調整(Glacial Isostatic Adjustment,GIA)的影響后，估計的南極和格陵蘭島質量變化速率分別為-171 Gt·a-1、-101 Gt·a-1.Velicogna(2009)發現2002—2009年期間南極和格陵蘭島的冰蓋消融速率并非均勻變化，它們分別以30 Gt·a-2、26 Gt·a-2的加速度項在加速消融.Nerem和Wahr(2011)通過比較去除GIA的J2項與南極、格陵蘭島冰蓋的質量變化，發現2002—2010年期間地球的扁率變化主要與南極和格陵蘭島冰蓋消融有關.Velicogna等(2014)詳細分析南極和格陵蘭島冰蓋子流域的質量變化后，發現格陵蘭島70%的冰川融化來自于其西南和東北區域，南極81%的冰川消融來自于阿蒙森海域地區和南極半島.Forsberg等(2017)通過分析2002—2015年期間南極和格陵蘭島冰蓋的質量變化，推斷出它們分別等效于全球平均海平面每年上升約0.72 mm和0.26 mm.Groh等(2019)采取多種GRACE數據后處理方法分別估計了南極和格陵蘭島冰蓋的質量變化，認為2003—2013年期間南極冰蓋的消融速率范圍為99～108 Gt·a-1，格陵蘭島的消融速率則為252～274 Gt·a-1.

1 數據和方法

1.1 Swarm衛星數據

(1)

(2)

(3)

1.2 GRACE和GRACE-FO數據

本文使用的GRACE/GRACE-FO level-2數據產品來自于CSR、GFZ和JPL三家公布的GSM 60階球諧系數(http:∥icgem.gfz-potsdam.de/series).其中GRACE數據的時間范圍選用2002年4月至2017年6月，GRACE-FO數據選用2018年6月到2020年3月.GRACE和GRACE-FO均為60階球諧系數，其后處理方法除采用300 km半徑的高斯濾波和nmax=60外，其余后處理方法與Swarm數據的處理步驟一致.本文將CSR、GFZ和JPL三家機構的算術平均值作為最終估計結果.此外，我們還選用了RL06 Version 2版本的CSR Mascon和 JPL Mascon數據來驗證本文的最終結果.

1.3 ARIMA-MC預測方法

為彌補GRACE/GRACE-FO衛星之間近一年的間斷期，本文除使用Swarm時變重力場數據外，還引入了ARIMA-MC(Autoregressive Integrated Moving Average Model-Monte Carlo)方法來預測間斷期兩極冰蓋的質量變化.ARIMA即自回歸移動平均模型，它是由Box等(2016)提出的一種基于統計學的時間序列預測分析方法.ARIMA包括自回歸模型(Autoregressive,AR)，差分模型(Integrated,I)和移動平均模型(Moving Average,MA)三個部分.首先，為保證樣本時間序列具有平穩性需要進行d階差分，差分后的時間序列采用自相關函數來判斷其是否具有平穩性.進行d階差分后，我們采用最小化信息量準則AIC(Akaike Information Criterion)來確定自回歸移動平均模型ARMA的階數p和q，最小化AIC公式為：

AIC=-2logL+2(p+q+1),(4)

這里的L為d階差分后時間序列的似然函數值.確定d、p和q之后，利用Matlab中的ARIMA函數獲取最優的擬合模型.ARIMA模型公式如下：

(5)

其中L為滯后算子，?i為自回歸模型系數，θi為移動平均模型系數，εt為白噪聲序列.確定最優的ARIMA模型后，依據ARIMA模型的概率分布隨機數進行蒙特卡羅模擬(Monte Carlo,MC)，并以此來預測2017年7月—2018年5月期間的南極和格陵蘭島冰蓋質量變化.本文的MC預測方法選用Matlab中的Simulate函數實現.

2 結果

2.1 兩極冰蓋質量變化的時空分布

圖1表示的是2002年4月至2020年3月期間，南極和格陵蘭島冰蓋質量變化趨勢的空間分布.南極半島(Antarctic Peninsula)和西南極阿蒙森海域(Amundsen Sea)的冰川消融信號較大，威爾克斯地區(Wilkes Land,WL)存在較小的冰川消融信號，而莫德女王地區(Queen Maud Land)和瑪麗伯德地區(Marie Byrd Land)則呈現出冰川累積.格陵蘭島冰蓋的主要消融區域位于其西北、西南、南部和東南的臨海區域，其北部和東北沿海地區也存在一定的冰川消融，但是格陵蘭島中央區域存在微弱的冰川累積.

圖1 2002年4月至2020年3月期間(a)南極和(b)格陵蘭島冰蓋質量變化趨勢的空間分布Fig.1 The spatial patterns of ice mass rate in (a)Antarctica and (b)Greenland from April 2002 to March 2020

2002年4月至2020年3月，南極和格陵蘭島冰蓋質量變化(等效海平面高變化)的時間序列如圖2所示.GRACE/GRACE-FO衛星探測到2002—2020年期間南極和格陵蘭島冰蓋質量損失速率分別為-119±23 Gt·a-1和-259±20 Gt·a-1，等效于全球平均海平面每年約上升0.33 mm和0.72 mm.這一結果與CSR RL06 Mascon v2和JPL RL06 Mascon v2結果十分接近.整個南極冰蓋質量變化的時間序列顯示2007年以前南極冰蓋質量處于凈平衡狀態(即該時間段內整個南極區域的冰蓋消融與累積維持收支相抵的狀態)，但2007年之后冰蓋消融的主導地位日益顯著，尤其是2017年和2020年暖季(當年11月份至次年3月份)出現較大的質量虧損.格陵蘭島冰蓋的時間序列顯示其冰蓋消融具有較強的季節性和趨勢項，且在2012年夏季和2019年夏季出現了年際尺度上的加速消融現象.

圖2 (a)南極和(b)格陵蘭島冰蓋質量變化(等效海平面高變化)的時間序列Fig.2 Time series of ice mass changes and equivalent sea level changes in (a)Antarctica and (b)Greenland,respectively

2.2 Swarm衛星和ARIMA-MC方法在間斷期的結果

如圖3所示，我們分別給出了GRACE/GRACE-FO、Swarm和ARIMA-MC在南極和格陵蘭島冰蓋質量變化的時間序列.2015年1月至2019年6月期間：Swarm與GRACE/GRACE-FO+ARIMA-MC(包含了GRACE/GRACE-FO與ARIMA-MC的預測結果)冰蓋質量變化時間序列在南極和格陵蘭區域的相關系數分別為0.58、0.72；Swarm和GRACE/GRACE-FO+ARIMA-MC觀測到南極冰蓋質量的趨勢變化分別為-70 Gt·a-1、-53 Gt·a-1，格陵蘭島冰蓋質量的趨勢變化分別為-92 Gt·a-1、-138 Gt·a-1.2017—2019年期間(主要包含間斷期)：Swarm與GRACE/GRACE-FO+ARIMA-MC結果在南極吻合較好，但在格陵蘭島2018年前后吻合較差.Swarm與GRACE/GRACE-FO結果在南極也存在一定的差異，Swarm結果在2015、2017、2019年的1—3月份相對于GRACE/GRACE-FO的冰蓋質量虧損的幅度更大,同樣在2016—2017年冰蓋質量增加的幅度較GRACE/GRACE-FO結果也更大.圖3b和3d的箱線圖顯示Swarm結果與GRACE/GRACE-FO+ARIMA-MC在2015年1月—2019年6月期間的數據分布情況：其中Swarm數據與GRACE/GRACE-FO+ARIMA-MC結果在南極的差異較大，GRACE/GRACE-FO+ARIMA-MC主要集中在-400～-1300 Gt，而Swarm數據波動的范圍較大；兩者的數據分布在格陵蘭島較為接近，主要集中在-1000～-1800 Gt.通過與GRACE/GRACE-FO+ARIMA-MC結果的對比可知，Swarm時變數據在格陵蘭島的表現更好.

圖3 (a)和(c)分別為南極和格陵蘭島冰蓋質量變化的時間序列，(b)和(d)分別為南極和格陵蘭島在2015年1月—2019年6月期間的不同數據分布Fig.3 Time series of ice mass changes in Antarctica (a)and Greenland (c),respectively.The boxplots (b)and (d)show the data distribution from Swarm and ARIMA-MC from January 2015 to June 2019

圖4a和4b顯示在2015年1月至2019年6月期間，Swarm與GRACE/GRACE-FO+ARIMA-MC兩種方法在南極和格陵蘭島冰蓋質量變化結果的空間相關性.兩者在南極吻合較好的區域主要位于西南極的大部分地區、北部的莫德女王地區和東部的威爾克斯地區.在格陵蘭島吻合較好的區域主要在整個南部地區和西北部地區.如圖4c—4f所示，兩者在南極和格陵蘭島質量變化趨勢的空間分布吻合較好.南極冰蓋的消融區域主要位于西南極的阿蒙森海域和威爾克斯地區，而南極的瑪麗伯德地區、北部的莫德女王地區和東部的威爾克斯地區呈現冰川累積信號.格陵蘭島的消融區域主要表現在其西部、東南部以及東北部.Swarm和GRACE/GRACE-FO+ARIMA-MC的空間趨勢分布也有一定的差異.GRACE/GRACE-FO+ARIMA-MC趨勢變化空間分布結果顯示，南極的瑪麗伯德地區的冰川累積較Swarm的結果更多，威爾克斯地區冰川質量虧損較Swarm結果更少.兩者在格陵蘭島北部以及東南部地區的空間分布也存在一定差異.因此在GRACE和GRACE-FO的觀測空白期，Swarm衛星和ARIMA-MC方法的結果均可有效填補南極和格陵蘭島冰蓋質量變化信息，但兩者之間也存在一定差異.

圖4 在2015年1月—2019年6月期間，Swarm和GRACE/GRACE-FO+ARIMA-MC結果分別在(a)南極和(b)格陵蘭島的相關性，白色區域表示相關系數在統計學上未通過95%的置信區間；(c)ARIMA-MCGRACE/GRACE-FO+ARIMA-MC和(d)Swarm 結果在南極的空間趨勢項；GRACE/GRACE-FO+ARIMA-MC (e)和Swarm (f)結果在格陵蘭島的空間趨勢項.圖中GFO為GRACE-FO的縮寫Fig.4 From January 2015 to June 2019,the correlation coefficients between the results from Swarm and GRACE/GRACE-FO+ARIMA-MC in (a)Antarctica and (b)Greenland,respectively.The white indicates that the correlation coefficients don′t statistically pass the 95% confidence interval.The Antarctic spatial trend is from GRACE/GRACE-FO+ARIMA-MC (c)and Swarm (d),respectively.The spatial trend of Greenland is from GRACE/GRACE-FO+ARIMA-MC (e)and Swarm (f),respectively.GFO is the abbreviation of GRACE-FO

2.3 Swarm和ARIMA-MC結果的不確定性

圖5a和5b顯示的是Swarm結果在南極和格陵蘭島的總誤差分布，該誤差包括了公式(3)中的Swarm觀測誤差、泄露誤差和尺度因子誤差.Swarm在南極和格陵蘭島格網點上的總誤差主要集中在0～500 mm區間，主要位于冰川質量穩定區域；總誤差的最大值約為2000 mm，且主要位于冰川消融和冰川累積區域.在圖5c和5d中，南極和格陵蘭島格網上的尺度因子主要集中在-5～5之間，其最大值約為15且主要集中在冰川消融和冰川累積區域.總而言之，Swarm數據的尺度因子和總誤差主要體現在冰川消融和冰川累積區域.

圖5 Swarm結果在(a)南極和(b)格陵蘭島的總誤差，用于恢復(c)南極和(d)格陵蘭島真實質量變化所采用的尺度因子Fig.5 The total error of Swarm′s results in (a)Antarctica and (b)Greenland,respectively.The scale factors are used to recover the true mass in (c)Antarctica and (d)Greenland,respectively

如圖6a和6b所示，利用ARIMA在南極和格陵蘭島建立質量變化模型對應的最小化AIC的空間分布較為均勻，這表明ARIMA模型和預測結果的穩健性較強.圖6c和6d表示的是利用MC方法預測間斷期南極和格陵蘭島冰蓋質量變化的不確定度.MC方法的不確定度范圍主要集中在0～500 mm區間，其最大值約為2000 mm且主要位于冰川消融和冰川累積區域.總之，ARIMA-MC預測方法具有較強的穩健性，且預測結果的不確定度與Swarm總誤差分布較為類似，均集中在冰川消融和冰川累積的區域.

圖6 ARIMA-MC方法在(a)南極和(b)格陵蘭島預測中的AIC值，ARIMA-MC方法在(c)南極和(d)格陵蘭島的不確定性Fig.6 The AIC of ARIMA-MC method in (a)Antarctic and (b)Greenland,respectively.The uncertainties of ARIMA-MC method in (c)Antarctica and (d)Greenland,respectively

3 討論

3.1 南極威爾克斯地區冰川的加速消融

近些年來，南極冰蓋特別是威爾克斯地區的冰川正在加速消融.圖7a和7b表示的是2002—2009年期間和2010—2020年期間，GRACE/GRACE-FO+ARIMA-MC南極冰蓋質量變化趨勢的空間分布.南極東部的威爾克斯地區和西南極的阿蒙森海域出現了明顯的加速消融，而西南極的瑪麗伯德地區和北部的莫德女王地區冰川出現了加速累積信號.圖7c顯示的是南極威爾克斯地區冰川質量變化的時間序列.威爾克斯冰川在2002—2009年期間的消融速率為3 Gt·a-1，但在2010—2020年期間的速率增加了10倍達到30 Gt·a-1.同樣Swarm結果也顯示，該地區在2015年1月—2019年6月期間的消融速率達到73 Gt·a-1.Shen等(2018)基于質量輸入輸出模型和Lansat影像數據發現威爾克斯地區冰川加速消融.Rignot等(2019)基于較新的冰流和冰厚度數據分析得到威爾克斯地區2009—2017年間的冰川質量損失為51 Gt·a-1.Miles等(2016)認為全球大量冰川的退縮與海冰減少可能會加速溫暖海水對冰川末端的侵襲，這可能會加速威爾克斯地區的冰川消融.Zhan等(2018)發現威爾克斯地區冰川的消融對厄爾尼諾影響下的海表溫度異常變化較為敏感，這可能是導致威爾克斯地區冰川加速消融的原因之一.Shen等(2018)認為威爾克斯地區附近的海洋熱通量會受到南極東部海域冰川動力學的影響，從而引發近年來威爾克斯地區冰川的加速消融.

圖7 (a)和(b)分別為南極在2002年4月—2009年12月和2010年1月—2020年3月期間冰川質量變化趨勢的空間分布，(c)是南極威爾克斯地區冰川質量變化的時間序列Fig.7 The spatial patterns of ice mass rate in Antarctic during the period of (a)April 2002—December 2009 and (b)January 2010—March 2020,respectively.(c)Time series of ice mass changes in Wilkes Land of Antarctica

3.2 夏季NAO對格陵蘭島冰蓋的影響

圖8a給出了格陵蘭島冰蓋質量在2002年4月至2020年3月的夏季期間的年際變化和北大西洋濤動指數(North Atlantic Oscillation,NAO).格陵蘭島質量在夏季的年際變化與北大西洋濤動夏季指數呈正相關，相關系數為0.78(通過95%置信區間檢驗).在2012年和2019年的北大西洋濤動夏季指數負相位期間，格陵蘭島冰蓋較平常夏季出現約300 Gt的質量損失加劇.圖8b顯示了NAO夏季指數與格陵蘭島夏季冰蓋質量年際變化相關性的空間分布，格陵蘭西北部和東部與NAO夏季指數呈現出較強的正相關，而格陵蘭島中央地區與NAO夏季指數呈較強的負相關.van Angelen等(2013)和Bevis等(2019)認為格陵蘭島冰蓋在2012年夏季明顯的加速消融與NAO有關.當NAO夏季指數出現負相位時，格陵蘭島上空會出現高壓和較好的能見度，這會增強太陽輻射的表面吸收、減少降雪和氣溫升高，從而導致冰川在更長時間段內加速消融(Fettweis et al.,2013;Van Angelen et al.,2014;Bevis et al.,2019).因此2019年夏季的明顯加速消融與2012年夏季的消融類似，是由于同時期負相位的NAO影響所致.在經歷2019年負相位的夏季NAO之后，我們推測2021年夏季格陵蘭島冰蓋的年際變化會類似2013年夏季，較正常年份的夏季出現消融減弱.

圖8 (a)NAO夏季指數與格陵蘭島冰蓋質量年際變化；(b)NAO夏季指數與格陵蘭島質量年際變化相關性的空間分布.JJAS表示選用6月、7月、8月和9月作為本文研究的夏季Fig.8 (a)The summertime NAO index (JJAS)and interannual variations of ice mass (JJAS)in Greenland.(b)The correlation coefficients between summertime NAO index and interannual variations of Greenland′s ice mass in spatial domain.JJAS denotes the months of June,July,August,and September in summer

3.3 南極和格陵蘭島冰蓋質量變化趨勢的估計

表1列舉了本文和不同學者關于南極、格陵蘭島冰蓋質量變化趨勢估計的研究結果.其中南極趨勢項誤差考慮了擬合誤差和GIA 20%的不確定性，格陵蘭島的趨勢項誤差參考Ran等(2018).在2003—2012、2003—2013和2003—2014年三個時間段內，本文基于球諧系數估計的趨勢項與CSRM RL06 v2、JPLM RL06 v2的結果吻合較好.通過對比Mu等(2017)結果發現，本文對于2003—2012年南極冰蓋質量變化趨勢的估計約有～10 Gt·a-1的差異，但在格陵蘭島變化趨勢的估計吻合較好，這可能是由于使用不同GIA模型對南極趨勢項的估計產生了一定差異.通過對比Schrama等(2014)、Velicogna等(2014)、Watkins等(2015)和Groh等(2019)結果發現，本文對于2003—2013年南極變化趨勢的估計處于所有估計值的中值，但不同學者由于采用不同的GIA模型使得南極的變化趨勢估計值具有很大的不同.不同研究對格陵蘭島冰蓋質量變化趨勢估計值的差異大約小于10 Gt·a-1，這說明不同方法和不同GIA模型對格陵蘭島冰蓋質量變化趨勢估計的影響較小.本文對2003—2014年南極、格陵蘭島冰蓋質量變化趨勢的估計均大于Harig和Simons(2015,2016)的結果，但與Mascon估計的結果吻合較好.基于以上不同時間段的對比結果發現，Harig和Simons(2015,2016)使用的Slepian方法可能會低估南極、格陵蘭島冰蓋質量變化的趨勢項，這可能是因為Slepian函數在外擴南極、格陵蘭島陸海邊界信號時存在信號損失等問題(高春春等,2019).總之，使用不同GIA模型對南極冰蓋質量變化趨勢的估計會造成較大的差異，而不同的后處理方法和GIA模型對格陵蘭島冰蓋質量變化趨勢估計的影響相對較小.

表1 南極和格陵蘭島冰蓋消融質量變化速率的研究結果(單位：Gt·a-1)Table 1 The ice mass loss rates in Antarctica and Greenland from recent studies (unit:Gt·a-1)

4 結論

本文結合Swarm數據和ARIMA-MC預測方法填補了GRACE/GRACE-FO間斷期南極和格陵蘭島冰蓋質量變化的時間序列，并詳細分析了其時空變化規律.結論如下：

(1)GRACE/GRACE-FO衛星探測到2002年4月至2020年3月期間南極和格陵蘭島冰蓋質量損失速率分別為-119±23 Gt·a-1和-259±20 Gt·a-1，等效于全球平均海平面高每年約上升0.33 mm和0.72 mm.南極地區的冰蓋消融主要在南極半島、西南極阿蒙森海域以及威爾克斯地區，南極莫德女王地區和瑪麗伯德地區則呈現出冰川累積.格陵蘭島冰蓋的主要消融區域位于西北、西南、南部和東南的沿海區域，其中央區域存在微弱的冰川累積.

(2)Swarm衛星數據和ARIMA-MC預測方法均能有效填補GRACE/GRACE-FO間隔期南極和格陵蘭島冰蓋質量變化的時間序列空白.Swarm衛星和ARIMA-MC方法可探測到空窗期南極、格陵蘭島冰蓋的消融信號，它們在主要的冰川消融區域呈現出較強的正相關，但在數據分布和質量變化趨勢的估計中存在一定差異.ARIMA-MC預測結果的不確定度與Swarm總誤差分布類似，它們的最大值分布在冰川消融和冰川累積區域.

(3)南極威爾克斯地區2010—2020年期間的冰川融化速率較2002—2009年期間增加了10倍，這可能與其附近的海洋溫度變化有關.夏季格陵蘭島冰蓋年際變化與NAO夏季指數相關系數為0.78，格陵蘭島2019年夏季明顯的加速消融與2012年夏季消融類似，均與負相位夏季北大西洋濤動有關.不同GIA模型會顯著影響南極冰蓋質量變化趨勢項的估計，但對格陵蘭島冰蓋質量變化趨勢項估計的影響則相對較小.

致謝感謝ICEGM (International Centre for Global Earth Models)網站提供的Swarm和GRACE/GRACE-FO衛星觀測數據(http:∥icgem.gfz-potsdam.de/home).感謝美國俄亥俄州立大學C.K.Shum教授和南方科技大學冉將軍助理教授對文章給出的建議.感謝審稿專家對本文提出的寶貴修改意見.