2.1 Geosat

The U.S. Navy launched the GEOdetic SATellite (Geosat) in March 1985, which operated until September 1989, providing limited Antarctic coverage between ±72^∘ latitude. The main goal of the mission was to provide the U.S Navy with detailed information about the marine gravity field. Geosat operations consisted of two separate missions, in which the initial 18 months were the classified “geodetic mission” (GM), in a 135 d repeat orbit ending in September 1986, and the “exact repeat mission” (ERM), in a 17 d repeat orbit lasting until the end of the mission. The mission carried a Ku-band (13.5 GHz) pulse-limited altimeter providing measurements every 670 m along track (10 Hz), with a pulse-limited diameter of ∼ 3 km. In this study we used “ice data record” (IDR) from the Radar Ice Altimetry Group at NASA Goddard Space Flight Center (GSFC) providing geolocated and corrected surface elevations. Only records with a valid retracking correction and waveforms containing a single return echo were used in the study to reduce noise in the derived surface elevations. We detected the presence of a bias in the automatic gain control (AGC) parameter of 1.23 dB between the geodetic mission (GM) and the exact repeat mission (ERM) phases. This is most likely due to the change in orbit, and it did not affect any of the other parameters, including the surface elevation change.

2.2 ERS-1 and ERS-2

The European Space Agency (ESA) launched the European Remote Sensing (ERS) satellites in 1991 (ERS-1) and 1995 (ERS-2), respectively. They operated continuously between ±81.5^∘ latitude until 1996 and 2003, respectively. Both missions carried conventional pulse-limited Ku-band (13.6 GHz) radar altimeters, with a pulse-limited footprint of ∼ 1.5 km and an along-track resolution of 370 m (20 Hz sampling rate). The two missions operated in a 35 d repeat orbit, though ERS-1 had several shorter mission phases early on that deviated from the standard repeat-track orbit. For this study the “REprocessing of Altimeter Products for ERS (GDR): 1991 to 2003” (REAPER), detailed in Brockley et al. (2017), is used to obtain surface elevation measurements. This product contains updated corrections and improved calibrations. For each satellite record we separated the data from the two operational modes, “ocean” and “ice”, excluding any data used for calibration. The product provides different retracking solutions from which we have chosen to use the ICE1 retracker, otherwise known as the “offset center of gravity” (OCOG) retracker (Wingham et al., 1986), using a 30 % threshold of the maximum waveform amplitude. The Ku-chirp and the ICE-1 20 Hz quality flags, available in the product, were used to exclude poor-quality observations from the analysis.

2.3 Envisat

The “Environmental Satellite” (Envisat) was launched by ESA in 2002 as a successor to the ERS mission and was officially decommissioned in 2012. Envisat was launched into a 35 d repeat orbit, operating with a pulse-limited radar altimeter with the same footprint, radar frequency, and sample frequency as the earlier ERS missions. For Envisat we used the “RA-2 Geophysical Data Record” (GDR) version 2.1. Only data collected during the period 2002 to 2010 were used due to changes in orbit initiated in October of 2010. The GDR product, as with the REAPER product, includes elevations determined using the ICE-1 retracker with a 30 % threshold of the maximum waveform amplitude, which we used for this analysis. We applied the same quality filter on the GDR records as with the ERS product, using the Ku-chirp and ICE-1 quality flags.

2.4 ICESat and ICESat-2

The National Aeronautics and Space Administration (NASA) launched the Ice, Cloud, and land Elevation Satellite (ICESat) in 2003, which operated from 2003 to 2009 in a 96 d repeat orbit. The mission carried a novel laser altimeter providing a 70 m beam-limited ground footprint, with 170 m along-track sampling (40 Hz). We used the latest version of the GLAS06 product (release 34), which has been corrected for the “Gaussian centroid offset” (Borsa et al., 2014) and detector saturation and converted to heights above the WGS84 ellipsoid. We did not apply any inter-campaign bias to the ICESat elevations as there is no consensus that these are required (Borsa et al., 2019). The records are further edited to remove poor-quality observations using the accompanying quality flags (elev_use_flh >0, sat_corr_flg >2, sigma_att_flg >0, i_numPk >0).

The ICESat-2 mission is a follow on mission to ICESat and was launched in October 2018 with the goal of continuing the long-term altimetry measurements of polar regions (Markus et al., 2017). It carries a new and novel photon counting laser altimeter that uses a 532 nm laser with a pulse repetition rate of 10 kHz and that operates in a repeat-track configuration over the continental ice sheets. In contrast to its predecessor's single beam, ICESat-2 collects ground measurements using six individual laser beams arranged in three pairs. Each of the beam pairs is separated by 3 km and each inter-pair beam by 90 m across track. This configuration allows for a direct estimate of the across-track surface slope that was not directly possible with ICESat's single beam configuration. The beam limited footprint for each beam is 12 m in diameter sampling every 0.7 m along track with a repeat frequency of 91 d. In this study surface elevation from the ATL06 product was used following the approach outlined in Smith et al. (2019, 2020). Here a segmentation filter (difference filter) was used to remove poor-quality observations if differences between consecutive points exceed a threshold of 2 m the point was rejected. Further editing was done using the ATL06 quality flag, keeping only data designated to be of good quality (“atl06_quality_summary  = 0”).

2.5 CryoSat-2

ESA's CryoSat-2 mission launched in 2010 with the primary purpose of monitoring changes in Earth's sea and land ice. This satellite carries a new type of delay Doppler radar altimeter (Raney, 1998) equipped with a dual antenna configuration allowing for interferometric measurements of surface elevations. The altimeter system, referred to as SIRAL, operates in two different modes over the ice sheets: a synthetic aperture radar interferometric (SARIn) mode over the marginal areas and a low-resolution mode (LRM) (a conventional Ku-band pulse-limited radar – identical to ERS and Envisat) over the ice sheet interiors. The Doppler/delay radar allows for increased along-track resolution compared to conventional pulse-limited altimetry. The SARIn mode has an effective resolution of 350 m along track and 1500 m across track, compared to the LRM modes' 1500 m along- and across-track resolution. Further, the dual antenna configuration allows for mapping of the exact position of the surface echo location by estimation of the across-track look angle from the difference in path length of the signals between the two antennas. In contrast to previous missions, CryoSat-2 operates in a drifting orbit, with a 369 d repeat and a 30 d sub-cycle. The drifting orbit offers improved spatial coverage compared to repeat-track orbits at the expense of larger across-track distances. We processed both the LRM and SARIn modes using the ESA L1b Baseline-C product for the time span 2010–2018 using a custom CryoSat-2 processor described in Nilsson et al. (2016). For the LRM mode we have chosen to use a 10 % threshold of the maximum waveform amplitude for retracking, similar to Schröder et al. (2019).

3.1 Slope-induced error correction

The largest source of error in radar altimetry over ice sheets is associated with the effects of surface slope inside the beam-limited radar footprint. This error stems from an inability to locate the surface from which most of the echo power originates (off nadir). Because of this, the echo is assigned the location of the sub-satellite point on the Earth surface. This introduces a slope-dependent measurement error on the order of 0–100 m (Brenner et al., 1983), which varies with the magnitude of the surface slope. There are a few ways of minimizing the slope-induced error (Bamber, 1994; Roemer et al., 2007). For this study we used the “relocation method” described in Nilsson et al. (2016). The relocation method corrects both the range and the coordinates to the echolocation (from nadir) using topographical information, such as surface slope, aspect, and curvature. This method has been shown to improve surface elevation retrievals compared to other approaches (e.g., Schröder et al., 2017). To compute the required surface slope, aspect and curvature, we used the “bedmap2” digital elevation model from Fretwell et al. (2013) resampled to 2 km horizontal resolution.

3.2 Elevation change estimation and algorithms

Surface elevation changes are determined as follows: the local mean topography within a specified search radius is removed from each mission and mode, leaving only the elevation anomalies that contain the time variable signal. Artificial trends and seasonal amplitudes in elevation anomalies that are introduced by changes in surface scattering characteristics are reduced proportionally to the correlation with the received radar waveform shape. Intermission biases in seasonal elevation anomalies are further minimized using a normalization scheme that references all seasonal elevation change amplitudes to those observed by CryoSat-2. A cross-calibration scheme is applied to adjust and merge elevation change from all missions and modes into a continuous monthly time series. Lastly, interpolation is used to generate a consistent gridded product with 1920 m horizonal resolution at monthly time steps from 1985 to 2020. The details of each step are provided in the following sub-sections.

3.2.2 Surface and volume scattering correction

The microwave pulses transmitted by spaceborne radar altimeters at Ku-band frequency (∼ 13.6 GHz) are sensitive to changes in the dielectric properties of the ice sheet surface (as determined by changes in the snow grain size, temperature, density, and water content, among others). Large-scale temporal and spatial changes in the scattering horizon induce changes in measured range, and thus surface elevation, and can introduce long-lived biases in the derived elevation change rates (Arthern et al., 2001; Davis and Ferguson, 2004; Khvorostovsky, 2012; Nilsson et al., 2015b; Wingham et al., 1998). To mitigate this effect, we use a retracking algorithm that tracks the leading edge of the return waveform (i.e., a maximum amplitude threshold between 10 % and 30 %). Such retrackers have been shown to be less sensitive to changes in ice sheet surface properties (Helm et al., 2014; Nilsson et al., 2016; Schröder et al., 2017). Another key step is removing elevation variability that is correlated with changes in the received radar waveform shape (Flament and Rémy, 2012; McMillan et al., 2014; Paolo et al., 2016; Simonsen and Sørensen, 2017; Zwally et al., 2005). The shape of the waveform is intricately linked to the medium in which it is propagated or reflected. Removing the correlation between changes in the shape of the radar waveform with elevation can largely reduce these artificial signals. For this study we approximated the shape of the radar waveform following the definition of Flament and Rémy (2012)​​​​​​​ and Simonsen and Sørensen (2017) using the backscatter (Bs), the leading-edge width (LeW), and the trailing edge slope (TeS).

The spatially variant scattering correction was estimated by computing the local sensitivity gradient (SG) between each waveform parameter and elevation residuals using a multi-variate least-squares inversion. The SG parameters were estimated for ascending and descending tracks separately. All waveform parameter time series were centered and normalized using the mean and standard deviation. Further, parameters were detrended by applying a difference operator, forming the following least-squares model:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathrm{\nabla }\left({\displaystyle \frac{h-\stackrel{\mathrm{‾}}{h}}{{\mathit{\sigma }}_{\mathrm{d}h}}}\right)& ={\mathrm{SG}}_{\mathrm{Bs}}\cdot \mathrm{\nabla }\left({\displaystyle \frac{\mathrm{Bs}-\stackrel{\mathrm{‾}}{\mathrm{Bs}}}{{\mathit{\sigma }}_{\mathrm{Bs}}}}\right)\\ & +{\mathrm{SG}}_{\mathrm{LeW}}\cdot \mathrm{\nabla }\left({\displaystyle \frac{\mathrm{LeW}-\stackrel{\mathrm{‾}}{\mathrm{LeW}}}{{\mathit{\sigma }}_{\mathrm{LeW}}}}\right)\\ & +{\mathrm{SG}}_{\mathrm{TeS}}\cdot \mathrm{\nabla }\left({\displaystyle \frac{\mathrm{TeS}-\stackrel{\mathrm{‾}}{\mathrm{TeS}}}{{\mathit{\sigma }}_{\mathrm{TeS}}}}\right),\end{array}\end{array}$$

where ∇ is the difference operator, h the elevation residual (elevation relative to time-invariant topography), σ the standard deviation, and the overbar the average value of the parameter.

The SG parameters were inverted for using the same adaptive search center approach as described in Sect. 3.2.1. The estimated SGs were then used to correct each observation within the search cap using the linear combination of the original waveform parameters and the estimated coefficients. Finally, we apply a linear space–time interpolation to estimate corrections at locations where the multi-variate fit did not provide a satisfactory solution.

Figure 2Original and scattering corrected area integrated time series for Lake Vostok in East Antarctica, which has been shown to have a height trend close to zero over recent decades (Richter et al., 2014). A discrepancy in uncorrected height trends is observed for the various missions due to differences in altimetry processing, orbit configuration, and the quality of the geophysical corrections. Envisat and ERS-2 (Ice) show the largest uncorrected magnitude in both trend and seasonal signal. Corrected height change records show significantly lower seasonal amplitudes and trends that are close to zero.

To determine the optimal search radius for generating the scattering correction, we performed a sensitivity study over Lake Vostok in East Antarctica (Fig. 2). Lake Vostok was selected due to its low surface slope, on average 0.03^∘, and highly stable surface (Richter et al., 2014), minimizing the impact of the static and time variable topography in the analysis. After varying the search radius from 1 to 5 km, we found that the 1 km solution provided the most accurate trend and seasonal amplitude for all missions and modes. We also found that the absolute magnitude of both the trend and amplitude increased linearly as the search radius increased. We interpret this result as a decrease in the efficiency of the correction possibly due to de-correlation with increasing spatial and temporal scales. The use of a 1 km search radius is also computationally efficient as fewer data are used in the inversion. Applying these lessons to the ice-sheet-wide processing, we found that the correction has a minor impact on the estimated trend for the CryoSat-2 SARIn mode and the Geosat missions. We also found that the application of the correction to the SARIn and Geosat data increased the seasonal amplitude of the local (single grid cell) time series (Sect. 3.2.4). Given that there is no physical justification for an increase in seasonal amplitude, we chose not to apply the correction to the Geosat mission and the SARIn-mode data. For the other missions, the magnitude of the correction varied across missions and modes of operation, in which the largest changes in trend and amplitude were found for Envisat and ERS-2 Ice mode and the lowest for CryoSat-2 LRM. By examining the changes in trend and amplitude we found significant spatial patterns, also varying across each mission and mode, as shown in Fig. 3. These patterns show strong correlations to both surface slope and roughness and signals of metrological origin (Armitage et al., 2014) and are mostly driven by katabatic winds and the re-distribution of snow. The wind effect can be observed in the RMSE plots for Envisat and for ERS-2 Ice as banded patterns in the East Antarctic sector following the main ridge lines.

Figure 3Change in elevation change rate and RMSE (seasonal amplitude) of the local time series after correction for temporal changes in scattering (penetration depth). Spatial patterns linked to surface conditions can be clearly observed. These effects are most prominent for Envisat and ERS-2.

3.2.3 Cross-calibration and integration

Removal of the time-invariant surface topography is done internally to each dataset such that elevation residuals are not aligned to the same surface (see Sect. 3.2.1). To align elevation anomalies to a common reference we first solve for intermission offsets. These offsets vary regionally (Khvorostovsky, 2012; Wingham et al., 2009; Zwally et al., 2005), depending on the underlying topography, physical interactions of the radar with the surface, and differing retracking methodologies. In contrast to previous studies (e.g., Davis and Ferguson, 2004; Davis, 2005; Khvorostovsky, 2012; Li and Davis, 2006; Schröder et al., 2019; Wingham et al., 2006, 2009; Zwally et al., 2005), we estimate these offsets using a least-squares adjustment. This approach allows for a simple, yet consistent, alignment of multiple relative elevation anomalies without requiring full overlap between missions to solve. The technique follows the approach of Bevis and Brown (2014), using the entire multi-mission record to constrain the solution while accounting for trend, seasonality, and intermission or mode offsets. The trend is represented by a polynomial, with a maximum order of six, a four-term Fourier series to account for seasonality, and Heaviside functions to solve for the intermission offset between missions and modes. The design matrix can be written as follows:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}h\left(t\right)& =\sum _{i=\mathrm{1}}^{n_p-\mathrm{1}}{p}_i(t-t_{\mathrm{r}}{)}^{i-\mathrm{1}}+\sum _{k=\mathrm{1}}^{n_{\mathrm{f}}}{s}_k\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }t}{T_k}}\right)\\ & +c_k\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }t}{T_k}}\right)+\sum _{j=\mathrm{1}}^{n_j}{b}_{j}h,\end{array}\end{array}$$

where n_p is the model order, t is the time in decimal years, t_r is the reference time in decimal years (t_r=2013.95), T_k is the seasonal period reference (T₁=1 and T₂=0.5), n_f is the number of Fourier series terms (n_f=4), and n_j is the number of missions and modes (n_j=10).

Here, we add offsets for 10 different missions and modes in the least-squares model (Geosat, ERS-1 Ocean and Ice, ERS-2 Ocean and Ice, Envisat, ICESat, CryoSat-2 LRM and SARIn, and ICESat-1) to all data falling within the search radius. To determine the order of the polynomial we use the Bayesian information criterion (BIC; Schwarz, 1978) to select the polynomial that produces the lowest BIC value estimated from monthly binned data.

The cross-calibration is performed on a 2 km polar-stereographic grid (EPSG: 3031) using a variable search radius of 1–10 km surrounding each grid cell. The radius is increased until 70 % of the time series is filled (monthly) or the maximum radius is reached. If the maximum search radius is reached and the 70 % criterion is not meet, we continue processing using all available data. In most cases the search radius is in the range of 2–10 km. Outliers in the original time series were initially removed using a 1-year running median filter in which values larger than 10 times the median absolute deviation (MAD) were rejected. The model is then fit to the time series using a robust least-squares inversion as in Sect. 3.2.1. Solutions are rejected if the absolute value of the linear rate is larger than 20 m a⁻¹ or if the RMS of the time series relative to the model is larger than 4 m. If any of the derived offsets are larger than 100 m, the offset is set to zero. The offsets estimated from the least-squares inversion are then subtracted from the time series providing an initial cross-calibrated record of elevation change. Further, a last outlier step is performed in which the model is used to filter the time series by omitting observations exceeding 10-times the MAD of the model residuals.

This approach has several advantages; it allows a first-order calibration of non-overlapping time series while also aligning overlapping missions and modes to their common mean. To account for time series that do not fully conform to our choice of a linear model, a secondary cross-calibration is performed for the four mission-specific offset coefficients (ERS-1 to ERS-2, ERS-2 to Envisat and ICESat, Envisat and ICESat to CryoSat-2, and CryoSat-2 to ICESat-2), using the post-fit model residuals. This approach was chosen as it facilitates the estimation of any residual offsets after removal of the majority of the trend and seasonality, making it simple to estimate the overall bias between the mission groups. The offsets for groups ERS-1 to ERS-2, ERS-2 to Envisat and ICESat, and CryoSat-2 to ICESat-2 were estimated by taking the median difference between the two datasets over their respective overlapping time periods. This approach was found to be suboptimal for the Envisat and ICESat to CryoSat-2 offsets due to the short period of overlap (less than 4 months) and large changes during the time period 2009–2011. To overcome this limitation, we applied three different methods, generating five different independent Envisat and ICESat to CryoSat-2 offsets at each search node. For method 1, we fit two second-order polynomials to the two residual time series and compute the median offset between the two functions over a 1-year overlap (2010–2011), as well as the difference between the two intercepts of the polynomials. For method 2, we applied a Kalman smoother with a state-space model consisting of a constant local level and a random walk trend (Kalman, 1960; Shumway and Stoffer, 1982) that better accommodate the variability in the time series. The filter was initialized with a variance rate of 1 mm² a⁻³ (Davis et al., 2012), with the observational noise given by the RMSE of each residual time series. Initial state values of the filter were set to zero for both the level and trend with large initial uncertainties (1×10⁶). The filter parameters were then optimized using the expectation__maximization (EM) algorithm (Shumway and Stoffer, 1982) with five iterations. The same approach as in Method 1 was used to generate the two estimates of the offset based on the 1-year overlap and the differencing of the two intercepts. For method 3, here the offsets were determined by computing the median difference between the two missions over the 2010–2011 time period. To determine which of the offsets produces the best cross-calibration, we apply each offset and compute linear rates of change from 2003 to 2019. These rates are then compared to rates estimated from unbiased ICESat and ICESat-2 measurements produced by Smith et al. (2020), and the offset with the smallest absolute difference was selected. Finally, the selected offset rate difference (radar minus laser) is checked against the difference computed without a residual cross-calibration. If the applied offset did not improve the rate compared to the ICESat and ICESat-2 record, then the residual offset was set to zero. Following Schröder et al. (2019), we remove outliers in the offsets using a 100×100 km 5-MAD moving spatial filter. The intermission offsets are then interpolated using a Gaussian kernel with a 20 km correlation length using the nine closest data points. This produces a spatially consistent field of offsets for the cross-calibration of the elevation residuals. Then the offsets estimated from the initial least-squares adjustment and the offsets estimated from the secondary residual calibration are then applied to create a fully calibrated local time series. Finally, the individually calibrated elevation time series for each mission or mode are averaged to monthly estimates of elevation change for each spatial grid cell, with an associated standard error. Once, the time series have been calibrated a seasonal amplitude correction is applied to the data to normalize amplitudes between missions. This is described in more detail in Sect. 3.2.4. Finally, the monthly normalized time series are then combined and integrated into a continuous record using the weighted average of the data within each overlapping temporal bin. Weights are specified as the inverse variance of each mission's accuracy and the random error estimated from the monthly averaging procedure (see Sect. 4.1).

Figure 4Maps of the residual cross-calibration offset and the corresponding error for the three main intermission transition periods. One should note that here ICESat has been grouped with Envisat in the initial calibration.

The initial least-squares adjustment provided good alignment between overlapping modes (ocean–ice mode) and missions (Envisat and ICESat), as well as a first-order correction for the three weakly overlapping missions that allows for better estimation of the residual biases from the detrended data. Initial offsets were determined to be as large as 10–15 m in areas of rapid change such as Pine Island Glacier. However, the least-squares adjustment was shown to be inadequate when large non-linear elevation changes are present. The magnitude of the estimated residual cross-calibration error (after least-squares adjustment) (Fig. 4) shows that most overlapping regions have a clear correlation with temporal coincident elevation change rates. This pattern is evident in the Envisat to CryoSat-2 transition (Fig. 4) for Dronning Maud Land (basins 5–8), Wilkes Land (basins 12–13), Bellingshausen Sea (basins 23–25), and the Amundsen Sea sector (basins 20–23) (Fig. 10: 2010–2012). For the ERS-2 to Envisat transition, we find a clear correlation between the magnitude of the offsets and the changes in elevation due to variations in surface mass balance in Wilkes Land (basins 12–13 seen in Fig. 1) over the 2001–2003 time period (Schröder et al., 2019).

3.2.4 Normalization of seasonal amplitude

The radar signal's interaction with the surface and sub-surface firn layers can create artificially large seasonal amplitudes and trends, as described in Sect. 3.2.2. We correct for these as best as possible using information contained in the waveform parameters. However, in many cases these corrections are unable to fully correct the artificial signals. This behavior can be seen in Schröder et al. (2019) and in our data, and even after the scattering correction has been applied there exist intermission variations in seasonal amplitude (Fig. 5). To further reduce this effect, we apply an amplitude correction (h_n) to each mission to normalize the seasonal signal over the entire record. We normalized the seasonal amplitudes of the ERS-1 and ERS-2 and Envisat records to match amplitudes computed from CryoSat-2. CryoSat-2, which is retracked with a much lower threshold of the maximum waveform amplitude (10 %) for LRM and a maximum gradient threshold for SARIn, has been shown to be less sensitive to changes in surface properties and produces seasonal amplitudes of the same magnitude as ICESat (Fig. 5) (Nilsson et al., 2016). After removal of the long-term trend, the amplitude normalization was computed for each mission, except for ICESat and CryoSat-2, according to

$$\begin{array}{}\text{(3)}& h_n=\left(\mathrm{1}-{\displaystyle \frac{a_i}{a_{\mathrm{r}}}}\right)\cdot \left[{\mathit{\alpha }}_{\mathrm{c}}\cos \left(\mathrm{2}\mathit{\pi }t\right)+{\mathit{\alpha }}_{\mathrm{s}}\sin \left(\mathrm{2}\mathit{\pi }t\right)\right],\end{array}$$

where (a_i) is the amplitude of the mission ($a^{\mathrm{2}}={\mathit{\alpha }}_{\mathrm{c}}^{\mathrm{2}}+{\mathit{\alpha }}_{\mathrm{s}}^{\mathrm{2}}$), (a_r) is the reference amplitude estimated from CryoSat-2 data, and ${\mathit{\alpha }}_{\mathrm{c},\mathrm{s}}^{\mathrm{2}}$ are the coefficients for the seasonal model. The correction is applied by subtracting it from each individual time series, and the normalization has the effect of producing more homogeneous amplitudes over the entire altimetry record. The application did not introduce any noticeable shift in the phase of the seasonal signal.

ICESat and CryoSat-2 LRM mode shows a similar magnitude in amplitude and supports the choice of using CryoSat-2 as reference in which the difference is most likely explained by the lower temporal sampling of ICESat. The slightly lower seasonal amplitude of ICESat-2 is mostly likely due to the short time span used to estimate the amplitude (2 years), as seen in Fig. 5.

Figure 5Median seasonal amplitude of the different missions and modes for the CryoSat-2 LRM (a) and SARIn (b) mode masks (south of 81.5^∘ S for LRM). The blue bars show the original seasonal amplitude with no corrections applied, the orange bars show the amplitude once the mission-dependent scattering correction has been applied, and the green bars show the normalized amplitude after adjustment using CryoSat-2 as reference.

Download

4.1 Uncertainties of elevation change time series and data

An internal crossover analysis was performed to determine the relative accuracy of each mission and mode in a similar manner as Brenner et al. (2007) and Schröder et al. (2019). We estimated the standard deviation of all crossovers with a time difference of less than 31 d. Crossovers were binned as a function of surface slope at intervals of 0.04^∘ (Fig. 6). The relative accuracy of each mission or mode was determined from the standard deviation at zero slope by fitting an error function (inside an interval of 0–0.4^∘) as shown in Table 1. To derive the uncertainty of each time series epoch we use the spatiotemporal variability inside each monthly time interval in the form of the standard deviation. This provides a random error for each monthly value that varies both in space and in time and encompasses measurement-related errors driven by topography, retracking, range corrections, etc. To quantify the total cross-calibration error for each time series we use the standard deviation of each grouped mission offsets (Sect. 3.2.3) and add them in quadrature to estimate the total cross-calibration error, similar to Schröder et al. (2019). We then have the total error (σ_m) for each month in each time series by summing the individual error sources as follows:

where (σ_m) the error due to the elevation change variability within each monthly interval for each time series, and (σ_c) is the total cross-calibration error for each time series. The estimated total error (${\mathit{\sigma }}_{\mathrm{m}}^{\mathrm{2}}$) is the provided RMS error in the product (varying both by location and time).

Figure 6Standard deviation (cm) of intra-mission and intra-mode crossovers for the Antarctic Ice Sheet as a function of surface slope (degrees). Precision decreases quasi-linearly as surface slope increases.

Download

Table 1Sensor and mode errors (σ_mission) as a function of the random (σ_noise) and slope-dependent (σ_slope) errors. Slope (α) is in degrees. Modeled error (σ_mission) is based on fitting the following function to the intra-sensor, intra-mode crossover data: ${\mathit{\sigma }}_{\mathrm{mission}}={\mathit{\sigma }}_{\mathrm{noise}}+{\mathit{\sigma }}_{\mathrm{slope}}{\mathit{\alpha }}^{\mathrm{2}}$.

Download Print Version | Download XLSX

4.2 Validation of rates of elevation change

To validate the data product, we computed elevation change rates and compared them to rates derived from near-coincident Operation IceBridge (OIB; MacGregor et al., 2021) and pre-OIB data spanning the period 2002 to 2019 using the Airborne Topographic Mapper (ATM; MacGregor et al., 2021) laser altimeter. Elevation change rates for ATM were derived following the approach of Nilsson et al. (2016), in which a linear model was solved at each measurement location using a search radius of 175 m. Following the approach of McMillan et al. (2014) and Wouters et al. (2015), the local slope was used to correct the measurements to the reference track, indicated as Track_Identifier = 0 in the product. Solutions were rejected if they contained fewer than two campaigns of ATM data, the magnitude of linear rate was larger than 10 m a⁻¹, the standard deviation of the solution exceeded 1 m a⁻¹, or if the solution contained fewer than 10 measurements, as well as if the time span was less than 2 years. The elevation accuracy of the ATM sensor family has an estimated error of less than 9 cm (Brunt et al., 2017), corresponding to an accuracy of roughly 0.5 cm a⁻¹ over the 18-year measurement period. Operation IceBridge coverage is concentrated to the western parts of the Antarctic Ice Sheet, providing very limited coverage in the east. To overcome this limitation, we also use elevation change rates estimated by Smith et al. (2020) that are based on crossover analyses of satellite laser altimetry (ICESat and ICESat-2: 2003–2019) that has an error of roughly 10 cm. This corresponds to an error in the rate of elevation change of about 0.6 cm a⁻¹, which is consistent with the error observed for ATM. These errors and their impact are discussed further in Sect. 5.

4.3 Area integrated error estimation

Area integrated error for each drainage region, based on the outlines from Zwally et al. (2012) (shown in Fig. 1), are estimated loosely following the approach of Nilsson et al. (2016). The total area integrated error is divided into three main components: the systematic bias, the random error, and the rate error estimated in the fitting procedure. These are then combined in quadrature to produce the total error according to

where σ_s is the systematic bias, σ_r the random error, ${\mathit{\sigma }}_{\dot{h}}$ the rate error, n the number of uncorrelated elevation change estimates (see below), and k the degrees of freedom in the least-squares model (k=2). The systematic bias and the random error are taken as the average and standard deviation of the difference in rate between the Jet Propulsion Laboratory (JPL; this study) and ICESat and ICESat-2 (Smith et al., 2020) products for the 2003–2019 period. We compute the error in the estimated rate using the variance–covariance matrix in the least square fitting procedure according to

where ${\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{m}}$ is the average monthly uncertainty from our product inside the time interval of interest, X is the design matrix of the linear model, X^T is the transpose of the design matrix, “diag” is the diagonal elements of the array, and ⁻¹ is the inverse of the dot products. The subscript $\dot{h}$ is the location of the rate error in the diagonal array. To account for spatial autocorrelation σ_r and ${\mathit{\sigma }}_{\dot{h}}$ are divided by n; n is estimated by dividing the total area of each drainage region with the correlation area: $n=A/\mathit{\pi }{\mathit{\rho }}^{\mathrm{2}}$, where A is the area of the region, and ρ is the correlation length. The errors for each drainage region are summarized in Table 2. The intrinsic quality of each mission was determined through internal crossover analysis (Sect. 4.1) of each mode and mission and is summarized in Table 1 and Fig. 6. Analyzing the correlation length of the laser-only versus JPL elevation change differences we find an ice-sheet-wide correlation length scale on the order of 20–100 km. To be conservative, a correlation length of 100 km was used to compute n.

Table 2Regionally averaged errors for the synthesized JPL record of elevation change, computed relative to the unbiased ICESat to ICESat-2 estimate of Smith et al. (2020). Errors were determined by differencing 2003–2019 linear rates of elevation change between products. The bias (mean: σ_s) and error (standard deviation: σ_r) are computed for each drainage basin (1–27; Fig. 1). Antarctic Ice Sheet (AIS), Antarctic Peninsula (AP), West Antarctic Ice Sheet (WAIS), and East Antarctic Ice Sheet (EAIS) statistics are determined using area weighted averages.

Download Print Version | Download XLSX

5.1 Accuracy of synthesis

Previous studies have relied on near-coincident airborne measurements to validate land ice elevation changes derived from multi-mission synthesis (McMillan et al., 2014; Nilsson et al., 2016; Simonsen and Sørensen, 2017; Wouters et al., 2015). This approach, however, is limited in both the spatial and temporal coverage. For Antarctica, airborne validation data have been collected during austral summer, mostly over rapidly thinning glaciers, such as Pine Island and Thwaites, in the western part of the ice sheet, with significant spatial coverage starting in 2002. The derived errors from these local comparisons are then extrapolated to the entire ice sheet into regions exhibiting very different surface and metrological conditions. With the launch of ICESat-2 in September 2018 we now have, for the first time, the ability to compare long-term unbiased laser-derived rates of elevation change on a continental scale. For this analysis we compare our synthesized rates of elevation change to those estimated by Smith et al. (2020) for the period 2003–2019 for each basin (Zwally et al., 2012) (Fig. 1). The results of this analysis are summarized in Table 2. We find an ice-sheet-wide error of $-\mathrm{0.8}\pm \mathrm{7.8}$ mm a⁻¹ (Fig. 7e) with a quadratic and linear increase as a function of surface slope in the systematic bias and random error, respectively (Fig. 7f and g). To determine the validity of this comparison we also compared ICESat and ICESat-2 rates with rates from ATM over the time period 2003–2018. Good agreement was found between the two datasets with an average difference of 2.3±22 cm a⁻¹ (Fig. 7c) over regions with an observed rate of elevation change from ATM ranging from −15 to 2 m a⁻¹. The main discrepancies between our product and the ICESat- and ICESat-2-derived elevation change are concentrated over areas of high relief and over regions with large magnitude changes, such as Pine Island and Thwaites glaciers (Fig. 7a). Here, differences larger than 10 cm a⁻¹ can be found, and for the main trunk of Pine Island Glacier we find a difference of 2±10 cm a⁻¹ (Fig. 7a). The magnitude of the ATM error compared to the ICESat and ICESat-2 product is larger. This is mostly due to the fact that the data comparison locations are in areas of rapid change. The correlation between the two laser datasets and our product is greater than 0.8 in both cases.

Figure 7Elevation change validation and comparison using rates derived from ICESat and ICESat-2 and airborne ATM data over the time period of 2003–2019 and 2001–2019, respectively. Panel (a) shows the spatial distribution of the elevation change differences from this study (JPL) differenced with rates derived from Smith et al. (2020). Panel (b) shows the comparison of rates derived from JPL with ATM at locations indicated in (d) with green flight lines. Panel (c) shows the comparison between ICESat- and ICESat-2-derived rates with ATM. Panel (e) depicts the ice-sheet-wide histogram of elevation change differences. Panels (b), (c), and (e) include the distribution mean (μ), standard deviation (σ), correlation (ρ), and number of observations (n). Panels (f) and (g) are the bias (mean) and error (standard deviation) as a function of surface slope for the JPL and ICESAT-1 and ICESat-2 validation.

The relative precision of the different satellite altimeters used in this study range from 5 to 40 cm over low slope surfaces (Table 1 and Fig. 6). Earlier missions such as Geosat, ERS-1, and ERS-2 are roughly 3 times less precise than later missions (Envisat, ICESat, ICESat-2, and CryoSat-2). However, it was also found that the ERS-1 and ERS-2 Ocean mode was ∼ 30 % more precise than ice mode data, bringing it closely in line with the later missions. Unfortunately, the data coverage of the ocean mode is far lower than the ice mode. For CryoSat-2, the lower relative precision of the SARIn mode can be attributed to the spatial coverage, with SARIn operating over rougher terrain compared to the LRM mode that operates over the interior of the ice sheet with a higher along-track resolution (i.e., smaller footprint). Similar effects were also seen in Schröder et al. (2019). The laser altimetry missions show the lowest noise levels, on the order of 5 cm over flat areas ranging up to 20 cm for slopes <0.8^∘, with ICESat-2 showing a factor-of-2 improvement in precision over its predecessor (ICESat) over all surface slopes.

5.2 Comparison to other studies and datasets

Previous long-term Antarctic Ice Sheet elevation change products have been produced by Dresden University of Technology (TUD; Schröder et al., 2019) and the Centre for Polar Observation and Modelling (CPOM; Shepherd et al., 2019). These products vary in both resolution and processing methodologies. The TUD product is provided at a spatial resolution of 10 km and as monthly elevation change estimates. In contrast, the CPOM product provides elevation change estimates every 5 years at 5 km resolution and basin-wide time series of mass change at quarterly resolution. The TUD dataset is comprised of Seasat, Geosat, ERS-1, ERS-2, Envisat, and CryoSat-2, while CPOM consists of data from ERS-1, ERS-2, Envisat, and CryoSat-2. To allow for a fair comparison between the different products we used our provided product without hypsometric extrapolation for the analysis.

The errors reported for our elevation change synthesis are slightly larger than those reported by TUD; this is due to the difference in retracking and the fitting procedure used to derive the error estimates. Comparing all three data products to the ATM validation data we find the best agreement with the JPL synthesis. (JPL 4±19 cm a⁻¹, TUD 6±20, and CPOM $+\mathrm{4}\pm \mathrm{53}$ to $-\mathrm{16}\pm \mathrm{61}$ cm a⁻¹). The JPL and TUD estimates were computed from the same ATM dataset and given the same editing criteria, while values from CPOM are the reported values from Shepherd et al. (2019). Applying the same analysis to the 2007–2011 and 2011–2016 elevation change solutions provided by CPOM, we found values of 29±41 cm a⁻¹ (2007–2011) and $-\mathrm{8}\pm \mathrm{30}$ cm a⁻¹ (2011–2016) for the comparison with ATM, as well as a weighted average of $-\mathrm{2.2}\pm \mathrm{33}$ cm a⁻¹ comparing data from overlapping locations. To further compare the noise level in the different datasets we use the elevation change from the common 1992–2016 time period (as CPOM only provides rates in 5-year intervals) of all products and compare against ICESat and ICESat-2 elevation change rate from 2003 to 2019. To reduce the impact of differences in time span, we initially compare only data between 81.5 and 90^∘ S (pole hole) as this spatial domain only contains ICESat and CryoSat-2 measurements and is thus the most closely aligned in time with the ICESat and ICESat-2 estimate. We also perform an ice-sheet-wide analysis, though the time spans are not identical. To compute the noise level, we simply differenced the three rate fields with the ICESat- and ICESat-2-derived rates and computed the average and standard deviation of the differences. This provided the following ice-sheet-wide results: $-\mathrm{0.32}\pm \mathrm{1.70}$ (JPL), $-\mathrm{0.45}\pm \mathrm{1.92}$ (TUD), and $-\mathrm{0.33}\pm \mathrm{2.59}$ cm a⁻¹ (CPOM). For the pole-hole region, 81.5 and 86^∘ S, the following results were obtained: $-\mathrm{0.33}\pm \mathrm{1.17}$ (JPL), $-\mathrm{1.37}\pm \mathrm{1.57}$ (TUD), and $-\mathrm{1.90}\pm \mathrm{3.15}$ cm a⁻¹ (CPOM).

Comparing the long-term rates for the overlapping time period 1992–2016, we find an overall good agreement for the three original products. Comparing only values north of 81.5^∘ S, we determine volume change rates of −58, −48, and −59 km³ a⁻¹ for JPL, TUD, and CPOM, respectively. Differences are well within the errors for all the three products. Studying the differences in spatial patterns (Fig. 8), using the JPL-derived rate as the reference, we find that the TUD and JPL products agree well over East Antarctica in basins 10–17, while a larger difference can be seen in basin 3 closer to the Weddell Sea. Larger differences between JPL and CPOM compared to JPL versus TUD can be observed in East Antarctica (East Antarctic Ice Sheet: EAIS). This is likely a result of different methodologies for correcting changes in the radar scattering. Dividing the estimates into different regions we find the following volume change estimates for the 1992–2017 period: WAIS (West Antarctic Ice Sheet; JPL −108, TUD −100, and CPOM −106 km³ a⁻¹), EAIS (JPL 61, TUD 48, and CPOM 43 km³ a⁻¹), and AP (JPL −11, TUD 4, and CPOM 5 km³ a⁻¹). The regional estimates agree well among products, with the largest discrepancies found in the Antarctic Peninsula. Here, both the TUD and CPOM products provide a positive volume change compared to the JPL product, highlighting the challenge in obtaining accurate estimates from this region. Comparing the JPL and TUD products with rates from Smith et al. (2020) (ICESat and ICESat-2) over the time period 2003–2017 (again using the original JPL product with no hypsometric extrapolation), we find that the two products agree well over WAIS (JPL −165, TUD −164, and LA −200 km³ a⁻¹) but lower in magnitude compared to ICESat and ICESat-2 due to the larger radar footprint. For EAIS (JPL 83, TUD 51, and LA 85 km³ a⁻¹) a disagreement of roughly 40 % is observed between the TUD and JPL products, in which LA and JPL values are practically identical. In the AP (JPL −19, TUD −7, and LA −39 km³ a⁻¹) both products are lower in magnitude compared to ICESat and ICESat-2, on the order of 50 %–80 % due to limitations in measuring over high relief topography.

Figure 8Comparison of overlapping long-term rates from the Technical University of Dresden (TUD) and Center for Polar Observation and Modelling (CPOM) altimetry product with rates from this study (JPL).

To understand how well these products can capture (and provide insight into) the change or variability in physical processes of the ice sheets, we compared our result with modeled changes in surface elevations (“zs”) from the IMAU (Institute for Marine and Atmospheric research Utrecht) firn densification model (FDM; Ligtenberg et al., 2012) that is forced by 6 h mass balance components (snowfall, rain, sublimation, and snowmelt), average surface temperature, and 10 m wind speed from the Regional Atmospheric Climate Model version 2.3p2 (van Wessem et al., 2018). The firn model only simulates changes in surface elevation due to changes in surface processes and does not account for thinning or thickening resulting from changes in ice dynamics (flow). To minimize dynamic signals, we mask areas with surface velocities larger than 30 m a⁻¹ using the velocity field provided by the ITS_LIVE project (Gardner et al., 2019) merged with phase-based estimates (Mouginot et al., 2019). The surface elevation long-term trend and acceleration fields (1992–2016), seen in Fig. 9, show that for Dronning Maud Land and Enderby Land (basins 4–11) there is generally good agreement in both the spatial pattern and the sign of the observed and modeled rate of elevation change. For these regions, the observed change can be attributed to an increase in accumulation (Boening et al., 2012). However, the magnitude between the modeled and measured rates of change differs by roughly 50 %. The altimetry-derived volume change for basins 4–11, over the time period 1992–2016, is estimated at 46 km³ a⁻¹ compared to a modeled change of 27 km³ a⁻¹. This disagreement becomes even more prominent for Wilkes Land (basins 12–14) where the difference between modeled and observed rates of change are larger and of opposite signs (Fig. 9). For these three basins, the estimated difference in volume change is on the order of 36 km³ a⁻¹ based on the difference in the modeled change of −25 km³ a⁻¹ compared to 11 km³ a⁻¹ from altimetry. The magnitude and sign of these results are consistent within all three altimetry products compared to the FDM. Further, comparing the differences in the magnitude of the seasonal amplitude for 1992–2016, we find that the TUD product has an annual amplitude that is ∼50 % larger than the JPL product (5.1±15 versus 2.7±4.9 cm). Our estimated value of 2.7±4.9 cm compares well with the 2.9±4.1 cm average FDM amplitude for the period 1992–2016. This analysis was not applied to the CPOM product as their provided basin time series are in units of mass after a firn correction has been applied.

Figure 9Spatial fields of rates (a, d), acceleration (b, e) and seasonal amplitudes (c, f) from our product (JPL; a–c) and modeled values from the IMAU firn densification model (FDM; d–f). Areas of fast flow (>30 m a⁻¹) have been masked out to minimize height changes caused by changes in ice flow. The altimetry data has been smoothed with a 50 km median filter to highlight large-scale spatial patterns.

5.3 Basin-scale time-evolving volume change

Analyzing the 1992–2020 record of surface elevation (Table 3 and Figs. 10 and 11), including the area between 81.5 and 90^∘ S, we determine an average rate of volume change of $-\mathrm{68}\pm \mathrm{11}$ km³ a⁻¹ over the entire ice sheet, with large losses from the West Antarctic Ice Sheet ($-\mathrm{113}\pm \mathrm{6}$ km³ a⁻¹) and gains for East Antarctic Ice Sheet ($+\mathrm{75}\pm \mathrm{5}$ km³ a⁻¹) which experienced large snowfall events in 2009 and 2011 (Boening et al., 2012). The Antarctic Peninsula (AP) is the most challenging region to measure elevation change due to its extreme surface relief and sparse data coverage. We anticipate that any estimate derived from conventional satellite radar altimetry will be biased positive due to the inability to measure low-elevation signals. That said, we estimate an overall negative trend for the AP of $-\mathrm{27}\pm \mathrm{8}$ km³ a⁻¹ for the 29-year record (1992–2020) (Figs. 10 and 11), which aligns closely with other estimates (Groh el al., 2021; Rignot et al., 2019; The IMBIE team, 2018; Zwally et al., 2021) but is highly dependent on the applied hypsometric extrapolation (Sect. 3.2.5). On decadal timescales we find that the large glaciers systems of Pine Island, Thwaites, Smith, and Kohler (basins 21 and 22) have shown relatively stable mass loss since the early parts of the satellite era, with signs of accelerated thinning since 2007–2009 (Fig. 11). WAIS has seen almost a doubling of its mass loss in the last decade (2011–2020) compared to the two previous decades (Fig. 11). EAIS has reverted back to its previous long-term decadal rate of ∼ +8 km³ a⁻¹, in line with the observed 5-year trend from Geosat over Dronning Maud Land (Figs. 11 and 12), down from +84 km³ a⁻¹ following the anomalous snowfall during the 2001–2011 period. AP was in balance and saw little observable change in the first decade (1991–2010) but increased its mass loss by a factor of 10 in the period of 2001–2011. The mass loss in the last decade was slowed by roughly 50 % due to a positive mass balance anomaly during the period 2016–2018. Over the Geosat time period from 1985 to 1989, and for latitudes <72^∘ S, a general stable and small positive rate of 6±16 km³ a⁻¹ was found for the EA1 region (basins 4–11; Fig. 12). This rate remained stable between 1985 and 2009 (∼ 10 km³ a⁻¹) until the onset of a precipitation event in 2009. For the EA2 region (basins 12–15; Fig. 1) a shift in both sign and magnitude was observed for the 1985–1989 period compared to the long-term positive rate for EA1. The mass loss over the 1985–1989 period was $-\mathrm{70}\pm \mathrm{22}$ km³ a⁻¹ and was found to be mostly driven by the Totten Glacier system in basin 13 (Fig. 10). This rate is based, however, on heavy extrapolation over the Totten region due to poor data coverage for the last 2 years of the mission and should be treated with caution. Trends for EA2 showed a stable negative rate (∼ 25–30 km³ a⁻¹) until 2001–2003 when a large positive change occurred due to an increase in SMB (surface mass balance; Fig. 10c). The region reverted back to a long-term negative trend after 2006 mostly modulated by changes in SMB.

Table 3Volume change rates spanning 1985 to 2020 for basins 1–27 (Fig. 1) and aggregate regions. Volume change errors are computed from the ICESat and ICESat-2 validation procedure, combined with the error in the estimated rate.

NA: not available

Download Print Version | Download XLSX

Figure 10Rates of Antarctic Ice Sheet elevation change. Elevation change rate (a) and acceleration (b) for the 1992–2020 period and average rates for (c) 3-year, (d) ICESat and ICESat-2, and (d) 10-year intervals. The * indicates a 5-year interval for Geosat.

Regionally, concentrated rates of thinning from accelerated glacier flow (Gardner et al., 2018; Rignot et al., 2019) are found to spread inland over time due to a regional dynamic imbalance (Shepherd et al., 2019). The marginal areas surrounding the Getz Ice Shelf (basin 20) also exhibit negative rates of elevation change but are more localized to the narrow glacier outlets due to inland topographic barriers and time since initiation of thinning (Figs. 10 and 11). This area saw a large break in the overall long-term trend around 2010 when rapid onset thinning was observed and attributed to short-term variations in both the surface mass balance and ice dynamics (Chuter et al., 2017; Schröder et al., 2019; Gardner et al., 2018). Basin 18, which contains the Kamb Ice Stream, experienced a relatively steady gain in volume over the last three decades resulting from the stagnation of the Kamb Ice Stream some 200 years prior (Catania et al., 2006) (Figs. 10 and 11). Totten Glacier (basin 13), part of the EAIS, has been losing mass since the late 1970s (Schröder et al., 2019) with the average trend mostly governed by ice dynamics and short-term variability and acceleration driven by changes in precipitation (Li et al., 2016). A major change in trend was observed in 2010 when a large-scale thinning of the entire glacier is observed, likely in response to a change in precipitation and possibly changes in ice dynamics driven by changes in ocean conditions (Khazendar et al., 2013; Li et al., 2016). The activation or reversal in trend of both the Totten and Denman glaciers in early 2009–2010 has disrupted the long-term equilibrium or gain that has been observed for most parts of Wilkes Land (basins 12 and 13; Fig. 1). A departure from the long-term trend can now be observed for large parts of Wilkes Land in the form of large-scale negative acceleration spreading inland (Fig. 10). In Dronning Maud Land and Enderby Land (basins 5–8), the previously mentioned snowfall events in 2009 and 2011 (Boening et al., 2012) are clearly observed in the regional elevation change trends. This pattern is most prominent along the Weddell Sea coast where the accumulation signal, in the form of precipitation, shows an earlier event in 2006 (basins 3 and 4) (Figs. 10 and 11). The glaciers flowing into the Bellingshausen Sea have shown a complex pattern of change over the last 29 years. Here, Palmer Land (basin 24) shows a steady increase in surface elevation over the initial 15 years of the record, following a long-term positive anomaly in precipitation from 1992. However, a reversal in this pattern was observed around 2007 when patterns of thinning (McMillan et al., 2014; Schröder et al., 2019; Shepherd et al., 2019; Wouters et al., 2015) (Fig. 10) can be observed localized to the major low-elevation outlet glaciers in the regions. The change can be largely attributed to a change in precipitation amount, with lesser contributions from changes in ice dynamics resulting from enhanced melting by the ocean (Gardner et al., 2018; Hogg et al., 2017). However, in the southern part of the Bellingshausen Sea, near Ferrigno Glacier in basin 23, we find a relatively stable trend during most of the record until 2009 when a large acceleration in ice loss can be observed. This acceleration can only be partially attributed to changes in ice dynamics (Gardner et al., 2018; Wouters et al., 2015), and it is likely that changes in precipitation are the major driver of change. Large changes in both spatial and temporal variability can be observed in the AP region in the last three decades, when large-scale reversals of signals can be observed over different time periods. Here, we find a large-scale positive elevation change anomaly in basins 23–26, superimposed on a long-term negative trend, over the time periods 1998–2000, 2004–2006, and 2016–2018. These changes are linked to changes in the short-term variability in SMB in the region due to increased precipitation. Examining the rates derived over the ICESat-2 time period (2018–2020) a large positive elevation change signal can be observed over the WAIS region, in contrast to the overall negative long-term trend. This anomaly is directly linked to large-scale snow accumulation resulting from an extreme precipitation event in the austral winter of 2019 which has been attributed to the landfall of atmospheric rivers (Adusumilli et al., 2021).

Figure 11Basin (Zwally et al., 2012) and ice sheet monthly elevation change time series for the period of 1992 to 2020.