Greenland Ice Sheet Mass Balance (1992–2020) From Calibrated Radar Altimetry

We present the first 1992–2020 record of Greenland Ice Sheet (GrIS) mass balance derived from multisatellite Ku‐band altimetry. We employ an empirical approach as an alternative detailed to radar‐propagation modeling, and instead convert elevation changes observed by radar altimetry into mass changes using spatiotemporal calibration fields. This calibration field is derived from a machine learning approach that optimizes the prediction of a previously published mass balance field as a function of ice sheet variables. Our mass balance record shows a GrIS contribution of 12.1 ± 2.3 mm sea‐level equivalent since 1992, with more than 80% of this contribution occurring after 2003. Our record also suggests that the 2017 hydrological year is the first year in the 21st century which, within uncertainties, the GrIS was in balance. Overall, the 28‐year radar‐derived mass balance record we present highlights the potential of the method to provide operational mass balance estimates derived from multisatellite Ku‐band altimetry.

The VMB method uses repeated measurements of ice-sheet surface elevation acquired by laser or radar altimetry. Most VMB estimates of GrIS mass balance to date are derived from laser altimetry using ICESat measurements (Sørensen et al., 2011;Zwally et al., 2011). Compared to both GMB and IOMB estimates, the ICESat-derived VMB provides unparalleled spatial resolution of ice-sheet mass balance, but this record only covers the 2003-2009 period (Sasgen et al., 2012). The launch of ICESat-2, in 2018, now extends the laser altimetry record. Using both ICESat and ICESat-2 measurements, Smith et al. (2020) estimated a VMB of −200 ± 12 Gt per year over the 2003-2019 period. Relatively poor understanding of the temporal changes in near-surface compaction rates, however, remains a challenge of the VMB.
The European Space Agency's (ESA) rich tradition of Ku-band radar altimeter missions provides a potential for radar altimetry to provide even longer VMB time series. The first Earth-observing satellite to enter a polar orbit was ERS-1 in 1991. Since ERS-1, ESA has ensured an unbroken series of Ku-band satellite altimeters with ERS-2, ENVISAT, and currently CryoSat-2 and Sentinel-3A/B. One challenge when estimating radar-VMB using Ku-band altimetry is that the relatively long Ku-wavelength results in appreciable surface penetration over snow/firn covered ice-sheet areas. This surface penetration depth can be dependent on local near-surface ice-sheet properties and climatology (Nilsson et al., 2015). In theory, this penetration can be as much as 15 m in cold, dry firn (Rémy et al., 2015). The issue for radar-VMB is that temporal variations in surface-penetration depth are directly imprinted on elevation measurements. Different strategies have been applied in the processing of Ku-band radar altimetry to suppress the influence of variable surface penetration, for example corrections for waveform parameters (Flament & Rémy, 2012;Simonsen & Sørensen, 2017), specific algorithms to track the surface height from the radar return waveform (retrackers) (Helm et al., 2014;Nilsson et al., 2015) or deconvolution of the waveform into surface and volume scattering components (Arthern & Wingham, 2001;Slater et al., 2019). Simonsen and Sørensen (2017) and Slater et al. (2019) both found approximately 5 cm per year biases between elevation changes observed by satellite radar altimetry and airborne laser altimetry. Over the area of the GrIS, such a bias translates into a mass change of 26-78 Gt per year, depending on the assumed density of the near-surface ice sheet. This has been a source of uncertainty when estimating GrIS mass balance from Ku-band measurements. The elevation change measured by the satellite altimeters is the combined result of SMB and dynamical imbalance, and the altimeters cannot provide information on the separation of the two components. This constitutes the major obstacle for computing mass balance from altimetry because the densities needed to perform the conversion from volume to mass depends entirely on the underlying glaciological process. This is handled in different ways in the literature (Mcmillan et al., 2016;Shepherd et al., 2019;Sørensen et al., 2011;Zwally et al., 2011). Furthermore, laser and radar altimeters are affected in different ways by the two processes because the laser tracks the snow-air interface and therefore is more sensitive to SMB changes, while the radar penetrates the snow and is more prone to track dynamical processes even in the accumulation zone.
Here, we employ an empirical approach to derive a mass balance record that incorporates measurements from multiple Ku-band altimetry satellites. The fundamental innovation of our approach is to develop a spatial cross-calibration field between the ICESat-VMB (Sørensen et al., 2011) and radar-VMB during their common 2003-2009 operating period. This permits incorporating apparent elevation changes that are not related to mass changes into the calibration field, including changes in radar penetration. This calibration field is input into a machine learning approach that finds correlations between the calibration field and independent a-priori data, such as ice velocity and changes in firn-air content. This enables us to estimate spatially and temporally varying volume-to-mass conversion fields for all radar observations during the 1992-2020 period. The strength of this empirical approach is that it give an alternative to radar propagation modeling, and is thus a candidate for operationalization within the ESA climate change initiative (CCI) or the European Commission Copernicus Climate Change Service (C3S). Sørensen et al. (2018) combined all available Ku-band altimetry observations (ERS-1, ERS-2, ENVISAT, and CryoSat-2) to provide the longest possible record of GrIS-wide elevation changes. This time series has since been extended to 1992-present within the C3S framework with the inclusion of the Copernicus Sentinel-3A mission. These data are now available in open access at the Copernicus Climate Data Store (www. cds.climate.copernicus.eu) and cover the Greenland ice sheet at 25 km grid resolution. Here, a modified C3S-processing chain is used to increase the spatial resolution of these data 5 km grid resolution (see Supporting Information for details). The mass change is given as rates, and in order to do so, we have used a temporal window to aggregate the radar observations. For the old satellites (ERS-1, ERS-2, and ENVISAT) this window is 5 years and for the present-day satellites (CryoSat-2 and Sentinel-3) it is 2 years. This results in a smoothening of the time series to resolve climatic changes and not seasonal weather. In addition to this elevation change record, we employ five mass balance time series: (1) ICESat-VMB, the average 2003-2009 GrIS mass balance derived by Sørensen et al. (2011). (2) IMBIE-MB, the monthly 1992-2018 IMBIE-2 record of reconciled GrIS mass balance (Shepherd et al., 2020). (3) GRACE-GMB, the monthly 2002-present GrIS gravitational mass-balance record from the combined processing of GRACE/GRACE-FO from the ESA GrIS CCI project based on Barletta et al. (2013). (4) Colgan-IOMB, the annual mass-balance record from 1995 to 2020 based on Colgan et al. (2019). (5) Mouginot-IOMB, the annual mass-balance record from 1972 to 2018 based on Mouginot et al. (2019). The first record, ICESat-VMB, is used as input in the cross-calibration of the VMB, similar to Xu et al. (2016), while the latter four records are used for independent validation of the resulting radar-VMB for the GrIS.

Cross-Calibrated ICESat and ENVISAT Mass Balance (2003-2009)
We establish an empirical volume-to-mass conversion between laser and radar altimetry measurements over the 2003-2009 period. We assume that the laser altimeter onboard ICESat observed the true elevation change of the ice sheet surface, that this elevation change can be converted into ice-sheet mass balance after correcting for processes that influence ice-sheet elevation but not mass: (1) vertical bedrock movement caused by both viscous deformation and elastic rebound of the Earth's crust, (2) time-varying biases in the ICESat laser altimeter, (3) changes in firn-air content. Subsequently applying a suitable near-surface density allows this corrected elevation change to be interpreted as VMB (Sørensen et al., 2011;Zwally et al., 2011). The ICESat-VMB provides a high spatial resolution, but sparse spatial sampling density, an estimate of the average mass balance during the period 2003-2009, and shows a net GrIS mass loss of 220 ± 28 Gt per year, not including the peripheral glaciers and ice caps.
For comparison, we initially estimate a radar-VMB based on ENVISAT-derived elevation change during 2003-2009 and assume firn density of 350 kg m −3 in the high-elevation accumulation area and ice density of 917 kg m −3 in the low-elevation ablation area (see Figure 1a for zonation). This initial radar-VMB, which does not include elevation corrections terms Sørensen et al., 2011) infers approximately half this mass loss (c. −108 Gt per year) during the same period. The majority of this difference may be attributed to the underestimation of the mass loss in high-sloping areas by the radar. However, we assume that the spatially variable ratio between the initial radar-VMB and the ICESat-VMB can be used as an empirical calibration field for the volume-to-mass conversion. By calibrating the coincident mass change fields and not elevation change fields, we will be able to circumvent the need for temporal modeling of surface densities which has been shown to be the main source of uncertainty in VMB results. This approach is analogous to previous studies enhancing the spatial resolution of GRACE-GMB by adding a-priori information from satellite altimetry and modeled surface mass balance (Colgan et al., 2015;Su et al., 2015;Xu et al., 2016). Figure 1b shows the calibration field derived for the temporally overlapping ICESat and ENVI-SAT VMB estimates during 2003-2009.

Time-Varying Calibration Field from Machine Learning
While the 2003-2009 calibration field applied to the full record of Ku-band radar satellite altimetry elevation changes (1992-2020) result in a radar-VMB record which agrees with annual IOMB and IMBIE assessments, it does not capture the high intraannual variability of the GRACE-GMB record (see supporting information, Figure S3). These inconsistencies between the records can likely be attributed to temporal variability in the relation between ICESat-VMB and radar-VMB for which we do not account. ICESat-2 is the obvious candidate to provide more information on the temporal evolution of the calibration field, but the ICESat-2 data record is not yet long enough to provide an independent calibration field. Instead, we seek to describe the calibration field as a function of other a-priori ice-sheet data, to build time-varying calibration fields. Using machine learning algorithms, we investigated ten a-priori datasets (location, surface elevation, surface slope, firn-cover, ice-sheet surface velocity, grid-interpolation distance, satellite sensing mode, radar elevation change, and firn air correction) to derive an empirical prediction model of a time-varying calibration field (see Supporting Information for a description of the individual datasets and tested machine-learning algorithms). These a-priori datasets introduce 768 plausible prediction models for the volume-to-mass calibration which are generated using machine learning. The preferred prediction model is found by optimizing agreement with independent IMBIE-MB and GRACE-GMB estimates, agreement with the static calibration field (2003)(2004)(2005)(2006)(2007)(2008)(2009), and an assessment of glaciological meaningfulness of spatial patterns in the 1990s. The preferred prediction model utilizes a-priori knowledge of latitude, surface slope, ice velocity, firn-air correction, satellite sensing mode, and observed radar elevation change. Given this a-priori knowledge, the monthly calibration fields differ as a function of the applied firn-air correction and observed radar elevation change itself, both of which vary at monthly time-steps. The series of calibration fields are then multiplied with the observed monthly radar elevation change given in pseudo-mass by the density assumption above. The uncertainty of the final radar-VMB estimate is found by combining the error estimate from the volume change observations and the error associated with the nonuniqueness of the machine learning prediction as uncorrelated errors. The volume error component is converted into mass SIMONSEN ET AL.  by the density assumption used for the initial radar-VMB, whereas the preferred prediction model has been run 180 times to derive the standard deviation of model predictions. It should be noted that shortening of the temporal window with the introduction of CryoSat-2, allows for a large imprint of weather variability on the stability of the plane-fitting algorithm. Hence, the VMB error estimate is larger post-2010.

Results
Applying the time-varying calibration fields to the 28-year Ku-band altimetry record produces the GrIS radar-VMB time series shown in Figure 2 (cyan curve), the underlying data is available tabulated in hydrological years in the supporting material, both for the full ice sheet and separated into major drainage basins. Annual maps of the GrIS radar-VMB are also available there. For intercomparison, Figure 2 also includes the 1972-2018 Mouginot-IOMB (red curve), the 1995-2020 Colgan-IOMB (magenta curve), and a 2-year running-mean mass-change rate based on the GRACE-GMB record (green curve) to match the temporal-averaging of the radar-VMB. Our radar-VMB agrees with each of these independent records within their individual uncertainties (see Figure S4 for all available mass balance time series). In terms of the GrIS contribution to global sea-level rise, the radar-VMB series infers a total of 12.1 ± 2.3 mm of eustatic sea-level contribution since 1992. Remarkably, 84% of this ice loss has occurred after c. 2003, with a record high of 3.5 ± 0.3 mm (28%) in the years 2010-2012. This high mass loss period is followed by a period of reduced mass loss (c. 2013-2017), before a sharp increase in mass loss after c. 2019.
The GrIS mass loss records differ slightly during the ICESat-era (2003ICESat-era ( -2009, with an average rate of mass loss of 220 ± 28, 206 ± 63, and 198 ± 48 Gt per year for ICESat-VMB, IMBIE-MB, and radar-VMB, respectively. These differences mainly originate before 2005, but might also incorporate spatially compensating errors. SIMONSEN ET AL.  Including GRACE-GMB in the analysis, we see agreement between the GRACE-GMB and ICESat-VMB on one hand and IMBIE-MB and radar-VMB on the other. After 2005, the midpoint of the ICESat-VMB record, the two independent datasets (GMB and radar-VMB) show an increasing correlation with time. This can be attributed to the introduction of CryoSat-2 data into the radar-VMB record in 2010, resulting in shorter temporal averaging that enables more interannual variability to be resolved within the radar-VMB record. Additionally, we observe changes in the agreement with the IMBIE time series (supporting material), which likely stem from the different sources for IMBIE values over the entire 1992-present period. Prior to 2003 IMBIE only consists of the two IOMB estimates and the radar-VMB has a better agreement with Mouginot et al. (2019Mouginot et al. ( ) prior to 2003Mouginot et al. ( , whereas in 2003Mouginot et al. ( -2005 the radar VMB agrees better with Colgan et al. (2019). This highlights the value of our novel altimetry-derived time-series, the first altimetry time series with consistent methodology over the entire IMBIE period. Including more homogeneous times-series in IMBIE will permit exploring the influence of changing data sources overtime on the IMBIE best-estimate time-series.
In addition to providing a long time series of integrated GrIS mass loss, the radar-VMB record can also infer spatial variations in GrIS mass change at high resolution. Figure 3 shows the spatial distribution and temporal evolution of mass loss, with total mass losses of 57 ± 43 Gt per year in the period 1992-2000, 163 ± 41 Gt per year in the period 2000-2010, and 241 ± 58 Gt per year in the period 2010-2020. In the 1990s, mass loss was focused at the main outlet glaciers, with a slight mass gain observed in the central and northwestern of the ice sheet. In the 2000s, mass loss accelerated and progressed inland from these outlet glaciers. The second stage of the acceleration occurs at the beginning of the 2010s, with a peak in the abnormal melt year of 2012 (Nghiem et al., 2012), after which the ice sheet seems to enter a new quiet phase (2016)(2017)(2018) resembling the pre-2000s.

Discussion
Our radar-VMB relies on spatiotemporal calibration fields to capture surface processes and their consequences on radar altimetry. Examining the static 2003-2009 calibration field (Figure 1b) we see a band of positive-signed calibration values inland from Melville Bay in Northwest Greenland. This positive signal is found in an area where no change in neither accumulation rate nor firn-air content was observed during the ICESat-era . Therefore, we believe that it is likely a consequence of a change in SIMONSEN ET AL.
10.1029/2020GL091216 6 of 10 the radar penetration depth in the area, possibly due to interannual variability in ice lens formation within the near-surface firn . In other areas, negative-signed calibration values are observed in the firn area (equivalent to negative pseudo densities in the conversion from volume to mass ). These suggest that radar altimetry sees a surface elevation increase, where ICESat observes the physical surface is decreasing. This can directly be interpreted as a result of a warming climate, in which the firn densification occurs faster. This lowers the surface elevation seen by ICESat, while the increasing near-surface firn density results in an upwards movement of the radar scattering horizon that is observed as an elevation increase in radar altimetry.
In addition to the glaciological interpretations of the calibration fields, insights into the physical nature of the radar retrieval are also observed. In the high-sloping regions at the margin of the ice sheet (most predominant in the southern and eastern ice sheet), the calibration field most often shows large positive values (>2). In these areas, the surface density is typically high and the high calibration has a relatively high impact on the overall mass loss of the final radar-VMB estimate. One should also note that the polar orbit of ICESat has difficulties in mapping elevation change in the southern part of the ice sheet. However, the large calibration values can mainly be explained as a consequence of the radar altimeter measuring the elevation of the point-of-closest-approach within the pulse limited footprint, possibly biasing the elevation change estimate toward higher elevation and less negative elevation change values. The Greenland south-east shows this most prominent where (as seen in the figure) steep topography of Kangerlussuaq, Helheim, and Køge Bugt glaciers result in the need for the largest calibration field. As the glacier flow here is restricted by nunataks, which again is the highest point within the radar footprint. Different approaches can be found in the literature to extrapolate the elevation change to the coast and counteract this positive bias by using the information on ice velocities (Hurkmans et al., 2012;Mcmillan et al., 2016). Our machine learning approach confirms this relation between ice velocities and elevation changes by favoring prediction models including the ice velocities as a priori information.
Ku-band radar missions provide continuous altimetry data since August 1991. These data are a valuable, albeit largely untapped, approach to bridge gaps in other satellite records, such as the GRACE/GRACE-FO and the ICESat/ICESat-2 missions. We find an average radar-VMB in the period 2003-2019 of −221 ± 51 Gt per year. Smith et al. (2020) provides the first mass loss estimate from combined ICESat/ICESat-2 observations, with an estimate of −200 ± 12 Gt per year. Although these two estimates agree within reported uncertainties, our slightly larger mass loss might be caused by the abnormal melt year of 2012, during which solid ice was melted and then subsequently replaced by snow, resulting in altered effective surface densities. Similarly, for the GRACE/GRACE-FO transition, we see the potential of radar-VMB to bridge these two satellites. There is no apparent bias in radar-VMB during these two gravimetry missions (insert Figure 2). This radar-VMB bridging of the GRACE/GRACE-FO gap agrees well with the Colgan-IOMB time series, which suggests that the machine learning approach has been able to capture much of the dynamics related to the rapidly evolving firn and radar penetration on the GrIS, but the ultimate answer to the long-term durability of the method will come when ICESat-2 has been operated long enough to provide an elevation change calibration field for validation. For now, the 28-year radar-VMB product that we present appears to have the potential to provide operational mass balance estimates for the GrIS.
The radar-VMB time series is explicitly calibrated to match the mean annual ICESat-VMB over the 2003-2009 period. In the ICESat-VMB processing algorithm, interannual variations in surface elevation are suppressed to derive climatological VMB. This hampers the radar-VMB's ability to resolve the intraannual variability seen in the GMB record, despite being a monthly product. With the inclusion of CryoSat-2 data after 2010, which greatly improves the spatiotemporal data coverage at the grid-cell level, we can solve for elevation changes at 2-year intervals. This improves the ability to capture interannual variability in the VMB time series. When compared with the original volume change time series (supporting material) we also observe the spatiotemporal calibration fields to increase the magnitude of the interannual variability for the newer satellites. Technical limitations of the older satellites (ERS-1, ERS-2, and ENVISAT) prevent this, and we are only capable of resolving 5-year intervals. The long-term radar-VMB record is therefore valuable for gaining insights into the climatological trends of the GrIS in the early years, while in the later years (and into the future) it also gives more insight into short-term mass balance variability. This improvement over time highlights how the novel orbit constellation and instrumentations of CryoSat-2 has transformed our ability to monitor GrIS mass balance. The CryoSat-2 satellite also introduces the novel SARIn sensing mode, and the machine learning informs us that the calibration should only be applied when the satellites operate in the Low-Resolution Mode (LRM) mode. This correction change in the middle of the 2010s gives rise to a discontinuity of the mass loss in Figure 3c, from the overestimation of the 2012 melt event in the LRM area of basin 6. The different time series of GrIS mass loss generally exhibit similar behavior at basin scales (see supporting information). However, the basin analysis suggests the radar-VMB to be better at estimating the ice sheet mass balance in areas of adjacent mass changes, which may leak into the GMB solution. Whereas it also points to the extreme year of 2012 not to be particularly well captured by the training period (2003)(2004)(2005)(2006)(2007)(2008)(2009), and thereby the mass loss is overestimated in specific basins (see supporting material). This might prove to be a major limitation to the method validity in a warming world, but performing a recalibration using ICESAT-2 data in the near future might help to also capture these high melt years.
A prominent feature in the radar-VMB record is the sharp increase in inferred ice-sheet mass balance between 2016 and 2017. This reflects an abrupt increase in surface mass balance, which is due to a combination of anomalously high cyclonic snowfall in fall 2016 and anomalously low surface melting in summer 2017. The HIRHAM5 regional climate model estimates that during October and November 2016 alone, cyclonic storms delivered c. 100 Gt more snowfall to the ice sheet, especially East Greenland, than average for the 1981-2010 period (Mottram et al., 2018). Positive North Atlantic Oscillation conditions during the subsequent 2017 melt season promoted cyclonic conditions that reduced incoming solar radiation, through increased cloudiness, and consequently decreased summer melting and increased summer snowfall. In the ablation area, summer ablation was average, or below average, at all 20 PROMICE sites (Tedesco et al., 2017). The high snowfall of fall 2016 likely also contributed to preconditioning the ice sheet for anomalously low surface melt in the summer of 2017. A thick spring snowpack can decrease subsequent summer melt by delaying albedo feedback processes, such as meltwater runoff and bare ice exposure, until later in the melt season (Box et al., 2012). Within uncertainties, the radar-VMB suggests that the high fall 2016 snowfall and low summer 2017 melt, likely made the 2017 hydrological year the first year in which the ice sheet was in balance during the 21st Century (Mattingly et al., 2018). In the 28-year radar-VMB record, only four other hydrological years are suspected to be in balance (1995, 1996, 1997, and 2018). Thus, 2017-2018 reflects a return to conditions not seen since the mid-1990s (The prediction intervals for all hydrological years are shown in Figure S4). Perhaps more remarkable is the rapid shift in ice sheet mass balance from record low to record high in only 5 years. This reflects a shift of c. 400 Gt per year, or more than 3σ when compared to the variability in the 28-year radar-VMB record, in 5 years. The extreme mass loss year of 2012 has been previously attributed to persistently negative NAO conditions. This results in atmospheric ridging across Greenland, with a stationary high-pressure cell that contributes to clear-sky conditions, elevated air temperatures, and surface melt across the ice sheet, with relatively little summer snowfall (Jeffries et al., 2012). The rapid shift between extreme low and high mass balance years, associated with the NAO, highlights the tremendous influence of short-term variability in atmospheric conditions on ice-sheet mass balance.

Conclusions
We present a new empirical approach for deriving operational mass balance estimates of the GrIS from multisatellite Ku-band altimetry. This approach relies on spatiotemporal calibration fields to convert changes in radar elevation to changes in mass. The resulting record shows a total eustatic sea-level rise contribution of 12.1 ± 2.3 mm from the ice sheet since 1992, with more than 80% of this contribution occurring after 2003. This 28-year record also shows the value of producing unbroken multisatellite records to bridge gaps between individual missions, such as ICESat/ICESat2 and GRACE/GRACE-FO, and improve the consistency of long-term ice-sheet monitoring. The 2017 hydrological year appears to be the first year in the 21st Century which, within uncertainties, the ice sheet was in balance. Contrasting this extreme high mass balance with the previously well-described 2012 extreme low mass balance highlights how changes in atmospheric conditions, especially the North Atlantic Oscillation, can rapidly shift the mass balance of the GrIS over 3σ of variability. Our empirical approach is an alternative to detailed radar-propagation modeling, yet reproduces multiple independent IOMB and GMB estimates within a reasonable uncertainty. This approach provides the first continuous 28-year ice-sheet mass balance record from Ku-band satellites and presents the potential for assessing operational mass balance estimates from satellite radar altimetry going forward.

Acknowledgments
The initial version of this study was conceived as a part of the ESA CCI Sea Level Budget Closure Project (ESA contract number 4000119910/17/I-NB). The long time-series radar altimetry was developed in the ESA CCI + Greenland Ice Sheet Project (ESA contract number 4000126523/10/I-NB) and operationalized within the C3S project (Framework agreement ECMWF/ COPERNICUS/2018/C3S_312b_Lot4_ EODC). William Colgan was supported by the Programme for Monitoring the Greenland Ice Sheet (PROMICE).