GPS Rates of Vertical Bedrock Motion Suggest Late Holocene Ice‐Sheet Readvance in a Critical Sector of East Antarctica

We investigate present‐day bedrock vertical motion using new Global Positioning System (GPS) timeseries from the Totten‐Denman glacier region of East Antarctica (∼77°–120°E) where models of glacial isostatic adjustment (GIA) disagree, glaciers are likely losing mass, and few data constraints on GIA exist. We show that varying surface mass balance loading (SMBL) is a dominant signal, contributing random‐walk‐like noise to GPS timeseries across Antarctica. In the study region, it induces site velocity biases of up to ∼+1 mm/yr over 2010–2020. After correcting for SMBL displacement and GPS common mode error, subsidence is evident at all sites aside from the Totten Glacier region where uplift is ∼1.5 mm/yr. Uplift near the Totten Glacier is consistent with late Holocene ice retreat while the widespread subsidence further west suggests possible late Holocene readvance of the region’s ice sheet, in broad agreement with limited glacial geological data and highlighting the need for sampling beneath the current ice sheet.


• We present new Global Positioning
System estimates of vertical bedrock motion in the Totten-Denman glacier region of East Antarctica • Surface mass balance loading displacements are a dominant signal across Antarctica and here contribute ∼+1 mm/yr of velocity bias (2010-2020) • A significant portion of this sector of East Antarctica is undergoing long-term subsidence, consistent with Holocene ice-sheet advance

Supporting Information:
Supporting Information may be found in the online version of this article.

Surface Mass Balance Loading Displacements (1979-2020)
Obtaining accurate and precise GPS velocities or timeseries that reflect long-term bedrock motion requires correction for elastic deformation of the solid Earth due to contemporaneous surface loading changes; in Antarctica largely atmospheric and snow/ice surface mass variations Thomas et al., 2011). Characterization of ice loading displacement is challenging due to uncertainties in input data sets (such as converting satellite altimeter ice elevation change to mass) or circularity in reasoning (using satellite gravimetry data which themselves depend on GIA models). With few exceptions (e.g., Koulali et al., 2022; these corrections have tended to be in the form of multi-year or multi-decadal linear rates of displacement (e.g., Argus et al., 2014;Thomas et al., 2011;Wolstencroft et al., 2015) with a focus on regions of substantial dynamic glacier change (Barletta et al., 2018;Nield et al., 2014;Samrat et al., 2020).
However, SMB time series exhibit a power-law like character in Antarctica over years to centuries , suggesting that subsequent displacements may produce biases in bedrock velocities and affect the interpretation of GIA (Koulali et al., 2022).
We model SMBL using SMB anomalies based on the regional model MAR v3.11 (35 km horizontal resolution; Kittel et al., 2021) and, for comparison, RACMO2.3p2 (27 km horizontal resolution; van Wessem et al., 2018). At the time of analysis MAR outputs spanned the period January 1979 to September 2020 and RACMO outputs spanned the period January 1979 to August 2018, and we work with monthly model outputs in both cases. We adopt MAR as the primary data set due to its extended period. To compute SMB mass timeseries, we first interpolated the input grids to a consistent 5 km resolution and computed SMB anomalies to the mean rate by subtracting the mean SMB per grid cell over the respective full data periods. Using SMB anomalies is necessary as SMB models, unlike hydrological models, only consider the mass change at the top surface of the ice sheet and do not consider mass discharge due to glacier flow. Using anomalies assumes the long-term inputs are entirely balanced by equivalent outputs; we address this further below. We then applied a consistent land mask, defined as the limit of the grounded ice sheet (Haran et al., 2005), cumulatively summed the timeseries, and converted to the mass to equivalent cylinders of 2.5 km radius.
We used these mass timeseries as inputs into the REAR v1.0 software (Melini et al., 2015) to compute elastic loading displacements in a center-of-solid Earth (CE) reference frame. We used the AK135 continental Earth model (Kennett et al., 1995) represented by 89-layers (2 km depth layers for the top 60 km then 50 km layers below), and computed elastic load Love numbers using the E-CL0V3RS software to degree 35,000 (Barletta et al., 2006). Using alternative Earth models produced differences <∼0.5 mm ( Figure S1 in Supporting Information S1); AK135 was chosen due to its similarity in this region with the recent three-dimensional model of Lloyd et al. (2020). We interpolated the model outputs to the GPS daily time series epochs.

Altimeter Timeseries (2003-2019)
For comparison, we also computed timeseries of elastic displacement based on satellite altimeter ice sheet volume change data (Shepherd et al., 2019). These data have the advantage of including the effects of both SMB and ice discharge but introduce systematic uncertainty, notably the conversion from volume to mass. We commenced with a merged CryoSat-2 and Envisat 5 km gridded elevation change product with 140-day sampling over January 2003 to February 2019 that covers the East Antarctic Ice Sheet (EAIS; Shepherd et al., 2019). Biases between Envisat and CryoSat-2 were already estimated and removed on a cell-by-cell basis to produce continuous timeseries. We applied the same land mask as for the SMB models and filled gaps using an efficient method based on discrete cosine transforms (Garcia, 2010;Wang et al., 2012). Volume changes were converted to mass changes using the simple approach of Shepherd et al. (2019), with elevation changes having the density of snow (350 kg/m 3 ) or ice (900 kg/m 3 ) depending on a pre-defined mask. In this region, only the Totten Glacier was within the ice mask, although we test applying a density of ice also for the Denman Glacier. We used the same approach for computing elastic displacements from the altimeter derived data set as for the SMBL.

GPS Timeseries (2006-2020)
We focus on six GPS sites in the sector of East Antarctica in the region of the Totten and Denman glaciers shown in Figure 1. Two of the sites are long-running International GNSS Service sites (DAV1 and CAS1) and one is an extended record of a site deployed since December 2006 (BHIL). Three of the sites are new and were deployed in December 2015 (CAD4, CAD5) or December 2016 (CAD6). We considered data spanning December 2006 to January 2019 (the start date was chosen to align with the start of the BHIL timeseries). Further site details may be found in Table S1 in Supporting Information S1. Opportunity for GPS deployment on bedrock is limited in this region, especially away from the coast, as indicated by the distribution of rock outcrops shown in brown in Figure 1. As such, the newly acquired data from CAD4, CAD5, and CAD6, and extended data from BHIL, CAS1, and DAV1, add rare and important constraints on GIA models.
We analyzed raw GPS data using GIPSY v6.4 software and NASA Jet Propulsion Laboratory orbit and clock products (repro3) in the IGS14 realization of ITRF2014 (Altamimi et al., 2016). We removed the effects of atmospheric pressure loading (ATML) from the time series using displacements in the center-of-figure (CF) frame based on the MERRA-2 atmospheric model available at massloading.net. GPS time series are shown in Figure 1. GPS site locations and bedrock uplift rates in the study region. GPS sites are labeled and estimated uplift rates are shown as circles. GPS velocities were determined following time series correction for common mode error and surface mass balance loading changes. The background shading is the predicted uplift rates from conventional forward (left) and empirical (right) glacial isostatic adjustment models. Bedrock locations are shown in brown and the ice grounding line and ice shelf fronts in black. All model predictions and GPS estimates are in a common mode frame (see text). Figure S2 in Supporting Information S1. To provide context, we also compute elastic loading due to surface loading (ATML, SMBL, and altimetry) at 259 existing or previous GPS sites across the entirety of Antarctica and its offshore islands.
We removed common mode error (CME) from the GPS timeseries using methods described in Text S1 in Supporting Information S1. In brief, we made use of the only long-running and stable GPS site in this region (DAV1) noting that CAS1 is subject to time-variable dynamical changes in the Totten Glacier. We decomposed the DAV1 timeseries using complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN; Torres et al., 2011) after removing SMBL and ATML and correcting for offsets related to hardware changes. We took the non-linear, low frequency components of the decomposition and subtracted that as CME from all timeseries. We return to the CME in the discussion. The improved linearity of BHIL timeseries after applying CME (and SMBL) at BHIL is clear in Figure S4 in Supporting Information S1.
To compute GPS site velocities and their uncertainties ( Figure S3 in Supporting Information S1) we used HECTOR v1.9 software (Bos et al., 2013) as described in Text S1 in Supporting Information S1.
The seven different predictions of present-day GIA uplift are shown in Figure 1, including the conventional forward models (left column) and empirical models (right column). Models not in a CM frame were shifted from a CE frame by adding −0.38 mm/yr based on the estimate of Sun and Riva (2020); this correction has a formal uncertainty (2-sigma) of 0.09 mm/yr. Of the models, IJ05R2, ICE6G_D and W12 were originally in CE. We note that there may be some inconsistency between the different realizations of CM, notably for G14 where the rates are relative to a region of low precipitation in central East Antarctica.
Comparison of the different GIA predictions ( Figure 1) reveals variation between models in terms of the mean uplift rate which exceeds 5 mm/yr, but also distinct differences in spatial patterns within this region. For the forward models this reflects a lack of Holocene ice history data from this region and uncertainty in the Earth rheology, while for the data-informed models it reflects, in part, differences in methodology and data corrections such as firn densification (Whitehouse et al., 2019). The REGINA model shows higher levels of variability, presumably associated with its basis in spherical harmonics and with little data to constrain that solution in this region; hence we focus mostly on the other models hereon.

Surface Mass Balance Loading Displacements
Timeseries of modeled SMBL displacement are plotted in Figure 2 (top) for illustrative sites CAS1 (green), DAV1 (blue), and BHIL (cyan) and they show ranges of up to 12 mm at CAS1 and BHIL but are smaller at DAV1. For context we also show modeled SMBL for a random subset of sites from across Antarctica (gray lines), which show substantial variation over all periods, ranging up to ∼±15 mm. Periods of rapid change in displacement are evident from time-to-time, and these are especially associated with large accumulation events; the most rapid changes will not be evenly distributed around zero as large negative SMB events are not yet common in Antarctica.
The respective power-spectra are shown in Figure 2 (middle) and reveal that the timeseries have a power-law like character, with an approximately linear increase in power with period in log-log space, although tapering off at long periods for DAV1. The power spectra for SMBL displacement at CAS1 and BHIL and the median power spectrum computed across all Antarctic sites (magenta) has a slope close to that of random walk noise (spectral index of −2). The median spectrum for the RACMO-based timeseries is like that of MAR but with slightly reduced power at the longest periods (closer to pure random walk).
Figure 2 (lower) shows the velocity biases that would be introduced if the SMBL displacements in Figure 2 (upper) were not corrected, as a function of time series duration, obtained by sub-sampling the modeled displacements with steps of 1 year. Biases at CAS1 and BHIL can exceed 1 mm/yr, and this includes the data period considered here where mean biases of 1 mm/yr are sustained over ∼2010-2020, and faster since ∼2016 (Figure 2, top).

GPS Velocities in the Denman-Totten Glacier Region
The GPS site vertical velocities, computed after correcting for ATML, SMBL (MAR), and CME, are shown in Figure 1, plotted over the different GIA predictions. Our estimates of bedrock motion are universally low, with highest uplift rates of ∼1.5 mm/yr near Totten Glacier (CAS1 and CAD6) or subsidence of ∼0.5-1.7 mm/yr around the Denman Glacier (CAD4, BHIL, and CAD5) and DAV1. We show these and their uncertainties in Figure 3 (green circles). We also show GPS velocities using time series corrected using one of the two different altimetry-based elastic corrections instead of the SMBL model: using the density of snow for Denman Glacier (gray circles) or the density of ice (light blue circles). While differences are evident between the different velocity solutions, the broad pattern remains the same. We adopt the SMBL-based GPS velocities as the primary solution hereon due to the presence of space-time noise in the altimeter data and uncertainties in converting altimeter volume change to mass, and with the knowledge that any unmodeled ice dynamic loss will shift the corrected GPS velocities toward the negative.
The weighted root-mean-square (WRMS) and weighted mean differences between the observations and models are given in Table S2 in Supporting Information S1. The best statistical agreement using WRMS is with W12, closely followed by IJ05R2. All models other than W12 have most of their variation orientated north-south and hence either match the signal in the eastern or western sectors but not both. The weighted-mean difference reveals the consistent bias in uplift rates between models and the velocities of ∼1 mm/ yr for IJ05R2 and W12 and up to ∼3 mm/yr for the other models. Including the C18 and RATES uncertainties noticeably reduces the WRMS and weighted-mean, with RATES values approaching those of the best forward models.
We also illustrate in Figure 3 the effect of not applying the CME correction or any form of ice mass loading model. The GPS velocities without either correction are shown as red circles with corresponding uncertainties (see also Figure S3 in Supporting Information S1). Aside from DAV1, these show rates of uplift of 2.5-3.5 mm/yr, generally above the prediction from all models but closest to G14. Removing CME from DAV1 (or any other correction) makes almost no difference to the estimated velocity there, with clear and consistent subsidence of ∼1.0 mm/yr that is well below all model predictions. However, at the other sites, the impact of removing the CME is to reduce uplift rates by ∼0.5-2.0 mm/yr (Figure 3, purple circles).
The effect of modeling and removing ice mass loading is to reduce the velocities further, reflecting the negative SMB anomaly over much of the data period (Figure 2, top). The underlying time series effects of the different modeling approaches is illustrated in Figure S4 in Supporting Information S1 for BHIL which shows, without these two corrections, substantial Velocity bias (mm/yr) 2 3 4 5 6 7 8 9 10 11 12 Span (Years) non-linearity. Both SMBL and altimetry solutions indicate BHIL subsidence after the application of these corrections, although the timeseries is more linear in the former case. Using the altimetry-derived loading that uses the density of ice for the Denman Glacier change further increases the rate of subsidence at BHIL and, especially, CAD4 (Figure 3) with negligible effect at other sites.
In all cases, the velocities based on timeseries with CME and SMBL (MAR) displacements modeled have either the smallest or second smallest velocity uncertainty (Table S2 in Supporting Information S1) suggesting that important timeseries noise is being removed in these solutions. Comparing modeled displacements based on MAR and RACMO shows good agreement since 2010 but differences before then ( Figure S4 in Supporting Information S1), and so this introduces some uncertainty into solutions based on SMBL corrections.

Discussion
Our finding of long-term subsidence in the sector of coastal East Antarctica west of Totten Glacier builds on previous studies that have analyzed earlier records from DAV1 and BHIL, and MAW1 650 km to the west of DAV1, indicating subsidence or marginal amounts of uplift along this coastline (Argus et al., 2014;Hammond et al., 2021;Thomas et al., 2011).  suggested uplift of ∼1.6 mm/ yr at DAV1 over the shorter period 2010-2014 (inclusive) but analysis of longer series (2006-2019) does not support this ( Figure S2 in Supporting Information S1), perhaps because of short-period SMBL variability. Our results provide clear evidence for subsidence along this coastline with both DAV1 and BHIL showing statistically significant subsidence, even before CME and ice mass loading corrections, which persist after these corrections are applied. While not the focus of our study, MAW1 shows clear subsidence of ∼1.2-1.6 mm/yr in our analysis (regardless of the correction choices; Figure S2 in Supporting Information S1). The addition of CAD4 and CAD5 provide further evidence that subsidence is occurring along substantial transects of this coastline once their series are corrected for CME and ice surface loading. The CAD4 finding is significant as it is the first inland site, ∼150 km from the coastline, for which a velocity is available in this sector of East Antarctica.
Our findings depend, to some extent, on the CME correction employed. CME is commonly removed from GPS timeseries (e.g., , including those from Antarctica (W. , but its computation typically uses a network of sites over a region spanning hundreds to thousands of kilometers across the region of interest. Applying a CME correction removes common error that is assumed to be systematic and of non-geophysical origin. Our analysis demonstrates coherence between DAV1 and BHIL over ∼1,000 km after removal of ATML and where very different SMBL displacements exist. This demonstrates that CME is a real and important signal to consider in precise point positioning (PPP) solutions in Antarctica, especially for short or patchy time series.
We note that while the GPS velocities are subject to reference frame uncertainty, it is unlikely to explain the observed subsidence, given the uncertainty of the origin of ITRF2014 maps to ±0.31 mm/yr (1-sigma) in the vertical component in this region (Riddell et al., 2017). Our comparison of GPS and model velocities also relies, in some cases, on the geocentre (CE-CM) translation we adopt. Here we use the values of Sun and Riva (2020) but geocentre estimates are uncertain and technique specific. Of the few published estimates, some of the largest are those of Métivier et al. (2020) which would reduce the forward-modeled uplift rates by a further 0.58 mm/yr from the value we used (the total correction would be 0.96 ± 0.12 mm/yr [1-sigma]). However, this would still not explain the observed differences entirely and would not affect the comparison to empirical solutions already in a realization of CM. As such, it is unlikely that differences in frame origins can resolve the differences we see.
Our results, therefore, suggest that present-day subsidence may not be unusual along this ∼1,800 km stretch of coastline in East Antarctica. In the Denman Glacier region this disagrees with the GIA models examined and runs contrary to an understanding of monotonic retreat of the ice sheet in this region from the Last Glacial Maximum (LGM) to present-day that is embedded within forward models. Regional paleoenvironmental records evidence Holocene ice margin fluctuations that are consistent with the GPS-derived pattern of modern-day subsidence and uplift. At the Windmill Islands, near Casey (CAS1), moraines record a retreat of nearby Law Dome by ∼1 km during the late Holocene (Goodwin, 1993). In contrast, at Bunger Hills (BHIL) and in the Davis (DAV1) region, Holocene moraines are restricted to at or near the modern ice margin (Colhoun & Adamson, 1992;White et al., 2009). Further, epishelf lake records at Bunger Hills suggest late Holocene ice re-advance (Berg et al., 2020) that may be driving present-day subsidence. In the Davis region, weathered clasts in marginal shear moraines at Vestfold Hills (Gore et al., 1994) indicate exposure and weathering of the landscape inboard of the present-day position followed by ice sheet re-advance. The timing of ice sheet re-advance is less clear, but may have been co-incident with a local late Holocene cool period recorded in sediments in marine inlets at Vestfold Hills (McMinn et al., 2001) and at Rauer Group (Berg et al., 2010). Geological evidence of relative sea level (RSL) at Vestfold and Bunger Hills has commonly been interpreted as indicating Holocene uplift, but is not inconsistent with late-Holocene subsidence. While there is a clear picture of relative sea level fall from an early Holocene RSL highstand of ∼10 m above sea level (e.g., Colhoun & Adamson, 1992;Zwartz et al., 1998), raised beach deposits and isolation basins only record RSL fall during the mid-Holocene, and there are no records that conclusively place RSL above the present day during the past 1000-2000 years. RSL records below sea level have not been thoroughly examined, but freshwater sediments deposited in Izvilistaya Inlet at Bunger Hills (sill height 3 m below sea level, Berg et al., 2020) are consistent with a RSL lowstand since ∼1.7 ka BP, indicative of limited late Holocene subsidence.
The late Holocene ice readvance advocated for here is not included in forward models of Antarctic GIA. This potentially missing ice history will not, however, contribute to the disagreement with purely empirical solutions of RATES, REGINA, and G14. That each of these solutions shows very different signal from one another and from the forward models suggests that the techniques themselves require further development in this region, although we note that in other places in Antarctica there is very good agreement between some of these solutions and GPS uplift rates Wolstencroft et al., 2015). This regional disparity suggests that localized signal is causing an issue in this region-perhaps because of SMB variability in this region and its impact on measurements of ice sheet volume change used in empirical solutions.
Lower rates of uplift than predicted by published models of GIA, or indeed subsidence, will alter estimated rates of mass change based on GRACE and GRACE Follow-On data. Notably, subsidence instead of uplift in this region will generally change the sign of the GIA correction for GRACE and lead to reduced estimates of mass loss (or greater gain), but not alter the non-linear parts of the timeseries. The hypothesis that increased ice loading during the late Holocene is the origin of the subsidence can be tested with new forward GIA models. Until those experiments are performed, GRACE analysis focusing on Totten or Denman glacier systems should use the W12 or IJ05R2 GIA models but with caution noting the biases present.

Conclusions
We report on a new GPS data set in a critical sector of East Antarctica, surrounding the Totten and Denman glacier systems, which have few geodetic or geological constraints on bedrock uplift or ice history, as are relevant to informing or validating models of glacial isostatic adjustment.
We first demonstrate the importance of correcting GPS time series for elastic deformation associated with fluctuations in surface mass balance, with effects in this region reaching +1 mm/yr over 10 years (see also Koulali et al., 2022). At GPS sites across Antarctica the modeled surface mass balance loading displacements show an approximately random-walk spectral signature, with displacements of several mm over a few months and sustained displacement rates over several years. This is an important effect to consider at all GPS sites in Antarctica (we discuss this further in Supporting Information S1).
We then show that, after correction for GPS common mode error and the effects of ice mass loading changes, the bedrock is most likely subsiding in the region from ∼77° to 107°E, including the Denman Glacier system. This is contrary to all current GIA models but is consistent with unmodeled Holocene readvance of the ice sheet in this region. Around the Totten Glacier we find modest uplift of ∼1.5 mm/yr, consistent with an understanding of late Holocene ice sheet retreat in this region.
Our analysis therefore provides new and important constraints on models of glacial isostatic adjustment, and we add new geodetic evidence to the hypothesis that much of this region may have experienced ice sheet readvance during the late Holocene.