Major Modes of Climate Variability Dominate Nonlinear Antarctic Ice‐Sheet Elevation Changes 2002–2020

We explore the links between elevation variability of the Antarctic Ice Sheet (AIS) and large‐scale climate modes. Using multiple linear regression, we quantify the time‐cumulative effects of El Niño Southern Oscillation (ENSO) and the Southern Annular Mode (SAM) on gridded AIS elevations. Cumulative ENSO and SAM explain a median of 29% of the partial variance and up to 85% in some coastal areas. After spatial smoothing, these signals have high spatial correlation with those from GRACE gravimetry (r∼ = 0.65 each). Much of the signal is removed by a firn densification model but inter‐model differences exist especially for ENSO. At the lower parts of the Thwaites and Pine Island glaciers, near their grounding line, we find the Amundsen Sea Low (ASL) explains ∼90% of the observed elevation variability. There, modeled firn effects explain only a small fraction of the variability, suggesting significant height changes could be a response to climatological ice‐dynamics.


Introduction
Observations of the changing volume of the Antarctic Ice Sheet play a major role in understanding ice-sheet change (e.g., Otosaka et al., 2023;Shepherd et al., 2012) from the whole-of-ice-sheet down to individual glaciers (e.g., Smith et al., 2020;Wingham et al., 2009).The now three-decade record of continuous ice volume change captures the variability and longer-term change of both surface mass balance (SMB), and related firn processes, and elevation effects of changing ice dynamics.These changes are, respectively, related to atmospheric and oceanic processes (Horwath et al., 2012;Smith et al., 2020).Several studies have examined the relationship between ice height changes and modes of climate variability, in particular linking them to both El Niño-Southern Oscillation (ENSO) and the Antarctic Circumpolar Wave (Kaitheri et al., 2021;Mémin et al., 2014Mémin et al., , 2015)).
Less studied in this context is the role of the dominant mode of climate variability in the Southern Hemisphere, the Southern Annular Mode (SAM).Despite SAM driving variability and trends in SMB across a wide range of timescales (Diener et al., 2021;Medley & Thomas, 2019;van den Broeke & van Lipzig, 2017), SAM has yet to be linked to observations of ice sheet elevation change, with one related study reporting no correlation to estimates of ice shelf elevation change (Paolo et al., 2018).By contrast, the cumulative sum of SAM (Diener et al., 2021;Kim et al., 2020) has recently been shown to be linearly related to the dominant signal in detrended surface mass time series derived from satellite gravimetry (King et al., 2023b), with large-scale spatially coherent signal across coastal regions at decadal timescales.
The ∼300 km spatial resolution of satellite gravimetry, combined with uncertainties in models of SMB (Mottram et al., 2021), meant that King et al. (2023aKing et al. ( , 2023b) ) were not able to separate the relative contributions of SMB and ice dynamical change forced respectively by the atmosphere and ocean (Hansen et al., 2021;Kim et al., 2020;Palóczy et al., 2018;Spence et al., 2017;Thomas et al., 2017;Verfaillie et al., 2022).In particular, ice dynamical change will have a distinct spatial pattern compared to SMB that is not detectable by GRACE but could be possible with altimetry (Smith et al., 2020).Detecting (or otherwise) a response of the grounded ice sheet to largescale climate variability via the oceans and ice shelves would provide important insights into ice-sheet sensitivity to climate change.
In this paper we analyze a recent gridded compilation of satellite altimeter data and compare these time series to cumulative climate indices.We compare the derived signals to those from space gravimetry and then, taking advantage of the high-resolution altimeter data, explore the signal over key ice streams: Thwaites, Pine Island, Totten, and Denman.

Altimeter Data Set
We make use of a gridded altimeter product (Nilsson et al., 2023) at 1,920 m horizontal resolution and covering the period from April 1985 to December 2020 (Nilsson et al., 2022).We spatially down-sample this to 5 km horizontal resolution.To facilitate comparison with space gravimetry data we only make use of data from 2002 to the end of the record.The data set contains monthly ice-sheet elevation-change data derived from a range of radar and laser altimeter missions; over the study period these are ERS-2, Envisat and CryoSat-2 and ICESat and ICESat-2.The approach to accounting for differences in reflection surfaces and other systematic effects is described by Nilsson et al. (2022).To reduce spatial noise, we apply a Gaussian smoother with widths specified in later sections.For these, width is defined at the half height of the Gaussian function, consistent with the definition commonly used in GRACE data smoothing (Wahr et al., 1998).

Space Gravimetry Data Set
We use the COST-G RL01 Level-3 50 km gridded ice-mass change per surface area GRACE and GRACE-FO V0002 data set obtained from http://gravis.gfz-potsdam.de/antarctica(Sasgen et al., 2020).We make use of data from March 2002 to December 2020, with the end point chosen to match the end of the altimetry data set.The data are spaced approximately monthly and with a data gap of ∼12 months between GRACE and GRACE-FO from mid-2017 to mid-2018.
We note that while this product is gridded at 50 km, the intrinsic GRACE resolution is ∼300-400 km.Postprocessing steps include replacement of low-degree GRACE coefficients and insertion of degree-1 terms using standard approaches (Dahle & Murböck, 2020;Sasgen et al., 2020).
Since we are interested in decadal variability and trends, we also lightly smooth the altimetry and GRACE data with a Gaussian filter with width 7 months (Wahr et al., 1998).

Climate Indices
We compare the altimeter and GRACE data primarily with SAM and ENSO indices, with additional comparison to Amundsen Sea Low (ASL) indices in the Amundsen Sea region.For the ASL indices, we make use of both the absolute center pressure (ASLP) and longitude (ASLλ) within the ASL Index version 3.20210820-era5 based on monthly ERA5 reanalysis data (Hosking et al., 2016).For SAM, we make use of the Marshall station index (Marshall, 2003).For ENSO, we make use of the Nino3.4index based on the HadISST1 data set (Rayner et al., 2003) and use a 6-month lag (King et al., 2023b;Paolo et al., 2018).We normalized each index with the mean and standard deviation computed over 1971-1999 inclusive, then cumulatively summed them, limited them to the data period, and then renormalized to produce SAM Σ , ENSO Σ , ASLP Σ , and ASLλ Σ .
The raw indices and their cumulative sums are shown in Figure S1 in Supporting Information S1.Correlations above 0.7 are evident between ASLP Σ and SAM Σ and between ASLλ Σ and ENSO Σ (Figures S1 and S2 in Supporting Information S1).This is due to the ASL being affected by larger-scale modes of climate variability, with SAM in particular modulating its absolute pressure and ENSO modulating the longitude of its center (Clem et al., 2017;Hosking et al., 2016;Turner et al., 2013).

Multi-Variate Empirical Orthogonal Functions
For a data-driven analysis we make use of Multi-variate Empirical Orthogonal Functions (MEOF) (Wang, 1992;Wu, 2023).MEOFs are an extension of conventional Empirical Orthogonal Functions but allow the dominant modes across multiple variables to be identified rather than treating each variable separately.We use MEOF to analyze the elevation and mass change gridded data sets after individual normalization.We first smooth the altimetry data set with a 50 km-wide Gaussian smoother and sub-sample the altimeter data set to match the 50 km horizontal grid resolution of GRACE.Given the limited sampling of altimetry in the northern Antarctic Peninsula we truncate that region from both data sets prior to computing MEOFs.

Regression
Using ordinary least squares, we solved the coefficients (a, b, c, d, and e) of the functional model describing timeevolving elevation (h) with time (t): Where f k = [1, 2] cycles per year.We adopted t 0 as the mid point of the altimeter series.

Data Uncertainty
For regression parameter uncertainties, we recognize the existence of temporal correlations in the altimeter time series (Ferguson et al., 2004), in part due to SMB variation (King & Watson, 2020), and take these into account.
Following King et al. (2023aKing et al. ( , 2023b)), we compared trend uncertainties from a linear regression using a Generalized Gauss Markov noise model to those generated using a white noise only (temporally uncorrelated) noise model using HECTOR v2.0 software (Bos et al., 2013).For regressions that included the SAM and ENSO terms, the white noise only model produced uncertainties a factor of 3 too small, taken as the median of the ratio of trend uncertainties, or factor 40 too small when not including the SAM and ENSO terms.We applied these scale factors to the uncertainties from the regression.For the GRACE uncertainties we used the scale factors of King et al. (2023aKing et al. ( , 2023b)).

Ice-Sheet Scale Analysis
Our data-driven MEOF analysis shows that ice elevation and mass time series are both dominated by decadalscale variability (Figures S3c and S3f in Supporting Information S1).Together, the two leading modes explain 65% of the non-linear variance of the combined and smoothed time series.Their corresponding principal components (PCs) correlate with detrended SAM Σ (r = 0.73) and 6-month lagged ENSO Σ (r = 0.89).The ASLP Σ and ASLλ Σ terms are not of direct relevance at the ice-sheet scale given the limited geographical footprint of influence of the ASL, but also have high correlations with the data.
GRACE and altimetry MEOFs have a high spatial correlation (Figures S3a and S3b, S3d-S3e in Supporting Information S1; r = 0.87 for MEOF1 and r = 0.75 for MEOF2) suggesting they are sensing the same signal and are both dominated by coastal changes.The potential in the high-resolution altimetry record is particularly evident in MEOF1 where the spatially diffuse signal in GRACE (Figure S3a in Supporting Information S1) is shown to be concentrated over small regions that coincide with the major ice streams of the Amundsen Sea Embayment and the coastline of the Bellingshausen Sea and Marie Byrd Land (Figure S3b in Supporting Information S1).We note that while MEOF3 (Figure S4 in Supporting Information S1) is partly affected by striping in the GRACE field, characteristic of GRACE systematic error, coherent signal is evident between GRACE and altimetry along the coastlines of the Bellingshausen Sea, Marie Byrd Land and Wilkes Land, suggesting the signal is robust in those regions, although the variance explained (5%) is much smaller than MEOFs 1 and 2. A similar signal to PC3, with periodicities of ∼4-7 years, has also been identified in analysis of GRACE (King et al., 2023b) or GRACE and altimetry data (Mémin et al., 2015).Beyond MEOF3, the modes explain little variance (<4%) and are dominated by noise, at least for GRACE (Figure S4d in Supporting Information S1).
To quantify the SAM and ENSO contribution to ice sheet elevation change we regress the altimetry time series against SAM Σ and ENSO Σ and the other parameters in Equation 1.Here we use the gridded data after applying a 10 km Gaussian spatial filter.The 5 km gridded altimeter regression analysis shown in Figures 1a and 1b reveals large-scale spatially coherent signal relating to SAM and ENSO around the coasts of Antarctica.Together, these two terms often explain more than 40% of the partial variance (R 2 partial ) of the timeseries around the coast and into the interior, with the partial variance the square of the partial correlation within which effects of the other regression terms are controlled.The median partial-variance explained across the ice sheet is 29% (Figure 1c).The SAM Σ coefficient is strongest in the Amundsen Sea Embayment where it centers on the Pine Island, Thwaites, Smith, and Pope Glaciers (Figure S5a in Supporting Information S1).The negative elevation signal in this region is linked to periods where positive SAM dominates negative SAM (positive SAM Σ ).Other strong signal exists along the coastal zone of the Bellingshausen Sea, Marie Byrd Land, and parts of coastal East Antarctica.A more diffuse signal is evident in the interior of West Antarctica and parts of East Antarctica (Figure S6a in Supporting Information S1).The ENSO Σ coefficient has particularly high positive values, indicating elevation increase associated with sustained El Niño, along the coast of the Bellingshausen Sea and well upstream into Pine Island Glacier (Figure S5b in Supporting Information S1).
Applying a 200 km Gaussian smoother to the altimeter data and rerunning the regression (Figures 1d-1e) produces coefficients with large-scale spatial coherence and larger partial variances explained, often exceeding 60% in key coastal regions but extending well into the interior of the ice sheet (Figure 1f).Comparing them to results of a regression with GRACE data (Figures 1g and 1h) (King et al., 2023a(King et al., , 2023b) ) shows high agreement in the signs and spatial distribution of the signal.We note that there are insufficient altimeter data in the Northern Antarctic Peninsula to analyze the signal in this region.Computing spatial correlations between the smoothed altimetry regression and the GRACE regression gives r = 0.65 for SAM Σ and r = 0.68 for ENSO Σ .
We next examine the role of SMB variability on the estimated coefficients from the altimetry regression.To do this we subtract the IMAU Firn Densification Model (IMAU FDM) v1.2A (Veldhuijsen et al., 2023) from the altimetry time series and repeat the regression.The results are shown in Figure 2. Comparing Figure 2a with Figure 1a shows that IMAU FDM effectively removes all the SAM-related signal in East Antarctic Ice Sheet

Geophysical Research Letters
10.1029/2024GL108844 (EAIS) but much of the SAM signal remains in West Antarctic Ice Sheet (WAIS).Much of the coastal EAIS ENSO-related signal is removed by IMAU FDM but with small over-correction evident for much of the ice sheet, including signal reversing sign in George V Land and WAIS.Repeating the regression but instead using the GSFC FDM v1.2.1 (Medley et al., 2022b) shows that there is significant sensitivity to the choice of FDM (Figures 2d-2f), with GSFC FDM apparently over-correcting ENSO-related signal in the Totten Glacier region but in much better agreement with the altimetry in WAIS.Given the decadal timescales of the signals, these intermodel differences are likely to have contributions from both the FDMs themselves and their underlying SMB models (Medley et al., 2022b;Verjans et al., 2021).
Next, we explore the origins of these signals further on a glacier-by-glacier basis.

Thwaites and Pine Island Glaciers
The partial variance explained by the SAM Σ and ENSO Σ terms (before subtracting an FDM) is above 60% for much of the Amundsen Sea Embayment (ASE; Figures 1c and 1f; Figure S5c in Supporting Information S1).Regardless of the FDM model adopted, much SAM Σ signal remains in the ASE broadly and ENSO Σ signal is evident in the Pine Island Glacier (PIG) region (Figure 2).Closer examination of these regions in Figure S5 in Supporting Information S1 (top row) indicates that the ASE signals are concentrated along low-elevation and fast flowing regions that correspond to PIG, Thwaites, and nearby glaciers.This is further evidenced through crosssections near to the front of these glaciers (Figure S6 in Supporting Information S1) along the yellow lines in Figure S5 in Supporting Information S1.It is notable that the phase of the SAM-related signal is switched in the low elevation and fast-flowing region of PIG.
Coefficient magnitudes generally decay upstream of the grounding line (Figure S7 in Supporting Information S1).Subtracting the IMAU FDM before performing the regression results in coefficients along the centerline and cross profiles that are shifted nearly uniformly but are not significantly altered in their spatial pattern (dashed lines Figures S6 and S7 in Supporting Information S1).
Along the coastal margin of the ASE the climatology is more directly controlled by the ASL than SAM and ENSO which modulate its depth and location (Clem et al., 2017;Turner et al., 2013).To explore this further we repeated the regression replacing SAM Σ and ENSO Σ in Equation 1 with ASLP Σ and ASLλ Σ .While the magnitude of the estimated coefficients differs between SAM Σ /-ASLP Σ and ENSO Σ /-ASLλ Σ the broader spatial pattern will be nearly identical due to the high correlations of these coefficient pairs over the data period (Figures S1 and S2 in Supporting Information S1) and so we just explore in detail the impact of estimating the ASL coefficients at one point location per glacier, at a centerline location about 20 km upstream of their respective grounding lines (Figure S5 in Supporting Information S1 yellow crosses; Table S1 in Supporting Information S1).
The detrended data are shown in Figure 3 (top row) where they reveal non-linear variability of several meters over the data period (blue plusses).Time series of estimated ASL coefficients sum to closely reproduce the data (black line).These two terms explain 84% (Thwaites) and 90% (Pine Island) of the partial variance of the altimeter time series.Interestingly, the phase of the ASLP Σ term is opposite between Thwaites and PIG, while the ASLλ Σ term is in phase.
Neither of the FDM models can explain the elevation variability at Thwaites or Pine Island glaciers (Figure S8 in Supporting Information S1, brown lines).This could be because the SMB models are unable to reproduce the precipitation in this region, especially in ∼2007 at Thwaites Glacier, but this would require a highly localized signal as this event does not occur at PIG.This is not implausible given existing SMB model limitations in lowaltitude coastal regions (Kappelsberger et al., 2023;Noël et al., 2023).The misfit could be caused by errors in background altimeter models, however we note we obtain nearly identical results using the alternative data set of Schröder et al. (2019) and misfit also exists in SMB-corrected GRACE fields (King et al., 2023b).A further possible source of the unexplained height signal is ice flow dynamics responding to large-scale climate variability.
The dynamic effect of ice flow and its influence on ice sheet mass and surface elevation at a given point can be estimated from satellite-derived glacier velocities and mass conservation (Text S1 in Supporting Information S1).The year-on-year changes in ice velocity since 2003 suggest several meters per year of dynamic elevation change.
Our analysis shows decadal variations in the lower parts of Pine Island and Thwaites, which could potentially be linked to climate variability and the way the glaciers respond dynamically to variations in ice shelf melt through a combination of advection and strain (Figure S9 in Supporting Information S1).

Totten and Denman Glaciers
The SAM and ENSO coefficients in the region of Totten and Denman glaciers have smaller magnitude and are much more diffuse than in the ASE (Figures S5d and S5e in Supporting Information S1).Nonetheless, these terms explain significant amounts of the partial variance (Figure S5f in Supporting Information S1) in this region.There is almost no non-linear signal to explain near the front of the Denman Glacier (Figure 3), with the largest SAM or ENSO signal in the Denman region located west of Denman.Nonetheless, SAM contributes about 30% of the partial variance at Denman.If the underlying surface lowering trend of Denman is affected by climate variability it is not obviously associated with SAM and ENSO over this period.
Despite the modest signal near Totten there is still evidence that significant SAM and ENSO signals exist in the fast-flowing region of Totten Glacier (Figure 3), at least in the 20-30 km above the grounding zone (Figures S6c and S7 in Supporting Information S1).Unlike the ASE glaciers, there is insufficient ice velocity time series for Totten Glacier to explore the cumulative impacts of time-varying ice dynamics on ice elevation.As noted above, the FDM-corrected results are model-dependent in this region and so the origin(s) of the Totten Glacier non-linear elevation change signal is unclear.

Discussion
Our analysis reveals the spatial fingerprints of SAM and ENSO on AIS elevation over 2002-2021, patterns which are confirmed by analysis of GRACE mass change data over the same period.These patterns may not be stationary with time.Indeed, circulation patterns associated with SAM are known to vary over decades (Marshall et al., 2013;Silvestri & Vera, 2009), with effects including variable precipitation in the Antarctic Peninsula (Goodwin et al., 2016).Within this context, it is therefore not unexpected that our pattern of SAM variability is different to the SMB-only SAM reconstruction of Medley and Thomas (2019) for the second half of the 20th century for instance.Differences with SMB-only reconstructions would also result if ice-dynamic effects on ice elevation and mass were non-negligible as hinted at by our data.
There are only a few previous studies exploring the relationship between ice dynamics, expressed as changes in ice mass, thickness, or elevation, and modes of climate variability, most notably in the Amundsen Sea Embayment region (Christie et al., 2023).Dutrieux et al. (2014) found reduced PIG ice shelf melt during a strong 2012 La Niña.Consistent with this, Paolo et al. (2018) found PIG ice shelf melting increased during El Niño, reducing ice shelf thickness, but that the ice shelf elevation increased overall due to increased accumulation.The way in which these decadal-scale changes impact upstream ice velocity and integrate with time to thickness and elevation variation requires investigation.
The rate of dynamic thinning or thickening depends on the relative contributions from strain and advection (Text S1 in Supporting Information S1).At PIG, both factors contributed to sustained thinning of >2 m per year in 2007-2008 in our analysis (Figure S9 in Supporting Information S1).In 2009, strain switched to thickening while thinning by advection also decreased.Hence, the glacier exhibited net thickening of >1 m per year during the 2012 La Niña.After 2015, observed velocities suggest the glacier thinned again, mostly due to strain.At Thwaites, strain thinning at rates of 4 m per year or greater outpaces the effect of advection, which is mostly positive (Figure S9 in Supporting Information S1).
The observed height anomalies (Figure 3) and dynamic elevation change (Figures S9a and S9b in Supporting Information S1) do not obviously correlate over the period velocities are available.This may have several explanations.To derive the latter, we assumed surface velocities are constant with depth, which is unlikely to be correct near the bed.Errors in the gridded velocity or ice thickness data also affect these estimates.Fundamentally, upstream ice elevation variability is unlikely to be a simple and immediate linear response to sub ice shelf climate variability (Snow et al., 2017).
The SAM/ASLP-related signal upstream of PIG, Thwaites, and other ASE glaciers is the largest unexplained signal in Antarctica, and it is especially strong at low elevations along the coast where atmospheric conditions may drive larger snowfall variations.While there is also a record of ice shelf geometry changes and grounding line retreat there, we suggest the large residual SAM/ASLP signal (Figure 2) in the ASE is related to SMB and/or firn densification because the ice-dynamic signal which is also present there has a different temporal pattern (compare Figure 3 with Figures S9c-S9f in Supporting Information S1).
We note that while the SAM Σ and ASLP Σ signals are correlated and our analysis cannot separate their different effects, they have different long-term implications for the ice sheet.As discussed by King et al. (2023aKing et al. ( , 2023b)), SAM Σ has a trend due to the positive phase of SAM that has emerged since the 1940s.ASLP Σ does not have a strong long-term trend, and so the extent to which the changes in coastal West Antarctica are related to the ASL rather than SAM will reduce the inferred contribution of SAM to ice-mass loss over recent decades (King et al., 2023b).

Geophysical Research Letters
10.1029/2024GL108844 Finally, our findings offer a simple way to remove decadal-scale variability from altimetry time series in the presence of imperfect firn densification models.This reduces correlated noise in the time series and will alter both the derived trends and, perhaps most significantly, the uncertainties of derived trends and other parameters if correlated noise is considered in the regression as it should (Ferguson et al., 2004;King & Watson, 2020;Williams et al., 2014;Wouters et al., 2013).

Conclusions
We analyzed gridded Antarctic ice elevation time series (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020) and show that much of the time series variance can be explained through a simple linear model based on the cumulative indices of the Southern Annular Mode and El Niño Southern Oscillation.The spatial pattern of this signal, once spatially smoothed, is in close agreement with the spatial pattern evident in GRACE data suggesting that observed ice elevation variability is robust and climatological.The Amundsen Sea Low is more directly relevant to the Amundsen Sea Embayment and we show that variations in its pressure and longitude linearly relate to ∼90% of the variance over Pine Island and Thwaites glaciers.
Subtracting the effects of modeled firn densification removes much, but not all, signal, with inter-model differences evident.Residual climatological signal is particularly large at the fronts of fast-flowing glaciers in the Amundsen Sea Embayment.Surface velocity and ice thickness data indicate that ice dynamics make a discernible contribution to decadal variability in upstream ice elevation.Further work is required to quantify the magnitude and response-times of upstream ice to changes in climatological variability in ice shelf melt.

Figure 1 .
Figure 1.Results of regression analysis of gridded data.Shown are the SAM Σ and ENSO Σ coefficients and variances explained for the altimetry (top row), altimetry after 200 km Gaussian smoothing (middle row), and GRACE (bottom row).The partial variances explained by SAM Σ and lagged ENSO Σ are in the right column.The hachuring indicates regions not significantly different to zero at the 95% confidence interval.

Figure 2 .
Figure 2. Results of regression analysis of FDM-corrected gridded altimeter data.Regression coefficients are shown (left and central columns) and the partial variances explained by SAM Σ and lagged ENSO Σ (right column).Shown are the coefficients and variances explained for the altimetry time series after subtracting of the IMAU FDM (top row) and GSFC FDM (bottom row).

Figure 3 .
Figure 3. Detrended elevation time series at glacier point locations.Time series are shown for sites ∼20 km upstream of the grounding line and along the centerline of flow (Figure S5 in Supporting Information S1 yellow crosses; Table S1 in Supporting Information S1).Shown are the altimeter time series after 10 km Gaussian smoothing and subtracting the estimated trend and harmonics (blue plusses), and the two components of the model (colored lines) and their sum (black line) for each glacier.For Thwaites and Pine Island glaciers (top row), ASL coefficients are shown, while for Totten and Denman glaciers (bottom row) SAM and ENSO terms are shown.The partial variances explained by the sum of the two coefficients are listed in each panel.Gray shading is the 1-sigma uncertainty of the model.Error bars represent the 2-sigma uncertainties of the data.