Return to rapid ice loss in Greenland and record loss in 2019 detected by the GRACE-FO satellites

Between 2003-2016, the Greenland ice sheet (GrIS) was one of the largest contributors to sea level rise, as it lost about 255 Gt of ice per year. This mass loss slowed in 2017 and 2018 to about 100 Gt yr−1. Here we examine further changes in rate of GrIS mass loss, by analyzing data from the GRACE-FO (Gravity Recovery and Climate Experiment – Follow On) satellite mission, launched in May 2018. Using simulations with regional climate models we show that the mass losses observed in 2017 and 2018 by the GRACE and GRACE-FO missions are lower than in any other two year period between 2003 and 2019, the combined period of the two missions. We find that this reduced ice loss results from two anomalous cold summers in western Greenland, compounded by snow-rich autumn and winter conditions in the east. For 2019, GRACE-FO reveals a return to high melt rates leading to a mass loss of 223 ± 12 Gt month−1 during the month of July alone, and a record annual mass loss of 532 ± 58 Gt yr−1. Mass loss from the Greenland ice sheet returned to record levels in 2019, following unusually small loss in 2017-18, according to an analysis of satellite data from GRACE and its follow-on mission GRACE-FO.

S tarting in the mid-1990s, the Greenland ice sheet (GrIS) transitioned from a modest contributor to global sealevel, adding 0.21 ± 0.1 mm year −1 in the 20th century   1 to causing 0.76 ± 0.1 mm year −1 of the total of 3.5 ± 0.2 mm year −1 global mean sea-level rise between 2005 and 2017 2 , almost equal to the contributions of all glaciers worldwide (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017). One-third of the ice loss is attributed to an 18 ± 1% increase of solid ice discharge into the ocean (2010-2018 relative to   3 and the other two-thirds is attributed to an increase in surface melting that reduced the surface-mass balance (SMB) by 48 ± 9% 3,4 . Melt production increases due to a variety of factors, including rising near-surface temperatures 5,6 , reductions in surface albedo 7 , migration of the snow line to higher altitudes 8 , an increase in the total melt area 9 and cloud-radiative effects 10 . Global climate models project Arctic annual temperatures to rise about two and a half times faster compared to the tropics, with a factor of one and a half estimated for the summer months only 11 . GRACE and GRACE Follow-On (GRACE-FO) enable us to quantify the ice sheet response to meteorological forcing on sub-annual time scales, helping to improve our understanding of feedbacks between surface processes and ice dynamics that are highly relevant for projecting mass changes of the GrIS.
Numerous studies have quantified the contemporary mass loss of the GrIS using different satellite-based Earth observations 3,12-14 . The consensus of these studies is that the GrIS has lost an average of −148 ± 13 Gt year −1 of ice from 1992 to 2018, with total losses (meltwater runoff and ice discharge) recently (2003-2016) exceeding total gains (snowfall) by more than one third 13 . The GRACE/ GRACE-FO time series provides evidence for an acceleration of mass loss from 2003 until 2012, caused by an increase in both the flowspeed of glaciers draining the ice sheet and melting at their surface 15 . Enhanced melt production correlates with more frequent anticyclonic circulation anomalies in mid-troposphere that are in-turn related to the presence of stable high-pressure systems over Greenland ('Greenland blocking') 16 . Although such changes are not projected by global models, anticyclonic conditions prevailed in about half of all summers since the end of the 1990s 17 . For example, the sizeable mass losses in 2010 and 2012, −462 ± 60 Gt year −1 and −464 ± 62 Gt year −1 respectively as measured by GRACE, were conditioned by an exceptional persistent high-pressure system over Greenland, forcing warm air from mid-latitudes along the west coast of the ice sheet 18 , along with cloud radiative effects 19,20 . Strong melt-albedo feedbacks amplified the response of the surface-mass balance to the atmospheric drivers, favouring anomalously dry and sunny summers over the ice sheet 21,22 .
In situ and remote sensing observations have previously documented the anomalous GrIS surface conditions in 2017 and 2018 13 . However, with end of science operations of the GRACE mission in June 2017, mass balance estimates from observations remained highly uncertain for that period. The GRACE-FO mission was launched on May 22, 2018. GRACE-FO mirrors the GRACE mission concept of two satellites measuring the variations in their along-track distance caused by differences in Earth's gravitational attraction using a microwave ranging system, recording redistributions of water and ice on the Earth surface 23,24 . As a technical demonstration, GRACE-FO is additionally equipped with a laser interferomer ranging system, currently being tested to facilitate more accurate tracking of the inter-satellite distance and, potentially, enable better-resolved gravity fields 25 . Both missions record nearly independent mass changes for distances greater than~300 km. For the GrIS 26 , GRACE/GRACE-FO mass changes are typically separated into six to eight drainage basins 13,27,28 . Because the GRACE-FO mission is able to track changes in gravity since the end of the GRACE mission, the most recent GRACE-FO measurements represent the cumulative mass change occurring between the last GRACE and first GRACE-FO observation. This allows quantification of the combined mass balance for 2017 and 2018.
Here, we show that early data from GRACE-FO compared to GRACE reveal a 58% slowdown in GrIS mass loss in 2017-2018 4 (−98 ± 29 Gt year −1 ; ±2σ), similarly indicated by satellite altimetry records 5 . Using regional climate modelling output 6 , we show that reduced losses in 2017 and 2018 result from two anomalous cold summers in western Greenland compounded by snow-rich autumn and winter conditions in the east. For 2019, GRACE-FO reveals a return to high rates of loss with a mass change of −223 ± 12 Gt month −1 during the month of July alone. We further explain how changes in synoptic conditions influence rates of loss; intensified summer cyclonic circulation over Greenland in 2017 and 2018 favoured southward movement of colder air masses along the melt-sensitive western flank of the GrIS. In summer 2019, conditions were largely reversed with a dominance of anti-cyclonic conditions over the ice sheet, similar to 2012, advecting warm mid-latitude air masses to northwestern Greenland. Combined with low snowfall, 2019 sets a new record GrIS annual mass balance of −532 ± 58 Gt year −1 (−464 ± 62 Gt year −1 in 2012), unprecedented in 1948-2019 and probably within the 20th century 7 . Despite the slowdown recorded in 2017 and 2018 and large year-to-year variability, the rapid mass loss in 2019 recorded by GRACE-FO indicates the GrIS has remained on a trajectory of increasing mass loss since the late 1990s in response to Arctic warming.

Results
Satellite gravimetry and regional climate models. Figure 1 presents the time series of mass change of the GrIS derived from 180 monthly GRACE and GRACE-FO gravity solutions (2002-2019) ("Methods" section). Using a piecewise linear model, we estimate biennial mass balances starting in January 1st of each year, centred on the strongest cumulative loss in July ("Methods" section). The estimate for 2017-2018 contains the five last months of GRACE data and the first five from GRACE-FO in 2018, even though gravity fields at the end of 2017 have higher uncertainties due to partial instrument failure (loss of accelerometer on GRACE-B satellite). The biennial mass balance averaged for the entire observation period is −235 ± 29 Gt year −1 .
The year 2012 marks a pronounced mass balance anomaly; mass loss increases to a peak biennial rate of −437 ± 26 Gt year −1 in 2011 and 2012 and slow to an average rate of −143 ± 17 Gt year −1 during the following 6 years (2013-2018). Rates of loss were especially low during the 2017-2018 period when GRACE and GRACE-FO detect an average rate of change more than four times less negative (−98 ± 26 Gt year −1 ) than the peak rates observed in 2011-2012. Nevertheless, all 2003-2018 biennial mass balances showed a statistically significant mass loss at the 5% confidence level.
To determine the cause of the 2017-2018 slowdown in rates of loss, we compare the GRACE/GRACE-FO time series with monthly mass budgets of the ice sheet. The SMB from the regional climate models MARv3.10 29,32 (20 km resolution; forced by NCEP-NCARv1 reanalysis 33 ) and RACMO2.3p2 30 (statistically downscaled from 5.5 km to 1 km; forced by ECWMF reanalysis 34 ) ("Methods" section), is contrasted against solid ice discharge (D) into the ocean, as estimated at sub-annual temporal resolution from feature tracking of optical and radar imagery 31 . The combination of SMB-D equals the net mass balance as observed by GRACE/GRACE-FO 4 . To account for possible biases in SMB and D 35,36 or the glacial-isostatic adjustment correction of the gravity fields ("Methods"), we apply a linear trend correction to the SMB-D data to reconcile with the To identify the processes underlying these anomalies, we separate the SMB-D* into its components of net accumulation (snowfall minus evaporation and sublimation), meltwater runoff (surface melt minus refreezing and retention), and solid ice discharge (D*). Figure 2 shows that the reduced mass loss in 2017-2018 is caused by both enhanced accumulation between 56 and 88 Gt year −1 (red bars in Fig. 2) and reduced meltwater runoff between +65 and 78 Gt year −1 (green bars in Fig. 2). The magnitudes of these anomalies are not exceptional over the study period-a comparable anomaly in snowfall and meltwater runoff occurred in 2003-2004 and 2013-2014 respectively. What is exceptional over the period of study is the occurrence of strongly positive snowfall and meltwater runoff anomalies in the same years. From 2003 to 2018, net runoff controls the largest variability in rates of biennial mass change (±61 to 69 Gt year −1 of ±83 to 98 Gt for SMB-D*), indicating a high sensitivity of the ice sheet to summer atmospheric forcing and changes in the surface radiation budget. The variability of biennial snowfall is only slightly lower at ±39 to 52 Gt year −1 .
The GRACE/GRACE-FO and SMB-D* mass anomalies in Fig. 2 show a striking regional difference in changes that are driven by snowfall in the east and melt in the west of Greenland. In fact, nearly all (72-88%) of the snowfall anomaly (+65 and 78 Gt year −1 ) in 2017-2018 occurred on the east and southeast sectors of the ice sheet, whereas 71-77% of the runoff anomaly occurred in the west and southwest (+50-55 Gt year −1 ). This spatial pattern of snow and melt anomalies are congruent with the climatic divide of the GrIS, with SMB variability dominated by snowfall in the east and melt in the west 29,37 . Ice discharge shows a slowdown of 7 Gt year −1 within the biennial period 2017-2018, largely due to the slowdown of Greenland's largest glacier 38 , but remains close to the 2003-2018 average.

Discussion
Warming of the Arctic over the past few decades has exceeded global average rates of warming, decreasing the near-surface air temperature gradient to the mid-latitudes by about 2°C 39,40 . Similar amplification, with factors ranging from one and a half in summer to three and a half in winter, is projected for the end of the century with contributions from various processes 11,41 . Today, circulation in the Arctic exhibits more stable highpressure systems and associated atmospheric blocking over Greenland in summer, as reflected by the Greenland Blocking Index 16 , enhancing anticyclonic circulation and advection of warm, continental mid-latitude air masses along the west coast of Greenland 42 . Figure 4 Fig. 1), circulation anomalies showed largely reverse conditions compared to 2012 and 2019; the summer months exhibited a strong low-pressure anomaly over Greenland, forcing cold Arctic air more southward along its west coast (Fig. 4). The modulation of the mass balance of the GrIS is a substantial response of accumulation and melt processes to changes in atmospheric forcing and a strong Arctic warming trend, further accentuated by feedbacks within the ice sheet/climate system, such as the ice-albedo and, on longer time scales, the ice-elevation feedback 45 . These changes are superimposed on substantial trends in solid ice discharge driven by climate-related (ocean and atmosphere) changes in frontal melting and terminus retreat 31,46,47 . Our study highlights that prolonged observations from GRACE/GRACE-FO and other satellite missions, together with modelled SMB are crucial to understand and quantify the processes controlling the ice sheet's response to changing meteorological conditions in the Arctic. Such insights will ultimately lead to improved model representation of the Arctic climate system in global projections of sea-level change.  Fig. 2).   52 , an improvement to the estimation method originally proposed by Swenson et al. 53 . These data are available as GRACE Technical Note 13 from https://podaac-tools.jpl. nasa.gov/drive/files/GeodeticsGravity/gracefo/docs/ for each SDS data set. (2) Replacement of highly uncertain C 20 coefficients from GRACE/GRACE-FO by more accurate estimates from Satellite Laser Ranging (SLR) 54 , accessible as GRACE Technical Note 14 from https://podaac-tools.jpl.nasa.gov/drive/files/ GeodeticsGravity/gracefo/docs/TN-14_C30_C20_GSFC_SLR.txt. Note that although the replacement of C 30 coefficients is recommended by the GRACE/GRACE-FO SDS centers for solutions starting August 2016, this issue is still under discussion from the user point of view since it introduces a discontinuity in the time series between months with nominal and anomalous GRACE/GRACE-FO accelerometer performance 55,56 . In this study, we do not adopt this procedure, as it moderately increases the difference with the modelled mass balance for the year 2019, while slightly decreasing the difference for 2016. For completeness, we state that replacing the GRACE/GRACE-FO C 30  The GRACE/GRACE-FO gravity fields are corrected for long-term mass trends caused by the Earth's viscoelastic relaxation to past changes in ice loading, the glacial-isostatic adjustment (GIA), using the GGG1.D model 57 . GGG1.D is tuned to fit measured GIA-induced uplift rates and represents a rate of gravity field change equivalent to +16 Gt year −1 of ice mass change over Greenland (Supplementary Fig. 3). The spread of ten GIA corrections presented in 13 ranges from −28 to +27 Gt year −1 with a mean of −5 Gt year −1 . As conservative uncertainty estimate, the standard deviation of ±17 Gt year −1 is included in the presented GRACE/GRACE-FO annual and long-term trends.

2017-2018
GRACE/GRACE-FO solution combination. To reduce the noise level of the solutions and mitigate the impact of possible outliers 58 , we estimate 180 combined GRACE/GRACE-FO monthly solutions (AV RL06). This combination is achieved by coefficient-wise weighting of the Stokes potential coefficients C jm from SDS centers (N sol = 3) for each time t, according to C AV;t jm ¼ P N sol n¼1 ðw n;t jm C n;t jm Þ= P N sol n¼1 w n;t jm , where w n;t jm represent weights corresponding to the inverse of the squared variance of the calibrated uncertainty of each coefficient,σ n;t À2 jm . Formal uncertainties provided with the GRACE/GRACE-FO coefficients result from different estimation procedures that prevent their direct use as quantitative weights. Therefore, we calibrate the uncertainties based on the noise level of each solution as follows; we estimate the temporal residual, .. after subtracting bias, trend, annual, semiannual and temporal variations longer than four months (using a moving average filter) from the GRACE/GRACE-FO coefficients' time series. We then determine the degree power in the noise-dominated spectral range (j min = 40 to j max = 60) acccording to P jmax j¼jmin P j m¼Àj ðC Res: jm Þ 2 , which is representative of the noise level in each solution. The formal uncertainties provided with the SDS centers are then calibrated by a single scaling factor to match the degree power of the residual estimated for GRACE/GRACE-FO coefficients, yieldingσ n jm . This approach adopts the error structure from the formal uncertainties with the error magnitude estimated from the residual.
We carry out the combination on detrended data, as the differences in the trends appear to be systematic and arising from different processing choices of the SDS centers (Supplementary Fig. 4). Therefore, the monthly weights derived from coefficients beyond degree and order 40 are not representative for the relative uncertainty of the trends. Applying these weights to the different temporal linear signals would introduce artificial monthly temporal variability of about 2 Gt. Therefore, we remove the trends before combination and restore them, assuming equal weights for each SDS center, in the combined solution, C AV jm . The GrIS mass change time series for the SDS solutions GFZ RL06, CSR RL06 and JPL RL06, as well as their combination AV RL06 are shown in Supplementary  Fig. 4.

Gravimetric inversion.
To retrieve the mass anomalies at a regional scale from the GRACE/GRACE-FO coefficients, we apply an inversion scheme based on forward modelling, using low-frequency damping and a non-uniform a priori mass distribution of known magnitude, and dividing the ice sheet into seven drainage basins. A summary of the inversion procedure is provided in the Technical Note of the GFZ RL06 Level 3 products ftp://isdcftp.gfz-potsdam.de/grace/GravIS/GFZ/ Level-3/ICE/GravIS_RL06_ICE_Technical_Note.pdf. Due to its proximity and possible signal overlaps, we additionally include Ellesmere Island in the estimation procedure. The gravity field signal arising from each basin and Ellesmere Island is calculated in the spectral domain, filtered identical to the GRACE/GRACE-FO data, and then adjusted in scale for each basin to least-squares minimize the difference between the modelled and observed spatial gravity field signal over the GrIS and Ellesmere Island. Using alternate GrIS mass change times series available at http://gravis.gfz-potsdam.de/home, which is based on the inversion schemes of 59 , as well as using AV RL06 with the inversion approach of 28 , we estimate that different assumptions underlying the gravimetric inversion introduce an uncertainty of ±25 Gt for the monthly estimate. These values are consistent with the ±26 Gt uncertainty presented in Table A Uncertainty of GrIS mass changes. To estimate the uncertainty of the GrIS monthly data arising from uncertainties in the GRACE/GRACE-FO coefficients, we first remove deterministic components from the time series of GrIS mass change, M G (t i ), consisting of bias, trend, annual, and semiannual and temporal variations longer than four months (using a moving average filter). The root-meansquared residual of the remaining mass variability M G Res: t i ð Þ is considered noise. This estimate is conservative, as episodic mass changes and fluctuations below four months are included in the noise estimate. The same detrending procedure is applied to the GRACE/GRACE-FO coefficients' time series. Then, the noise level of each monthly solution is estimated by the cumulative degree-power of the residual GRACE/GRACE-FO coefficients C Res jm t i ð Þ in the spectral range j = 40-60, P 60 j¼40 P j m¼Àj C Res jm t i ð Þ, reflecting the noise level of the global gravity field solution. To get monthly uncertainties of mass change, we divide the global noise level of each monthly gravity field solution by the mean noise level for all gravity field months, and multiply this by the mean residual mass variability for the GrIS. Note that this noise level depends on latitude as the ground-track density of GRACE/ GRACE-FO measurements increases towards the poles (Supplementary Fig. 2). The uncertainty estimates for each SDS data set (GFZ RL06, CSR RL06 and JPL RL06) and the combination (AV RL06) are shown in Supplementary Fig. 5.
SMB-D. The SMB and its components over the GrIS are obtained from daily output of the Modèle Atmosphérique Régional (MAR) Version 3.10, run at 20 km resolution and 6-hourly forced by the NCEP-NCARv1 climate reanalyses 29 . The forcing with NCEP-NCARv1 enables near-real-time simulations of MARv3.10 used in this paper. Daily values are aggregated to monthly values by summation and integrated within the drainage basins of the contiguous ice sheet. The temporal coverage is January 1948 to December 2019. As an alternative SMB model, we employ RACMO2.3p2 at a resolution 1 km downscaled from 5.5 km 30 . The model is forced by ERA-40 (1958ERA-40 ( -1978 34 and ERA-Interim (1979-2019) 60 and covers the time period January 1958 to August 2019.
The main SMB components shown are net accumulation (snowfall minus evaporation and sublimation) and meltwater runoff (meltwater production minus meltwater refreezing and deposition and retention). Uncertainties in the SMB 36 are not presented in this paper, however, amount to roughly ±22 Gt for the total de-trended SMB (Supplementary Fig. 6). Note however, that uncertainties in the individual components may compensate each other and that both models may show biases with regard to the unknown true mean or long-term trend values of SMB. If possible, values for both MARv3.10 and RACMO2.3p2 are given in the text and figures, to provide an estimate of the effect of model setup, lateral forcing and grid resolution. We use a temporally extended version of monthly-resolved ice discharge measurements from King et al. 31 , which includes glacier-scale discharge time series at all large (width > 1 km, N > 200) outlet glaciers. Changes in ice discharge are measured through flux gates drawn orthogonal to the mean direction of ice flow and are positioned within 5 km of the most retreated calving front positions, where valid data exist. Time-varying ice discharge records are derived by integrating changes in ice thickness, from differencing ASTER, SPOT-5 and ArcticDEM elevation data with BedMachine v3 61 bed topography, and velocity, from Landsat 7 and 8, RADARSAT 1 and 2, and TerraSAR/TanDEM-X. Velocity data gaps are filled at the individual glacier scale by applying a Kalman filter that uses the mean seasonal variability observed at each glacier. The 2-sigma uncertainty of D ranges between 10 and 20 Gt year −1 .
The estimated correction to reconcile GRACE/GRACE-FO and SMB-D trends for 2003-2019 amounts to −7 Gt year −1 (RACMO2.3p2 forced by ECMWF reanalysis) and 52 Gt year −1 (MARv3.10 forced by NCEP re-analysis). For completeness we state that the correction is 27 Gt year −1 and 26 Gt year −1 for MARv3.10 and MARv3.9 forced by ECMWF re-analysis, and 41 Gt year −1 for MARv3.9 forced by NCEP re-analysis. Note that the correction for the previous version RACMO2.3p1 amounts to 59 Gt year [−138 .
Time series analysis. To enable a consistent comparison, we process the MARv3.10 output of mass changes similar to the GRACE/GRACE-FO time series. First, we calculate the cumulative sum of the MARv3.10 mass changes, representing storage variations of the ice sheet due to surface-mass balance and its components, M MAR (t), and subsample the time series by linear interpolation at the mid-point, t i , of the GRACE measurement epochs (typically 30 days, but varying due to data quality), M MAR (t i ). To obtain mass fluxes, dM t i ð Þ dt , presented in Fig. 3, we calculate the temporal derivative of the GRACE/GRACE-FO and interpolated MARv3.10 time series, using an unweighted central-difference scheme around t i , involving the proceeding (t i−1 ) and succeeding (t i+1 ) time step. Note that this procedure acts as a moderate smoothing filter, reducing the peak values compared to the initial monthly mass fluxes from MARv3.10. Biennial mass balances are determined by unweighted fitting of piecewise linear functions to the mass time series with breakpoints in January 1st of every second year, without prior removal of annualor semi-annual components. Time intervals are centred on summer peak losses to encompass the complete melt seasons, and to allow us to bridge the gap between GRACE and GRACE-FO in 2017-2018. The climatological mean of the seasonal mass fluxes from GRACE (Fig. 3) is obtained by linearly interpolating the GRACE time series dM t i ð Þ dt to regular monthly time intervals, excluding data gaps, followed by unweighted averaging.
RACMO2.3p2 output of the surface-mass balance and its components are freely available from the authors upon request and without conditions. To submit a request, please contact Brice Noël: b.p.y.noel@uu.nl.
Dynamic ice discharge data are available upon request from Michalea D. King: michaleaking@gmail.com.