Self‐Consistent Ice Mass Balance and Regional Sea Level From Time‐Variable Gravity

Measurements of time‐variable gravity from the Gravity Recovery and Climate Experiment (GRACE) and the GRACE Follow‐on (GRACE‐FO) missions are an invaluable tool for monitoring changes in the mass of the Earth's glaciated regions. We improve upon estimates of glacier and ice sheet mass balance from time‐variable gravity by including instantaneous spatiotemporal variations in sea level. Here, a least squares mascon technique is combined with solutions to the sea level equation to iteratively correct the GRACE/GRACE‐FO data for the induced sea level response on a monthly basis. We find that variations in regional sea level affect ice sheet mass balance estimates in Greenland by approximately 4% and in Antarctic by approximately 5%. Since 2002, the Greenland ice sheet has been losing mass at an average rate of 263 ± 23 Gt/yr, and the Antarctic ice sheet has been losing mass at average rates between 90 ± 52 and 122 ± 53 Gt/yr depending on the rate of glacial isostatic adjustment. The mass losses from both ice sheets represent an increase of 15.6 ± 2.0 mm to global mean sea levels since 2002.


Introduction
Time-variable gravity measurements from the National Aeronautics and Space Administration (NASA)/DLR Gravity Recovery and Climate Experiment (GRACE) and the NASA/GFZ GRACE Follow-on (GRACE-FO) missions are an essential method for estimating the mass change of the Earth's glaciated regions. The gravity measurements from GRACE and GRACE-FO are global, near monthly and are directly related to changes in mass. This measurement technique has distinct advantages over traditional measurements of the ice mass balance, which are difficult to scale over broad regions due to sparse sampling and the complex behavior of glacial systems (Gardner et al., 2013). However, because several mass transport processes can occur concurrently for a given region, the total time-dependent geopotential from GRACE/GRACE-FO can relate to multiple time-varying components (Forsberg et al., 2017;Wahr et al., 1998). Isolating the mass change of a glacial system from time-variable gravity data requires the removal of all other temporally varying Earth system components (Velicogna, 2006;Wahr et al., 1998). Uncertainties in the components removed from the GRACE/GRACE-FO data will directly impact the precision of the final mass balance estimate.
One specific component of regional mass change that may affect estimates of glacier and ice sheet mass change is due to regional variations in sea level. A variation in mass of the Earth's glaciated regions leads to both an elastic response of the solid Earth and a change in the Earth's gravitational field (Tamisiea & Mitrovica, 2011). The combination of the solid Earth deformation and gravity change induces a sea level variation that deviates strongly from a uniform barystatic change (Clark & Lingle, 1977;Gregory et al., 2013; Figure 1). This nonuniform redistribution of sea level is often referred to as the effects of self-attraction and loading (SAL) for the gravitational and deformation effects respectively or as a sea level fingerprint (SLF) of the land mass change (Tamisiea & Mitrovica, 2011). The SLF signal has been unambiguously detected apart from changes in barystatic sea level using time series from ocean bottom pressure recorders (Hsu & Velicogna, 2017). The detection of the SLF signal in oceanographic data is direct evidence that the redistribution of sea level is occurring on rapid timescales with current land-sea fluxes (Hsu & Velicogna, 2017). Modern-day changes in ice sheet mass are some of the largest contributors to present-day sea level rise (Church et al., 2011) and will likely be one of the key drivers of future sea level rise .  (Tamisiea et al., 2010). A reduction in land ice results in both a gravitational change and an elastic uplift of the solid Earth, which results in a reduction of the height of the water column near the coast. The initial state from (a) is indicated in (b) by the dashed red lines.
The sea level equation derived in Farrell and Clark (1976) provides a method for calculating the spatial pattern of sea level variation induced by a change in the Earth's mass distribution Clark & Lingle, 1977;Kendall et al., 2005;Mitrovica & Milne, 2003). Nonuniform variations in local sea level in turn impact estimates of ice mass change from time-variable gravity due to the limited spatial resolution of the GRACE/GRACE-FO data (Adhikari & Ivins, 2016;Sterenborg et al., 2013;Velicogna & Wahr, 2013). In Velicogna and Wahr (2013), SAL effects were calculated using a synthetic mass distribution that redistributed the estimated mass change over a defined spatial pattern. The spatial pattern was defined from auxiliary data sets of ice sheet mass change as the pattern is not estimated using the optimized averaging kernel technique . In Velicogna et al. (2014) and Sutterley et al. (2014), SAL effects were corrected over the long-term by calculating the SLF of mass changes calculated using the mass budget method . The sea level estimates in Velicogna et al. (2014) and Sutterley et al. (2014) were scaled to take into account differences between the GRACE and mass budget method solution over the long term. These previous corrections did not account for accelerations in sea level rise that are evident with longer observations (Nerem et al., 2018). Here, we create a self-consistent monthly record of ice mass change and regional sea level rise using an iterated least squares mascon technique combined with solutions to the sea level equation to correct for leakage from land-sea fluxes. In this study, we (i) calculate glacier and ice sheet mass balance estimates that are fully compatible with the induced sea level response, (ii) investigate the sensitivity of mass balance estimates to different sea level redistributions, and (iii) contextualize the ice sheet mass balance results with other source of uncertainty. In the following sections, we discuss how the method is implemented, the mass balance and sea level estimates resulting from the method, and the overall implications of the results for gravimetric ice sheet mass balance studies.

Materials and Methods
We use 177 monthly GRACE/GRACE-FO Release-6 (RL06) gravity solutions provided by the University of Texas Center for Space Research for the period April 2002 to August 2019 (Bettadpur, 2012;Tapley, 2004) provided by the NASA Physical Oceanography Distributed Active Archive Center. Each CSR gravity field solution consists of fully normalized spherical harmonic coefficients (C lm , S lm ) up to degree, l, and order, m, 60. We substitute the C 20 coefficients with monthly estimates from Satellite Laser Ranging provided by NASA Goddard Space Flight Center due to anomalously large variability in the original coefficients (Loomis et al., 2019;Velicogna & Wahr, 2013). For periods when the satellite pairs had a single fully operating accelerometer, we replace the C 30 coefficients derived from GRACE/GRACE-FO data with Satellite Laser Ranging-derived estimates provided by Goddard Space Flight Center (Loomis et al., 2019; supporting information section S5). Variations of the Earth's geocenter are inferred using GRACE/GRACE-FO data of degree 2 and higher using an iterated self-consistent solution with self-attraction and loading effects (Sutterley & Velicogna, 2019;Swenson et al., 2008). We correct the GRACE/GRACE-FO mass changes for glacial isostatic adjustment (GIA), the Earth's viscoelastic response to the glacial unloading since the last glacial maximum using coefficients from Ivins et al. (2013), Whitehouse et al. (2012), and Simpson et al. (2009). Here, we calculate the mean GIA rate for each deglaciation model considering a range of rheological parameters and calculate the GIA rate uncertainty as the root mean square (RMS) off of the mean rate (supporting information Table S1). We remove the variability from terrestrial water storage using outputs from the Global Land Data Assimilation Systems (GLDAS) NOAH land surface model (Version 2.1) (Rodell et al., 2004). The monthly GLDAS terrestrial water storage estimates are calculated by combining the snow water equivalent, canopy water storage, and soil moisture variables for nonglaciated regions (Syed et al., 2008). The total land-sea flux from terrestrial land water storage variability is distributed over the ocean taking into account self-attraction and loading effects, and converted to spherical harmonics to be comparable to the GRACE/GRACE-FO data.
To isolate each ice mass time series, we use the least squares mass concentration (mascon) approach described in Velicogna et al. (2014). We updated the mascon configuration to include each ice sheet and the glaciated regions of Alaska, Canada, Iceland, Svalbard, Franz Josef Land, Svernaya Zemlya, Novaya Zemlya, and Patagonia (supporting information Figure S1). Each mascon is an equal-area 3 • diameter spherical cap with a mass equal to a uniformly distributed centimeter of water (Farrell, 1972;Velicogna et al., 2014). The mascons are constructed in the spherical harmonic domain and are converted to be in terms of mass Tiwari et al., 2009). Each mascon is freely evolving and has its total mass evenly redistributed over the ocean. We smooth the mass coefficients using a 250 km radius Gaussian averaging function to reduce the impact of random spherical harmonic errors (Jekeli, 1981;Wahr et al., 1998). The mass coefficients are simultaneously fit to the geophysically corrected GRACE/GRACE-FO coefficients to obtain estimates of the monthly mass variability for each mascon.
GRACE/GRACE-FO spherical harmonic data is spectrally limited, which limits its effective spatial resolution. Changes in sea level can affect estimates of ice mass change, particularly when the effects of SAL are incorporated (Sterenborg et al., 2013;Velicogna & Wahr, 2013). The MPIOM model used to correct GRACE/GRACE-FO data for ocean circulation does not include the effects of land-sea exchange and the corresponding sea level response . Here, we use our monthly mascon estimates as an input distribution for solving the sea level equation using a pseudo-spectral approach (Farrell & Clark, 1976;Mitrovica & Milne, 2003;Mound & Mitrovica, 1998; supporting information Section S3 ). The resultant sea level fingerprint is used to correct the GRACE/GRACE-FO data for the next iteration of the mascon fit ( Figure 2). This process is repeated until the solutions converge, which usually takes three to four iterations with most of the correction occurring in the first iteration. After the solutions are converged, we calculate the final mass anomaly time series, M(t), for glaciated regions through the summation of each region's local mascons. We follow the same procedure using only the mass change from both ice sheets and then each ice sheet individually to help determine the indirect geophysical leakage components from glaciers, ice caps, and the other ice sheet. Here, we define indirect geophysical leakage as leakage of mass from a region through a third process, such as the redistribution of sea level by glacial mass change, and direct geophysical leakage as the leakage of mass from an adjacent region or by a concurrent process in the region of interest.
Uncertainty in the GRACE/GRACE-FO regional ice mass estimates is a combination of GRACE/GRACE-FO measurement error, geophysical signal leakage, statistical uncertainty, and GIA uncertainty. The contributions of each of these error sources were evaluated as described below . For every mascon, we calculate the corresponding sensitivity kernel that describes how mass at any given point within the region contributes to total time series (Tiwari et al., 2009). The effects of monthly errors in the GRACE/GRACE-FO measurements are calculated by convolving each mascon's sensitivity kernel with estimates of the GRACE/GRACE-FO spherical harmonic uncertainty calculated as described in Wahr et al. (2006). The combined effects of spherical harmonic truncation, Gaussian smoothing, and least squares misfit cause the leakage of signal both in and out of each mascon. For ice sheets, we evaluate the leakage effect for each mascon as described in Velicogna et al. (2014) by convolving its sensitivity kernel with a synthetic signal constructed using a simulated, realistic ice mass change rate over all mascons except for the mascon of the iteration, which is left empty. We calculate the least squares fit uncertainty by convolving the mascon sensitivity kernels with the residual time series obtained by removing our iterated mascon solution from the GRACE/GRACE-FO data . The contribution of GIA uncertainty on our GRACE/GRACE-FO estimate is calculated considering a range of rheological parameters ; supporting information Table S1). This estimated GIA uncertainty will not take into account any deficiencies for a given GIA model, such as uncertainties in the model deglaciation history or assumptions of mantle lateral homogeneity (Caron et al., 2018;Geruo et al., 2013). The uncertainty in atmospheric circulation is estimated using the difference between an ensemble of reanalysis outputs  supporting information Figure S3). We estimate the contribution from uncertainty in ocean circulation using the difference between MPIOM outputs provided by the GRACE/GRACE-FO project and ocean bottom pressure outputs from ECCO-JPL near-real-time simulations (Fukumori, 2002; supporting information Figure S4). The uncertainty in terrestrial water storage is estimated using the RMS between GLDAS land surface models (Rodell et al., 2004; supporting information Figure S5).

Results
We calculate monthly maps of ice mass change and associated sea level response from the GRACE/GRACE-FO spherical harmonic data after correcting for GIA, terrestrial water storage, and sea level feedback. From our monthly maps, we calculate the average ice mass and sea level change rates using a least squares model following Velicogna et al. (2014). Figure 3 shows the GRACE/GRACE-FO-derived ice sheet change rates, dM∕dt, for 2002-2019 after correcting for the induced sea level response. We determine the significance of each point in the GRACE/GRACE-FO-derived map using a variant of the Akaike Information Criterion (AIC c ) for use with small sample-sized data sets (Burnham & Anderson, 2002;Velicogna et al., 2014). In Antarctica, the land ice mass change is concentrated most heavily over the Amundsen Sea Figure 2. Flowchart of our processing scheme for estimating self-consistent ice mass change and regional sea level estimates from time-variable gravity data. Purple nodes denote input data sets, orange nodes denote calculations, and the green node denotes the final converged solution.   Figure 3b). The interpretation of the spatial pattern of ice mass change is limited by the spatial resolution of the GRACE/GRACE-FO data (∼350 km) compared to the size of individual glaciers . Figure 4 shows the rate of change in induced sea level, dS∕dt, from changes in mass of glaciers and ice sheets over the same period. If left uncorrected, the redistribution of sea level would bias our GRACE/GRACE-FO ice mass change estimates, particularly in regions adjacent to large mass change. The pattern of sea level change differs strongly from the global average sea level change rate if the mass change was uniformly redistributed, which is delineated in blue in Figure 4 and supporting information Figure S7.
We calculate the mascon-derived time series of Greenland and Antarctic ice sheet mass change and equivalent sea level contribution both before and after correcting for the induced sea level response (Figures 5a-5d). We calculate the average mass change over the entire period by simultaneously fitting a least squares model with constant, linear, and quadratic terms with annual, semiannual, and 161 day oscillating terms to the mass time series to test the effect of different sea level redistributions . The 161 day oscillating terms are included to account for aliasing of the S2 tidal constituent in the monthly GRACE/GRACE-FO fields (Chen et al., 2009). For each resultant estimate, the measurement errors, localization uncertainties and geophysical correction uncertainties are added in quadrature following Velicogna et al. (2014). In Greenland, we find an average ice mass change, after correcting for GIA, terrestrial water storage, and sea level feedback, of −263 ± 23 Gt/yr for the 2002-2019 period (Figure 5a). The Antarctic ice sheet is in a state of overall mass loss during the GRACE/GRACE-FO period (Figures 5b-5d). The average ice mass change from Antarctica is −90 ± 48 and −122 ± 50 Gt/yr for the entire GRACE/GRACE-FO period using GIA coefficients from the Ivins et al. (2013) and Whitehouse et al. (2012) models, respectively. The acceleration in Antarctic ice sheet mass change is not dependent on the GIA correction because the signal is secular over time scales of years to decades . We find that the mass loss from Antarctica is accelerating at 8 ± 5 Gt/yr 2 from 2002-2019. Without correcting for the induced sea level response, our GRACE/GRACE-FO estimates of mass change would be biased high, or more negatively in the case of mass loss (Table 1 and Figures 5a-5d). For Greenland, the mass difference is 9 Gt/yr between the uncorrected and the SLF corrected estimates with an acceleration difference of 0.3 Gt/yr 2 . In Antarctica, the mass difference is 3-5 Gt/yr depending on the GIA correction.

Discussion
We use solutions to the sea level equation to quantify the regional sea level variation for glaciated regions from GRACE/GRACE-FO solutions (Hill et al., 2011;Riva et al., 2010). Here, we use the induced sea level distribution to help improve the original mass balance estimate to reach a self-consistent solution. The SLF correction to the ice mass balance is small over the entirety of the ice sheet (∼4% in Greenland and ∼5% in Antarctica). In Antarctica, locally positive and locally negative sea level change rates compensate to partially negate the overall SAL effect (Figure 4a). The SAL effect can be more significant at smaller spatial scales and for regions surrounded by regions of large mass change (supporting information Figures S8 and S9). More in-depth investigations of regional mass change will be possible as the time series of time-variable gravity extends and data processing methods improve (Mohajerani et al., 2018). Investigating smaller glaciated regions with less signal strength requires the geophysical corrections to be applied in a rigorous manner. There is both month-to-month and long-term variability in the SLF signal (Figures 5e-5h). While the SLF signal is presently small compared to the overall uncertainty (Table 2), the effect is cumulative and acts as a bias. Leakage of the SLF signal will increase with increasing glacier and ice sheet mass change (Moon et al., 2018). The effects of sea level redistribution will impact the mass balance time series from GRACE/GRACE-FO regardless of the processing method, which includes the new mascon solutions provided by GRACE/GRACE-FO processing centers as they are not presently corrected for land-sea fluxes. The magnitude of the SLF correction for a GRACE/GRACE-FO estimate will depend on the GRACE/GRACE-FO processing center, and the processing methods used to derive ice mass change (Velicogna & Wahr, 2013).
Modifications of the sea level correction technique can easily be incorporated into other spherical harmonic methods for regional mass change, and to help partition the new mascon solutions between the ocean and land components. The method described here will account for both climatic variations and the proper acceleration signal of glacial mass loss . In addition, deriving the SLF signal from GRACE/GRACE-FO data will take into account regional and rapid variations in mass change, and include regions where auxiliary data sets may not be available.
The sea level change induced from glaciers impacts estimates of Greenland ice sheet mass balance from GRACE/GRACE-FO (Table 1 and supporting information Table S2). In Table 1, the ice sheet SLF includes the mass change of just Greenland and Antarctica, and the total SLF includes both glacier and ice sheet mass contributions (supporting information Figure S7). The sea level change from peripheral glaciers affects our Greenland estimate by 2 Gt/yr. This indirect geophysical leakage component is currently small compared to uncertainty in GIA but could increase with time within continuing glacial mass losses (Moon et al., 2018). For Greenland and Antarctica, leakage from far-field signals is currently negligible compared to near-field sea level signals. Indirect leakage from the sea level change induced by the Antarctic ice sheet affects Greenland estimates by 0.2%. The sea level response from glaciers and ice caps affects Antarctic estimates by 0.4% and the sea level response from Greenland ice sheet mass change affects Antarctic estimates by 0.5%.
In this study, we estimate the effects of terrestrial water storage (TWS) fluctuations using outputs from the GLDAS NOAH land surface model. We redistribute the TWS mass change over the ocean taking into account the SAL effects. The indirect geophysical leakage from terrestrial hydrology predominantly affects the seasonality of our mass change estimates, particularly when taking into account SAL effects (Velicogna & Wahr, 2013; supporting information Figure S6). Reager et al. (2016) report that increases in non-glacial TWS withdraw 0.71 ± 0.20 mm/yr from global average sea levels over the GRACE/GRACE-FO period. The increase in water storage is comparable to the non-ice sheet glacier mass losses over the period but could be due to fluctuations in climate rather than being reflective of a long-term trend (Reager et al., 2016). This accumulation in TWS could also be partially due to changes in groundwater withdrawal and other anthropogenic changes in land water storage, which would not be estimated in GLDAS (Reager et al., 2016). Using GRACE/GRACE-FO-derived estimates of TWS change or placing mascons over the major TWS change regions would improve the method by including the anthropogenic TWS effects .
While the GRACE mission provided data for over 15 years and GRACE-FO has been providing data for over a year, the period of observation is still short compared to some climatic oscillations . While the GRACE-FO mission is not expected to last as long as the original GRACE mission, extending the time series of gravity measurements will help improve estimates of ice sheet mass change, help improve projections of sea level rise, and help validate ice sheet model estimates. The longer time series of time-variable gravity is necessary to help separate climate oscillations from accelerating gains or losses Velicogna et al., 2014). However, the cumulative effect of the sea level signal will be increasingly important for gravimetric studies with a longer time series of observations, particularly with sustained and increasing mass losses (Shepherd et al., 2012(Shepherd et al., , 2018.

Conclusions
The Greenland and Antarctic ice sheets are presently contributing to sea level and that contribution has magnified since the start of the GRACE mission. The effects of self-attraction and loading are a potential bias between different gravimetric ice mass balance estimates, particularly at smaller spatial scales. Leakage from sea level change will be an increasingly important factor in mass balance estimates from time-variable gravity with continued glacier and ice sheet mass change. We enumerate a method for quantifying and resolving this potential bias in gravimetric estimates that results in fully compatible mass change and sea level calculations. We investigate the effects of different treatments of the sea level response on regional and total ice sheet mass balance. We show that the sea level response is a minor correction compared to the effects of GIA but is likely correctable using time-variable gravity data alone. In addition, the sea level adjustment correction will likely become increasingly important with continued glacier and ice sheet mass losses. Lastly, we find that estimates of ice mass variations using data from the GRACE-FO mission are largely consistent with estimates from the GRACE mission, which enables the extension of the ice mass balance time series going forward.