Separating Long‐Term and Short‐Term Mass Changes of Antarctic Ice Drainage Basins: A Coupled State Space Analysis of Satellite Observations and Model Products

Satellite gravimetry and altimetry measurements record gravity and elevation changes, respectively, which are useful for determining mass and volume change of the Antarctic Ice Sheet. Common methods employ products from regional climate modeling and firn modeling to aid interpretation and to link volume changes to mass changes. Estimating deterministic parameters over defined time periods is a conventional way to evaluate those changes. To overcome limitations of deterministic analyses with respect to time‐variable signals, we have developed a state‐space model framework. Therein, we jointly evaluate four mass and volume data sets by coupling of temporal signal variations. We identify long‐term signals of ice drainage basins that are observed by the satellite gravimetry mission GRACE and several satellite altimetry missions from April 2002 until August 2016. The degree to which we can separate long‐term and short‐term variations strongly depends on the similarity of the mass and volume change time series. For the drainage system of the Pine Island Glacier (West Antarctica), our results show noticeable variations of the long‐term trend with an acceleration of the contribution of ice dynamics to the mass balance from −11 ± 8 to −58 ± 8 Gt a−1. Our results in Dronning Maud Land (East Antarctica) show a positive long‐term contribution to the mass balance at almost a constant rate. The presented approach can fit time‐variable changes without artificial selection of periods of interest. Furthermore, because we only enforce common long‐term time variations between mass and volume data, differences in mean trend rates help to uncover model discrepancies.

changes of the ice surface geometry. Reconciled estimates of the mass balance of the AIS are −73 ± 52 and −219 ± 43 Gt a −1 over the period 2002-2007 and 2012-2017, respectively (Shepherd et al., 2018).
The mass balance of an ice sheet has two (major) components. One component is the surface mass balance (SMB) involving precipitation, surface and drifting snow sublimation, drifting snow erosion, and meltwater runoff (van Wessem et al., 2018). Those processes take place mainly in the snow and firn layer of the ice sheet, hereinafter referred to as SMB-driven signals. The other component is ice dynamics (ID) taking place in the ice layer, hereinafter referred to as ID-driven signals. We do not account for changes of basal (or bottom) melting of the grounded ice sheet separately. There are very few direct observations of basal melting of the grounded ice sheet (Seroussi et al., 2017). Pattyn (2010) found that the mean of the basal melt rate of the grounded AIS is ∼3% of the total surface accumulation, but it remains an open question how this contribution changes temporally. If the ice sheet is in a state of equilibrium, SMB-driven and ID-driven mass changes are in balance. Due to changing boundary conditions, for example, as a result of a changing climate, we expect mass and volume changes on long-term and short-term time scales. Long-term changes are characterized by a rate that varies only slowly at decadal scales. Short-term changes oscillate around zero at sub-decadal scales and tend to average out over decadal scales (Zwally et al., 2015). In Antarctica, the dominant ID-driven changes are on long time scales (Rignot et al., 2019). Here we consider long-term ice dynamic changes as dynamic thinning or dynamic thickening if the ice discharge is larger or smaller than the long-term SMB, respectively.
Satellite gravimetry and altimetry cannot distinguish between the two components of the mass balance, but only observe the sum of all mass and volume changes. SMB-driven mass and volume changes can be obtained from the SMB product from a regional climate model and from the firn thickness change product from a firn densification model (FDM). On drainage basin scale GRACE and SMB as well as altimetry and FDM show strong correlations of the nonlinear signal components (Figure 1), which motivates us to further investigate the temporal variations of mass and volume changes. Horwath et al. (2012) and Mémin et al. (2015) found common interannual variability in GRACE and altimetry observations and could make links to their geophysical origin. They point out the limitations of treating long-term signals as a simple linear trend with a constant rate. Davis et al. (2012) and Didova et al. (2016) applied state space modeling to geodetic observations of ice mass changes. These studies showed the advantages of state space modeling for estimating trends with time-variable rates and more realistic trend uncertainties from GRACE and global navigation satellite system (GNSS) time series, compared to deterministic results from conventional least-squares adjustment. Zammit-Mangion et al. (2015) and Martín-Español et al. (2016) used a Bayesian hierarchical model to estimate annual ID-driven, annual SMB-driven, and linear glacial isostatic adjustment (GIA) mass changes from satellite altimetry, satellite gravimetry, and GNSS. In this approach, spatio-temporal parameters are taken from auxiliary observations and models, and were adjusted to GRACE, altimetry, and GNSS observations using a statistical inversion scheme.
In this study, we separate the long-term changes and the short-term changes on the drainage basin scale in Antarctica. For this purpose, we examine the GIA-corrected residual GRACE and altimetry time series after subtracting modeled SMB-driven mass change and firn-related volume change, respectively, to account for the bulk of SMB and firn signals. In our approach, we assume that (a) short-term variations of those residual time series consist of regional climate model and firn model errors next to observational errors; (b) gravimetry and altimetry are sensitive to the same long-term ID-driven variations; and (c) decadal and centennial SMB-driven signals and firn-thickness trends (e.g., Medley & Thomas, 2019) unaccounted by the model products are not predominant in the regions under investigation. We use a coupled state space model to jointly evaluate these four data sets. We estimate the time-variable signals and simultaneously relate common temporal variations from gravimetry and altimetry data.
The basins in this study refer to drainage systems defined by Zwally et al. (2012). In the main text we focus on three drainage basins (Figure 1), namely a part of Dronning Maud Land (basin 6), Wilkes Land (basin 13), and the drainage basin of the Pine Island Glacier (basin 22). We have chosen these basins because they differ in the dominance of long-term and short-term changes (Rignot et al., 2019) and should illustrate the feasibility of our approach. In the supporting information (SI), we provide material of all investigated drainage systems.

Model Products
To model the SMB-driven mass change we use the RACMO2.3p2 SMB product (van Wessem et al., 2018). The continent-wide version of this product has a spatial resolution of 27 km and is sampled monthly from January 1979 to December 2016 for the whole AIS. SMB anomalies are derived by assuming a representative climatology over a long-term period, that is, more than 30 years (van den Broeke et al., 2009). SMB anomalies are the residuals to this long-term mean SMB. In our case we define the whole modeling period from January 1979 to December 2016 as the reference period. This is consistent with the period considered during the spin-up of the FDM that we use. All our mass change estimates refer to this reference period. Table S1 summarizes the mean SMB over the reference period with an uncertainty estimate analogous to Wouters et al. (2015). For this, we have calculated the standard deviation of mean SMB values over all 25-, 30-, and 35-year time periods within the reference period (Table S1). Because satellite gravimetry and satellite altimetry observe cumulated mass balance anomalies over time, we temporally integrate SMB anomalies, which we hereinafter refer to as cumulated surface mass balance anomalies (cSMBA).
The IMAU-FDM (Ligtenberg et al., 2011)   (a) Basins that we investigate in this study are colored, the focus basins (6, 13, 22) are indicated with dark yellow. The polar gap of ERS-2 and Envisat is indicated with a red circle. Basin borders and numbers from Zwally et al. (2012). (b-d) GRACE-derived mass time series, cumulated surface mass balance anomalies (cSMBA), and the difference between both (GRACE-cSMBA) for basins 6, 13, and 22. (e-g) Altimetry-derived volume time series and the volume change from the firn densification model (FDM), and the difference between both (Altimetry-FDM).The correlation coefficient shows the correlation between detrendend GRACE and cSMBA (ρ(GRACE*,cSMBA*)) as well as detrended ALT and FDM (ρ(ALT*,FDM*)). Table S1 summarizes the deterministic linear trends for GRACE, cSMBA, altimetry, and FDM-derived time series, the GIA correction applied to GRACE, correlation coefficients, and the basin area of all investigated basins. per ice sheet boundary. It is available with the same spatial resolution. During model creation the spin-up of the firn layer is achieved by looping over the time series from January 1979 to December 2016 until an equilibrium firn layer is established. This assumes that the 38-years reference period is representative for the centuries before 1979 during which the firn layer was formed. Ligtenberg et al. (2014) found that the present-day RACMO2 forcing is a realistic reference climate. We also assume that the firn does not exhibit a trend during the reference period itself, because we consider it an equilibrated firn layer due to the spin-up run (Ligtenberg et al., 2011).

GRACE Time Series
We use CSR RL06 monthly gravity field solutions (Bettadpur, 2018) from April 2002 to August 2016. Because of accelerometer issues after August 2016 (Loomis et al., 2020), we only include solutions up to this time. We compute degree-1 coefficients following Swenson et al. (2008) and Sun et al. (2016) to complement the monthly gravity field solutions. Further, we replace the C 20 coefficients with a product derived by satellite laser ranging (Loomis et al., 2019). We apply error and leakage optimized sensitivity kernels on the complemented gravity field time series to derive basin-wide mass changes . We correct for GIA by using the model product from Ivins et al. (2013). A residual GIA signal can still be present in the time series, because GIA models are subject to large uncertainties. The net GIA signal predicted by current models differ over a range of ∼40-80 Gt a −1 (Whitehouse et al., 2019). We do not attempt to take GIA errors into account in the state-space filtering approach. They will be reflected in the temporal mean of our results (Sections 2.4 and 2.5). In Table S1, we provide the GIA uncertainty and the spread for every drainage basin that is estimated from a GIA model ensemble. We propagate the GIA uncertainty to uncertainty estimates involving GRACE. The SI provides technical details.

Altimetry Time Series
Schröder et al. (2019a) provide a monthly sampled altimetry-derived elevation change product from a multi-mission analysis using a repeat altimetry approach. Temporal and spatial smoothing is applied during processing using a moving average of 3 months and a 10-km-sigma Gaussian smoother, respectively. Within the radius of 10 km there are at least four satellite ground tracks (depending on orbit geometry), each is measured at least twice in 3 months. This means that the monthly surface elevation change in a grid cell is based on an average of at least eight satellite passes.
For our investigation, we use data over the period from April 2002 to August 2016 corresponding to the availability of GRACE observations. The altimetry-derived elevation changes during this time period include observations from ERS-2, Envisat, ICESat, and CryoSat-2. We selected all drainage basins that are covered at monthly resolution over that period. This criterion requires coverage by ERS-2 and Envisat and excludes basins 1, 2, 3, 17, and 18 that extend beyond the 81.5°S limit of coverage of these two missions. Furthermore, we do not include basins where the altimetry data are of lower quality due to strong topographic gradients. This is the case for the basins in the region of the Antarctic Peninsula (25,26,27) and Victoria Land (15,16). For further details on uncertainties and altimetry quality limitations, see Schröder et al. (2019a) and Strößenreuther et al. (2020).
We correct for the elastic deformation of the solid earth and use the model product from Ivins et al. (2013) to correct for GIA. GIA and elastic-induced elevation changes are very small compared to elevation changes induced by SMB-driven and ID-driven changes. To illustrate the order of magnitude in case of basin 13, the integrated elastic-induced volume change and the GIA-related volume change are ∼0.4 and ∼0.4 km 3 a −1 , respectively, whereas the corrected altimetry observed volume change (Table S1) is −27.1 ± 1.3 km 3 a −1 from April 2002 to August 2016.

Long-Term Uncertainties in Model Products
If we examine long-term signals, we have to accept three limitations imposed by the data that we use: (a) Any long-term SMB-driven signal would violate the assumption that the mean SMB over the chosen reference period is a long-term mean SMB and would thereby introduce a SMB contribution to the mean rate in our analysis. (b) Any long-term firn-thickness change over the reference period will affect the mean trend of firn-induced elevation change. Thomas et al. (2017) showed that, especially, in the Antarctic Peninsula, there seems to have been a long-term increase in snowfall, by about 10%-15% over the past century. Ligtenberg et al. (2014) predicted an increase in firn air content of 150 km 3 a −1 over the 21st century for the AIS. Thus, the assumption that the firn is in steady-state at the start of the integration may not be valid everywhere. Furthermore, the firn densification model output is affected by errors in the input data, which include SMB variations apart from temperature variations. We expect that errors of firn-thickness change will be correlated with errors of SMB because the IMAU-FDM is forced by RACMO2 outputs. However, uncertainty estimates are not available. Data from firn cores could quantify this long-term uncertainty, but this is beyond the scope of this study. (c) Any uncertainty in the GIA trend directly propagates to the temporal mean of the estimated trend. These three long-term uncertainties are superimposed on the mean rate of mass and volume time series, but gravimetry and altimetry are affected by them differently. First, long-term trends in SMB induce volume trends that involve some effective firn density, which is generally smaller than ice density. Therefore, the volume effects of long-term trends in SMB are more pronounced than the mass effects. Second, any long-term trend of the firn air content is solely part of altimetry observations. Finally, an unaccounted GIA-signal predominantly affects the gravimetry observations due to the large effective density of GIA-induced deformation of the solid earth (Wahr et al., 2000).

Signal Separation Strategy of Basin Time Series
GRACE and altimetry detect long-term and short-term mass and volume changes, which we define as a rate varying on decadal scales and variability at sub-decadal scales, respectively. The model products account for a large part of the short-term variability. The differences between (a) GRACE and cSMBA (grace-csmba) and (b) altimetry and FDM (alt-fdm) include the ID-driven variability (Zwally et al., 2015) in addition to long-term errors (Section 2.4). Furthermore, the differences include a remaining short-term signal, due to modeling or observational errors.
We assume that the same ID-driven signal is present in grace-csmba and alt-fdm and that this ID-driven signal is most likely only long-term in Antarctica (Rignot et al., 2019). Thereby, we do not presuppose that the rate of the long-term signal is constant over time. Rather, we suppose that changes of this rate are slow. We parameterize the long-term signal through a trend with a time-variable rate. This is in contrast to deterministic approaches where, for example, a trend is used with a constant rate. To enable a coupled evaluation of the four data sets (1) we link the temporal variations of the rate of grace-csmba and altfdm with the density of ice (917 kg m −3 ) in the state space model framework. This density is a conceptual assumption and we assume it is free of errors.
(2) We model the long-term signal for both time series using the seemingly unrelated time series model (Commandeur & Koopman, 2007). This implies that the changes in the rate of the long-term parts of both time series are fully correlated. However, the mean rate may differ over the entire time series. In this way, we allow potential long-term errors (Section 2.4) to be reflected differently in both results.
We approach the remaining short-term signal of the grace-csmba and alt-fdm time series in three ways: (1) We use annual and semi-annual cycle components with time-variable phase offset and amplitude. (2) We use a residual component to model uncorrelated noise. (3) We assume that noncyclic short-term signals can be described as an integrated random walk starting at zero at the first observation epoch. We assume that unmodeled SMB can be represented with white noise. Therefore, the errors of cumulated SMB anomalies are cumulated white noise, that is, an integrated random walk. It can be modeled as an autoregressive process of the order one (AR(1) process) in which each time step is the previous one (AR-coefficient equals 1) plus a disturbance term (Section 2.6). King and Watson (2020) showed that Generalized Gauss Markov models are better suited than an AR(1) process to explain the full SMB-driven variability. However, here we presume an AR(1) process is suited for modeling the remaining short-term variability of differential time series, which contains observational errors beside the SMB-driven variability. We test which AR-coefficient between 0 and 1 best describes the auto-correlation length of the noncyclic short-term signals. This implies that any error that is not an unmodeled SMB signal but behaves like an autoregressive process is interpreted with the same AR(1) process. That is, the AR(1) process also absorbs correlated short-term GRACE and altimetry errors. Consequently, we do not necessarily expect that the same remaining short-term signal is present in grace-csmba and alt-fdm time series and we do not enforce a coupling of the remaining short-term variations.

State Space Setup
A state space model describes a time series, for every time t i , i = 1, …, n with [p × 1] observation vector y i , using the [m × 1] time-variable state vector α i of m model coefficients. We extend the state space model approach from Frederikse et al. (2016) to the simultaneous filtering of multiple time series, where in our case the dimension is p = 2 (bivariate): the first dimension is alt-fdm and the second is grace-csmba. We set up a state space model as (Durbin & Koopman, 2012): with trend μ i , cycle (seasonal harmonic) terms c i,j , AR(1) process ζ i and irregular term (residual) ϵ i , each a vector of size [p × 1], in total m = p ⋅ 7 = 14 model coefficients. Each of these terms is modeled as a time variable process, where at each epoch a stochastic disturbance term allows for time variability.
First, the trend is defined as: that is, dt i is generally close to 1 for an approximately regularly spaced time. As the trend μ is updated at each epoch with a time-variable rate ν-due to the disturbance term ξ i -a smoothly changing trend is possible. Second, cycle terms can be iteratively defined as (Harvey, 1990): where c i,j is the cycle component, and * , i j c an auxiliary cycle term that allows for the recursive description of the cycle term, Third, the autoregressive AR(1) process is modeled as (Harvey, 1990;Laine et al., 2013): where the [p × m] design matrix Z i relates the current state α i to the observations y i . Subsequent states are related by [m × m] transition matrix T i : and [m × m] variance matrix Q i that contains all disturbance variance-covariances. The SI contains an explicit description of all relevant vectors and matrices.

Estimation of the State and Disturbance Variance-Covariance
We use a Kalman filter and smoother to estimate the state α, irregular term  and their error variances for both time series simultaneously (Koopman, 1993). We co-estimate the uncertainties of the time-variable model coefficients such as of the trend and its rate. We describe technical details in section B of the SI.
As we do not have reliable knowledge about the disturbance and irregular variances and covariances, as well as the AR-coefficients ϕ, we estimate these parameters a priori. A common approach in state space modeling is to statistically optimize these parameters by maximizing the likelihood: the goodness of fit of the observations given a choice of (co)variance and AR parameters. Using the expectation maximization (EM)-algorithm, we iteratively estimate variance-covariance matrices (Koopman, 1993) while we do a two-dimensional grid search for an optimal ϕ v -ϕ m -combination. We enforce a ratio of the rate disturbance variances between grace-csmba In this way, we couple the common long-term variability of grace-csmba and alt-fdm with ice density, because we assume that this is predominately ID-driven (Section 2.5).

Mean Trend Rates
We compute the mean rate from the trend and propagate estimated time variable uncertainties for a mean rate uncertainty. The SI provides the formal mathematical description. For comparison to the results from state-space filtering, we compute deterministic results by estimating trends with constant rates of gracecsmba and alt-fdm. The deterministic rates are co-estimated with bias, annual and semi-annual cycle components using least-squares adjustment.

Results
Primarily, we focus on the long-term contributions to the mass balances of the drainage systems. We illustrate the basin-integrated time series of the three focus basins in Figures 1b-1g. Figure 2 shows the estimated components together with the original observation, and model time series. As it shows, the AR(1) and cycle terms absorb the largest part of the differential time series that is not explained by the trend. The residual, irregular term absorbs the uncorrelated noise of grace-csmba, but is negligible in case of altfdm. Figure 3 shows the grace-csmba and alt-fdm time series (observation minus model), the trend with time-variable rates and its uncertainty, and the deterministic trend of the focus basins. We observe that the uncertainties of the trend with time-variable rates vary among basins. High uncertainties are accompanied with a high auto-correlation of the estimated AR(1) component ( Figures S2 and S6) and reveal that the trend is less independent of the AR(1) process (Section 4.1). Figure 4 visualizes that the estimated time-variable rate from grace-csmba and alt-fdm is fully correlated as the time variability for both is identical, even though the mean rate may differ (Section 2.6).
When we compare between the mean rates of grace-csmba and alt-fdm for the focus basins 6, 13, and 22, we find differences of ∼2.1, ∼1.5, and ∼0.9 Gt a −1 , respectively, (we use a density of 917 kg m −3 to convert alt-fdm volume changes to mass changes). Figures 5a and 5d illustrate mean rates of grace-csmba and alt-fdm. Further we compare the mean rates to the deterministic results in Figure 5, and include the mean rate uncertainties (values are provided in Table S2). The absolute values of the mean rate are highest for the Amundsen Sea Embayment (basins 20, 21, and 22). In basins 19 and 24 the mean rate and the deterministic rate from grace-csmba have opposite signs to the mean rate and the deterministic rate from alt-fdm.
The root mean square (RMS) of the AR(1) component, RMS AR(1) , and the irregular component, RMS irr , reflect the magnitude of the signal that is not explained by the trend or cycle terms. Similarly-in the deterministic case-the RMS of post-fit residuals, RMS resid , represents the unexplained signal. Table S2 compares these values. In the case of alt-fdm a large part of the remaining short-term signal can be explained by an AR(1) process. In the case of grace-csmba there is a significant irregular component, whereas this is negligible in alt-fdm. The RMS of the post-fit residuals from the deterministic fit is highest in basins 22 and 14 in case of grace-csmba and alt-fdm, respectively.

Interpretation of the Results
The state space model is able to separate long-term variability shared by both grace-csmba and alt-fdm, and remaining short-term signals in these differential time series. We argue that the estimated long-term WILLEN ET AL. signal is approximately equivalent to the ID-driven time variable mass and volume changes. This has the limitation that the long-term signal may be affected by errors of the model products and the assumptions involved therein (Section 2.4). Particularly, the SMB anomalies and firn thickness anomalies considered in grace-csmba and alt-fdm imply assumptions on a steady-state mean SMB and firn structure as a reference. Consequently, a negative rate (positive rate, respectively) of the ID-driven contribution means that ice discharge is larger (smaller) than the mean SMB thus defined. Based on the trend uncertainties, we conclude that the separation performs best when the trend rates are either large (basins 21 and 22) or the unexplained, auto-correlated, signals have short correlation lengths (e.g., basins 7 and 11-14 in Figure S6). When the unexplained AR signals contain substantial long-term signals (large correlation length) either in grace-csmba or alt-fdm, it will be increasingly difficult to discriminate between the AR process and the trend. This is reflected by larger estimated uncertainties of the trend. The cross-correlation of the AR(1) between time series for grace-csmba and alt-fdm is small or negative ( Figure S6, except for basins 6 and 8), affirming that the common signals have been absorbed in the trend. Because the alt-fdm time series are more sensitive than grace-csmba to variability in SMB, due to the relatively low firn density, the grace-csmba is dominant in the definition of the trend.
Our estimated trends vary substantially with the region. We find a clear difference between the WAIS and the EAIS ( Figure 5). In particular, long-term signals are dominant over basins 21 and 22. An accelerated decrease of volume and mass is already visible in the basin integrals of the satellite observations (Figures 1d  and 1g). The accelerating long-term change can be fitted with the trend with time-variable rates and interpreted as an accelerating ID-driven signal.
Basin 6 shows strong inter-annual variations, with a significant step in the year 2009. This step can be attributed to an accumulation anomaly (as already present in the SMB model) (Boening et al., 2012;Lenaerts et al., 2013). The rate of the trend is almost constant over time. From this, we conclude that there is a positive long-term contribution to the mass balance, which only slightly changes over time. Similarly to basin 6, a positive trend with a low variability is visible in all investigated basins of the EAIS except for basins 10, 11, WILLEN ET AL.
10.1029/2020JF005966 9 of 16 13, and 14 ( Figures S2 and S3). In these basins, we interpret the positive trend as ice dynamic thickening. However, this positive trend may also have causes in the violation of the assumptions we made in Section 2.4, for example, that there is a long-term SMB-driven signal that is not part of the model products.
In basin 13, the trend is small. Here, the altimetry-derived and the FDM-derived volume time series are very similar, whereas the mass time series show a small long-term difference compared to the cSMBA values (Figure 4b). Explanations for this discrepancy include: altimetry missions have insufficiently sampled the mass change of Totten Glacier; there are unaccounted long-term climate signals that affect mass and volume changes differently (Section 2.4); the GIA correction of GRACE is underestimated. The spread of the GIA signal predicted by models is 3 Gt a −1 in basin 13, which is higher than the difference of the mean rates (1.5 Gt a −1 ).

Differences in Mean Rates
For both grace-csmba or alt-fdm the mean rates using the state space model agree well with the deterministic trends using least-squares adjustment. Mean rates of grace-csmba and alt-fdm agree within the indicated 2-σ-uncertainties for basins: 4-13 and 22 (Figure 6) using the state space model. For the deterministic approach the grace-csmba and alt-fdm rates match within the 2-σ-uncertainties for basins: 9-12 and 22. The state space model allows for a more realistic uncertainty estimation compared to unrealistic small standard uncertainties derived by least-squares adjustment. Further biases remain unconsidered. For example, GRACE detects mass changes only at a low spatial resolution and signal leakage may distort the results (Horwath & Dietrich, 2009). On the other hand, the altimetry product provides a high resolution but biases are present due to topographic characteristics or interpolation (Strößenreuther et al., 2020). Other factors that may contribute to the differences between the mean rate of grace-csmba and alt-fdm are due to GIA uncertainties and assumptions about long-term equilibrium in the SMB and firn models (Section 2.4). Our approach allows for different mean values to reflect those potential long-term uncertainties. For example, if we use a different reference period to calculate cSMBA, this would lead to a different mean rate-or in other words it would lead to a shift of the time-variable rates of the trend-and therefore create a different discrepancy between the mean rates from grace-csmba and alt-fdm. Due to limitations of the input datasets, we cannot conclude whether grace-csmba or alt-fdm results have a smaller (systematic) bias. Because we assume that the AR(1) component represents an accumulating error, we force it to start from zero (from which it may deviate by incorporating disturbance variance). The deterministic fit does not include this assumption and therefore it may lead to an underestimation of the residual and its correlation length. It cannot capture the auto-correlated signal around the trend (e.g., basins 4, 5, 6, and 9 in Figure S2).

Unexplained Short-Term Signals
Except for basins 9-11, the estimated short-term signals are small compared to the cSMBA and FDM time series ( Figure S1 and Table S2), which indicates that the cSMBA and FDM time series already explain a large part of the short-term variability of the GRACE and altimetry data. In basins 9-11, we find residual short-term components in the same order of magnitude as the time series of total mass and volume variations, especially in the volume time series ( Figure S1). Furthermore, the discrepancy between satellite observations and model products is evident by the low correlation between the time series here (Table S1). We find the lowest correlation of 0.12 in basin 10 between the altimetry and the FDM time series. The rates of the trends over time, their uncertainties (purple and orange), and the mean rates (blue and red) from grace-csmba and alt-fdm, respectively. The rate from alt-fdm is converted to mass change with ice density. grace-csmba uncertainties include GIA uncertainties. We provide results of all investigated basins in the SI.
The estimated AR(1) process absorbs two main parts. (1) It absorbs short-term SMB-driven signals that are not included in the model products. In this sense, a part of the AR(1) process can be understood as modeled short-term error of cSMBA and FDM. (2) Errors of the satellite data, or errors due to different temporal sampling of observations and model products, are absorbed in the AR(1) process, too. We find a larger AR(1) process in alt-fdm than in grace-csmba. A major reason for this is the temporal smoothing applied during the processing of altimetry (Schröder et al., 2019a). This artificially correlates temporally uncorrelated errors. The AR(1) process captures these correlated errors and this also explains the small irregular component in case of alt-fdm. We observe that the AR(1) process from alt-fdm generally contains relatively large signals on timescales of more than a year in basins 4-11. We also find larger residual seasonal signals in alt-fdm than grace-csmba ( Figures S1 and S4). If the seasonal signal and the AR(1) process are an unmodeled SMB signal, a low firn density could explain the comparatively large volume change and relatively small mass change. In order to enable the attribution of the error to either the satellite data or the model products, these errors would have to be parameterized separately. So far we have not identified a suitable parameterization for either errors in the model products, or errors in the satellite data. We have tried to constrain the errors of the model products by means of an error assessment analogous to Willen et al. (2020). Doing so, we made use of the SMB product of another regional climate model, namely MAR (Agosta et al., 2019). But the error estimation based on only two SMB products does not provide sufficient error information, resulting in an unrealistic error budget. Moreover, no FDM based on MAR is available.
The filtering results of basin 6 show that most of the short-term variations in grace-csmba can be explained by the irregular component and can thus be described as temporally uncorrelated Gaussian noise. The AR(1) process is very small and accounts for a minor trend (Figures 2a and S4). This can be prevented by choosing a slightly reduced irregular disturbance variance. This would result in an WILLEN ET AL.
10.1029/2020JF005966 11 of 16 Figure 5. The mean rates, the 1-σ-uncertainties, and differences to deterministic rate are color coded. Results from alt-fdm are converted to mass with ice density (917 kg m −3 ). grace-csmba uncertainties include GIA uncertainties. Table S2 provides the underlying numbers.  . Table S2). (b) Comparison of mean trend rates to ID-driven estimates published by Zwally et al. (2015). (c) Comparison of mean trend rates to the mass balance estimates published by Martín-Español et al. (2016). To enable the comparison, we removed the mean SMB anomalies during the indicated time period from the published mass balance estimates. We used basin numbers from Zwally et al. (2012) (first row of labels), basin numbers from Martín-Español et al. (2016) indicated in the second row of labels. The sum of basins 5 and 6 approximately equals the sum of basins 305 and 306 and the sum of basins 9 and 11 approximately equals basin 310. The time periods are indicated in the subfigure titles. The mean rates from alt-fdm are converted to mass change with ice density (917 kg m −3 ). 2-σ-uncertainties are indicated with error bars. gracecsmba uncertainties include GIA uncertainties.
increased part of the short-term variations being explained by the AR(1) process, at the expense of a lower likelihood.
Furthermore, we have investigated the sensitivity of the approach to an alternative GRACE product based on monthly GRACE solutions from Mayer-Gürr et al. (2016). A clear difference in the estimated AR(1) process is visible for most basins ( Figure S7). From this we conclude that a large amount of the AR(1) process in the grace-csmba time series can be explained by errors due to GRACE data processing, which needs further investigation (Section 4.5). Our results suggest that considerable signals are present in the grace-csmba and alt-fdm time series that are not correlated for the two time series and cannot be attributed to long-term changes. -Español et al. (2016) estimate annual ID-driven and SMB-driven mass changes along with linear GIA from altimetry, GRACE, and GNSS. A priori information is included from models (among others RACMO2.3p2 outputs) and further observations to constrain the spatio-temporal characteristics of ID, SMB, and GIA. Table 2 2016) within 2-σ-uncertainties. This is also the case comparing grace-csmba results, except for basins 4 and 22 where our mass rate from grace-csmba is slightly more positive and more negative, respectively. The temporal evolution of the estimated long-term signal (Figures 4b and 4c) and the ID-driven signal from Martín-Español et al. (2016) (Figures 4c and 4a therein, respectively) are very similar. Zwally et al. (2015) use ICESat data and global reanalysis products to separate SMB-driven and ID-driven signals from October 2003 to December 2008 in altimetry data. Our mean rates of the trends during this time period agree well with the ID-driven mass changes published by Zwally et al. (2015) (Figure 6b). If we compare mean rates over this time period from alt-fdm and Zwally et al. (2015) they agree within the uncertainties except for basin 19. We find the largest differences (>10 but <15 Gt a −1 ) for basins 10, 12, 14, and 23. The smallest discrepancies (<1 Gt a −1 ) arise for basins 5, 11, and 20. It is remarkable that even the results with the largest differences agree within the uncertainties. Zwally et al. (2015) use a different methodological framework and use laser altimetry observations (ICESat) only, whereas we use an altimetry product (Section 2.3) that includes radar altimetry (Envisat) in addition to ICESat laser altimetry and involves a calibration of ICESat laser operation period biases (Schröder et al., 2017) that differs from that by Zwally et al. (2015). Note that , we do not investigate all drainage basins of the EAIS (Figure 1), where Zwally et al. (2015) identified large positive ID-driven mass changes, which are however under debate (Martín-Español et al., 2017;Richter et al., 2016).

Martín
High loss rates and an ID-driven acceleration are known in the Amundsen Sea Embayment of the WAIS (Shepherd et al., 2018). Therefore, we discuss this region in more detail. Basin 21 (includes Thwaites Glacier) and basin 22 (includes Pine Island Glacier) belong to this region. Martín-Español et al. (2016) published mass changes of −64.8 ± 4.5 and −37.4 ± 3.5 Gt a −1 for basins 21 and 22 (equal to 321 and 322 from Martín-Español et al. (2016)), respectively during January 2003 until December 2013. SMB anomalies during this time period contribute −2 and −1 Gt a −1 to the mass balance, respectively (removed in Figure 6c). The results from Martín-Español et al. (2016) and Zwally et al. (2015) are similar to the estimated mean rates from grace-csmba and alt-fdm ( Figure 6). Furthermore, our results in the Amundsen Sea Embayment agree well with ice discharge estimates from Rignot et al. (2019), which they estimated from ice velocity and ice thickness data. The drainage systems defined by Rignot et al. (2019) differ to those that we use. (1) Basin 21 corresponds approximately to the aggregation of Thwaites Glacier, Haynes Glacier, Crosson Glacier, and Dotson Glacier; (2) basin 22 corresponds to Pine Island Glacier (Table 1 in Rignot et al., 2019). For the time period 2009-2017, they quantified the ice discharge of (1) and (2) at 190.7 ± 4.7 and 133.2 ± 5.8 Gt a −1 , respectively. The absolute ice discharge values from our results (sum of the mean SMB and the estimated mean rate) between January 2009 and August 2016 from grace-csmba for basin 21 and basin 22 are 176 ± 2 and 136 ± 2 Gt a −1 , respectively. alt-fdm results are 186 ± 2 and 135 ± 2 Gt a −1 , respectively.

Outlook
During the selected time period of 14 years and 4 months covering the available GRACE observations, we have not investigated all Antarctic drainage basins due to the limited quality of the altimetry data (Section 2.3). High-quality altimetry data of the almost entire AIS is available from the CryoSat-2 mission and ICESat-2 mission since July 2010 and October 2018, respectively. Gravity fields from the GRACE-FO mission are available since June 2018. Through the continuous observation of the AIS it will become possible to investigate further areas of Antarctica within the near future. Moreover, an error model for the satellite data and the model products is needed to improve the estimation of the time-variable signals. An error model should account for the temporal heteroscedasticity (nonconstant variability of errors), especially of the observation products. To properly understand the mechanisms that are responsible for ice volume and mass changes we argue that more attention should be focused on the temporal behavior of the observed mass and volume changes. The state space method is a likely candidate for future temporal analyses of ice mass changes. In particular, the approach can unravel interannual signals and aid to overcome limitations due to deterministic methods pointed out, for example, in Horwath et al. (2012) and Mémin et al. (2015). Our estimated AR process and time-variable (seasonal) cycle time series can be instrumental in assessing unmodeled SMB-driven signals.

Conclusions
Our data-driven approach is able to estimate common time-variable signals in both geodetic data sets with the aid of products from regional climate modeling and firn modeling. We interpret these common signals as most likely ID driven. We find residual auto-correlated and seasonal signals in these time series, with often significant variance if compared to the climate and firn model products. However, we cannot yet attribute the short-term signals to a specific source.
Furthermore, our results confirm the accelerating ice loss of the WAIS in the Amundsen Sea Embayment. Our approach allows us to fit the acceleration without the artificial selection of time periods. The results for the investigated basins of the EAIS show small long-term signals with a low temporal variability. The temporal variability of mass and volume changes of the EAIS can mainly be attributed to the SMB component. Residual short-term signals are most likely not ID driven because these signals are not positively correlated between the GRACE and altimetry time series.
Limitations of the presented approach are due to the treatment of short-term and long-term errors. So far we have not been able to assign uncertainties to the input data sets. The estimated mean rate is sensitive to long-term uncertainties of the SMB and the FDM product, for example, by the chosen reference period. However, we do not expect large errors in the mean rates of the model products for the Antarctic drainage basins.