The influence of Antarctic ice loss on polar motion: an assessment based on GRACE and multi-mission satellite altimetry

Increasing ice loss of the Antarctic Ice Sheet (AIS) due to global climate change affects the orientation of the Earth’s spin axis with respect to an Earth-fixed reference system (polar motion). Here the contribution of the decreasing AIS to the excitation of polar motion is quantified from precise time variable gravity field observations of the Gravity Recovery and Climate Experiment (GRACE) and from measurements of the changing ice sheet elevation from altimeter satellites. While the GRACE gravity field models need to be reduced by noise and leakage effects from neighboring subsystems, the ice volume changes observed by satellite altimetry have to be converted into ice mass changes. In this study we investigate how much individual gravimetry and altimetry solutions differ from each other. We show that due to combination of individual solutions systematic and random errors of the data processing can be reduced and the robustness of the geodetic derived AIS polar motion excitations can be increased. We investigate the interannual variability of the Antarctic polar motion excitation functions by means of piecewise linear trends. We find that the long-term behavior of the three ice sheet subregions: EAIS (East Antarctic Ice Sheet), WAIS (West Antarctic Ice Sheet) and APIS (Antarctic Peninsula Ice Sheet) is quite different. While APIS polar motion excitations show no significant interannual variations during the study period 2003-2015\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2003-2015$$\end{document}, the trend of the WAIS and EAIS polar motion excitations increased in 2006 and again in 2009 while it started slightly to decline in 2013. AIS mass changes explain about 45%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$45\%$$\end{document} of the observed magnitude of the polar motion vector (excluding glacial isosatic adjustment). They cause the pole position vector to drift along 59∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$59^{\circ }$$\end{document} East longitude with an amplitude of 2.7 mas/yr. Thus the contribution of the AIS has to be considered to close the budget of the geophysical excitation functions of polar motion.


Introduction
In recent years, climate change has led to increasing ice loss of the Antarctic Ice Sheet (AIS) which has an significant impact on polar motion. Due to redistribution and motion of masses within the Earth system the Earth's pole is continuously in motion. While the circular motion -mainly described by the Chandler and the annual oscillation -is caused by atmosphere and water variations, the non-linear drift results from changes in the solid Earth and ongoing ice melting as well as related sea level changes. The last two climate related mass relocations (ice and water changes) are principally responsible for the abrupt eastward turn of the mean pole position around the year 2005 (Chen et al. 2013). Whereas the reversal of polar motion after 2012 towards west may result from regional differences in terrestrial water storage changes (Adhikari and Ivins 2016). There exists a large number of studies focussing on atmospheric, oceanic and hydrological polar motion excitation mechanisms derived from geophysical models and space geodetic observations (e.g., Chen et al. 2017;Göttl et al. 2015;Jin et al. 2010; Open Access *Correspondence: franziska.goettl@tum.de 1 Technische Universität München, Deutsches Geodätisches Forschungsinstitut, Arcisstraße 21, 80333 Munich, Germany Full list of author information is available at the end of the article Göttl et al. Earth, Planets and Space (2021) 73:99 Meyrath and van Dam 2016; Nastula et al. 2007;Seoane et al. 2011), while there exist only a few studies regarding the contributions of the cryosphere (e.g., Chen et al. 2013;Adhikari and Ivins 2016). In contrast to the other subsystems of the Earth no adequate geophysical fluid models for the cryosphere are available up to now. Therefore in this study several gravimetric and altimetric data products are used to investigate the impact of AIS mass changes on polar motion in more detail.
Since 2002 ice mass changes in Antarctica can be observed by the Gravity Recovery and Climate Experiment (GRACE) in terms of time variable gravity field changes (e.g., Barletta et al. 2013;Luthcke et al. 2013;Velicogna and Wahr 2013;Groh and Horwath 2016;Su et al. 2018). GRACE is the only space-borne sensor which allows to observe mass redistribution directly with a spatial resolution of 200 -500 km of the entire AIS with monthly resolution. Besides the limited spatial resolution the accuracy of GRACE ice mass change estimates is limited by noise (meridional error stripes), leakage effects and uncertainties of glacial isosatic adjustment (GIA) models (e.g., Horwath and Dietrich 2009;Barletta et al. 2013), used to correct for the solid Earth's still ongoing response to past ice mass changes. Moreover, significant differences may arise from the algorithm applied to estimate the mass changes (Groh et al. 2019).
Since 1992 the temporal evolution of ice mass changes on the AIS can be derived from surface elevation changes measured by satellite radar and laser altimeter missions such as ERS-1, ERS-2, Envisat, CryoSat-2 and ICESat (e.g., Wingham et al. 2006;McMillan et al. 2014;Zwally et al. 2015;Schröder et al. 2019;Shepherd et al. 2019). In contrast to GRACE, satellite altimetry observes elevation changes which need to be converted to mass changes, causing other sources of errors. The spatial resolution of the altimetry data is significantly higher ( ∼ 2 km or less along the tracks according to the footprint size, ∼ 20 km or less according to the separation of neighboring satellite ground tracks) whereas the temporal resolution is comparable to GRACE (monthly repeat orbits for most of the missions). Due to a polar gap, reaching from 81.5 • S for ERS and Envisat down to 88 • S for CryoSat-2, 21% to 1% of the AIS near the South Pole cannot be observed by altimetry. Furthermore, the tracking of the satellite altimeters fails in sloping and rugged terrain near the coast, which is less than 2% of the total area . Nevertheless, these coastal areas are the location of the largest dynamic mass changes, which are not observable by altimetry. The accuracy of satellite altimetry derived ice mass change estimates is limited by waveform retracking, slope related relocation errors (especially in rugged terrain) and the density assumption in the volume-to-mass conversion. Here uncertainties due to GIA models are almost negligible because GIA induced height changes are significantly smaller than satellite altimetry observed ice surface height changes.
Both space geodetic observation techniques -satellite gravimetry and altimetry -have different strengths and weaknesses with respect to ice mass change estimations. To access the accuracy of the gravimetry and altimetry derived polar motion excitation functions we use different GRACE gravity field models and multi-mission satellite altimetry data. In order to reduce the systematic errors and to increase the robustness of the geodetic derived Antarctic polar motion excitation functions we combine the individual solutions. Based on these new improved time series we investigate not only the impact of the entire AIS on polar motion but also the contributions of the subregions: East Antarctic Ice Sheet (EAIS), West Antarctic Ice Sheet (WAIS) and Antarctic Peninsula Ice Sheet (APIS). Our focus lies in particular on the long-term behaviour and interannual variability.

Antarctic excitation of polar motion
In this section we explain briefly how Antarctic massrelated polar motion excitations can be determined from gravimetry and altimetry derived equivalent water height (EWH) anomalies of the AIS. Individual geophysical mass-related excitations of polar motion are mathematically described by equatorial angular momentum functions χ e 1 and χ e 2 ( Barnes et al. 1983;Gross 2015;Wahr 2005), where the index e denotes what kind of excitation mechanism is described by the equatorial angular momentum functions (e.g. AIS =: mass effect of the entire Antarctica, EAIS =: mass effect of the East Antarctic Ice Sheet, WAIS =: mass effect of the West Antarctic Ice Sheet, APIS =: mass effect of Antarctic Peninsula Ice Sheet, see Fig. 1). These so-called excitation functions are directly related to the fully normalized subsystem specific degree-2 potential coefficients C e 21 and S e 21 which are derived from equivalent water heights ewh using the global spherical harmonic analysis (GSHA) where θ and denote the spherical co-latitude and longitude of the computation point, n and m are the spherical harmonic degree and order, a is the Earth's equatorial radius, ρ e = 5517 kgm -3 is the mean density of the Earth, ρ w = 1025 kgm -3 is the density of sea water, k ′ n are the degree n load Love numbers, P nm (cos θ) are the normalized associated Legendre functions and σ denotes the unit sphere (Wahr et al. 1998). Note, that due to the application of the load Love numbers the elastic deformation of the solid Earth caused by surface load is taken into account in the potential coefficients. Thus the real parts of the mass-related polar motion excitation functions can be written as where G is the gravitational constant, GM is the geocentric gravitational constant, C is the axial moment of inertia of the Earth and A ′ = (A + B)/2 is the average of the equatorial principal moments of inertia. The constant α takes into account the rotational deformation of the Earth. It is calculated via where is the Earth's mean angular velocity, A c is the equatorial moment of inertia of the Earth's core, ǫ c is the ellipticity of the Earth's core and σ 0 = 2π/T is the real part of the Chandler frequency derived from the Chandler period T. The imaginary parts of the polar motion excitation functions can be neglected because they are smaller than 1% of the real parts of the polar motion excitation functions. In this study, we determine Antarctic polar motion excitations, therefore we apply the constants of a tide-free Earth model listed in Table 1.
Equation 2 corresponds to that of Gross (2015) and Dobslaw and Dill (2019) within 1% , the differences are based on dissimilar values of the numerical constants. The dimensionless polar motion excitation functions are converted to milliarcseconds (mas) by the factor f = 360 · 60 · 60 · 1000/(2π) . They describe the position of the excitation axis with respect to the Earth-fixed terrestrial reference frame.

Data and data processing
This section provides an overview of the GRACE gravity field models and multi-mission satellite altimetry solutions which are used within this study to determine AIS mass changes and their impact on polar motion.

Gravimetry derived EWH anomalies of the AIS
In this study we use EWH anomalies of the AIS derived from four monthly GRACE gravity field solutions: CSR RL06M (Save et al. 2016;Save 2019), JPL RL06M (Watkins et al. 2015;Wiese et al. 2016Wiese et al. , 2019, ITSG-Grace2018 (Mayer-Gürr et al. 2018) and LDCmgm90 (Chen et al. 2019(Chen et al. , 2020. As investigated by Chao (2016) one should keep in mind that EWH and surface mascon solutions cannot represent internal processes in a physically meaningful way but they are a appropriate representation of surficial gravitational processes such as ice mass changes. Important for Earth rotation studies is that these gravity field models are all based on the new linear mean pole model, because this has a significant impact on the potential coefficients C 21 and S 21 as shown in Göttl et al. (2018). The degree-1 coefficients of all gravity field solutions are replaced by estimates listed in the GRACE Technical Note 13 (Swenson et al. 2008;Sun et al. 2016) to take into account that a redistribution of masses is referred to a coordinate system attached to the Earth's Table 1 Parameters of the Earth used for the determination of Antarctic polar motion excitation functions. Sources: (a) Petit and Luzum (2010), (b) Wilson and Vicente (1990), (c) Göttl (2013), (d) Seitz et al. (2012), (e) Mathews et al. (1991)  crust which moves relatively to the Earth's center-ofmass frame applied within the GRACE data processing. Whereat the motion of the center-of-mass (CM) with respect to the center-of-figure (CF) of the solid Earth surface is defined as geocenter motion. Furthermore the inaccurate C 20 coefficient is replaced by an external satellite laser ranging (SLR) solution from Loomis et al. (2020) (GRACE Technical Note 14). GIA induced mass changes are removed from the total GRACE signal to identify ice mass changes in Antarctica. Martin-Español et al. (2016) compared eight forward and inverse GIA models for Antarctica. They found that GIA models induce uncertainties on GRACE derived present-day ice mass variations of about 60 Gt/yr. In terms of GRACE derived AIS polar motion excitation functions the uncertainties are about 0.2 mas/yr. For consistency reasons we decide to apply the GIA model ICE6G_D (Peltier 2015) to all gravity field solutions because it is already used as the standard model for the GRACE mascon solutions (CSR RL06M, JPL RL06M) and the applied degree-1 coefficients. All gravimetry derived EWH anomalies are interpolated to a regular 1 • × 1 • grid and masks for the entire AIS as well as for the subregions EAIS, WAIS and APIS are applied. In general, EWH anomalies for a particular subsystem of the Earth do not conserve mass globally.
To ensure mass conservation, the EWH anomalies may be complemented by an additional EWH layer over the ocean. This layer represents the ocean's gravitational selfconsistent response to the changing surface loads (Clarke et al. 2005) over the ice sheet region under investigation, including the rotational feedback of the additional ocean load (Rietbroek et al. 2012). By using the Eqs. 1 and 2 the polar motion excitation functions for the specific regions are determined. The focus of this study is on the time period 2003 to 2015, i.e., excluding the "GRACE single ACC" period where the accelerometer (ACC) on board GRACE-B was turned off. The quality of the "GRACE single ACC" gravity field solutions is less and a much stronger decorrelation and smoothing is required to identify mass changes (Dahle et al. 2019). Therefore we do not include these GRACE gravity field models in our investigations.

GRACE mascon solutions
The CSR RL06M and JPL RL06M gravity field models are based on geolocated spherical cap mass concentration functions. One advantage of the mass concentration (mascon) parameters over spherical harmonics is a higher reduction of noise and errors in the gravity field solutions due to the limited longitudinal sampling and orbital configuration of the satellite mission GRACE. Thus no post-processing in terms of destriping or smoothing is required -such as filtering and forward modeling -to reduce the north-south stripes and to estimate meaningful mass anomalies (Andrews et al. 2015). Besides the improved signal-to-noise ratio the GRACE mascon solutions have a higher spatial resolution and lower leakage errors from neighboring mass processes. The native resolution on an equal-area geodesic grid is 1 • for the CSR RL06M solutions versus 3 • for the JPL RL06M solutions, but the provided mass anomaly grid files are given on a regular 0.25 • × 0.25 • and 0.5 • × 0.5 • grid, respectively, to represent the mascon grids at coastlines properly. While at CSR (Center for Space Research, Austin) the hexagonal grid tiles along the coastlines are split into two tiles to minimize the leakage between land and ocean signals, at JPL (Jet Propulsion Laboratory, Pasadena) the Coastal Resolution Improvement (CRI) filter has been developed to reduce leakage errors across coastlines in a post-processing step. A further difference between the two mascon solutions is that CSR uses the Tikhonov regularization along with the L-ribbon approach (Save et al. 2012) to derive time dependent regularization parameters from GRACE information only, whereas at JPL the correlated errors are removed by introducing realistic geophysical data or model during the solution inversion.

GRACE spherical harmonic solutions
The ITSG-Grace2018 gravity field solution is based on the global spherical harmonic potential coefficients. As mentioned before unconstrained spherical harmonic gravity field solutions suffer from erroneous meridional stripes. In this study we use two data sets of EWH anomalies derived from ITSG-Grace2018 by using two different post-processing methods: (1) The method of tailored sensitivity kernels developed at TU Dresden (TUD) in the frame of the Antarctic Ice Sheet project of the European Space Agency's (ESA) Climate Change Initiative (CCI) (Groh and Horwath 2016) and (2) the filter effect reduction approach on global grid point scale developed in the frame of the German Research Foundation (DFG) project CIEROT (Combination of geodetic space observations for estimating cryospheric mass changes and their impact on Earth rotation) at Deutsches Geodätisches Forschungsinstitut der TU München (DGFI-TUM) (Göttl et al. 2019). Therefore we abbreviate these gravimetry EWH anomalies solutions as ITSG-Grace2018/TUD and ITSG-Grace2018/TUM, respectively, according to the applied post-processing approaches. The first method is a dedicated variant of the regional integration approach that implicitly applies tailored sensitivity kernels. The tailored sensitivity kernels are designed to minimize both GRACE error effects and signal leakage. Hence, they are based on information about the variances and covariances of GRACE monthly solution errors as well as on geophysical signals that induce leakage. The ITSG-Grace2018/TUD mass anomalies are derived on a polar stereographic grid with a resolution of 50km × 50km . The second method has been developed especially for Earth rotation studies. It is independent from geophysical model information and works for several subsystems of the Earth. Here the filter effects (attenuation and leakage) are reduced due to application of global grid point gain factors estimated from once and twice filtered GRACE gravity field solutions. In a second step on the level of polar motion excitation functions scaling factors are applied to counteract the damping of the gain factors due to the filtering processes.

Multiple-data spherical harmonic solutions
In contrast to the gravity field solutions mentioned above, LDCmgm90 is a combined gravity field model based not only on GRACE observations but on multiple data. The potential coefficients C 20 , C 21 and S 21 are based on GRACE and SLR gravity field models as well as information from geophysical models for the atmosphere, oceans and continental hydrosphere, whereas the other potential coefficients are a combination of the CSR and JPL GRACE mascon solutions; for more details see Chen et al. (2017) and Yu et al. (2018). Like the mascon solutions the multiple-data-based gravity field model does not suffer from meridional error stripes and the leakage effects are significantly reduced. EWH anomalies are obtained by using the global spherical harmonic synthesis (GSHS) No post-processing in terms of destriping or smoothing is required.

Altimetry derived EWH anomalies of the AIS
In this study we use EWH anomalies of the AIS derived from two multi-mission satellite altimetry solutions, one for monthly gridded surface elevation changes (SEC) in Antarctica determined at TUD (Schröder et al. 2019) and one for monthly integrated mass changes in the AIS drainage basins provided by the University of Leeds (UL) . To ensure meaningful comparisons with gravimetry derived EWH anomalies GIA induced height and mass changes respectively are removed with the same GIA model ICE6G_D (Peltier 2015). Furthermore, the altimetry derived mass anomalies are interpolated to a 1 • × 1 • grid like the GRACE data and masks for the subregions EAIS, WAIS and APIS are applied. To guarantee mass conversation the EWH anomalies are complemented by an additional EWH layer over the ocean as for the GRACE data. Finally the polar motion excitation functions for the specific regions are determined via Eqs. 1 and 2.

Surface elevation changes on grid scale
The monthly TUD SEC data of the AIS are based on altimeter observations from the satellite missions ERS-1, ERS-2, Envisat, ICESat and CryoSat-2 covering the time period 1992 to 2017. In order to obtain comparable results the focus of this study is on the time period 2003 to 2015 like the GRACE data, where mainly the altimeter missions Envisat, ICESat and CryoSat-2 were operating covering the globe up to 81.5 • , 86 • and 88 • latitude, respectively. While Envisat carries a classical pulse-limited radar altimeter on board (2 km pulse-limited footprint size), CryoSat-2 is equipped with a Delay-Doppler/ Synthetic Aperture Radar (SAR) altimeter with a significantly improved along-track footprint size of about 0.3 km. ICESat carries a laser altimeter with high spatial resolution (70 m footprint size) on board. For a consistent combination of the different satellite missions TUD applied a refined reprocessing of the radar altimetry data, which also contributed to significantly improved accuracy of the observation data; for details see Schröder et al. (2019). The SEC data are given on a 10km×10 km polar stereographic grid. They originate from changes in the ice flow, changes in the surface firn (as a lack or excess in snowfall) or even elevation changes in the underlying bedrock. The latter are the combined effect of the Earth's elastic response on present-day ice mass changes and GIA induced by past ice mass changes. We use the volume-to-mass conversion according to Schröder et al. (2019) to estimate EWH anomalies of the AIS. In a first step, the uplift rates due to GIA have to be removed. To maintain consistency, we use the GIA model ICE6G_D (Peltier 2015). In the next step, the reduced SEC data are multiplied by the scaling factor s ela = 1.0205 to account for the elastic solid Earth rebound effect (Groh et al. 2012). The corrected SEC data represent the ice sheet thickness changes and are multiplied with a time-invariant density mask ρ , adapted after McMillan et al. (2014), with varying firn density from Ligtenberg et al. (2012). We interpolate the retrieved EWH anomalies to a regular 1 • × 1 • grid like the GRACE data and fill the data gaps beyond the southern limit of the satellite orbits ( 1 to 21% ) and in (5) �ewh = (SEC − GIA) · s ela · ρ rugged terrain ( < 2% ) where the satellite altimeters failed to track the elevation changes with data from the ITSG-Grace2018/TUM gravity field solutions. In regions with large scale mass variations due to precipitation changes this method seems to be correct, whereas in regions with small scale mass variations due to changes in the ice dynamic one has to keep in mind that GRACE data underestimate the effect due to the fact that the signal is smeared over a larger region.

Mass anomalies on basin scale
The monthly integrated mass change time series for the individual AIS drainage basins provided by the University of Leeds are based only on radar altimetry observations from the satellite missions Envisat and CryoSat-2 in the investigation period 2003 to 2015 of this study. Therefore this altimeter solution relies on different input data sets (no ICESat data) as well as data processing, outlier detection and merging strategy. However, the main difference between both satellite altimetry data sets is the applied volume-to-mass conversion. While Schröder et al. . Thus the UL mass anomalies contain no detailed information about the exact location of the mass changes within the basin. Therefore we distribute the mass changes equally within each basin to obtain gridded mass anomalies. In order to quantify the consequences of this lack of spatial information in the UL mass data, we produced a comparable data set from the TUD data and compared the results. Therefore, we formed cumulated basin time series from the gridded data, redistributed the signal equally over the whole basin and calculated the resulting polar motion excitation functions. Compared to the results of the spatially distributed data set (Fig. 2), we see that the temporal structure of the signals is very similar but the trend of equally distributed mass changes on polar motion is about 4% for χ 1 and 7% for χ 2 smaller. Hence, the impact of AIS mass changes on polar motion derived from the UL solutions is slightly underestimated. Equivalent comparisons for APIS, WAIS and EAIS show similar results. The relative standard deviations (RSD) due to equally distributed mass changes are about 5% for APIS, 2 or 5% for WAIS ( χ 1 , χ 2 ) and 4 or 6% for EAIS ( χ 1 , χ 2 ). We conclude that this approximation is adequate to study AIS polar motion excitations.

Results, combination and validation
In this section we show Antarctic polar motion excitations derived from GRACE and satellite altimetry data. We combine the individual gravimetry and altimetry solutions to reduce systematic errors and to improve the reliability of the geodetic derived AIS polar motion excitations. To quantify the quality of the individual and combined results comparisons with solutions from the multiple-data-based gravity field model LDCmgm90 are performed.

Gravimetry and altimetry results
In Fig. 3 all gravimetry and altimetry solutions of this study are shown, not only for the entire AIS but also for the subregions APIS, WAIS and EAIS (according to Rignot et al. (2011), see Fig. 1). In this way the individual contributions of the subregions to Earth rotation variations can be investigated. The GRACE derived excitation functions show higher agreement among themselves than the satellite altimetry derived excitation functions. This follows from the fact that the altimetry solutions are based on different approaches for the volume-to-mass conversion which is one of the largest error sources. Sasgen et al. (2019) have shown that the uncertainties of altimetry derived ice mass changes for Antarctica due to the density assumption, retracking and adjustment method are significantly higher than the uncertainties of GRACE derived ice mass changes due to GIA corrections and uncertainties of the potential coefficients. In this study the GRACE solutions are homogenized by an uniform GIA model to eliminate one major source of discrepancy. Therefore the uncertainties of the GRACE derived Antarctic excitation functions due to different post-processing strategies are 0.1 mas for χ 1 and 0.7 mas for χ 2 which is about 1 and 6% of the magnitude of the AIS excitation functions, respectively. Uncertainties due to the GIA model amount 0.2 mas for χ 1 and 0.3 mas for χ 2 (RSD: 4% and 3% ). The uncertainties of altimetry derived AIS excitation functions are significantly higher about 0.8 mas for χ 1 and 1.1 mas for χ 2 (RSD: 6% and 15% ). In Fig. 4 the mean correlation coefficients and relative standard deviations between the individual gravimetry solutions (10 combinations) and the both altimetry solutions for AIS, APIS, WAIS and EAIS are shown, respectively. The higher the correlation value the higher is the accordance of the signal structures. The highest differences can be seen for the APIS. Reasons therefore could be that the area of APIS (232, 000 km 2 ) is significantly smaller than the area of WAIS (2038, 000 km 2 ) and EAIS (9620, 000 km 2 ) and that the portion of coastal mass variations of the narrow peninsula is very high. It holds, that the smaller the catchment area the larger the leakage effect and the lower the accuracy of the GRACE derived mass variations. Due to the fact that the signal strength is quite different for the four regions AIS ( −13 to 19 mas), APIS ( −2 to 3 mas), WAIS ( −9 to 15 mas), EAIS ( −4 to 6 mas) we need to have a look at the RSD instead of the root mean square (RMS) differences. The RSD is defined as the ratio of the RMS to the maximum value of the arithmetic mean of all gravimetry or altimetry solutions, respectively. Depending on the region, the accordance of the solutions is quite different (see Fig. 4). While polar motion excitations for the WAIS are well determined by both geodetic space observation techniques, the results for APIS and EAIS show significantly higher differences. The altimetry solutions for the EAIS show the largest differences. We assume that these discrepancies result from uncertainties of the volume-tomass conversion caused by changes of the ice dynamics and in the firn pack.
Comparing the GRACE and satellite altimetry solutions for the EAIS reveals that the trends of the altimeter time series are significantly lower than the trends of the GRACE solutions. Figure 5 shows that these differences result mainly from the basins: Jpp-K, K-A, A-Ap, Ap-B, B-C. Whereat the basins A-Ap (Dronning Maud Land) and Ap-B (Enderby Land) make the largest contributions.
Here the trend of the GRACE-derived polar motion excitation functions is significantly larger than the trend of the satellite altimetry-derived polar motion excitation functions, especially after 2011. One reason for this could be that the TUD altimetry solutions due not take into account that parts of the mass changes in Dronning Maud Land and Enderby Land could have been caused by ice dynamics (e.g. dynamic thickening) (Schröder et al. 2019). The UL altimetry solutions consider ice dynamic effects but they cannot be fully compensated because of uncertainties in the surface mass balance (SMB) models. Furthermore the UL solutions suffer from uncertainties due to equally distributed mass changes in the individual drainage basins (see RSD values in Fig. 5). Another reason could be that GRACE overestimate the mass signal in EAIS, for example due to uncertainties in GIA models. Martin-Español et al. (2016) show that for EAIS the uncertainties of the GIA models are generally lower than for WAIS and APIS but due to the large area of the EAIS basins these uncertainties have a larger impact on GRACE-derived integrated mass changes than it is the case in WAIS and APIS. Furthermore they showed that within EAIS for the basins: Jpp-K, K-A, Ap-B, B-C, C-Cp and E-Ep the uncertainties of the GIA models are higher. This coincide with our investigations on basis of the polar motion excitation functions.

Combination of GRACE and satellite altimetry data
In this section gravimetry and altimetry solutions for Antarctic polar motion excitation functions are combined in order to cancel systematic errors of the different processing strategies and to gather the strengths of the different space geodetic observation techniques. Fig. 2 Comparison of the Antarctic polar motion excitation functions χ AIS 1 and χ AIS 2 derived from TUD mass changes given on a 1 • × 1 • grid (blue) and derived from TUD integrated mass changes in the AIS drainage basins (red). The corresponding correlation coefficients and relative standard deviations (RSD) are provided Advantages of the satellite mission GRACE are that ice mass changes can be observed directly for the entire AIS whereas advantages of the altimeter satellite missions are the significantly higher spatial resolution and no smearing of spatially detailed information like it is the case by GRACE observations.
The individual input solutions are consistent regarding the determination of the equivalent water heights (temporal and spatial resolution, GIA, elastic deformation) and the Antarctic polar motion excitation functions (Earth's parameter), but the geophysical background models of the input solutions are partly different. This has the advantage that systematic errors of the background models can be adjusted within the combination. We determine four combined solutions for Antarctic polar motion excitations: (1) arithmetic mean of the gravimetry solutions ITSG-Grace2018/ TUM, ITSG-Grace2018/TUD, CSR RL06M and JPL RL06M, (2) arithmetic mean of the altimetry solutions TUD and UL, (3) arithmetic mean of the gravimetry  (4) weighted combination of the gravimetry and altimetry solutions based on the least squares adjustment approach described in Göttl et al. (2012) to show that via a combination an improved excitation time series can be retrieved. Hereinafter we refer to the results from the latter combination approach as "adjusted solutions". The least squares adjustment approach will only be described briefly: The stochastic model takes into account the empirical variances of the observations as well as the time dependency of the noise of the unknown parameters. The empirical variances of the gravimetry and altimetry estimated excitation functions are calculated via where χ e j (t k ) is the average of all Antarctic polar motion excitation functions χ e j,p with j ∈ 1, 2 derived from different processing centers and observation data with different processing strategies ( p ∈ 1, ..., P ) at the discrete time moment t k with k = 1, ..., K (total number of months). In Fig. 6 the empirical standard deviations of the GRACE and satellite altimetry solutions are shown. The higher the empirical standard deviation the lower the weight within the least squares adjustment. It becomes evident that depending on the region and the polar motion excitation functions the weighting of the gravimetry and altimetry solutions differs significantly. Except for the subregion APIS the altimetry solutions of UL have the lowest weights whereas the weights of the TUD satellite altimetry solutions are in the same order like the weights of the gravimetry solutions. Via the empirical auto-covariances the temporal dependency of the noise of the unknown parameters can be considered in the stochastic model. As shown by Göttl (2013) in this way the reliability of the formal errors of the adjusted time series can be improved. No correlations between the time series based on observations from the same measurement technique are taken into account because they are not precisely known. Due to the application of different processing strategies it can be assumed that they are small. Simulations have shown that the estimated variances of the unknown parameters differ only by about 2 to 16% from the true variances. The formal errors of the adjusted GRACE and satellite altimetry solutions are given in Table 2. For EAIS they are larger than for APIS and WAIS.

Comparison with LDCmgm90
Due to the lack of accurate geophysical model results for Antarctic polar motion excitation functions we compare Fig. 4 The mean relative standard deviations (RSD) are displayed along the horizontal axis and the mean correlation coefficients are displayed along the vertical axis. The higher the correlation and the lower the RSD (upper left corner) the higher is the agreement of the gravimetry solutions (blue) among themselves and the altimetry solutions (red), respectively. The agreement between the gravimetry and altimetry solutions is shown in green our individual and combined satellite gravimetry and altimetry solutions with estimates from the multiple-data spherical harmonic solution LDCmgm90. As mentioned in the data description section the potential coefficients C 20 , C 21 and S 21 are based besides several GRACE gravity field solutions on SLR gravity field solutions, geophysical model information and precise Earth orientation parameters whereas the other potential coefficients are only based on GRACE mascon solutions. Altimetry solutions are not considered in the LDCmgm90 solution as well as spherical harmonic GRACE solutions for potential coefficients beyond degree and order 2. Therefore the gravimetry and altimetry solutions of this study are as far as possible independent from the LDCmgm90 solution. Figure 7 shows the differences of the combined solutions: DGFI-TUM and LDCmgm90. The correlation coefficients and relative standard deviations for all individual and combined gravimetry and altimetry solutions are given in Table 3. We expect the higher the correlation and the lower the RSD value with our gravimetry, altimetry and combined solutions the higher is the reliability of the resulting time series. We found that except for the polar motion excitations χ APIS 1 , χ APIS 2 and χ EAIS 2 the accordance with the GRACE solutions is very high (RSD: 2 to 12% ). The best agreement is shown with the ITSG-Grace2018/TUD and ITSG-Grace2018/TUM solutions. For APIS the accordance with the altimetry solutions is partly higher (RSD: 15 to 27% ) than with the gravimetry solutions (RSD: 16 to 43% ) while for EAIS the accordance with the altimetry solutions is lower (RSD: 8 to 60% versus 2 to 12% ). It seems to be that the leakage errors in GRACE mass estimates for the small and narrow Fig. 5 Comparison of polar motion excitation functions for the single EAIS basins derived from the combined GRACE gravity field model LDCmgm90 (black) and derived from satellite altimetry data provided by TUD (red) and UL (brown). The relative standard deviations (RSD) of the UL solutions due to equally distributed mass signals are provided Antarctic peninsula are higher than the errors in satellite altimetry mass estimates due to rugged terrain and the density assumption. Except for APIS the agreement of the LDCmgm90 and TUD altimeter solutions is higher than with the UL altimeter solutions. By determining the arithmetic mean of the individual time series we want to prove if the reliability and consistency of the merged solutions can be improved. By considering all regions and both components of the polar motion excitations, the combination of all gravimetry time series reveals a larger overall agreement (RSD: 4 to 21% ) compared to the individual gravimetry time series whereas a combination of the altimetry solutions only slightly improves the agreement with the LDCmgm90 results (RSD: 2 to 56% ). By merging the two different space geodetic observation techniques the accordance can be further improved especially for the contributions of the APIS and WAIS. Reasons therefore might be that the uncertainties of the GRACE solutions for these regions are larger due to the greater leakage effects and therefore can be adjusted by the altimetry solutions. These improvements confirm that even due to a simple combination (arithmetic mean) of GRACE and satellite altimetry data systematic and random errors can be reduced and the robustness of the geodetic derived Antarctic polar motion excitation functions can be increased. By applying the weighted least squares adjustment approach described by Göttl et al. (2012) further improvements can be achieved (RSD: 2 to 19% ). The polar motion excitation function χ EAIS 2 show the highest uncertainties, this phenomenon requires further investigation. The best results for the entire AIS can be obtained through summation of the weighted adjusted GRACE and altimetry solutions for the subregions APIS, WAIS and EAIS (RSD: 2% ). Thus the results of the two different combination approaches show a high concordance. Figure 8 shows the adjusted AIS polar motion excitation functions. While χ 1 is dominated by mass changes close the meridians 0 • and 180 • , χ 2 is dominated by mass changes close to the meridians ±90 • . Thus, ice mass changes in EAIS have the largest impact on χ AIS 1 whereas ice mass changes in WAIS have the largest impact on χ AIS 2 . Due to the fact that most of the WAIS sits on rock below sea level, which is not the case for EAIS, it is more sensitive to global warming. Thus WAIS is dominated by a strong ice loss whereas EAIS is still dominated by a  smaller ice gain mainly due to snow accumulation. Therefore the magnitude of χ AIS 2 is about 5 mas higher than the magnitude of χ AIS 1 . Furthermore, ice loss in APIS slightly increase the trend of χ AIS 2 and decrease the trend of χ AIS 1 . In order to study the interannual variability of the AIS polar motion excitation functions Fig. 8  is almost constant during the study period and amounts to −0.1 mas/yr. For χ 2 the trends feature different characteristics. For χ WAIS 2 the trend began to rise already in 2006 from 0.6 mas/yr to 1.2 mas/yr and again in 2009 to 2.3 mas/yr. In 2013 the trend started slightly to drop to 1.9 mas/yr. One reason for this is that due to global warming and shifts in wind patterns warmer ocean water especially in the Amundsen Sea sector induce a stronger deglaciation (Holland et al. 2019). Table 3 Correlation coefficient/relative standard deviation (RSD) between GRACE and satellite altimetry solutions for Antarctic polar motion excitation functions and results from the combined gravity field solution LDCmgm90. The best results for the entire AIS can be maintained through summation of the weighted adjusted GRACE and altimetry solutions for the subregions APIS, WAIS and EAIS. The best individual and combined results are highlighted by bold values  is nearly constant during the study period it amounts to 0.3 mas/yr. In Fig. 9 the piecewise linear trends of the geodetic observed polar motion excitation functions are shown. They are derived from the precise pole coordinates x and y provided by the International Earth Rotation and Reference System Service (IERS) via the EOP (Earth Orientation Parameters) 14 C04 time series. To ensure consistency we remove the GIA induced linear drift signal from the observed polar motion excitations with the same GIA model ICE6G_D (Peltier 2015) as used within the determination of the AIS polar motion excitation functions. The trend increased significantly in 2006 and started to decrease around 2013. We found that these turning points coincide well with the turning points of the AIS polar motion excitation functions.

Impact on polar motion
It holds the drift of the excitation pole equals the drift of the celestial intermediate pole (CIP) (Gross 2015). According to the adjusted AIS polar motion excitation functions the pole vector drift along 59 • East longitude of amplitude 2.7 mas/yr due to AIS ice mass changes during the study period 2003 − 2015 . These findings are close to the results in Adhikari and Ivins (2016). The observed pole position vector drift along 46 • East longitude at a rate of 6 mas/yr after removing the GIA signal. Thus AIS polar motion excitations alone explain about 45% of the observed magnitude of the polar motion vector (excluding GIA) and the AIS induced pole drift deviates only about 13 • from the GIA reduced observed pole drift.

Summary and conclusions
In this study we determined for the first time Antarctic polar motion excitation functions. They were computed from two independent data sources, namely from GRACE time variable gravity fields and from multimission satellite altimetry surface elevation changes. To assess the accuracy of the gravimetry derived polar motion excitations we use different GRACE gravity field models based on mass concentration parameters (CSR RL06M, JPL RL06M) or on spherical harmonics (ITSG-Grace2018). The combined gravity field model LDCmgm90, which is based on multiple-data such as GRACE, SLR, geophyiscal models and EOP, is used as reference time series. While the differences of the gravimetry solutions for WAIS are small ( 3% ), the differences for APIS are significantly higher ( 25% ). Reasons for this are that the area of the Antarctic peninsula is significantly smaller and the portion of coastal mass changes is very high. To assess the accuracy of the altimetry derived Antarctic polar motion excitations we use monthly gridded SEC data (TUD) and integrated mass change time series for the individual AIS basins (UL). For WAIS the accordance is high ( RSD = 6% ) whereas for EAIS and APIS the uncertainties are significantly higher ( 20 to 30% ). One reason for this is that the narrow peninsula is extremely rugged and another reason might be that parts of the mass changes in EAIS (Dronning Maud Land, Enderby Land) may have been caused by ice dynamic effects which are difficult to be taken into account. We show that due to the combination of the gravimetry and altimetry solutions the systematic errors of the data processing can be reduced and the robustness of the geodetic derived Antarctic polar motion excitation functions can be increased. The agreement with the combined gravity field solution LDCmgm90 can be significantly improved. The RSD values for the AIS amount 2% only although the combination approaches are based on different input data: SLR, geophysical models and EOP versus satellite altimetry. Thus the impact of ice mass changes in the Antarctica on polar motion can be determined accurately due to combination of different geodetic space observation techniques. Based on these investigations we found that ice mass changes in EAIS have a larger impact on the x pole coordinate whereas ice mass changes in WAIS have a larger impact on the y pole coordinate. The trend behavior of the three ice sheet subregions EAIS, WAIS and APIS is quite different. While APIS polar motion excitations show no significant interannual variations due to global climate change during the study period, the trend of the WAIS and EAIS polar motion excitations increased in 2006 and again in 2009 while it started slightly to decline in 2013. Within this study we found that AIS mass changes induce the pole position vector to drift along 59 • E by 2.7 mas/yr during the study period 2003 − 2015 , which explain about 45% of the observed magnitude of the polar motion vector (excluding GIA). In comparison, mass variations of Greenland and the continental hydrosphere cause the pole position vector to drift approximately along 38 • W by 3.6 mas/yr and 80.5 • E by 2.4 mas/yr (Adhikari and Ivins 2016), respectively. While ice mass changes of inland glaciers and collective mass redistribution of the atmosphere and oceans have no significant impact onto the drift of the polar motion vector. Therefore it is important to reduce not only atmospheric and oceanic signals from observed polar motion, as it is usually done, but also the contributions of Antarctica and Greenland to identify for example hydrological signals in observed polar motion. While there exist adequate geophysical models for the atmosphere and the oceans there exist up to now no adequate models for the cryosphere. The current study contributes to overcome this deficiency. Furthermore these improved AIS polar motion excitation time series can be used to improve modelling of ice mass redistribution in the Antarctica.