Elevation Changes of the Antarctic Ice Sheet from Joint Envisat and CryoSat-2 Radar Altimetry

: The elevation changes of ice sheets have been recognized as an essential climate variable. Long-term time series of these changes are an important parameter to understand climate change, and the longest time-series of ice sheet elevation changes can be derived from combining multiple Ku-band satellite altimetry missions. However, unresolved intermission biases obscure the record. Here, we revise the mathematical model commonly used in the literature to simultaneously correct for intermission bias and ascending–descending bias to ensure the self-consistency and cohesion of the elevation time series across missions. This updated approach is applied to combine Envisat and CryoSat-2 radar altimetry in the period of 2002–2019. We tested this approach by validating it against airborne and satellite laser altimetry. Combining the detailed temporal and spatial evolution of elevation changes with ﬁrn densiﬁcation-modeled volume changes due to surface processes, we found that the Amundsen Sea sector accounts for most of the total volume loss of the Antarctic Ice Sheet (AIS), mainly from ice dynamics. However, surface processes dominate the volume changes in the key regions, such as the Totten Glacier sector, Dronning Maud Land, Princess Elizabeth Land, and the Bellingshausen Sea sector. Overall, accelerated volume loss in the West Antarctic continues to outpace the gains observed in the East Antarctic. The total volume change during 2002–2019 for the AIS was − 68.7 ± 8.1 km 3 / y, with an acceleration of − 5.5 ± 0.9 km 3 / y 2 .


Introduction
The elevation changes of ice sheets are the result of changes in ice sheet ice dynamics and weather-driven changes on the surface. The long-term evolution of ice sheet volume changes is an essential climate variable defined by the Global Climate Observing System (GCOS) to assess the impact of climate change and its contribution to sea level rise. Satellite altimetry can accurately measure the elevation of ice sheets and thereby provide valuable data for monitoring ice sheet changes at both local and continental scales. Since Zwally et al. (1983) [1] first used satellite altimetry measurements in ice sheet surveys, altimetry observations have been essential for studying ice sheet evolution [2][3][4][5][6][7][8]. However, limited by the effective lifetime of the satellite, a single mission cannot characterize the long-term evolution of the ice sheet [9,10]. Consequently, combining observations from multiple argued that such a threshold retracker is more precise in terms of ice sheet elevation measurements. The erroneous height records are culled using standard quality flags.
CryoSat-2 was launched in April 2010 using a 369-day full repeat orbit with a subcycle of 30 days and is still in operation. It operates at an altitude of 717 km with an inclination of 92 • and provides coverage up to 88 • S. CryoSat-2 measures the interior parts of the AIS in the Low-Resolution Mode (LRM), while using the novel Synthetic Aperture Radar Interferometric (SARIn) mode around the AIS margins. The LRM is equivalent to conventional Ku-band radar altimeters. SARIn is very suitable for measuring steep and complex topography across ice sheet margins because of the use of dual antennae operating in an interferometric mode [29]. The retracked LRM observations with the OCOG retracker from July 2010 to April 2019 were used in the following study to match the Envisat ICE-1 retracker elevations. As for the SARIn observations, the only retracked observations by the Wingham/Wallis model fit retracker provided by the product were used. The erroneous height records and SARIn records where the interferometric across-track location failed were removed. Therefore, the period from July 2010 to April 2012 could be used to calibrate the intermission biases.

Laser Altimeter Elevation Datasets from the Operation IceBridge (OIB) Project and ICESat
In this study, based on previous studies, we presented an updated regression model to correct the intermission bias and the A-D bias simultaneously. To evaluate the performance of the model, we needed to validate its results. Here, we used laser altimeter elevation datasets obtained from NASA's OIB project and ICESat mission.
The OIB project operated since 2009 to bridge the gap between the ICESat and ICESat-2 missions and maintain a continuous observational record over the Arctic and the Antarctic. The Airborne Topographic Mapper (ATM) is an important airborne laser altimeter payload for OIB projects that can provide surface elevation measurements with an accuracy of 10 cm or better [30]. In addition, some pre-IceBridge ATM elevation observations were also collected during 1993-2008.
From 2003 to 2009, ICESat was the first satellite laser altimeter mission, and, by operating in the near inferred spectrum, it provided precise observations of the surface elevation. Here, the latest release of ICESat data (GLA12, Version 34) was used to derive surface elevation change estimates in a 20 km grid via a repeat-track analysis [31]. Intercampaign bias correction and saturation correction were also applied. The final elevation change rate for each grid was generated from the median of surface height change rates to minimize the impact of the remaining outliers, thereby providing a good dataset for benchmarking our elevation change estimates from Envisat.

The Updated Least-Squares Regression Model
Previous studies have demonstrated that elevation time series' can be obtained by performing a least-squares regression on all elevation measurements within a grid cell, also known as "plane-fitting" [3, 10,17,23]. Following these studies, the regression parameters can be expressed as h(lon i , lat i , t i ) = h 0 + a 0 (lon i − lon 0 ) + a 1 (lat i − lat 0 ) + a 2 (lon i − lon 0 ) 2 + a 3 (lat i − lat 0 ) 2 + a 4 (lon i − lon 0 )(lat i − lat 0 ) + a 5 (t i − t 0 ) + a 6 cos(2π(t i − t 0 )) + a 7 sin(2π(t i − t 0 )) + a 8 bs i − bs + b AD (−1) AD + b m (−1) m + res(lon i , lat i , t i ), (1) where h(lon i , lat i , t i ) is the elevation measured by the satellite altimetry missions at longitude (lon i ), latitude (lat i ), and time (t i ); h 0 is the mean altitude of the bin with the center coordinates (lon 0 , lat 0 ) in which the least-squares regression is performed; t 0 is the reference epoch set to May 2010; a 0~a4 are the coefficients of a biquadratic surface model accounting for topography; a 5 is the elevation change rate Remote Sens. 2020, 12, 3746 4 of 20 constant; a 6~a7 describe the seasonal variations in elevation; a 8 is for the anomaly of backscattered power (bs i − bs) to account for variations in the penetration depth of the radar signal; b AD is for the intramission A-D bias, and AD is assigned a value of 1 for ascending tracks or 0 for descending tracks; b m is only used in CryoSat-2 data processing to account for the possible bias caused by the different retracking algorithms used for LRM and SARIn, and m is 1 for LRM or 0 for SARIn; and res(lon i , lat i , t i ) denotes the residuals.
However, the regression model above is designed for CryoSat-2. It is not suitable for multimission processing, unless observations from multiple radar altimeter missions are also cross-calibrated. Cross-calibration is used to correct the bias between different satellite radar altimeter missions. As mentioned in Section 1, the bias between two radar altimeter missions, including the intermission bias and the intermission A-D bias, can be regarded as invariants in the local area [10,14,15,18] and can be modeled like A-D in Equation (1). Hence, the intermission extension of Equation (1) is given by where b E a E d and b C a C d are for the intramission A-D bias of Envisat or CryoSat-2, and b C a E a , b C d E a , b C a E d and b C d E d are for the total static systematic biases of the intermission bias and the intermission A-D bias between Envisat and CryoSat-2. In this study, we referred to these total static systematic biases as the total intermission biases. The internal consistency between these biases (assuming they are constant over time) would give additional constraints, which were handled by the last three equations in Equation (2). Notably, here we discarded b m (−1) m , as this bias could be modeled as a part of the intramission A-D bias.

Calculation of the Combined Elevation Time Series
We first performed the least-squares regression, as shown in Equation (2), at the intersections of the ascending and descending satellite ground tracks of Envisat to estimate all the parameters. For each of the intersections, an independent least-squares regression was estimated based on all Envisat and Cryosat-2 observations within 2.5 km of the center of an intersection. The first-round estimated parameters require that the number of observations of each satellite in an ascending or descending orbit shall not be less than 50. Then, an iterative procedure continues until parameters are estimated that pass the 3σ outlier rejection criteria, as in [6], or until the remaining observations in the bin cannot meet the thresholds. The thresholds require that the total observations be more than 100 and that the number of ascending or descending observations of each satellite be no less than 20. The estimations of the systematic biases Then, we computed the least-squares regression shown in Equation (2) at a 2.5 km × 2.5 km grid to ensure the resolution of the results. For each of the grid nodes, all Envisat and Cryosat-2 observations within 2.5 km of the center of a grid node were used in the independent least-squares regression estimator. The same 3σ outlier rejection criteria and thresholds used in the least squares regression at the intersections above were also applied to ensure robust least square estimations. Due to a lack of observations, some of the 6 systematic bias parameters were not able to be estimated at some grid points. Here, interpolation of the corresponding parameters calculated at the intersections mentioned above were used. Then, the corrected elevation was calculated as Here, the mean altitude of the bin h 0 is removed to facilitate the study of elevation changes; then, the elevation for a given month t j is derived from where n is the number of corrected elevations in month t j . With ∆h lon 0 , lat 0 , t j computed for all months in a given bin, a time series of elevation can be formed. A floating median is a robust low-pass filter and was used in this study to identify and remove remaining outliers, similar to that in Schröder et al. (2019) [10]. Then, the monthly elevation estimations were smoothed using a spatiotemporal median filter to generate the final 5 km grid of the elevation time series. In the temporal domain, median filtering was performed on the monthly elevation estimations over 3 months (i.e., the given month, the previous month, and the next month) to eliminate high frequency errors. In the spatial domain, median filtering was operated over the monthly elevation estimations within a radius of 10 km from a grid node to minimize the impact of possible outliers. If the number of monthly elevation estimations within the 10 km radius was less than 9, the filtering radius would be expanded in steps of 10km until the radius was 30km, and the scaled median absolute deviation was used as the realistic error uncertainty of the elevation of the grid node as in Ewert et al. (2012) [32]. Then, the trend, acceleration, and other spatial-temporal characteristics of elevation and volume over the AIS can be derived from the elevation time series [4].

Calculation of Volume Changes for the AIS from the Combined Elevation Time Series
Long-term time series of volume changes are of great scientific value for the assessment of changes in the AIS and its response to climate change. After correcting for Glacial Isostatic Adjustment (GIA)-induced vertical crustal deformation and elastic solid earth rebound effects, the altimetric volume time series can be calculated from the elevation time series. In this study, the ICE-6G_C model [33] was used to correct for the uplift due to GIA. Elastic contributions were corrected by multiplying the elevation change by a scale factor α = 1.0205, which was inferred in Groh et al. (2012) [34]. i.e., by applying a scaling factor derived from a 50 km autocorrelation radius to the squared uncertainties. To analyze and interpret volume changes of the AIS, we also calculated the surface volume changes due to surface processes (surface mass balance variations, liquid water processes, and firn compaction) using the Institute for Marine and Atmospheric Research Firn Densification Model (IMAU-FDM) [35].

The Total Intermission Biases and the Intramission A-D Bias
A larger amount of observations from different missions are used simultaneously in the least-squares model to account for the total intermission biases and the intramission A-D bias is the highlight of our proposed algorithm. Figure 1 shows estimates of these biases derived at the intersections of the ascending and descending satellite ground tracks of Envisat. In this study, the mode mask boundary between CryoSat-2's LRM and SARIn modes is used to simply divide flat areas and steep areas: CryoSat-2 uses the SARIn mode to measure the steep areas at the margins of the ice sheets, and the LRM mode over the flat interior of the ice sheets [2,29]. Corresponding statistics can be found in Table 1.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 21 intersections of the ascending and descending satellite ground tracks of Envisat. In this study, the mode mask boundary between CryoSat-2's LRM and SARIn modes is used to simply divide flat areas and steep areas: CryoSat-2 uses the SARIn mode to measure the steep areas at the margins of the ice sheets, and the LRM mode over the flat interior of the ice sheets [2,29]. Corresponding statistics can be found in Table 1. In terms of the intramission A-D bias, the spatial distributions of CryoSat-2 and Envisat exhibit similar patterns. Over the interior of the AIS, the bias is close to zero, while in the steeper SARIn area, the amplitude of bias exceeds 1 m (Figure 1c (Table 1). Whether using Cryosat-2 or Envisat, the mean and RMS of the intramission A-D bias over the SARIn area are much greater than those over the LRM area. This is mainly because the azimuthally anisotropic effects of polarization on altimetry measurements can be negligible over flat and homogeneous areas due to the similar dielectric conditions and low mispointing, but the opposite is true in steep and very active areas [27]. In general, the spatial distribution patterns and magnitude of the intramission A-D bias calculated in this study are consistent with those in Khvorostovsky (2012) [12], Davis (1997) at the intersections of Envisat's ascending tracks and descending tracks. The mode mask boundary between CryoSat-2's LRM and SARIn modes [29] is shown in orange.
In terms of the intramission A-D bias, the spatial distributions of CryoSat-2 and Envisat exhibit similar patterns. Over the interior of the AIS, the bias is close to zero, while in the steeper SARIn area, the amplitude of bias exceeds 1 m (Figure 1c,f). The median, mean, RMS, and the 10th and 90th percentile range (P 90 -P 10 ) over the SARIn mode area are 0.00 m, −0.14 m, 4.71 m, and 1.81 m for CryoSat-2 and 0.00 m, −1.38 m, 13.86 m, and 0.94 m for Envisat, while those over the LRM mode area are relatively close (Table 1). Whether using Cryosat-2 or Envisat, the mean and RMS of the intramission A-D bias over the SARIn area are much greater than those over the LRM area. This is mainly because the azimuthally anisotropic effects of polarization on altimetry measurements can be negligible over flat and homogeneous areas due to the similar dielectric conditions and low mispointing, but the opposite is true in steep and very active areas [27]. In general, the spatial distribution patterns and magnitude of the intramission A-D bias calculated in this study are consistent with those in Khvorostovsky (2012) [12], Davis (1997) [27], and McMillan et al. (2014) [5].
Given the total intermission biases between CryoSat-2 and Envisat, it can be seen from Figure 1 that all of them share the similar spatial distribution patterns. These biases are highly correlated to the topography slope of the ice sheets. In steep and complex areas, these biases can reach up to several meters, while in the flat interior of the AIS, these biases appear to be relatively uniform at a magnitude of about 0.2 m. According to Frappart et al. (2016) [15], the intermission bias in the flat interior of the AIS is mainly caused by orbital errors and differences in the position of the centers of gravity of the two satellites and the phase of the antennae, while the bias over the steep area originates from the geographical small-scale bias. Thus, the total intermission bias in the flat interior of the AIS shows a uniform distribution, but a higher amplitude of bias is found where the topography is steeper. Overall, the total intermission biases estimated in this study agree with those found in Frappart et al. (2016) [15]. This agreement gives us confidence in our updated least-squares model.

Comparison with OIB and ICESat
The results of our least-squares model are validated against the laser altimeter elevation datasets obtained from the OIB project and ICESat. Both datasets are taken from near-inferred laser altimeters, and their accuracy has been validated [36,37]. The results of our validation are shown in Figure 2.  [40] and the OIB measurements, which are comparable with our combined time series in their processing strategy and resolution.  However, OIB ATM elevation observations are strongly concentrated along the outlet glaciers of West Antarctica and the Antarctic Peninsula. As a supplement, we also used ICESat data to validate our merged time series. For ICESat, the average elevation change rates for the period coinciding with ICESat over the AIS are first derived from our combined elevation time series by fitting a first-order polynomial and a sine with a 1-year period to account for the average elevation Surface elevation changes can be obtained wherever there are coincident OIB ATM elevation measurements. To validate our results, we used the surface elevation changes from OIB ATM L4 [38]. These changes were computed from the OIB ATM elevation measurements observed during 2002-2016. The differences between elevation changes derived from OIB and the combined elevation time series (∆h = ∆h OIB − ∆h ALT ) are shown in Figure 2a. To gain insight into the results of the intercomparison, we separated the evaluation results according to where or when the OIB ATM elevation measurements were observed to derived the OIB surface elevation changes for a grid cell, hereby referred to as the steeper SARIn area and the flat LRM area, and the periods before 2010, after 2010, and from 2002 to 2016. The OIB surface elevation change from 2002 to 2016 here refers to that is computed from a pair of coincident measurements, with one obtained before 2010 and the other after 2010. This allowed us to evaluate our time series in terms of their temporal and spatial distribution. The mean and RMS of these differences are given in Table 2. Meanwhile, statistics of the comparison between the OIB product and time series of Schröder et al.'s (2019) [10] elevation time series are also given in Table 2 to demonstrate the effects of different approaches to account for the intermission biases on the combined elevation time series. However, as shown in Figure 2a, the validation results are mainly concentrated along the outlet glaciers of the West Antarctic Ice Sheet (WAIS) and the Antarctic Peninsula Ice Sheet (APIS). Thus, intercomparison over the LRM area may be less persuasive. Over the SARIn area, the RMS of our results is nearly 30% smaller than that in Schröder et al. (2019) [10]. This suggests that our method is effective in correcting the systematic biases in the SARIn area. The decrease of the RMS by nearly 40% of the comparison of from 2002 to 2016 between ours results and that in Schröder et al. (2019) [10] further demonstrates the effectiveness of our method in intermission bias correction. Table 2. Statistics of the comparison between elevation changes derived from the combined elevation time series and those from Operation IceBridge (OIB) observations. The mean and RMS of the differences between elevation changes over the AIS and CryoSat-2 LRM and SARIn mask areas, as well as those acquired before 2010, after 2010, and from 2002 to 2016, are given. Note that the intercomparison results are mainly concentrated along the outlet glaciers of West Antarctica and the Antarctic Peninsula. To make a consistent intercomparison, we used all surface elevations during 2009-2018 from OIB ATM L2 [39] within a 2.5-kilometer radius and a 4-day interval of our radar altimetry grid nodes. When comparing with the OIB surface elevations, the mean altitude of the bin (i.e., the h 0 in Equation (2)) is added back into the combined elevation time series. Figure 2b shows the map of the differences between the combined elevation time series and the OIB measurements. The differences are again the greatest in steep regions of the AIS, which is likely because that the radar altimeter has larger measurement errors in steeper areas due to the higher surface slope and roughness. Over the AIS, the median and RMS of the differences are 5.7 m and 28.7 m, respectively. They are better than the median of 25.4 m and RMS of 138.6 m of the differences between the interpolated grid cells of DEM in Slater et al. (2018) [40] and the OIB measurements, which are comparable with our combined time series in their processing strategy and resolution.

Different Regions
However, OIB ATM elevation observations are strongly concentrated along the outlet glaciers of West Antarctica and the Antarctic Peninsula. As a supplement, we also used ICESat data to validate our merged time series. For ICESat, the average elevation change rates for the period coinciding with ICESat over the AIS are first derived from our combined elevation time series by fitting a first-order polynomial and a sine with a 1-year period to account for the average elevation change rate and seasonal variations. Subsequently, the average elevation change rates from our combined results within a 10 km radius of an ICESat grid node are compared with the corresponding average ICESat rate. A map of the differences is shown in Figure 2c. Again, the largest differences are observed at steep areas. In the flat interior of the AIS, the differences are small. The higher surface slope and roughness are still the main sources of these differences. Changes in the penetration depth of the microwave signal into the snow due to changes in the snow surface characteristics is another important source of difference [41]. The median and RMS of the differences over the AIS are −1.1 cm/y and 10 cm/y, respectively, while the median and RMS are −1 cm/y and 11 cm/y for the differences between the average ICESat elevation change rates and those in Schröder et al.

Elevation Time Series
As mentioned above, the average rates of elevation change over any subinterval can be obtained from the elevation time series using a fitting method with a first-order polynomial and a sine with a 1-year period. Unlike the information obtained directly from elevation observations, this method provides more accurate average rates for specific periods.

Elevation Time Series
As mentioned above, the average rates of elevation change over any subinterval can be obtained from the elevation time series using a fitting method with a first-order polynomial and a sine with a 1-year period. Unlike the information obtained directly from elevation observations, this method provides more accurate average rates for specific periods.    [10]. Many previous studies have attributed the thickening of the ice sheet over the period in Dronning Maud Land to increased snowfall [5,10,15,43,44].
Lake Vostok is the largest of the subglacial lakes located far inland of the EAIS and is used to represent the slow-flowing interior. The corresponding elevation time series are shown in Figure 4d. The elevation change rates at the selected points of 3 ± 0.7 mm/y (Point A), 1.4 ± 0.5 mm/y (Point B), and 0.7 ± 0.5 mm/y (Point C) are in good agreement with the estimates of −3 ± 9 mm/y ~ 5 ± 9 mm/y from altimetry data in Schröder

Volume Time Series for the AIS
To quantitatively examine the changes in volume at a regional scale, the volume time series for six typical regions selected based on elevation change features are shown in Figure 5. It is very challenging to use radar altimetry data to estimate the surface elevation changes of ice sheets due to signal penetration into the snow. Fortunately, Legrésy and Rémy (1998) [45], Davis and Ferguson (2004) [46], and Slater et al. (2018) [40] have shown that the annual and subannual changes in the snowpack volume/surface scattering ratio can produce spurious changes in altimetric elevation. Therefore, in Figure 5, we also show the time series without seasonal oscillations. To clearly illustrate the relationship between surface volume changes and altimetric volume changes, time series with linear trends and seasonal oscillations removed are shown in Figure 6.
challenging to use radar altimetry data to estimate the surface elevation changes of ice sheets due to signal penetration into the snow. Fortunately, Legrésy and Rémy (1998) [45], Davis and Ferguson (2004) [46], and Slater et al. (2018) [40] have shown that the annual and subannual changes in the snowpack volume/surface scattering ratio can produce spurious changes in altimetric elevation. Therefore, in Figure 5, we also show the time series without seasonal oscillations. To clearly illustrate the relationship between surface volume changes and altimetric volume changes, time series with linear trends and seasonal oscillations removed are shown in Figure 6. It can be clearly seen that there are strong interregional variations in the altimetric volume curves for smaller regions within Antarctica during 2002-2019. The Amundsen Sea sector (Figure 5a) displays a large coherent trend featuring a volume loss of −82.2 ± 2.4 km 3 /y with a significant acceleration of −5.4 ± 0.3 km 3 /y2. The time series reveals that this region experienced significant acceleration around 2006. Surface processes cannot cover the large coherent trends of accelerating volume loss. This means that sustained increases in ice discharge are the main reason for accelerated thinning, which agrees with many studies using observations from the Gravity Recovery and Climate Experiment (GRACE) [41][42][43][44]. There is a good consistency among the interannual changes in both volume time series (Figure 6a), with altimetric volume changes lagging behind those of the surface volume by several months.
The Getz sector and the Bellingshausen Sea sector are the other two regions of volume loss in the WAIS. The Getz sector has experienced a rapid volume loss of 21.1 ± 2.1 km 3 /y and acceleration of −1.4 ± 0.2 km 3 /y2 (Figure 5c). Figure 5e suggests that the Bellingshausen Sea sector has progressed to a less negative volume change at a rate of −3.4 ± 2.4 km 3 /y and an acceleration of −0.9 ± 0.3 km 3 /y2. The time series reveals that these two regions both experienced small losses in 2002-2006, followed by an increase in volume loss in 2007-2015 and regained volume in 2016. However, in contrast with the Bellingshausen Sea sector, discrepancies between altimetric and surface volume changes are very large in some subintervals in the Getz sector. This can be clearly seen from the time series of interannual changes (Figure 6c,e). Moreover, the interannual changes in altimetric volume in the Bellingshausen Sea sector are in good agreement with those of the surface volume change, except for a slight lag. This means that surface processes control the volume change in the Bellingshausen Sea sector, while ice discharge sometimes dominates the volume changes in the Getz sector.
In Dronning Maud Land, the volume change and acceleration for the entire period are 7.8 ± 3.1 km 3 /y and −4.9 ± 0.4 km 3 /y 2 , respectively. However, the relevant time series reveals that there is a significant change in its volume trend partway through the period of 2002-2019 (Figure 5b). Prior to the end of 2008, this area showed a slow increase in volume. Subsequently, the Dronning Maud Land experienced a continuous increase in volume gain. It can be seen from Figure 5b that surface processes dominate the volume change trends in Dronning Maud Land. The impact of two extreme accumulation events in 2009 and 2012 [48,49] can also be clearly seen in the volume time series, especially in Figure 6b. Several months after the surface volume changes, the altimetric volume responds. Glacier sector (f)) from our combined altimetric time series (green dots) and IMAU-FDM (purple triangles). After removing the long-term trends, seasonal oscillations were removed using a 13-month moving average. Due to the defects of the moving average algorithm, the 6-month values at the beginning and end of the time series are not displayed. The regions covered by each localization are shaded red in the inset. The boundaries refer to the Antarctic drainage systems in [47].
In Dronning Maud Land, the volume change and acceleration for the entire period are 7.8 ± 3.1 km 3 /y and −4.9 ± 0.4 km 3 /y 2 , respectively. However, the relevant time series reveals that there is a significant change in its volume trend partway through the period of 2002-2019 (Figure 5b). Prior to the end of 2008, this area showed a slow increase in volume. Subsequently, the Dronning Maud Land experienced a continuous increase in volume gain. It can be seen from Figure 5b that surface processes dominate the volume change trends in Dronning Maud Land. The impact of two extreme accumulation events in 2009 and 2012 [48,49] can also be clearly seen in the volume time series, especially in Figure 6b. Several months after the surface volume changes, the altimetric volume responds.
In contrast with Dronning Maud Land, Princess Elizabeth Land, and the Totten Glacier sector show a gradual change from initial volume accumulation to volume loss (Figure 5d,f). The accelerations in volume change of these sectors are −3.0 ± 0.2 km 3 /y 2 and −3.0 ± 0.4 km 3 /y 2 , respectively. The fluctuations of volume changes in 2008-2009 caused by surface process anomalies are most obvious in the Totten Glacier sector, as previously reported in Li et al. (2016) [42]. The trends of volume changes, especially accelerations, are dominated by surface processes in the two regions. The rapid volume loss of the Totten Glacier sector began in 2010. The difference is that the interannual oscillations of the Totten Glacier sector are consistent with its surface volume changes for the whole period, while those in Princess Elizabeth Land are not. Figure 7 shows the derived volume time series for the entire AIS and the subregions of EAIS, WAIS, and APIS north of 81.5 • S over 2002-2019. An overall volume loss of −68.7 ± 8.1 km 3 /y and an acceleration in loss of −5.5 ± 0.9 km 3 /y 2 during 2002-2019 were detected. Overall, increasing volume loss is mainly a result of volume changes of the WAIS. The volume loss of the WAIS during 2002-2019 was −88.4 ± 3.9 km 3 /y with an acceleration of −8.1 ± 0.4 km 3 /y 2 . The APIS contributes to a small part of volume loss at a lower rate of −18.8 ± 2.2 km 3 /y, and its loss gradually slows, with an acceleration of +1.9 ± 0.2 km 3 /y 2 . The EAIS has experienced a steady volume gain at a rate of 38.5 ± 6.7 km 3 /y. The acceleration of the EAIS is 0.8 ± 0.8 km 3 /y 2 .
accelerations in volume change of these sectors are −3.0 ± 0.2 km 3 /y 2 and −3.0 ± 0.4 km 3 /y 2 , respectively. The fluctuations of volume changes in 2008-2009 caused by surface process anomalies are most obvious in the Totten Glacier sector, as previously reported in Li et al. (2016) [42]. The trends of volume changes, especially accelerations, are dominated by surface processes in the two regions. The rapid volume loss of the Totten Glacier sector began in 2010. The difference is that the interannual oscillations of the Totten Glacier sector are consistent with its surface volume changes for the whole period, while those in Princess Elizabeth Land are not. Figure 7 shows the derived volume time series for the entire AIS and the subregions of EAIS, WAIS, and APIS north of 81.5° S over 2002-2019. An overall volume loss of −68.7 ± 8.1 km 3 /y and an acceleration in loss of −5.5 ± 0.9 km 3 /y 2 during 2002-2019 were detected. Overall, increasing volume loss is mainly a result of volume changes of the WAIS. The volume loss of the WAIS during 2002-2019 was −88.4 ± 3.9 km 3 /y with an acceleration of −8.1 ± 0.4 km 3 /y 2 . The APIS contributes to a small part of volume loss at a lower rate of −18.8 ± 2.2 km 3 /y, and its loss gradually slows, with an acceleration of +1.9 ± 0.2 km 3 /y 2 . The EAIS has experienced a steady volume gain at a rate of 38.5 ± 6.7 km 3 /y. The acceleration of the EAIS is 0.8 ± 0.8 km 3 /y 2 .     Figures 7 and 8 show that surface volume changes can only explain a small proportion of the overall volume changes in the AIS. The biggest difference in long-term trends comes from the WAIS, although its interannual variations in surface and altimetric volume changes show good consistency. According to previous studies, the differences in long-term trends in the WAIS are possibly due to ice dynamics [4,5,44,50]. Similarly, because of dynamic thinning, surface processes cannot fully explain the altimetric volume changes in the APIS, especially during 2002-2006. Following the disintegration of the Larson B Ice Shelf in 2002, increased ice discharge occurred [51,52]. Although surface volume changes in some subregions are in good agreement with their altimetric volume changes, there are differences in both the long-term trends and the interannual variations between the two integrated volume time series for the EAIS. These differences can also be seen in Schröder et al. (2019) [10]. dynamic thinning, surface processes cannot fully explain the altimetric volume changes in the APIS, especially during 2002-2006. Following the disintegration of the Larson B Ice Shelf in 2002, increased ice discharge occurred [51,52]. Although surface volume changes in some subregions are in good agreement with their altimetric volume changes, there are differences in both the long-term trends and the interannual variations between the two integrated volume time series for the EAIS. These differences can also be seen in Schröder et al. (2019) [10].  [47].

Discussion
Accurately correcting intermission biases is crucial for combining observations from different radar altimetry missions to study the long-term changes of ice sheets. Imperfect methods may lead to non-self-consistent combined results. This is why there are leap or step signals at the junctions between ERS-2 and Envisat and between Envisat and CryoSat-2 in the elevation time series shown in Figure 10 of Schröder et al. (2019) [10]. This can be attributed to the uncertainties in the intermission cross-calibration due to the short overlapping period. In this study, we constructed an updated planefitting least-squares regression model by introducing the comprehensive parameters of the intermission bias and the A-D bias and their internal consistency into the existing least-squares regression model. In contrast with previous methods, it allowed us to use a larger amount of observations in the same regression model for the integrated least-squares adjustment to improve the accuracy of cross-calibration of the data from different missions, and thereby to ensure a selfconsistent and accurate elevation time series. We used Envisat and Cryosat-2 observations to test our model and successfully established a 5 km grid of monthly elevation time series over the AIS. The estimations of the total intermission biases and the intramission A-D bias are consistent with those in previous studies [5,15,18]. Moreover, the self-consistent elevation time series at the joint of the two satellites ( Figure 4) show that our algorithm is effective in cross-calibration. The intercomparison between the surface elevation changes from OIB ATM L4 and from the combined elevation time series also demonstrate the effectiveness of our method in creating time series using data from

Discussion
Accurately correcting intermission biases is crucial for combining observations from different radar altimetry missions to study the long-term changes of ice sheets. Imperfect methods may lead to non-self-consistent combined results. This is why there are leap or step signals at the junctions between ERS-2 and Envisat and between Envisat and CryoSat-2 in the elevation time series shown in Figure 10 of Schröder et al. (2019) [10]. This can be attributed to the uncertainties in the intermission cross-calibration due to the short overlapping period. In this study, we constructed an updated plane-fitting least-squares regression model by introducing the comprehensive parameters of the intermission bias and the A-D bias and their internal consistency into the existing least-squares regression model. In contrast with previous methods, it allowed us to use a larger amount of observations in the same regression model for the integrated least-squares adjustment to improve the accuracy of cross-calibration of the data from different missions, and thereby to ensure a self-consistent and accurate elevation time series. We used Envisat and Cryosat-2 observations to test our model and successfully established a 5 km grid of monthly elevation time series over the AIS. The estimations of the total intermission biases and the intramission A-D bias are consistent with those in previous studies [5,15,18]. Moreover, the self-consistent elevation time series at the joint of the two satellites ( Figure 4) show that our algorithm is effective in cross-calibration. The intercomparison between the surface elevation changes from OIB ATM L4 and from the combined elevation time series also demonstrate the effectiveness of our method in creating time series using data from different missions. In addition, comparisons with OIB-derived surface elevations and ICESat surface elevation changes also increased confidence in our merged time series. It should be noted that the difference in footprint size between Envisat data and CryoSat-2 SARIn-mode data may be a source of uncertainty in the estimates of intermission bias and in creating of the merged time series. Following the study by Adusumilli et al. (2018) [16], we believe that the point of closest approach (POCA) applied to generate Envisat and CryoSat-2 Level 2 products reduces the difference between topography sampled with the different satellite radar altimeter techniques. Certainly, it still needs further study to accurately calibrate the different techniques against each other in the future.
Our gridded monthly elevation time series allow for the analysis of temporal and spatial evolutions over 2002-2019. Well-known patterns of thinning and thickening were also confirmed in this study. Volume losses are mainly concentrated in the coastal sectors of the WAIS, especially in the Amundsen Sea sector, while the volume gains are mainly observable in the Dronning Maud Land of the EAIS [2,4,5,10,14]. A simple comparison of the elevation change rates of the two periods (2002-2010 vs. 2010-2019) reveals an expansion of thinning in the Pine Island Glacier and interannual fluctuations of elevation changes in Princess Elizabeth Land and the Bellingshausen Sea sector. Although these findings were previously reported [10,53,54], they indicate that our time series are reliable.
Due to signal penetration into the snow, the elevation changes derived from radar altimeters might not reflect the actual surface elevation changes [17,55]. However, including the backscatter coefficient in the least-squares regression model can only account for a part of the effect of surface penetration, as shown in Simonsen and Sorensen (2017) [3] and Sørensen et al. (2018) [17]. Indeed, there were spurious seasonal fluctuations caused by surface penetration in our combined results. However, this spurious signal also exists in Schröder et al. (2019) [10], even though the authors reprocessed their radar altimetry observations with a low threshold retracker able to produce spurious changes in elevation due to changes in surface penetration [27,28]. Therefore, radar altimetry data is generally considered to be very difficult to use for studying the changes caused by surface processes directly, especially for surface accumulation events. Nevertheless, the Greenland 2012 melt event provides us with an extreme example for which radar altimeters can detect the changes caused by surface processes. The formation of a refrozen melt layer can increase the reflective surface for radar altimeters instantaneously [55]. Comparably, other surface processes, such as firn compaction, can also form ice lenses to change the reflecting surface, allowing them to be observed by radar altimeters. That is to say, after a time-consuming process of firn compaction, accumulation events can be observed with time lags by radar altimeters. After removing seasonal variations to minimize the effects of surface penetration, more reliable estimates of the interannual and long-term variations caused by surface processes can be clearly seen from our combined altimetric time series, such as snowfall in Dronning Maud Land in 2009 and 2012 and snow accumulation in 2008-2009 in the Totten Glacier sector.
Therefore, together with surface volume changes, the integrated volume time series over larger regions enabled us to further improve our understanding of elevation change processes. The agreement between altimetric and surface volume changes in Dronning Maud Land indicates that surface processes dominate the variability of its volume changes. In other words, changes in ice dynamics do not significantly affect the time variability of its volume signal. This confirms that they are distal to areas with fast ice flow [56]. Modest thickening in Dronning Maud Land was associated with increases in snowfall due to its extreme accumulation events in 2019 and 2012 [54,57]. The long-term trend of altimetric volume changes of Princess Elizabeth Land are consistent with those of the surface volume. This result agrees with that in Schröder et al. (2019) [10], showing that thinning here is related to less snowfall. In the Totten Glacier sector and the Bellingshausen Sea sector, the altimetric and surface volume changes also agree well. The difference is that excessive ice discharge and grounding line retreats have made Totten Glacier and several Bellingshausen Sea glaciers thin [42,56]. However, the surface volume changes are so large that they still dominate the volume changes in these regions, such as the accelerated thinning of Totten Glacier over 2010-2015 [42].
As reported in previous studies [4,5,10,44,54], changes in the Amundsen Sea sector ice sheet are mainly driven by rapid ice discharge. Mouginot et al. (2014) [53] and Harig et al. (2015) [57] reported that flow acceleration and widening of Pine Island Glacier and Thwaites Glacier led to rapid acceleration loss in the Amundsen Sea sector. These factors agree with the pronounced differences that we observed between altimetric and surface volume changes here. In the Getz sector, we found that neither ice dynamics nor surface processes dominate volume changes over the entire period of 2002-2019. As reported in Shepherd et al. (2019) [54], the patterns of elevation changes in this sector reflect a complex mixture of surface processes and dynamic ice imbalances.
Overall, the accelerated ice volume loss in the coastal regions of WAIS continues to outpace the gains made in the EAIS. The volume loss of the AIS is mainly caused by ice discharge in the Amundsen Sea sector. Meanwhile, in Princess Elizabeth Land, the Totten Glacier sector of the EAIS, and the Bellingshausen Sea sector and Getz sector of the WAIS, surface processes also make a considerable contribution to volume loss. The contribution of surface processes in Princess Elizabeth Land and the Totten Glacier sector of the EAIS to the acceleration of volume loss over the AIS also cannot be ignored. Volume gains due to surface processes were mainly observed in the EAIS, especially in Dronning Maud Land.

Conclusions
In this study, we presented an updated plane-fitting least-squares regression model that combines different radar altimetry missions. Using this approach, we constructed a monthly elevation time series using 5 km grid cells over the AIS. Further analysis of the output estimates for the systematic biases and elevation time series shows that the intermission bias and the A-D bias were successfully corrected. Validation with airborne and ICESat laser altimetry observations indicates that our combined gridded monthly elevation time series are reliable.
Further analysis shows that our resulting elevation time series can provide detailed insight into temporal and spatial evolution of elevation for the AIS over 2002-2019. In addition to the well-known patterns of thinning and thickening, acceleration and increases in thinning in Pine Island Glacier and Thwaites Glacier, as well as interannual fluctuations in Princess Elizabeth Land and the Bellingshausen Sea sector, were shown very clearly. Overall, the volume loss in the WAIS continues to outpace the gains made in the EAIS. From 2002 to 2019, the total volume change estimates for the AIS were −68.7 ± 8.1 km 3 /y, with an acceleration of −5.5 ± 0.9 km 3 /y 2 . Most of the volume loss was from the Amundsen Sea sector and the Getz sector, and most of the acceleration was from the Amundsen Sea sector, Princess Elizabeth Land, and the Totten Glacier sector.
After removing the seasonal variations, we obtained the interannual and long-term variations caused by surface processes in our combined time series. With the help of additional surface volume changes derived from IMAU-FDM, we determined that ice discharge in the Amundsen Sea sector provides the largest accelerated volume loss of the AIS. However, surface processes still play an important role in the volume changes of the AIS. In addition to thickening in Dronning Maud Land, surface processes in regions such as the Totten Glacier sector, Princess Elizabeth Land, and the Bellingshausen Sea sector are the main sources of acceleration in volume loss for the AIS.