Altimetry , gravimetry , GPS and viscoelastic modeling data for the joint inversion for glacial isostatic adjustment in Antarctica ( ESA STSE Project REGINA )

The poorly known correction for the ongoing deformation of the solid Earth caused by glacial isostatic adjustment (GIA) is a major uncertainty in determining the mass balance of the Antarctic ice sheet from measurements of satellite gravimetry and to a lesser extent satellite altimetry. In the past decade, much progress has been made in consistently modeling ice sheet and solid Earth interactions; however, forward-modeling solutions of GIA in Antarctica remain uncertain due to the sparsity of constraints on the ice sheet evolution, as well as the Earth’s rheological properties. An alternative approach towards estimating GIA is the joint inversion of multiple satellite data – namely, satellite gravimetry, satellite altimetry and GPS, which reflect, with different sensitivities, trends in recent glacial changes and GIA. Crucial to the success of this approach is the accuracy of the spacegeodetic data sets. Here, we present reprocessed rates of surface-ice elevation change (Envisat/Ice, Cloud,and land Elevation Satellite, ICESat; 2003–2009), gravity field change (Gravity Recovery and Climate Experiment, GRACE; 2003–2009) and bedrock uplift (GPS; 1995–2013). The data analysis is complemented by the forward modeling of viscoelastic response functions to disc load forcing, allowing us to relate GIA-induced surface displacements with gravity changes for different rheological parameters of the solid Earth. The data and modeling results presented here are available in the PANGAEA database (https://doi.org/10.1594/PANGAEA.875745). The data sets are the input streams for the joint inversion estimate of present-day ice-mass change and GIA, focusing on Antarctica. However, the methods, code and data provided in this paper can be used to solve other problems, such as volume balances of the Antarctic ice sheet, or can be applied to other geographical regions in Published by Copernicus Publications. 494 I. Sasgen et al.: Altimetry, gravimetry, GPS and viscoelastic modeling data the case of the viscoelastic response functions. This paper presents the first of two contributions summarizing the work carried out within a European Space Agency funded study: Regional glacial isostatic adjustment and CryoSat elevation rate corrections in Antarctica (REGINA).


Introduction
Glacial isostatic adjustment (GIA), the viscoelastic deformation of the solid Earth in response to climate-driven ice and water mass redistribution on its surface, is poorly constrained in Antarctica. The primary reason is the sparseness of geological evidence for the past ice sheet geometry and local relative sea-level change. These are important constraints on the glacial forcing exerted and on the viscoelastic structure of the lithosphere and of the mantle, which together determine the signature of GIA (e.g., Peltier, 2004;Ivins and James, 2005;Whitehouse et al., 2012;van der Wal et al., 2015). The predictions of GIA in Antarctica remain ambiguous (Shepherd et al., 2012) and cause a large uncertainty in gravimetric mass balance estimates of the ice sheet of the order of the estimate itself (Martín-Español et al., 2016b). Measurements of bedrock uplift by GPS are inconsistent with the predictions of existing GIA. In many regions, uplift rates and thus mass increase due to GIA is overpredicted (Bevis et al., 2009), biasing estimates of present-day Antarctic icemass loss from Gravity Recovery and Climate Experiment (GRACE) to more negative values. However, for regions with a weak Earth structure, large uplift signals are recorded by GPS (e.g., Groh et al., 2012), which are likely caused by load changes within the past few thousand year, and are often not accurately represented in GIA predictions (Wolstencroft et al., 2015).
Much progress has been made in reconstructing the ice sheet evolution from geomorphological evidence (Bentley et al., 2014) and inferring the underlying Earth structure from seismic observations Heeszel et al., 2016). However, an independent approach to constraining GIA is to make use of the different sensitivities of the various types of satellite data to recent glacial changes and GIA, respectively. Separating both signals in a joint inversion approach has been pursued by, e.g., Wahr et al. (2000), Riva et al. (2009), Wu et al. (2010, Gunter et al. (2014) and Martín-Español et al. (2016a). Another approach used regional patterns of GIA from forward modeling and adjusted them to GIA uplift rates in Antarctica .
In this paper, we present methods and data inputs in preparation to solve the joint inversion problem for GIA in Antarctica. As the GIA process is gradual, causing an approximately constant rate of change within a decade, we first process the satellite data to recover optimal temporal linear trends. We focus on the trends derived for the time period [2003][2004][2005][2006][2007][2008][2009] in which GRACE and Ice, Cloud,and land Elevation Satellite (ICESat) operated simultaneously. Note that the stationarity of the trend is a key assumption underlying our approach when including GPS rates covering a longer time span (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013). However, limiting the GPS data to the time span 2003-2009 leads to a significant reduction in the number of stations for which reliable trends can be estimated and, hence, to a loss of spatial coverage. For comparison, the reader is referred to the data archive, in which GPS uplift rates for the time periods 2003-2009 and 2003-2013 are made available.
In this paper, we present refined procedures for estimating trends in the data sets on surface-ice elevation changes, surface displacement and gravity field changes. The rates of surface-ice elevation change from Envisat and ICESat satellite altimetry are improved by combining both data sets based on their respective uncertainties, increasing the spatial coverage and accuracy of the elevation rates (Sect. 2). Bedrock displacement from in situ networks of GPS stations in Antarctica are improved in coverage by allowing for campaignbased data and carefully assessing the uncertainty in the trend with a noise model (Sect. 3). Compared to the rates in Thomas et al. (2011) more stations and longer time series are also included The gravity field changes from GRACE are refined compared to previous work by optimizing the destriping filtering for the region of Antarctica (Sect. 4). The processing aims at fulfilling the requirement of the joint inversion to combine input data based on the same time period (not possible for GPS without having to ignore a large number of stations) and covering all of Antarctica, accompanied by a realistic description of the uncertainties.
We also present forward-modeling results of viscoelastic response functions to disc load forcing for the range of Earth structures likely to prevail in Antarctica (Sect. 5). The viscoelastic response functions allow us to combine the surface displacement and gravity changes based on the physical description of the Earth's viscoelastic response for a specified Earth structure. In addition, the response functions enable us to combine data sets of different spatial resolutions, as is the case for GPS, GRACE and altimetry.
The determination of viscoelastic response functions is a classic topic in solid Earth modeling (e.g., Peltier and Andrews, 1976), though uncommon in the application to joint inversion studies of satellite data. Although this paper focusses on Antarctica, the response functions and data processing techniques presented here are applicable to other regions. The response kernels represent a wide range of Earth structures and can be used for the separation of superimposed present-day (elastic) and past (viscoelastic) signatures of mass change in other regions with a similar Earth struc-ture, for example hydrological storage changes and GIA in North America and Alaska. The response functions give insight into the temporal and spatial scales of deformation expected for Antarctica and are crucial when combining the input data streams.
The data sets and modeling results presented in this paper are accessible in the PANGAEA archive (https://doi.org/10.1594/PANGAEA.87574); subsections provide user guidance and point to data and code stored in the archive. As mentioned above, the data sets and modeling results are of value to address other research questions as well. For example, the GPS rates provided are useful for the validation of forward-modeling GIA solutions, the GRACE gravity rates can be used for mass balance studies, and altimetry data [2003][2004][2005][2006][2007][2008][2009] can be extended with the ongoing CryoSat-2 mission to infer volumetric mass balances, also over the ice shelves. The viscoelastic response functions are based on Earth model parameters potentially suitable to other geographical regions, as well; they are useful for similar studies combining different data sets of geodetic observables, surface deformation, gravity field change and topographic change in glaciated areas.
The actual method of the joint inversion is described in a second contribution of the Regional glacial isostatic adjustment and CryoSat elevation rate corrections in Antarctica (REGINA) project team (Sasgen et al., 2017a). In this second paper, the resulting GIA estimate is also compared to previous studies.

ICESat elevation rate determination
We use along-track altimetry measurements from ICESat 633 Level 2, providing high-resolution elevation change observations for the period February 2003 to October 2009. Two corrections are applied to this data set: the range determination from transmit-pulse reference-point selection (centroid vs. Gaussian; Borsa et al., 2014) available from the National Snow and Ice Data Center (NSIDC) and the intercampaign correction (Hofton et al., 2013). The centroid-Gaussian correction is a well-established correction and has been incorporated into the latest ICESat release (634). Concerning the ICESat inter-campaign bias (ICB) correction, uncertainties are available at Hofton et al. (2013). Furthermore, several studies have determined this correction from different methodologies. For a summary of published ICESat ICB corrections, see Scambos and Shumman (2016). Note that these corrections are part of a widely accepted procedure, and their effect on the elevation rates and uncertainties caused by varying the processing choices have not been evaluated. Because the ICESat tracks do not usually overlap, a regression approach is used in which topographic slope (both across-track and along-track) and the rate of surfaceelevation change y h ICESat are simultaneously estimated using the "plane" method (Howat et al., 2008) over areas that span 700 m and are a few hundred meters wide. A regression is only performed if a plane has at least 10 points from four different tracks that span at least 1 year. Regression was carried out twice. First, individual elevation measurements with corresponding residuals outside the range of two standard deviations were detected; then, the regression was repeated omitting these outliers. The standard deviation of the regression coefficient, here taken as the uncertainty in the elevation rate, σ h (here, ICESat) is calculated by the propagation of the residuals of the input data and the estimated topographic heights, to the trend parameter (see Eq. 1 in Hurkmans et al., 2012), where e is the vector of residuals, n is the sample size (i = 1, 2, . . ., n) and x is the vector of input elevations with mean x. This standard deviation (σ ICESat ) takes into account the sample size and the variance in both input data and residuals of the regression (Hurkmans et al., 2012). The residuals of the regression are used as they quantify the approximation of fitting the data with a plane. The exact ICESat observation periods are shown in the Appendix (A1 , Table A1). Then, the elevation rate and its uncertainty are interpolated (bilinearly) to a common 10 km × 10 km grid in polar-stereographic projection (central latitude 71 • S; central longitude 0 • W; origin at the South Pole; WGS-84 reference ellipsoid).

Envisat elevation rate determination
We use a time series of elevation changes derived from alongtrack Envisat radar altimetry data for the interval January 2003 to October 2009 (coeval to ICESat time span). Elevation rates y h Envisat are obtained at points every 1 km along track by binning all the echoes within a 500 m radius. Then, a 10-parameter least squares model is fitted in order to correct for the across-track topography and changes in snowpack properties. The least square model is defined in Flament and Remy (2012). The estimated parameters include parameters determined for the backscatter, leading-edge width and tailing-edge slope, the mean altitude, quadratic surface slope parameters to define surface curvature, and a linear time trend. A digital elevation model was not used for the correction of the topographic slope. For processing reasons, the temporal resolution is resampled from 35 days to monthly periods for each grid cell before estimating the elevation rates. This has a minor effect on the elevation rate estimate (smaller than ±1 cm yr −1 ) and reduces the standard deviation by about 14 %. As with ICESat, the elevation rate is interpolated to a common 10 km × 10 km polar-stereographic grid (and 20 km × 20 km for download in the archive), and the standard deviations of the rates within each grid cell are taken as an estimate of the measurement uncertainty, σ Envisat , according to Eq. (1).

Combination of Envisat and ICESat
We produce a combined rate of surface-elevation change product from the ICESat and Envisat data sets for the Antarctic ice sheet, y h . The aim is to take advantage of the high spatial resolution of ICESat data and the high temporal resolution and high-track density of the Envisat data.
We combine the two altimetry data sets based on their common 10 km × 10 km polar-stereographic grid. At each location, the elevation rate with the smallest standard deviation is chosen from either Envisat or ICESat data sets. We prefer this masking procedure instead of a weighted average, in order to avoid introducing possible biases associated with gridded elevation rates of very high uncertainty. Figure 1 shows the resulting mask underlying the combination. It is evident that some grid points are only represented by either ICESat or Envisat. Most prominent is the narrowing of the polar gap with ICESat data, resulting from the 81.5 • S latitude limit for Envisat compared to 86 • S for ICESat due to satellite orbit inclination. On the Antarctic Peninsula, Envisat picks up some points that are not present due to a sparser track coverage in the ICESat data set. As expected, ICESat outperforms Envisat in terms of the uncertainty in the elevation rate over steep topographic slopes and along the ice sheet margins. This is due to the smaller footprint of the laser altimeter, its higher accuracy and lower slope-dependent uncertainty (e.g., Brenner et al., 2007). On some flat areas and over some faulty ground tracks, where ICESat data measurements are scarce, however, Envisat provides better temporal and spatial coverage, leading to better accuracy of the resulting elevation rates. The resulting combined data set of surface-elevation rates and its uncertainties are shown in Fig. 2.

Firn correction
The elevation rates derived from ICESat and Envisat are corrected for changes in the firn layer thickness using the firn compaction model of Ligtenberg et al. (2011), which is driven by the regional atmosphere and climate model RACMO2/ANT (Lenaerts et al., 2012). We determine the firn compaction for January 2003 to October 2009, with respect to the mean of the years 1979 to 2002, and estimate a temporal linear trend, h comp . The model output is re-gridded onto the 10 km × 10 km common grid using nearest-neighbor interpolation. The standard deviation of the re-gridding is less than 1 cm yr −1 , causing a maximum change of 2 % of the firn compaction rate. Note that the firn compaction model has a spatial resolution of 27 km, potentially neglecting finerscale processes relevant for the altimetry data. Clearly, the regridding uncertainty stated above is merely a minimum estimate, neglecting, for example, uncertainties in the calibration or the atmospheric forcing of the firn compaction model.
The data were resampled from every 2 days to monthly mean time periods for every grid cell before estimating el- Figure 1. Mask for the combination of Envisat and ICESat. ICE-Sat but not Envisat available -yellow; σ ICESat ≤ σ Envisat -green; σ ICESat > σ Envisat -turquoise; Envisat but not ICESat availableorange; and no data -blue. No interpolation is used. evation rates. As with the Envisat and ICESat data, no seasonal terms were co-estimated and removed (i.e., annual and semiannual). We do not apply an a priori correction for surface mass balance (SMB) trends, in accordance with GRACE processing (Sect. 5), which requires defining a climatological reference period. Note that applying the commonly used reference period (1979 to present) leads to spurious accumulation anomalies in the altimetry data (see Appendix A2, Fig. A1). The derivation of an adequate climatological reference epoch in the RACMO2/ANT simulations is in itself challenging and beyond the scope of this paper.
The total uncertainty in the rate of elevation change from satellite altimetry is calculated by where the standard deviation of the firn correction, σ Firn , is the formal regression uncertainty (neglecting model uncertainties, as these are not available), and we assume the error sources to be uncorrelated. It is recognized that neglecting uncertainties in the firn model leads to underestimated values of σ h . However, the magnitude of the firn correction itself is small (see Appendix A2) compared to the observational uncertainties, and the associated underestimation of σ h is likely to be small.

Data availability
Annual elevation trends from a combination of Envisat and ICESat data are provided for the time period between February 2003 and October 2009. Trends have been corrected for firn densification processes using RACMO2/ANT. Elevation trends are provided in a 20 km polar-stereographic grid (central meridian 0 • , standard parallel 71 • S) with respect to the WGS84 geoid. X and Y are given in kilometers, and the elevation rate and its standard deviation are given in meters per year.

ICESat elevation trend for the time period between February 2003 and October 2009
The data set is provided in a 10 km grid in polarstereographic projection (central meridian 0 • standard parallel 71 • S) with respect to the WGS84 geoid. X and Y are given in kilometers, and the elevation rate and its standard deviation are given in meters per year.

Envisat elevation trend for the time period between February 2003 and October 2009
The data set is provided in a 10 km grid in polarstereographic projection (central meridian 0 • , standard parallel 71 • S) with respect to the WGS84 geoid. X and Y are given in kilometers, and the elevation rate and its standard deviation are given in meters per year.

ICESat and Envisat combination for time period between February 2003 and October 2009
Elevation changes have been corrected for firn densification processes using a firn densification model. The data set is provided in a 10 km grid in polar-stereographic projection (central meridian 0 • standard parallel 71 • S) with respect to the WGS84 geoid. X and Y are given in kilometers, and the elevation rate and its standard deviation are given in meters per year.
2.5.4 Annual elevation trends from CryoSat-2 derived from a single trend covering the time period 2010-2013 An acceleration term in areas with dynamic thinning was added to the linear trend to obtain annual rates. Elevation trends are provided at 10 km resolution in a polarstereographic grid (central meridian 0 • , standard parallel 71 • S) with respect to the WGS84 geoid. X and Y are given in kilometers and the elevation rate and its standard deviation are given in meters per year.

Elevation changes from firn model
Annual firn densification rates over 2003-2013 rates obtained from RACMO2.3. Data are provided in a 27 km polarstereographic grid (central meridian 0 • , standard parallel 71 • S) with respect to the WGS84 geoid. X and Y are given in kilometers and the annual firn densification rates in meters per year.

Snow/ice density map
To perform the conversion of volume change to mass change, a density map is provided at 20 km resolution in a polarstereographic grid (central meridian 0 • , standard parallel 71 • S) with respect to the WGS84 geoid. X and Y are given in kilometers and density in kg m −3 . We provide the data set at this point for completeness; more details on the generation of this density map are given in Sasgen et al. (2017a).

ICESat-Envisat combination mask
The mask used for combining ICESat and Envisat is provided with a 10 km resolution in polar-stereographic coordinates..

498
I. Sasgen et al.: Altimetry, gravimetry, GPS and viscoelastic modeling data X and Y are coordinates in kilometers, and the ID indicates whether ICESat or Envisat has been used to construct the elevation change combination: -4: only Envisat was available -3: only ICESat was available -2: ICESat lower errors -1: Envisat lower errors.

GPS uplift rate estimation and clustering
The aim of the GPS time series analysis is to derive uplift rates, y u that represent the vertical ground motion at the sites as accurately and robustly as possible. We derive uplift rates based on GPS records from a total of 118 Antarctic sites. Data were processed from 1995 day of year (doy) 002 to 2013 doy 257 (decimal year 1995.0-2013.7), but data at individual sites are of varying length and quality. The processing and uplift rate and uncertainty estimation resemble that of Thomas et al. (2011), but with more recent processing software (GIPSY 6.2) and model updates (including second-order ionospheric and earth radiation models): an initial satellite orbit and clock estimation step is performed, using a carefully selected balanced stable global network of GPS sites (at the time of processing, Jet Propulsion Laboratory (JPL)-reprocessed orbits for these state-of-the-art options were not available). The orbits and clocks are then used to perform precise point positioning (PPP) processing of all the available Antarctic sites of interest. A mini-ensemble was created to investigate systematic processing uncertainties and the effects of possible systematic errors in the time series on uplift rates. The mini-ensemble investigation showed that decisions taken when analyzing time series tended to have larger effects on uplift rates and uncertainties than the effects of small processing strategy changes. Outliers and systematic errors, such as offsets due to equipment changes or other causes, were removed where possible. Due to the varying characteristics of the time series, it was not possible to use the same approach at all sites. The strategy was as follows (and is summarized in Appendix A3, Fig. A3). For sites with over 2000 days of data, uplift rates and associated uncertainties were estimated using the CATS software (Williams, 2008). We co-estimated a white-noise scale factor for the formal uncertainties and a power-law noise amplitude with the index fixed to −1 (flicker noise), along with the temporal linear trend (rate), seasonal (annual and semiannual) parameters, and sizes of the offsets (at the specified epochs).
The median values of the white-noise scale (1.6) factor and the power-law noise amplitude (13.4 mm), derived from these long time series, were then used to propagate rates and uncertainties for the shorter time series, for which CATS cannot produce reliable estimates of the error model. For the propagation, the time series with fewer than 2000 epochs are additionally subdivided into two categories; continuous sites (≥ 2.5 years), for which periodic parameters are estimated in the propagation of uncertainties, and very short continuous sites (< 2.5 years) and campaign sites for which periodic parameters are not estimated. For each campaign, 1 mm of noise was added when propagating the uncertainties, to allow for tiny differences when setting up equipment again.
Finally, for each site, the uplift rate y u and its uncertainty σ u are assessed by manually removing portions of the time series (for example deleting campaigns in turn). If the rate changes by an amount larger than the propagated uncertainty for the site, the uncertainty is determined as ± the maximum difference in rate, and the rate is adjusted, if necessary, to the values of the most likely part of the range. Sites with only two campaigns were assigned an uncertainty of ±100 mm yr −1 , unless there was further evidence for or against the existence of systematic errors. Table 1 summarizes the rate estimation methods and the number of sites for each. For further details and full information on individual rates and time series, see Petrie et al. (2018a) for a full description of the processing and ensemble evaluation and Petrie et al. (2018b) for details of time series analysis and rate and uncertainty estimation. Further work was undertaken to combine or "cluster" the rates regionally for inclusion in the estimation process -see Section 3.2 of the second REGINA paper (Sasgen et al. 2017a) for more details..

Comparison with existing results
Next, we briefly compare the uplift rates at individual sites (data span: decimal year 1995.0-2013.7) derived from the GPS processing described above with those available from three previous studies: Thomas et al. (2011) Thomas et al. (2011) rates are in ITRF2005 (which has negligible scale or translation differences to ITRF2008) and the Argus et al. (2014) rates are in a reference frame specific to their paper, which they note yields 0.5 mm yr −1 more uplift than ITRF2008 at high southern latitudes. All rates from Argus et al. (2014) in Tables 2, 3 and 4 are shown as given in the original paper.
Due to the large number of Antarctic sites (in total 118), we focus the comparison on the uplift rates and uncertainties derived by the methods "cats, cats" (Table 2) and "prop, prop" (Table 3). Uplift rates resulting from our study are provided in Appendix A4 for all sites (Table A2). Table A3 shows comparisons with the values of Thomas et al. (2011) and Argus et al. (2014) for "prop, eman" sites not shown in the main text. All uplift rates, y u , are in mm yr −1 , with un- CATS rate and uncertainty ("cats, cats") 18 CATS rate, manually increased uncertainty ("cats, eman") 2 Propagated rate and uncertainty ("prop, prop") 28 Propagated rate and manually increased uncertainty ("prop, eman") 50 Manually adjusted rate and manually increased uncertainty ("rman, eman") 20 Table 2. Uplift rates y u and associated uncertainties σ u (mm yr −1 ) for selected GPS sites with more than 2000 epochs of data, compared to data published by Thomas et al. (2011) and Argus et al. (2014). Temporal components and noise characteristics are derived using the CATS software (Williams, 2008), i.e., the "cats, cats" method.
Site REGINA Thomas et Argus et al. (2011) al. (2014) certainties reflecting 1σ standard deviations, σ u . Sites with particularly complex nonlinear time series, such as those at O'Higgins (ohi2, ohig) and Palmer (palm) in the Antarctic Peninsula, are omitted here, as comparison with different studies is potentially misleading due to the effects of different measurement time periods. Table 2 shows data for selected sites with long time series, where uplift rate and uncertainty were derived using the CATS software (Williams, 2008). Uplift rates at the majority of the GPS sites agree within uncertainty, except syog (Syowa), where the REGINA value is between that from the other two studies. The uncertainty limits for the REGINA value and the Argus et al. (2014) just meet at 0.9 mm yr −1 , even when allowing for the ∼ 0.5 mm difference in reference frames, but the Thomas et al. (2011) value does not. This may be due to the fact that Thomas et al. (2011) estimate two offsets in the series. Table 3 shows uplift rate comparisons for sites where the "prop, prop" method was used; the noise characteristics are derived from median values from CATS software results for longer site records and then propagated in the parameter estimation Table 3. Uplift rates y u and associated uncertainties σ u (mm yr −1 ) for selected GPS sites with fewer than 2000 epochs for data, compared to data published by Thomas et al. (2011) and Argus et al. (2014). Noise characteristics are derived median values from CATS software results for longer station records and propagated in the parameter estimation ("prop, prop" method). See Appendix A4, Table A2, for a full list of rates from this study.
Site REGINA Thomas et Argus et al. (2011) al. (2014) in which annual and semiannual parameters were also estimated along with the trend. Again, the rates agree within uncertainty, except for site belg where there is disagreement with Thomas et al. (2011). This may be due to their shorter data span. Table 4 shows comparisons for sites where the REGINA rates and uncertainties have been manually evaluated based on the spread of rates obtained by subsampling the time series ("rman" method). There is a large difference (over 10 mm yr −1 ) in the values at capf (Cape Framnes) between the REGINA value (4.0 ± 1.4 mm yr −1 ) and the Argus et al. (2014) value (15.0 ± 4.2 mm yr −1 ). Interestingly, the Wolstencroft et al. (2015) rate values for bean, gmez, lntk, mkib and trve are all systematically higher than the REGINA values, by an average of just over 3 mm yr −1 , and the uncertainties we assigned are also several times larger. The systematic differences between Wolstencroft et al. (2015) and the REGINA values for Palmer Land are currently unexplained and a matter of ongoing investigation.

Data availability
The GPS data and related code are directly accessible in the PANGAEA repository: http://hs.pangaea.de/model/ Sasgen-etal_2017/In_situ_GPS_uplift_rates.zip.

Bedrock uplift rates
Bedrock uplift rates derived for the REGINA project are available in the text file "REGINA_rates_full.txt", as presented in Tables A2 and A3 of the Appendix A4. The files "REGINA_rates_03-13.txt" and "REGINA_rates_03-09.txt" contain subsets of the data, with the temporal coverage limited to 2003-2013 and 2003-2009, respectively. The files are organized as follows: longitude ( • ), latitude ( • ), uplift rate (mm yr −1 ), uncertainty in the uplift rate (mm yr −1 )s and GPS site ID. These *.txt files are the input to the clustering script described below. No elastic correction has been applied.

Clustering script
In addition to the uplift rates for individual GPS sites, we provide a bash script "cluster.sh" for clustering the heterogeneous data according to their geographic locations, for a predefined threshold value. The idea is to reduce stochastic and geophysical noise of neighboring stations in order to obtain uplift rates that are better regional representations of the length scale recovered with GRACE (ca. 200 km). In an iterative procedure, the script selects neighboring sites within a threshold ranging from 10 to 220 km and calculates the weighted average of the uplift rates and a simple average of the station locations. Input to the script is the REGINA rate files, specified in the previous Sect. 3.2.1. Further details and the application to the GPS data set can be found in the second REGINA paper (Sasgen et al., 2017a). Note that the script relies on the open-source program suite Generic Mapping Tools (http://gmt.soest.hawaii.edu/; Wessel et al., 2013). Similar clustering can be achieved with the function kmeans in Matlab ® or its open-source alternative GNU Octave.

GPS time series
The GPS time series were created as part of the RATES project (UK NERC grant NE/I027401/1), not solely the REGINA study. The data can be obtained upon request from coauthor Elizabeth Petrie.

Gravimetry data analysis
We investigate the Release 5 (RL05) GRACE coefficients of the Centre for Space Research (CSR; Bettadpur, 2012) and the German Research Centre for Geosciences (GFZ; Dahle, 2013), provided up to spherical harmonic degree and order j max = 96 and 90, respectively, in the Science Data System (SDS). For reasons of comparison, we adopt j max = 90 for both GRACE solutions. A temporal linear trend in the ocean bottom pressure variations modeled by the atmospheric and oceanic background models (GAD; Bettadpur, 2012) was re-added to the monthly solutions, according the GRACE Science and Data System recommendation (Dobslaw et al., 2013). The GRACE coefficients C 20 were replaced by estimates from satellite laser ranging (SLR) provided by Cheng et al. (2013). In our analysis we apply the cutoff degrees j max = 50, which has been commonly used, and j max = 90, which is considered experimental in terms of the remaining signal content.
The determination of the rate of the gravity field change over Antarctica follows the scheme sketched in Fig. 3. The rate of the gravity field change, expressed as equivalent water height variations, is estimated in the spatial domain by adjusting a six-parameter function consisting of a constant, a temporal linear trend, and annual and semiannual harmonic amplitudes. A quadratic term was not co-estimated due to the project's focus on the rates (i.e., temporal linear trends). It should be stated that including a quadratic term would slightly reduce the residual uncertainties, particularly in the Amundsen Sea sector, where an ice-dynamic acceleration of mass balance rates occurs that is not accounted for by nonlinear SMB variations in the ice sheet (see Sect. 4.2).
The post-processing of the GRACE coefficients follows three main steps below.
Step 1: Optimization of de-striping filter Due to effects like the propagation of measurement noise and temporal aliasing, a large proportion of the variations contained in the monthly solutions is related to noise. The noise of the monthly solutions is lowest close to the pole and exhibits a characteristic north-south-oriented stripe pattern. This is visible in the gravity field rate and the propagated root-mean-square (rms) uncertainties shown in Fig. 3. In order to remove the stripe pattern, we apply the de-correlation filter of Swenson and Wahr (2006) (hereinafter, "Swenson filter") specifically tuned to optimize the recovery of the gravity field rate over the region of Antarctica, which is detailed in Sect. 4.1. Figure 3 shows that the de-striping procedure reduces the rms uncertainty in the rate by approximately 1 order of magnitude.
Step 2: Reduction in nonlinear mass variations For isolating gravity field rates, the second step in the processing is the reduction in de-trended variations in the surface mass balance, caused by accumulation events. The data set used for this purposes is the RACMO2/ANT (Lenaerts et al., 2012) converted into monthly sets of spherical harmonic coefficients. The reduction in these accumulation variations does not change the temporal linear trend, but it reduces rms uncertainties especially in coastal regions (Fig. 3). Details are provided in Sect. 4.2. Table 4. Uplift rates y u and associated uncertainties σ u (mm yr −1 ) for selected sites where uplift rates are manually evaluated based on the spread of rates obtained by subsampling the time series ("rman" method), compared to data published by Thomas et al. (2011), Argus et al. (2014) and Wolstencroft et al. (2015). See also "rman" sites in Table Appendix A4, Table A2 Step 3: Month-dependent weighting The performance of the GRACE satellite system was weaker in the early mission phase due to issues with the star cameras of the satellites (Christoph Dahle, GFZ German Research Centre for Geosciences, personal communication, 2015; Fig. 5). A rate estimate with uniform weighting of all months does not account for these variations. Therefore, in the last step, month-dependent uncertainties are estimated and applied as weights during the linear regression of the temporal linear trend. This slightly changes both the resulting rate estimate, as well as its rms uncertainties. Details are provided in Sect. 4.3. Finally, after post-processing and evaluation of the gravity field rate (Sect. 4.4), we select the GRACE release and cutoff degree providing the lowest uncertainty level (Sect. 4.5) as a reference input for our joint inversion for present-day icemass change detailed in the second REGINA paper (Sasgen et al., 2017a).

Optimization of de-striping filter
The Swenson filter has been proven to effectively reduce the typical north-south-correlated error structures of GRACE monthly solutions. The filter is based on the observation that these structures correspond to correlated patterns in the spherical harmonic domain, namely correlations within the coefficients of the same order and even or odd degree (Swenson and Wahr, 2006). The standard way of fitting and removing these patterns is by adjusting polynomials to the respective sequences of spherical harmonic coefficients, independently for individual months. Parameters to choose are the degree of the polynomial n pol and the minimum order m start , starting from which this procedure is applied. In principle, a higher-degree polynomial reduces the variability in coefficients of even/odd degree and results in stronger filtering also at a lower minimum order -however, the behavior of the filter may differ for regional applications, as discussed below. Note that the tuning of other parameters has been presented previously, e.g., the window width (Duan et al., 2009) or the degree range to which the filter is applied. Chambers and Bonin (2012) have assessed these parameter options with regard to the new GRACE RL05 solutions and global oceanic signals. Here, we perform a detailed analysis of the choice of the Swenson filter parameters in order to optimize the signalto-noise characteristics of the rate of the gravity field change over Antarctica.
We assess signal corruption by applying the filter to a synthetic test signal, which is based on high-resolution elevation rates from satellite altimetry and reflects the prevailing signatures of present-day ice change with sufficient realism. For each choice of filter parameters, the signal corruption is assessed as the rms difference between the original and the filtered synthetic signal (rms signal ). The rms is evaluated in terms of water-equivalent height per year for the signal components within the region south of 60 • S latitude.
For assessing the noise and noise reduction in the filtered fields, we face the task of separating the noise from the geophysical signals in the gravity field rates derived from GRACE. Here we attempt such a separation by reducing a priori information on the rate of ice-mass change from the GRACE fields and considering the residual as an upper bound representation of noise. The a priori information is, again, based on elevation rates. For the noise assessment we then take the rms of the residual rates in terms of water-equivalent height per year, rms noise , again for the region south of 60 • S latitude. Since the residual gravity field rates may still contain some geophysical signal, we consider this noise estimate as an upper bound for the true GRACE uncertainties. It should be stated that, after the Swenson filtering, an additional Gaussian filtering is applied to the sig- Post-processing steps applied to the GRACE gravity fields; shown is the impact on the gravity field rate y g (left) and the associated rms uncertainty σ g (right). Small maps show change in the gravity field rate between two subsequent steps. Color scale is mm w.e. yr −1 . GRACE data are GFZ RL05a. nal and noise models with a 200 km filter width, which was determined to be the optimal smoothing half width for the signal-to-noise ratio in the GRACE spectra by Wiener optimal filtering (Sasgen et al., 2006) as reflected in the degreeamplitude spectrum. Figure 4 shows the assessed signal corruption and noise reduction as a function of the two Swenson filter parameter choices, the polynomial degree n pol and the minimum order m start . The results are shown for the gravity field expanded to degree and order j max = 90 of the GFZ RL05a coefficients, even though using j max = 50 and CSR RL05 yields similar results. As expected, the signal corruption (rms signal ) increases with increasing strength of the Swenson filter, that is with increasing n pol and the decreasing minimum order m start . In terms of noise reduction, we see as expected that stronger filtering (increasing n pol ; decreasing m start ) decreases the rms noise (Fig. 4). However, we find that for filter parameters m start < 10, this pattern is reversed, and rms noise increases again. A closer analysis indicates that the consideration of the low orders in the Swenson filtering transfers energy (both from signal and noise) from low latitudes to midlatitudes to the polar regions, which leads to a considerable signal corruption over the region of interest. We avoid this degradation by limiting the range of filter parameters in the subsequent optimization of the gravity trends to m start ≥ 10.
To define the optimal filter parameters a quadratic sum of the signal corruption and noise reduction is computed, allowing us to balance both effects. The optimal values are m start = 12 and n pol = 7, as indicated in Fig. 4c. These filter parameters are subsequently used. For comparison it is stated that Chambers and Bonin (2012) found m start = 15 and n pol = 4 to be optimal for oceanic applications. Note that the signal corruption is assessed only to optimize the destriping filter. Possible signal degradation due to de-striping is not included in the uncertainty estimate of the optimally filtered GRACE trends. However, signal loss due to the additional smoothing with a 200 km Gaussian filter is accounted for by applying the same filter to the viscoelastic response functions, as well as the altimetry-based input fields (Appendix A5).

Reduction in nonlinear mass variations
The temporal variations in the Antarctic gravity field show a strong year-to-year fluctuation, apart from the linear trend (Wouters et al., 2014). A large portion of the nonlinear signal in geodetic mass and volume time series is well explained by modeled SMB fluctuations (Sasgen et al., 2010;Horwath et al., 2012). Towards the ultimate goal of isolating the linear GIA signal from time series of mass change, we removed nonlinear effects of modeled SMB variations from the GRACE time series; for this we calculate the monthly cumulative SMB anomalies with respect to the time period 1979 to 2012 obtained from RACMO2/ANT (Lenaerts et al., 2012).
We then transfer the monthly cumulative SMB anomalies in terms of their water-equivalent height change into the spherical harmonic domain and subtract them from the monthly GRACE coefficients. In principle, the reduction in the SMB variations from the GRACE time interval has two effects: first, it may change the overall gravity field rate derived from GRACE, depending on the assumption of the SMB reference period. Ideally, the reference period reflects a state of the ice sheet in which input by SMB equals the outflow by ice discharge, and SMB anomalies estimated for today reflect the SMB component of the mass imbalance. However, any bias in the SMB in the reference period leads to an artificial trend in the ice sheet mass balance attributed to SMB. This is an undesired effect, and to avoid it, we de-trend the cumulative SMB time series for the time interval coeval  to the GRACE analysis (February 2003 to October 2009) before subtracting it from the GRACE gravity fields, yielding zero difference in the gravity field rates before and after processing Step 2 (Fig. 3). The second effect is the reduction in the post-fit rms residual for this known temporal signal variation. After reducing the SMB variations, the propagated rms uncertainty in the derived gravity field rate becomes closer 504 I. Sasgen et al.: Altimetry, gravimetry, GPS and viscoelastic modeling data to the uncertainty level of the GRACE monthly solutions (Fig. 3).

Month-dependent weighting
The quality of GRACE monthly solutions changes with time, for example due to changing orbital sampling patterns (Swenson and Wahr, 2006). Figure 5 shows the temporal evolution of rms uncertainties in the monthly GRACE gravity fields in the Antarctic region; shown are residual mass anomalies, integrated over Antarctica, after the grid-based removal of the temporal linear trend and annual oscillation components and after the application of the filtering described in Step 1 and the removal of the SMB fluctuations in Step 2. Note than an annual oscillation component is included to remove possible seasonal fluctuations in SMB not captured by the regional climate model. However, omitting the annual oscillation component yields similar results. The residual monthly mass anomalies are attributed to noise and are used as weights in our least-squares linear regression, applied as Step 3 of the GRACE processing. Figure 5 shows that these uncertainties are higher during early 2003. Applying the monthly dependent weighting has the effect of reducing the influence of the first months of the year 2003 on the estimated gravity field rate, which is similar to shortening the time series, given the relatively large uncertainties. As expected, the post-fit rms uncertainty associated with the rate reduces if the early months of the year 2003 are excluded. Altogether, the month-dependent weighting reduces the magnitude of stripe patterns characteristic of the uncertainty in GRACE monthly solutions and yields a more realistic estimate of the uncertainty associated with the gravity field rates (Fig. 3). Figure 6 shows the estimated rms uncertainty in the gravity field rate over Antarctica, after post-processing. It is evident that the largest uncertainties are located in a ring south of −80 • S latitude. This is explained by the design of the Swenson filter; little or no noise reduction is achieved close to the poles, as the gravity field is represented by near-zonal coefficients, which pass the filter mostly unchanged (m start = 12). Extending the kernel of the Swenson filter to these near-zonal coefficients (m start ≤ 10) creates high signal corruption and is not suitable for the optimal rate estimate over Antarctica (see Sect. 4.1). Larger uncertainties are also estimated for the Ronne and Ross ice shelf areas, which are most likely a consequence of incomplete removal of the ocean tide signal during the GRACE de-aliasing procedure (Dobslaw et al., 2013). It should also be stated that the rms uncertainty estimate does not include possible systematic errors in the GRACE solutions, e.g., due to long-term drift behavior of the observing system.

Selection of GRACE release
Our evaluation of the monthly GRACE uncertainties (Fig. 5), as well as the propagated rms uncertainty in the temporal linear trend (Fig. 6) indicates that the lowest noise level for the Antarctic gravity field rate (February 2003 to October 2009) is currently achieved with GRACE coefficients of CSR RL05, expanded to j max = 50. We therefore refrain from including coefficients with j max > 50 in order not to compromise the rate estimates by unnecessarily increasing the noise level (see Appendix A5, Fig. A3). We adopt CSR RL05 with j max = 50 as our preferred solutions for the representation of the gravity field rates over Antarctica, even though GFZ RL05 with j max = 50 yields very similar rates (Fig. 6). This choice is supported by the joint inversion, as CSR RL05 with j max = 50 provides the highest level of consistency (lowest residual misfit) with the altimetry and GPS data sets (see the second REGINA paper, Sasgen et al., 2017a, Supplement Sect. S3), which we interpret to indicate a minimum of spurious signals in the trends. To account for the uncertainty related to our choice of the solution, we consider not only rms uncertainties in the GRACE rates but also solution differences, in the uncertainty in the final GIA estimate (Fig. 6). The solution differences represent the absolute deviation between trends from GFZ RL05 and CSR RL05 (February 2003 to October 2009; cutoff degree j max = 50). These are then summed up and squared with the propagated rms uncertainties. It is acknowledged that the solution differences contain systematic noise arising from the GRACE processing; the pattern and magnitude may change over time. However, they provide a measure of how much the results will change if a GRACE release alternative to CSR RL05 is considered. The difference between GRACE rates filtered with Gaussian smoothing of 200 km and the optimized Swenson filter together with Gaussian smoothing of 200 km is shown in the Appendix A5, Fig. A4

Stokes coefficients of gravity field change
The monthly GRACE gravity field solutions from the data system centers GFZ and CSR are available atftp://podaac.jpl. nasa.gov/allData/grace/L2/ or http://isdc.gfz-potsdam.de/ as spherical harmonic (SH) expansion coefficients of the gravitation potential (Stokes ccoefficients). More information is available in Bettadpur (2012). The data archive contains temporal linear trends in the fully normalized Stokes coefficients in the "geodetic norm" (Heiskanen and Moritz, 1967), complete to degree and order 90, inferred from these time series according to Sect. 4. We provide data for GFZ RL05 and CSR RL05 for the time period 2003-2009 and 2003-2013 and for various combinations of filtering. The coefficients are organized as Degree j , Order m, c_jm, s_jm.

Code for de-striping filtering
The Matlab ® function "KFF_filt" performs decorrelation filtering for sets of spherical harmonic coefficients, typically from GRACE gravity field solutions, following the idea of Swenson and Wahr (2006). An open-source alternative to Matlab ® is GNU Octave: https://www.gnu.org/ software/octave/. The function is called KFF_filt = swen-son_filter_2(KFF, ord_min, deg_poly, factorvec, maxdeg), where variables ord_min and deg_poly equal m start and n pol , respectively, in Sect. 4. KFF contains the sets of spherical harmonic coefficients in the "triangular" format (not memory-efficient but intuitive). For example, for a set of coefficients with maximum degree j max = 3 and maximum order m max = 3, the set of coefficients is stored in a j max × m max matrix in the following way: % 00s_11c_10c_1100; % 0s_22s_21c_20c_21c_220; % s_33s_32s_31c_30c_31c_32c_33]

Viscoelastic modeling
The Earth structure of Antarctica is characterized by a strong dichotomy between east and west, separated along the Transantarctic Mountains (e.g., Morelli and Danesi, 2004). Recent seismic studies have produced refined maps of crustal thicknesses also showing slower upper-mantle seismic velocities in West Antarctica, indicating a thin elastic lithosphere and reduced mantle viscosity Heeszel et al., 2016). Moreover, yield strength envelopes of the Earth's crust and mantle suggest the possibility of a viscously deforming, ductile layer (DL) in the lower part of the crustal lithosphere (Ranalli and Murphy, 1987), a few tens of kilometers thick and with viscosities as low as 10 17 Pa s (Schotman et al., 2008). High geothermal heat flux is in agreement with the seismic inferences of a thin elastic lithosphere and low mantle viscosity and would favor the presence of such a DL also in West Antarctica (Shapiro and Ritzwoller, 2004;Schroeder et al., 2014). The choice of the viscoelastic modeling approach used to determine load-induced surface displacements and gravitational perturbations is governed by three main requirements: (i) to accommodate lateral variations in Earth viscosity, (ii) to allow for Earth structures with a thin elastic lithosphere and low viscosity layers, in particular including a DL, and (iii) to provide viscoelastic response functions for the joint inversion of the satellite data described in the second REGINA paper (Sasgen et al., 2017a). With regard to point (iii), it should be mentioned that the viscoelastic response functions provide a geophysically meaningful way to relate surface displacement and gravity field changes, considering also dynamic density changes within the Earth's interior . Moreover, it allows us to consider the changes in the ratio of surface-displacement and gravity field changes caused by the Earth structure, in particular, the lithosphere thickness. Another advantage is that different filtering can be applied to the viscoelastic response functions in order to match the filtering of the input data set, avoiding the introduction related biases (Appendix A5).
To meet these requirements, we adopt the time-domain approach (Martinec, 2000) for calculating viscoelastic response functions of a Maxwell continuum to the forcing exerted by normalized disc loads of constant radius. Then, the magnitudes and spatial distribution of the surface loads are adjusted according to the satellite data to obtain the full GIA signal for Antarctica. The forward modeling of viscoelastic response functions is a classic topic in solid Earth modeling (e.g., Peltier and Andrews, 1976); however, their application to inverting multiple-satellite observations for present and past ice sheet mass changes is new and applicable to other regions, such as Greenland or Alaska.
The viscoelastic response function approach allows for high spatial resolution at a low computational cost in the numerical discretization of the Earth structure as well as in the representation of the load and the response. In addition, we can accommodate a high temporal resolution, which is required when considering low viscosities and associated relaxation times of only a few decades. The spherical harmonic cutoff degree for the simulations shown in the following is j max = 2048 (ca. 10 km).

Load model parameters
The load function σ (t, ϑ) is disc shaped with a constant radius of ca. 63 km. The radius of 63 km matches the mean radius of the discs south of 60 • S of the geodesic grid (here, the Icosahedral Nonhydrostatic (ICON) Model 1.2 grid, status 2007, e.g., Wan et al., 2013), which underlies the joint inversion of the altimetry, gravimetry and GPS observations (see the second REGINA paper, Sasgen et al., 2017a). The resolution of the geodesic grid is chosen to allow for an adequate representation of the load and viscoelastic response with regard to the input data sets while minimizing the computational cost. The disc load experiment consists of a linear increase in the ice thickness at a rate of 0.5 m yr −1 continuing until a new dynamic equilibrium state between load and response is reached. After the application of the constant loading rate, two extra time steps are done with no loading change to give the purely viscoelastic response. For West Antarctica, the loading rate is held constant for 2000 years, for East Antarctica it is 15 000 years, which are longer times than needed to reach dynamic equilibrium (see Appendix A8). With reference to the assumed ice density of 910 kg m −3 , this thickness increase corresponds to a mass gain of ca. 5.6 Gt yr −1 . Then, to obtain the signal component of the viscous Earth response only, the elastic response and the direct gravitational attraction of the load are subtracted. The experiment is designed as an increasing load, for example representative of the ceasing motion of the Kamb Ice Stream (Ice Stream C; Retzlaff and Bentley, 1993), West Antarctica. Due to the linearity of the viscoelastic field equations, it is not necessary to calculate separately the equivalent unloading experiment, −σ (tϑ), for example corresponding to the past and present glacier retreat of the Amundsen Sea sector, West Antarctica (Bentley et al., 2014 andRignot et al., 2014, respectively). Among others, the combined inversion of the altimetry, gravity and GPS data (the second REGINA paper, Sasgen et al., 2017a) solves for the magnitude and the sign of the load, allowing for ice advance as well as ice retreat.

Earth model parameters
We set up an ensemble of 58 simulations representing different parameterizations of the viscosity structure (Table 5), split into West Antarctica (56 simulations) and East Antarctica (2 simulations). The ensemble approximately covers the range of values of the viscosity and lithosphere thickness inferred from Priestley and McKenzie (2013). For West Antarctica, varied parameters are the lithosphere thickness, h L (30 to 90 km in steps of 10 km), the asthenosphere viscosity (1 × 10 18 to 3 × 10 19 Pa s in four steps) and the presence of a DL, with 10 18 Pa s. For East Antarctica, we employ parameter combinations appropriate for its cratonic origin with h L of 150 and 200 km and an asthenosphere viscosity equivalent to the upper-mantle viscosity of 5 × 10 20 Pa s. These values lie in the range of previously applied viscosity values in Antarctica (Nield et al., 2012;Whitehouse et al., 2012;van der Wal et al., 2015). For the radial layering of the elastic properties, we adopt the Preliminary Reference Earth Model (PREM; Dziewonski and Anderson, 1981). Later, in the joint inversion, the distribution of viscoelastic response functions is based on the Earth structure model of Priestley and McKenzie (2013). Priestley and McKenzie (2013) provide a global distribution of viscosity values up to a depth of 400 km, which is sampled at the location of the geodesic grid. We then define a threshold value for the viscosity (here, 10 22 Pa s) above which the Earth response is considered purely elastic and infer the associated thickness of the elastic lithosphere. The impact on the final joint inversion estimate of changing the threshold value of 10 22 Pa s is presented in the second REGINA paper (Fig. S4 in Sasgen et al., 2017a). Note that the Earth response in the equilibrium state only depends on the lithosphere thickness (independent of viscosity), which is therefore considered as the main Earth model parameter in the joint inversion. Further details are presented in the second REGINA paper, Sasgen et al. (2017a).

Gravity and displacement rate response functions
The calculated response functions for surface deformation (radial displacement) and gravity (geoid-height change) are discretized along 1507 latitudinal points within the range 0 ≤ ϑ ≤ 90. Simulations are typically run over 2 kyr with a temporal resolution of t = 10 years (plus two time steps with constant load thickness). For East Antarctic parameterizations, the simulation period was extended to 20 kyr due to the higher upper-mantle viscosities and associated slower relaxation. However, note that the ratio of geoid-height change versus radial displacement falls off to a factor of 1/e 2 relative to the initial value after ca. 2 kyr of simulation (Appendix A6, Fig. A5). The forcing expected in central East Antarctica is an increase in accumulation towards presentday conditions after ca. 7 ka BP (van Ommen et al., 2004), also justifying the use of equilibrium kernels for East Antarctica. The time derivatives of the radial displacement y u and of the geoid-height change y g are calculated with a central difference scheme.
Examples of response functions to the loading detailed in Sect. 5.1 for the rate of radial displacement, y u , and rate of geoid-height change, y g , are shown in Figs. 7 and 8, respectively. Instantaneously, the increasing load,σ (t) = const., induces an elastic response that is characterized by subsidence and an increase in the direct gravitational potential (dashed lines in Figs. 7 and 8, respectively). This is the elastic response function adopted in the joint inversion. Note that the elastic response function will not differ between East and West Antarctica, as it is entirely based on the distribution of densities and elastic parameters provided by the PREM. As the load buildup continues, the instantaneous response is followed by the viscoelastic response, which depends in timing and magnitude on the underlying lithosphere and viscosity structure, further increasing the displacement rates, y u (blue to red lines in Fig. 7). The compensation by solid Earth deformation is reflected in the decreasing geoid rate, y g (Fig. 8). After a certain time, which depends on the value of the asthenosphere viscosity, a new dynamic equilibrium state is reached at which y u and y g do not change in time any more. In the last two time steps, the load is kept constant (y g (t) =) and the responses in y u and y g are only caused by the relaxation of the Earth's viscoelastic deformation (solid black line in Figs. 7 and 8), which is the viscoelastic response function adopted in the joint inversion.

Discussion of effects of selected earth model
parameterizations on GIA response Figure 9 shows the response of y u for four end-member sets of Earth model parameters with a thick lithosphere, weak asthenosphere (TkWk: h L = 90 km; η AS = 1 × 10 18 Pa s); a thick lithosphere, strong asthenosphere (TkSg: h L = 90 km; η AS = 3 × 10 19 Pa s); a thin lithosphere, weak asthenosphere (TnWk: h L = 30 km; η AS = 1 × 10 18 Pa s); and a thick lithosphere, strong asthenosphere (TnSg: h L = 30 km; η AS = 3 × 10 19 Pa s), without a ductile layer, DL. In this context, thick/thin and strong/weak refer to values in comparison to the "average" value of the ensemble for West Antarctica; an elastic lithosphere of thickness 90 km (here, "Tk") is in the range of the global average continental lithosphere usually applied in GIA studies (e.g., Peltier, 2004) or that of East Antarctica (150 to 200 km). Figure 10 shows the response in y u for the same end-member set of Earth model parameters with a DL included. It should be stated that the Earth structure with h L = 30 km and a DL is considered very extreme because in this case the ductile layer extends down to the asthenosphere and an elastic mantle lithosphere is missing. Figures 9 and 10 show that for the weak asthenosphere (η AS = 1 × 10 18 Pa s), viscoelastic deformation is visible already after 1 decade of loading (or unloading), leading to considerably larger subsidence rates compared to the purely elastic case even on very short timescales. For these Earth model parameters, a new dynamic equilibrium state is achieved within a few centuries. The rates of subsidence in this equilibrium then primarily depend on the support provided by the flexure of the elastic lithosphere.
For the extreme TnWk case, equilibrium rates of −45 mm yr −1 are achieved at the load center, and considerable subsidence of −20 mm yr −1 already occurs after 10 years of loading (Fig. 9). Increase in asthenosphere viscosity (TnSg case) reduces the viscous material transport and leads to a slower adjustment towards the dynamic equilibrium state, which takes more than 1 kyr. It should be stated that in our definition of the ensemble parameters, reducing the lithosphere thickness in turn increases the thickness of the asthenosphere (bottom depth of asthenosphere is fixed), which facilitates lateral material transport inside the asthenosphere.
The consideration of the DL in the Earth structure causes a thinning of the effective elastic lithosphere. As a consequence, greater and more localized subsidence rates are produced for all sets of parameters (Fig. 10). Interestingly, in the case of a thick elastic lithosphere (90 km), the radial displacement exhibits a local minimum at around 120 and 160 km distance from the load center (Fig. 10), which is a consequence of the viscous material transport inside the ductile layer. The maximum equilibrium rate of −76 mm yr −1 is achieved for the ThWk case with DL, where the viscous deformation leads to rates of −25 mm yr −1 already after 10 years of loading.

Assumptions and limitations
Although the approach of modeling response functions to axisymmetric disc loads and subsequently superposing them is very efficient in terms of the computational cost, this simplification introduces some limitations. First, the superposition of response functions representing different Earth struc-tures neglects the transmission of stresses between these regions -a problem that can only be resolved with fully threedimensional solid Earth modeling (e.g., van der Wal et al., 2015). The largest impact for the displacement rates is expected in regions with lateral contrasts in lithosphere thickness and mantle viscosity such as the Transantarctic Mountains. Second, the constant disc radius of about 63 km implies that finer-scale deformation cannot be resolved. Although this resolution is adequate for interpreting GRACE data (spatial half wavelength of ca. 200 km) smaller-scale loading excitement may be necessary for interpreting local GPS measurements near the loading, particularly for the elastic response to present-day glacial changes. Furthermore, the viscoelastic response functions describe the Earth response in an equilibrium state for a constant rate of load change; if the load exhibits more complex temporal variations, this assumption is violated. Finally, it is assumed that the lithosphere thickness, upper-and lower-mantle viscosities are approximately known.

Data availability
The viscoelastic response functions and related ancillary data are directly accessible in the PANGAEA repository: http://hs.pangaea.de/model/Sasgen-etal_2017/Viscoelastic_ response_functions.zip.

Viscoelastic kernels
Output files contain 1507 latitudinal points (0≤ ϑ ≤ 90) covering a region greater than the size of the Antarctic domain, as well 203 time steps of West Antarctica (213 time steps for East Antarctica because of extending the simulation period to 15 kyr). The time derivative of the radial displacement, u, is calculated with a central difference scheme (y u = [u(t + t/2) − u(t − t/2)]/ t). The difference between two time steps is t = 10 years. The same applies to the rate of geoid-height change, y g . Note that the load is constant during the last two time steps (no rate of change); therefore, the kernels represent the viscoelastic relaxations only, without the instantaneous elastic deformation or the direct gravitational attraction of the load.
The results are stored independently for the rheology of East and West Antarctica, the latter with and without a ductile layer in the elastic part of the lithosphere. The data are stored in a Matlab ® file format, which is also readable with GNU Octave https://www.gnu.org/software/octave/ .
The response kernels for East Antarctica are organized in an analogous way, [HL, LAT, TIME, VAR]. HL is the lithosphere thickness at 150 and 200 km (two entries) Note that the asthenosphere and upper mantle viscosity is constant at 5 × 10 20 Pa s and therefore has no entry.
The spectral resolution underlying these fields is spherical harmonic cutoff degree 2048. The user should apply an adequate smoothing filter when using for inverting GRACE gravity fields. Filtered kernels are available upon request by the author.

Geodesic grid
The computation of the geodesic grid is not an original contribution of the authors but based on the grid generator of the ICON GCM project (http://icon-downloads.zmaw.de/). For completeness, we provide the data set with disc locations based. An alternative resource for downloading geodesic grids at different resolutions in netCDF format can be found here: http://kiwi.atmos.colostate.edu/BUGS/geodesic/.
The files format is as follows: vert-7.mask.cont_and_shelf.re.dat: longitude ( • ), latitude ( • ) vert-7.mask.cont_and_shelf.re.proj.dat: The thickness of the elastic lithosphere at the locations of the geodesic grid for different values of the viscosity threshold is applied to the data set of Priestley and McKenzie (2013).

Open-source code for viscoelastic modeling
The open-source software package SELEN allows the computation of the Maxwell viscoelastic Earth response to userdefined ice sheet evolutions, in particular also a simplified disc load forcing as presented in this paper. The program is downloadable at: https://geodynamics.org/cig/software/ selen/.

Data availability
The altimetry, gravimetry, GPS and viscoelastic modeling data used in this project are available at https://doi.org/10.1594/PANGAEA.875745 in the www.pangea.de archive. The data description and user documentation are given for each data type within the respective subsection of this paper (Sects. 2 to 5).

Conclusions
In this paper, we have presented refined temporal linear trends in surface elevation, gravity field change and bedrock displacement based on Envisat-ICESat (2003), GRACE (2003) and GPS (1995-2013, respectively. In addition, we have performed forward modeling of the viscoelastic response of the solid Earth to a disc load forcing. These response functions are particularly suited to represent the distinct geological regimes of East and West Antarctica in the joint inversion of multiple satellite data. Similarly, the functions can be applied to the other geographical regions as well. The data and code necessary to reproduce our results, or apply our approach to a different problem, are provided at www.pangea.de (https://doi.org/10.1594/PANGAEA.875745; Sasgen et al., 2017b; follow link "View dataset as HTML").
We have refined surface-elevation rates for the Antarctic ice sheet for the time interval 2003-2009 by combining Envisat and ICESat altimetry data. The straightforward approach performs a grid-based comparison of the noise in the elevation rates obtained from Envisat and ICESat. For large parts of the ice sheet, the elevation rate is based on ICE-Sat data, particularly for the rough terrain along the coast as well as close to the South Pole (polar gap of Envisat). Envisat contributes in some low-relief areas in East Antarctica and along the Antarctic Peninsula, as well as along single spurious ICESat tracks. Thus, the composite elevation rates are maximized in terms of spatial coverage and minimized in terms of the uncertainties.
The GPS processing carried out as part of the RATES and REGINA projects has produced a comprehensive data set of 118 Antarctic GPS records, which, for continuous sites, span a longer time interval (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) than those of previous studies (Thomas et al., 2011(Thomas et al., : 1995(Thomas et al., -2011Argus et al., 2014Argus et al., : 1994Argus et al., -2012Martín-Español et al., 2016b. The ensemble processing done for the REGINA project has allowed us to assess the contribution of systematic error sources. In addition, for sites where there is potential doubt about the quality of the metadata or the behavior of the site, we have adopted a "conservative but realistic" approach to assigning new confidence limits. The screening of GPS data for outliers involved careful manual assessment, encompassing the review of measurement logs and notes on problems in the field. The data quality is reflected in the uncertainty estimates for the GPS rates, which therefore represents more reliable input data than GPS rates based on processing without manual intervention. Note, however, that SMB variations might also contribute to the GPS uplift rates given that the time spans of these data vary. We have optimized the post-processing sequence for estimating the temporal linear trend and its uncertainty in the GRACE gravity field solutions for the region of Antarctica. In particular, we have derived optimal parameters for destriping the monthly gravity fields over Antarctica according to Swenson and Wahr (2006). In addition, we have removed de-trended SMB fluctuations from the GRACE time series, to obtain a more representative uncertainty estimate based on the post-fit rms residual. We have included month-dependent weighting in the least-squares estimate of the gravity field rates to account for the varying quality of the monthly GRACE solutions. The optimization of the de-correlation filter of Swenson and Wahr (2006) to the signals expected in Antarctica reduced the residual uncertainty and improved the reliability of inferred mass anomalies.
With the aim of joining the multiple satellite data using the knowledge of the geophysical processes involved, we have calculated elastic and viscoelastic response functions of the solid Earth. The viscoelastic response functions represent the gravity field change and surface displacement to a disc load forcing for a variety of Earth model parameters; particularly, values of mantle viscosity and lithosphere thickness strongly varying between the distinct geological regimes of West and East Antarctica.
In particular, we have investigated the effect of a ductile layer in the crustal lithosphere on the viscoelastic re-bound signature. We show that for moderate load changes of 0.45 m yr −1 water equivalent (here, applied as disc load with a radius of ca. 63 km), uplift rates reach the cm yr −1 level within decades assuming asthenosphere viscosities < 10 19 Pa s and a lithosphere thickness < 50 km; both as plausible values for parts of West Antarctica. Including a ductile layer in the crustal lithosphere further attenuates the uplift rates and localizes the deformational response. This suggests that GIA in West Antarctica may locally be a result of more recent, centennial load changes, most notably in the Amundsen Sea Embayment and in part of the Antarctic Peninsula (Nield et al., 2012). Similar conclusions were reached by Ivins and James (2005) and Nield et al. (2014), even though it is not possible to constrain the exact timing of the load from our approach.
The advantage of the viscoelastic response kernels is that a meaningful ratio of the rate of the gravity disturbance versus the rate of the surface displacement is calculated for each choice of the Earth model parameters, avoiding the approximation with an average rock density (e.g., Riva et al., 2009;Gunter et al., 2014). Using the response functions allows us to reconcile GIA signatures with measurements of large bedrock uplift and small gravity field increase in the Amundsen Sea Embayment, associated with weak Earth structures. Clearly, the response functions adopted here represent only the viscoelastic equilibrium state and, thus, are considered only an intermediate step to full dynamic modeling of the GIA response. Nevertheless, this approximation represents a significant improvement of other joint inversion methods, as it bases the joint inversion on physically meaningful response kernels. With extra data on the past ice evolution, such as paleo-thickness rates, our approach can be expanded to address the temporal evolution as well.
In the second REGINA paper (Sasgen et al., 2017a), we perform the joint inversion for present-day ice-mass changes and GIA in Antarctica, based on the input data sets and viscoelastic response functions presented here. We validate our results using forward-modeling results and other empirical models and show the impact on CryoSat-2 volume and GRACE mass balances, respectively. Note, however, that the post-processing methods and viscoelastic functions presented here are applicable also to other geographical regions with superimposed present-day mass change and GIA signatures.

A2 Firn compaction and SMB corrections
We apply rates of firn compaction, h comp , using the output of the firn compaction model provided by Ligtenberg et al. (2011), which is driven by RACMO2/ANT (Lenaerts, 2010). However, we do not apply a correction for anomalies in the SMB, δh SMB , as, e.g., undertaken by Gunter et al. (2014)  A4 Uplift rates at all GPS site used in this study Table A2. GPS uplift rates for this study. The columns are as follows: site name, estimated uplift rate σ u (mm yr −1 ), estimated uncertainty σ u (mm yr −1 ), rate method, uncertainty method, approximate latitude (dec. degrees) and approximate longitude (dec. degrees). Methods are as follows: values estimated by the CATS noise analysis software ("cats"), propagated median uncertainty from CATS sites ("prop"), manual intervention in rate due to potential systematic uncertainties ("rman") and manual intervention in uncertainties due to potential systematic uncertainties.
Site name y u σ u Method y u Method σ u Lat (    In this study, we identify GRACE coefficients of CSR RL05 up to degree and order 50 appropriate to yield the most robust gravity field rates over Antarctica. Figure A3 provides another indication based on the degree-power spectrum of the geoid rates. It is visible that GFZ RL05 and CSR RL05 are very similar up to degree and order 50, where the power spectra show minima. For higher degrees, however, the power of the gravity field recovered with GRACE increases due to increasing noise; for the unfiltered coefficients, this increase is faster for GFZ RL05 than for CSR RL05. Figure A3. Degree-amplitude spectrum of the rate of geoid-height change (mm yr −1 ) for unfiltered (diamond-dashed lines) and for Swensonfiltered (solid lines with circles) solutions. Red: GFZ; green: CSR; black combination of GFZ and CSR with equal weights. Figure A4. Spatial rate of geoid-height change (a) and rate of equivalent water-height change (b) (mm yr −1 ) for the difference between the GRACE trends processed by a Gaussian smoothing of 200 km and the optimal Swenson filter and Gaussian smoothing. The solutions are CSR RL05; the spherical harmonic cutoff degree is 50.
The filtering of the GRACE gravity fields was optimized for reducing noise over Antarctica. The effect on the rms uncertainties is shown in Fig. 3. Additionally, Fig. A4 presents the difference between the GRACE rates filtered only with a Gaussian smoothing filter of 200 km and with the optimized Swenson filter. It is visible that the differences in the rate of geoid-height change and the associated rate of equivalent water-height change show a stripe-like noise pattern. This suggests that the de-striping is superior to conventional Gaussian smoothing, even at high latitudes, where GRACE ground-track spacing is very dense. It is also important to note that the filter does not introduce any magnitude bias or change the spectral content of the gravity field rates, which is important when applying only a Gaussian smoothing of 200 km (without Swenson filtering) to the altimetry data set and response kernels.

A6 Evaluation of the assumption of a viscoelastic equilibrium state
The viscoelastic response kernels employed (Sect. 5) describe the viscoelastic equilibrium state for the forcing with a disc load of constant radius and constant rate of mass increase (likewise mass loss). We neglect transitional changes in the solid Earth for load changes that have not reached the equilibrium state in terms of geoid-height change and surface displacement. Although the deformation and gravity signature in equilibrium eventually only depends on the lithosphere thickness, the time to reach the equilibrium is controlled by the viscosity parameters chosen. Figure A5 shows the evolution of the standardized ratio of the geoid-height change vs. surface displacement over time, calculated as r = [r (t) − r (t = t tmax )] /max [r (t) − r (t = t tmax )]), where r = y g (t)/y u (t) is evaluated at the load center. It is visible that for the weaker West Antarctic rheology (asthenosphere viscosity between 1 × 10 18 and 3 × 10 19 Pa s), r falls to 1/e 3 within 500 years. For East Antarctica (1 × 10 20 Pa s), r = e −2 is reached within 2 kyr. With this quasi-stationary solution approach, the inference about the timing of the past ice-mass change is limited to an upper limit in terms of magnitude and a lower limit in terms of load duration; a similar ratio is achieved by a thinner lithosphere thickness, which has not reached viscoelastic equilibrium state, and earlier load changes are fully relaxed.

A7 The assessment of SMB fluctuations in GPS uplift rates
We assess the impact of SMB fluctuations on the uplift rate at the GPS station locations using the modeled SMB of RACMO2 for the years 1979-2010. We compute the elastic deformation related to cumulative monthly SMB, de-trended for the entire simulation period 1979-2010. We then estimate the temporal linear trends at the GPS station locations for a moving window of varying width from 3 to 16 years. Figure A5. Standardized ratio of the rate of geoid-height change versus the rate of radial displacement for different values of the asthenosphere viscosity. Note that the ratio is calculated at the load center. Figure A6. Standard deviation (2σ ) of the uplift rates caused by accumulation variability for different GPS stations and time periods. It is visible that the uncertainty decreases with record length; for most regions, trend uncertainties are below 0.4 mm yr −1 for the actual GPS record length.
Then, for each window width, we estimate the standard deviation of the apparent trend induced by SMB for selected stations (Fig. A7). Typically, the uncertainty in uplift rate due to SMB variability is below 0.4 mm yr −1 for the actual GPS record length. An exception is PALM, which is located on the Antarctica Peninsula -a region with annual accumulation of up to 4 m yr −1 equivalent water height. Here, even after 12 years of measurements, GPS uplift rates are likely to contain accumulation signals of 4 mm yr −1 . A similar effect of the SMB fluctuations is expected at vesl. The load increases, with a fixed radius, at a constant rate of ca. 5.6 Gt yr −1 until an approximate equilibrium state is reached: 2 kyr for West Antarctica and 15 kyr for East Antarctica (Fig. A7). Then the load is applied without a change to obtain the purely viscoelastic response of the Earth model, i.e., without direct gravitational attraction of the load and the instantaneous elastic response. The associated Earth response constitutes the viscoelastic response functions adopted in the joint inversion. Figure A7. Load function applied to obtain the viscoelastic response functions.