Response of Pacific-sector Antarctic ice shelves to the El Niño/Southern Oscillation.

Satellite observations over the past two decades have revealed increasing loss of grounded ice in West Antarctica, associated with floating ice shelves that have been thinning. Thinning reduces an ice-shelf's ability to restrain grounded-ice discharge, yet our understanding of the climate processes that drive mass changes is limited. Here, we use ice-shelf height data from four satellite altimeter missions (1994-2017) to show a direct link between ice-shelf-height variability in the Antarctic Pacific sector and changes in regional atmospheric circulation driven by the El Niño-Southern Oscillation. This link is strongest from Dotson to Ross ice shelves and weaker elsewhere. During intense El Niño years, height increase by accumulation exceeds the height decrease by basal melting, but net ice-shelf mass declines as basal ice loss exceeds lower-density snow gain. Our results demonstrate a substantial response of Amundsen Sea ice shelves to global and regional climate variability, with rates of change in height and mass on interannual timescales that can be comparable to the longer-term trend, and with mass changes from surface accumulation offsetting a significant fraction of the changes in basal melting. This implies that ice-shelf height and mass variability will increase as interannual atmospheric variability increases in a warming climate.

Projections of global sea-level change over the next century are highly uncertain due to insufficient understanding of the processes causing mass loss from the Antarctic and Greenland ice sheets[1]. The West Antarctic Ice Sheet (WAIS), which is mostly grounded below sea level[2] (a marine ice sheet), contains sufficient ice above floatation to increase global sea-level by over 3 m[3]. Ice-sheet models[4-6] suggest that mass loss will accelerate as glaciers and ice streams respond dynamically to internal instability mechanisms[5,7]. The predicted ice loss on timescales of decades to centuries from WAIS is about 1 m of global sea-level equivalent, with full ice-sheet loss within a few millennia[5,6].

The acceleration of grounded ice loss in the Amundsen Sea (AS) sector of WAIS has been attributed to reduced backstress as the fringing ice shelves thin and their grounding lines retreat[5,8]. The rapid, sustained thinning of the AS ice shelves[9,10] appears to be caused by increasing wind-driven flow of warm Circumpolar Deep Water (CDW) into the ocean cavities beneath ice shelves, enhancing basal melting[10,11]. Models[12,13] and limited observations[14,15] suggest that ice shelves might respond to changes in CDW circulation on interannual timescales. However, a paucity of time series of ocean observations on the continental shelf offshore of AS ice shelves and in the sub-ice cavities limits our ability to confirm this hypothesis, leading us to seek indirect measures of the sensitivity of ice-shelf mass change to large-scale climate variability.

Climate variability in the Antarctic Pacific sector

The El Niño-Southern Oscillation (ENSO), the Southern Annular Mode (SAM, or Antarctic Oscillation), and variability of the Amundsen Sea Low (ASL) are well-known climate drivers of interannual changes in the Antarctic Pacific sector[16-19]. ENSO is the leading mode of ocean-atmosphere variability on timescales of 2–7 years in the tropical Pacific, and is the strongest interannual climate fluctuation at the global scale[20]. ENSO causes much of the observed variability of the atmosphere, ocean, and sea ice in the Amundsen-Bellingshausen Sea sector[12,13], which exhibits the largest climate fluctuations around Antarctica[17,21]. Observed regional responses to ENSO include changes in snowfall[22-24], surface air temperature[25], sea ice extent[16,26,27], upwelling of CDW near the front of Pine Island Glacier’s ice shelf[14], and variation in basal melting under Getz Ice Shelf[15].

The SAM is a major driver of climate variability in the Southern Hemisphere, strongly influencing precipitation and temperature patterns from the subtropics to Antarctica[28,29]. SAM is usually strongest in austral spring and summer[30]. The phase of SAM influences the effect of ENSO in Antarctica, with the strongest Pacific sector response to ENSO when SAM is weak or in opposite phase[31] (i.e. with the combinations La Niña/SAM+ and El Niño/SAM− strengthening the atmospheric circulation anomalies in the mid-to-high latitudes).

The ASL is a persistent atmospheric low-pressure system located within the Amundsen and Bellingshausen seas, and plays the dominant role in determining the regional-scale pattern of atmospheric circulation across West Antarctica[17-19]. Through changes in strength (i.e. central pressure) and position, ASL determines the wind anomalies in all seasons, strongly influencing snowfall, temperature distribution and sea ice conditions near AS ice shelves[19]. Variations in the ASL position and strength are driven by tropical Pacific ocean-atmosphere variability (ENSO) and fluctuations in southern hemisphere pressure[17-19]. ASL central pressure tends to be lower during positive SAM conditions and La Niña years, and higher during El Niño years[17,31,32].

In this study, we show that ice-shelf height and mass changes in the Pacific sector, in particular the AS sector (Fig. 1), are correlated with interannual variability in regional atmospheric and oceanic circulation driven by ENSO.

Figure 1

Relationship between ice-shelf height anomalies and ENSO index

a, AS averaged ice-shelf height anomaly (12-month running mean; blue curve; top horizontal bars denote the time period of each satellite mission), 1-sigma bounds from 2000 bootstrap samples are plotted as dashed lines. Red line is the best fit between the height record and a combination of (12-month running integral) Oceanic Niño Index (ONI) and Amundsen Sea Low relative central pressure (ASL), both lagged by 4–6 months (preceding height); b, ONI; colored areas denote moderate-to-very-strong El Niños (red) and La Niñas (blue) as defined by NOAA (http://ggweather.com/enso/oni.htm).

Interannual changes in ice-shelf height

We merged data from four European Space Agency satellite radar-altimetry missions (ERS-1, ERS-2, Envisat, and CryoSat-2) to create individual records of ice-shelf height in the Pacific sector for the period 1994 to 2017. Each 23-year record samples every three months and represents an area of about 30 km × 30 km[9,33]. We detrended each height record using locally weighted regression (LOWESS)[34]. We then calculated interannual anomalies in ice-shelf height, δh(t), using 12-month running means.

We used 12-month running integrals of atmospheric and climate-index anomaly records (Methods) to represent climate forcing driving changes in δh. Correlating time-integrated forcing variables with δh(t) is analogous to correlating the climate drivers with dh/dt anomalies, while avoiding the high-frequency noise introduced by explicit differentiation of δh(t). We assessed the statistical significance of correlations from the 95% confidence intervals of 2000 bootstrap samples, taking into account autocorrelation in the time series[35,36] (Methods).

We evaluated an averaged δh(t) from 77 height records covering the AS sector (Fig. 1), and correlated this time series with time-integrated climate indices (Table 1). AS-averaged δh(t) is moderately-to-strongly correlated (r = 0.47 to 0.61) with the Oceanic Niño Index (ONI)[37], a measure of ENSO that tracks sea-surface temperatures in the east-central tropical Pacific. We also found a weak-to-moderate negative correlation between δh(t) and the ASL relative central pressure (r = −0.25 to −0.39). We found no significant correlation between δh(t) and SAM. The best fit between δh(t) and a linear combination of ONI and ASL indices (Fig. 1a) showed correlation values up to r = 0.67. Combining the time-integrated ONI and ASL central pressure into a single index lagged by 4–6 months with respect to δh(t) (Fig. 1a) slightly increased correlation values and their significance (Table 1); that is, some large changes in ice-shelf height that are unexplained by ONI (such as the height increase circa 2006) can be explained by the contribution from the ASL.

Table 1

Correlation between Amundsen Sea ice shelves and climate indices

Correlation between interannual anomalies in ice-shelf height (Fig. 1) with the Oceanic Nino Index (ONI), Southern Annular Mode (SAM), Amundsen Sea Low relative central pressure (ASL), and linear combination of ONI and ASL: δhfit = a ONI + b ASL (‘Best Fit’; red curve in Fig. 1). ONI, SAM and ASL are 12-month running integrals (Methods). Values in parentheses are correlations with ONI and ASL lagged by 4–6 months (preceding ice-shelf height).

	ONI	SAM	ASL	AS Ice Shelf
ONI	—
SAM	−0.04	—
ASL	−0.22*	+0.01	—
AS Ice Shelf	+0.47* (+0.61*)	−0.01	−0.25*(−0.39*)	—
Best Fit	+0.95* (+0.92*)	−0.04 (−0.09)	−0.51* (−0.59*)	+0.49* (+0.67*)

Correlations significant at the 95% level are marked with an asterisk, otherwise they are not significantly different from zero.

Correlation patterns between AS-averaged ice-shelf height anomaly and the time-integrated zonal (westerly) and meridional (southerly) components of the wind from ERA-Interim reanalysis[38] lagged by 4–6 months (Fig. 2), reveal a strong link between ice-shelf height variability and local environmental controls. Zonal wind anomalies modulate ice-shelf basal melting through upwelling and downwelling of coastal waters[13,14], and meridional wind anomalies modulate the transport of warm moist air (and associated precipitation) from the ocean to the AS ice shelves. The same fields (zonal and meridional wind, plus sea-ice concentration) showed approximately the same spatial pattern and correlation values when correlated with ONI (Extended Data Fig. 2).

Figure 2

Spatial pattern of correlation between wind and ice-shelf height anomalies

Correlation between 12-month running means of ice-shelf height anomalies in the AS (ice shelves in dark red) for the period 1994–2017 and 12-month running integrals of wind anomalies from ERA-Interim reanalysis for: a, the zonal (westerly) component,; b, the meridional (southerly) component. Wind anomalies (preceding ice-shelf height changes) were lagged by 4–6 months. The arrows indicate the direction of the wind anomaly that correlates with increase in ice-shelf height (e.g. in panel b, negative correlation between meridional wind and ice-shelf height means that onshore wind correlates with increase in height).

Extended Data Figure 2

Spatial pattern of correlation between wind, sea ice and ONI

Correlation between the Oceanic Nino Index and anomalies in: a, the zonal (westerly) component of the wind, associated with upwelling/downwelling of CDW underneath the ice shelves; b, the meridional (southerly) component of the wind, associated with transport of warm moist/cold dry air between the open ocean and West Antarctica; c, sea-ice concentration. Correlations are between 12-month running integrals of ONI and wind records (without lag; a and b) and 12-month running mean of sea ice (lagged by 4–6 months; c). See Extended Data Fig. 6 for similar (panel b) “reversed” pattern in SMB between the AS and BS-AP sectors.

Ice-shelf height and mass changes during 1997–2000

The most rapid changes in δh(t) generally coincided with large changes in ONI (Fig. 1b), with some modulation of the response to ENSO by changes in the ASL (Fig. 1a). The largest change in ice-shelf-height anomalies in the AS sector occurred between 1998 and 2001 (Fig. 1a), with more than 0.5 m of averaged surface lowering. This period encompassed one of the two strongest El Niños on record since 1950 (1997–1998), followed by a strong and long-lasting La Niña period between 1999 and 2001 (Fig. 1b).

ERA-Interim reanalysis revealed large differences in environmental conditions between the strong El Niño of 1997–1998 (EN1997) and the subsequent strong La Niña of 1999–2000 (LN1999) (Fig. 3). During EN1997 the westerly wind was stronger than normal (defined here as the 1979–2017 climatological mean), increasing upwelling of ocean waters on the continental shelf (Fig. 3a). The regionally-averaged air temperature was also warmer than normal, although with a cold anomaly present along the coast from Getz Ice Shelf to the western side of Pine Island Bay (Fig. 3b). Precipitation was higher than normal during EN1997, with the largest anomaly centered on Getz (Fig. 3c). There was an overall reduction in sea ice during EN1997 relative to LN1999, although its spatial structure was complex (Fig. 3d). Concurrently, AS ice-shelf height increased overall during EN1997 and decreased during LN1999, with respect to the longer-term trend (Fig. 3d).

Figure 3

Average oceanic, atmospheric and ice-shelf conditions during two distinct ENSO phases: El Niño (1997–1998) and La Niña (1999–2000)

Annual-average anomalies in local oceanic, atmospheric and ice-shelf conditions for 1997–1998 (left) and 1999–2000 (right). From top to bottom: a, wind velocity vectors, and vertical transport of coastal waters (upwelling/downwelling) derived from the wind-stress curl; b, air temperature; c, precipitation rate; d, observed sea-ice concentration, and ice-shelf height anomaly (six months after peak ENSO activity; see Methods). Fields in panels a–c are from ERA-Interim reanalysis. Surface mass balance and sea-level pressure during these conditions are shown in Extended Data Figs. 6, 7 and 8.

We assessed the relative contribution of anomalies in atmospheric-driven accumulation and ocean-driven melt by partitioning the observed (tide-corrected) ice-shelf height change (ΔhOBS) between two epochs into the sum of height changes due to atmospheric pressure (the inverse barometer effect; ΔhIBE)[39], buoyancy-compensated surface accumulation (ΔhSMB), basal melting (ΔhBMB) and ice dynamics (ΔhDYN), such that: We used Eq. (1) to estimate the relative magnitudes of accumulation- and melt-driven changes between EN1997 and LN1999 averaged for the AS ice-shelf area. We assume that ΔhDYN is zero: within this short time interval (annual averages centered at 1998.2 and 1999.5), ΔhDYN is about one-to-two orders of magnitude smaller than ΔhOBS (see Methods), and regional averaging smooths out non-coherent changes in ice dynamics across different ice shelves. We estimated ΔhIBE from the ERA-Interim pressure field and ΔhSMB from modeled surface mass balance (SMB)[40,41] and surface (firn) density. We obtained SMB from the regional climate model RACMO2.3[41] (which is forced by ERA-Interim), and from the ERA-Interim precipitation field directly. Estimates of SMB from both approaches are consistent (Extended Data Fig. 3). We estimated ΔhSMB for a range of firn densities (a Normal distribution with μ = 495 and σ = 25 kg/m3; Methods) that captures the likely densities of the top one-meter firn layer on AS ice shelves (Extended Data Fig. 4). We then estimated ΔhBMB as the residual term in Eq. (1) while including uncertainties for processes such as compaction rate and ice dynamics through Monte-Carlo methods (Methods). The spatial variability of these terms is shown in Fig. 4.

Extended Data Figure 3

Oceanic and atmospheric contribution to ice-shelf height change between El Niño and La Niña

(Same as Fig. 4 but now using ERA-Interim precipitation directly instead of RACMO2.3-based SMB) a, Measured ice-shelf height (anomaly) from altimetry; b, atmospheric pressure (IBE) from ERA-Interim (pressure decrease: height increase); c, SMB derived from ERA-Interim (reduction in accumulation: height decrease); d, melt rate inferred from the previous three: d = a – b – c (reduction in melting: height increase). SMB was converted to buoyancy-compensated height-equivalent using firn density of 490 kg/m3 for visualization purposes. The color map depicts the sign of each contribution to the measured height change (anomaly during La Niña minus that of El Niño).

Extended Data Figure 4

Modeled surface density profiles

Mean density (1-year average) of the firn layer (0 to 1 m depth). Two different epochs: El Niño of 1997–1998 (red) and La Niña of 1999–2000 (blue), on Pine Island Glacier ice-shelf front (a) and grounding line (b). These profiles are from a firn densification model [58], and are representative of the AS ice shelves (i.e. mean firn density varies smoothly from ice shelf to ice shelf in the AS). The density of fresh snow is constant with time in the model, although it does vary spatially as a function of temperature, accumulation rate and wind speed [58].

Figure 4

Oceanic and atmospheric contributions to ice-shelf height anomalies between El Niño (1997–1998) and La Niña (1999–2000)

a Altimetry-derived ice-shelf height changes (La Niña minus El Niño); b height change due to atmospheric pressure, derived from ERA-Interim (pressure decrease: height increase); c surface mass balance from RACMO2.3, for visualization purposes, surface mass balance was converted to buoyancy-compensated height-equivalent using surface density of 490 kg/m3 (accumulation decrease: height decrease); d basal mass change inferred from the previous three: d = a – b – c (melting decrease: height increase).

Anomalies in AS-averaged ice-shelf height between EN1997 and LN1999 were: ΔhOBS = −0.31 ± 0.09 m, ΔhIBE = 0.05 ± 0.01 m, ΔhSMB = −0.79 ± 0.21 m, and ΔhBMB = 0.43 ± 0.23 m (values are the mean and standard deviation from 2000 realizations of a Monte Carlo test; Methods). The surface-lowering rate anomaly (relative to the trend) from EN1997 to LN1999 (ΔhOBS/1.3 yr ≈ 0.24 m/yr) is comparable to the AS-averaged 23-year trend of ~0.21 m/yr (Extended Data Fig. 1).

Extended Data Figure 1

Average trend and variability in Amundsen Sea ice-shelf height

a, Average of all time-series of ice-shelf height change, relative to zero mean, in the Amundsen Sea (AS) sector. Green line shows trend component. b, Decomposition of average AS time series into its trend, interannual, and annual-plus-noise components. The ‘Annual + Noise’ component includes the residual of seasonal variation in backscatter after standard correction (for variation in ice-surface properties) was applied [33]. The interannual component has higher energy content than the annual component (see power spectrum, Extended Data Fig. 5).

We calculated the surface-mass anomaly from RACMO2.3, integrated over the AS ice-shelf area for the 1.3-year EN1997–LN1999 period, to be −42 ± 11 Gt. For the range of assumed surface densities, the basal mass-balance anomaly required to explain ΔhOBS in Eq. (1) is about 207 ± 111 Gt. This positive basal mass balance (i.e. less mass during EN1997 relative to LN1999) is consistent with previous reports of increased melt rate during El Niño and decreases during La Niña[12-14]. The surface-mass anomaly (accumulation) is about 20% of the basal mass anomaly (melting) in this case. Mass-anomaly change rates during EN1997–LN1999 of −32 ± 8 Gt/yr for SMB and 159 ± 85 Gt/yr for BMB are comparable in magnitude to the reported[42,43] longer-term changes in ice-shelf mass balance for the AS sector, with a mean SMB of ~+55 Gt/yr and a total mass-loss trend of about 156–173 Gt/yr (Extended Data Table 1).

Extended Data Table 1

Estimates of ENSO-driven changes in ice-shelf mass

Rate of change in ice-shelf mass from 1997–1998 (El Niño) to 1999–2000 (La Niña) in the AS sector. ΔSMB is change in accumulation; ΔBMB is change in basal melting; (for comparison) ‘SMB Trend’, ‘Thickness Trend’ and ‘BMB Total’ are the longer-term change (>5-year average) from accumulation, thinning (mass loss) and total meltwater production, respectively. These values are the average of previous estimates [42,43]. Densities used are 495 ± 25 kg/m3 (firn), 917 kg/m3 (solid ice), and 1028 kg/m3 (ocean water). ‘# Points’ are the co-located cells between height-change (~30 km) and ERA-Interim (~80 km) grid products.

Ice shelf	Area (km2)	# Points	Rate of mass change (Gt/year)
			ΔSMB		ΔBMB		SBMTrend	ThicknessTrend	BMBTotal
			ERA-lnterim	RACMO2.3	ERA-lnterim	RACMO2.3
Cosgrove	3000	2	−0.9 ± 0.2	−1.1 ± 0.3	+9.5 ± 2.4	+10.9 ± 2.8	+2	−3	−10
Pine Island	6000	2	−1.9 ± 0.5	−2.5 ± 0.6	+8.2 ± 4.8	+13.3 ± 6.1	+4	−32	−98
Thwaites	4600	4	−1.8 ± 0.5	−3.7 ± 0.9	+2.6 ± 5.6	+19.7 ± 9.5	+4	−30	−83
Crosson	2700	2	−1.3 ± 0.3	−3.3 ± 0.8	−2.3 ± 4.9	+16.1 ± 8.6	+4	−19	−39
Dotson	5400	5	−2.8 ± 0.7	−3.7 ± 1.0	+3.7 ± 8.9	+11.7 ± 10.9	+6	−17	−45
Getz	33200	49	−28.9 ± 7.1	−26.8 ± 6.7	+145.7 ± 76.4	+126.2 ± 70.5	+35	−65	−140
Amundsen	56000	70	−40.7 ± 10.2	−42.1 ± 10.5	+194.0 ± 106.9	+206.8 ± 110.5	+55	−165	−460

Coherence and extent of ENSO influence across the region

We assessed the spatial coherence of height changes in the AS by accounting for the covariation between the different locations, decomposing the set of 77 individual time series across the AS sector into principal modes of variability (i.e. empirical orthogonal functions; EOFs) using multivariate Singular Spectrum Analysis (SSA) (Methods)[44,45]. This spectral decomposition allowed us to detect small, coherent signals (with time-variable amplitude and phase) in short and noisy records while testing the significance of extracted modes[46]. This is a standard procedure for identifying the dominant frequencies of ENSO in time series of climate indices[44].

The spectrum of changes in ice-shelf height (Extended Data Fig. 5), identified two statistically significant (at the 95% level) spectral regions: a band with characteristic periodicity of 4–5 years (1st and 2nd EOFs); and the less energetic annual cycle (3rd and 4th EOFs). We obtained the same dominant interannual mode (~4–5 year timescales) when we applied univariate SSA to the full ONI record (1950 to 2017), also represented by the leading pair of EOFs[44] (not shown). This agreement between dominant timescales for ONI and ΔhOBS supports the presence of a causal link between ENSO and the variability of the AS sector ice shelves.

Extended Data Figure 5

Power spectrum of changes in AS ice-shelf height

Spectrum of AS ice-shelf height changes estimated from 77 time series, using multivariate Singular Spectrum Analysis (SSA). a, Spectral energy distributed along regularly sampled frequencies. b, Significance of each eigenvalue (mode of variability, or EOF) as identified by SSA. Error bars are 95% confidence intervals for the estimated background noise. Eigenvalues above the error bars are statistically significant. Pairs of EOFs with the same frequency represent an oscillation identified by SSA. The characteristic frequency of each EOF is obtained by least-square fitting the corresponding EOF to a sine function. The noise spectra were produced with a Monte-Carlo approach [46].

To quantify the extent of the ENSO influence across the region, we derived a “similarity index” (i.e. a distance metric) between the variance-normalized time-integrated ONI and δh(t) averaged over each ice shelf in the Pacific sector (Methods). Dotson and Getz ice shelves in the AS sector, and Nickerson, Sulzberger and Ross, have the greatest similarity with ENSO (Fig. 5). Getz Ice Shelf dominates the AS sector-averaged variability because it accounts for about 60% of the total AS ice-shelf area, is centered on the ENSO precipitation anomaly (Fig. 3c and Fig. 4c), and has a large inferred ENSO-driven basal melt-rate anomaly (Fig. 4d and Extended Data Table 1). These results suggest that Getz Ice Shelf is the largest contributor to ENSO-driven variability of freshwater flux to the AS.

Figure 5

Relative influence of ENSO along the Antarctic Pacific margin

a, Regional variation of the “similarity index” (size and color of squares) between ice-shelf height-anomaly records and the time-integrated Oceanic Niño Index (ONI). b, 12-month running integral of ONI (i.e. ENSO) lagged by 4–6 months (top plot) and 12-month running means of ice-shelf height anomalies for the AS (AMU) and six individual ice shelves, the shaded area highlights the large height change resulting from the 1997–2001 El Niño-to-La Niña transition. Ice shelves are: Pine Island (PIG), Dotson (DOT), Getz (GET), Nickerson (NIC), Sulzberger (SUL), and Ross (ROS).

Global and regional environmental controls on ice-shelf change

Besides ENSO and ASL, other regional modes of interannual variability (e.g. the Zonal Wave Three[47]) may also contribute to observed changes in the Pacific sector ice shelves, even though their correlations with δh(t) might be negligible; e.g. SAM (Table 1). These other modes interact with ENSO and with each other to drive complex patterns of ice-shelf height change that cannot be solely explained by a linear combination of climate modes. Further, studies[17,19,27,32] correlating ENSO tropical forcing with Pacific sector climate indicators, such as the ASL strength, sea-ice extent and Antarctic Peninsula temperature, have found that correlations with ENSO are significant for some seasons while not for others, with reversals of the sign of the correlation from season to season in some cases. Changes in the ASL strength, location and extent, which all exhibit large interannual variability (Extended Data Fig. 8), means that each distinct ENSO event yields different responses in individual ice shelves, depending on their location. Our ice-shelf height and mass assessments for the AS sector (Fig. 4; Extended Data Table 1) reflect the unique environmental conditions of the strong 1997–2000 ENSO period, during which the change in ASL between El Niño and La Niña cycles was exceptionally large. In contrast, the strong 2015–2016 El Niño, characterized by remarkably large changes in tropical sea-surface temperature, did not impact AS ice shelves in the same way (Fig 1a and Fig. 5). This was likely a result of El Niño interaction with a positive SAM during the development phase, a northward displacement of the ASL relative to previous strong El Niños[48] (Extended Data Fig. 8), and absence of a following strong La Niña.

Extended Data Figure 8

Amundsen Sea Low during El Niño and La Niña years

Changes in intensity of the Amundsen Sea Low (ASL) pressure system during El Niño and La Niña years within 1994–2017 (derived from ERA-Interim). The ASL changes consistently between El Niño years (positive anomaly / weaker ASL) and La Niña years (negative anomaly / deeper ASL). Each panel matches (approximately) a moderate-to-strong El Niño and La Niña year as defined by NOAA [and depicted in the Oceanic Niño Index, Fig. 1: 1997.7 (EN), 1999.0 (LN), 2002.5 (EN), 2008.5 (LN), 2009.8 (EN), 2011.2 (LN), and 2015.9 (EN)].

The dominant effect of El Niño on AS ice-shelf mass is likely the increased basal melting associated with onshore flow of CDW and coastal upwelling as westerly wind stress intensifies (Fig. 3a)[12-14]. On interannual timescales, this basal mass loss anomaly, relative to the longer-term mass loss trend[9], is partially offset by increased snowfall. This precipitation increase is consistent with the northerly wind anomaly during El Niño events (Fig. 3a), possibly including increased local moisture uptake from the coastal ocean due to a reduction in regional sea-ice concentration (Fig. 3d).

The reversed pattern of meridional-wind anomalies between the AS and Bellingshausen Sea–Antarctic Peninsula (BS/AP) sector (Extended Data Fig. 2) is consistent with modeled surface mass balance for WAIS (Extended Data Fig. 6) and with satellite-gravity measurements over the AS and AP sectors[22]. Precipitation increases during El Niño in the AS sector while it decreases in the BS/AP region, and vice versa during La Niña. The principal mechanism affecting local winds is the response of the ASL pressure system to ENSO and SAM[17,19,31], with ASL consistently changing its intensity within El Niño and La Niña years (Extended Data Fig. 8). The rapidly thinning Pine Island Glacier ice shelf showed a weaker response to ENSO (Fig. 5), consistent with the weaker, and in some cases reversed, atmospheric and sea-ice anomalies in Pine Island Bay relative to the broader AS (Fig. 2 and 3).

Extended Data Figure 6

Antarctic surface mass balance during El Niño and La Niña

a, Anomaly in surface mass balance (SMB) during El Nino of 1997–1998 and, b, during La Nina of 1999–2000 (reversed pattern), from a regional climate model (RACMO2.3) [41]. The AS sector shows the strongest ENSO-driven anomaly in SMB of the entire Antarctic continent, consistent with the ENSO-driven anomaly in precipitation (from ERA-Interim; Extended Data Fig. 7). c, The largest precipitation (and surface density) values in Antarctica occur along the West Antarctic coast as evidenced by the long-term mean SMB (and Fig. 7 in Ligtenberg et al.) [58].

We have shown how modes of tropical Pacific variability affect the height and mass of ice shelves in the AS sector of WAIS. The response in height is the combined effect of two opposing processes, which are both intensified during El Niño events: surface snow accumulation and ocean-driven basal melting. The result is an overall height increase, but net mass loss, since the ice lost from the base has higher density than the fresh snow being gained at the surface. Ice-shelf response to ENSO variability is strongest between Dotson and Ross ice shelves, with a weak response in Pine Island Bay, the Bellingshausen Sea, and west of the Ross Sea. Given expected increases in total precipitation[49] and frequency of extreme ENSO events[50] as Earth’s atmosphere warms, our results imply that interannual variability of ice-shelf height and mass will also increase, stressing the need to quantify surface accumulation relative to basal melting to project future changes in Antarctic ice shelves.

Methods

Ice-shelf height

Ice-shelf height change time series are derived from measurements acquired by four overlapping satellite radar altimetry (RA) missions (ERS-1, 1991–1996; ERS-2, 1995–2003; Envisat, 2002–2012; CryoSat-2, 2010–2017). Our records span 23 years (1994–2017), have a temporal resolution of ~3 months, and represent grid cells of ~30×30 km.

The processing steps for the conventional altimeters (ERS-1, ERS-2 and Envisat) are fully described by Paolo et al.[9,33]. We performed the following processing steps for CryoSat-2, applied to ESA’s SARIn L2 Baseline C product over the Antarctic ice shelves: we corrected for a 60-meter range offset for data with surface types ‘land’ or ‘closed sea’; removed points with anomalous backscatter (>30 dB)[51]; estimated heights with a modified (from McMillan et al.[52]) surface-fit approach, with a variable rather than constant search radius to account for CryoSat-2’s heterogeneous spatial sampling; removed height estimates < 2 m above the Eigen-6C4 geoid[53] to account for ice-shelf mask imperfections near the calving front; applied all the standard corrections to altimeter data over ice shelves (e.g. removed gross outliers, and residual heights (with respect to mean topography) >15 m; ran an iterative three-sigma filter; and minimized the effect of variations in backscatter[9,33]).

When merging the records, we used only grid cells that were sampled by all four satellites, only accepted time series that overlapped by at least 0.75 years to ensure proper cross-calibration, and removed all merged time series with standard deviation of residuals (with respect to the trend) >1 m. This removed records with, for example: satellite mispointing, anomalous backscatter fluctuations, grounded-ice contamination, high surface slopes and geolocation errors.

Sea-ice concentration

Daily time series of sea-ice concentration, averaged to a 25-km polar stereographic grid, are from the National Snow and Ice Data Center’s Daily Polar Gridded Sea Ice Concentrations Version 1 (Near-Real-Time Defense Meteorology Satellite Program Special Sensor Microwave Imager/Sounder)[54], processed with the NASA Team Sea Ice Algorithm for sea-ice concentration retrieval[55]. We averaged the daily records to monthly anomaly time series, i.e., the anomaly with respect to the climatological average (1979–2017) for each month.

Climate indices

The Oceanic Niño Index (ONI), which tracks the running three-month mean sea-surface temperatures in the east-central tropical Pacific with limits 5°N–5°S and 120°W–170°W (the Niño 3.4 region), is our primary metric of ENSO. This standardized series is derived from the Extended Reconstructed Sea Surface Temperature Version 4 (ERSST.v4)[37,56]. We also reproduced the analysis using an alternative metric of ENSO, the Southern Oscillation Index (SOI), a standardized time series of observed sea-level pressure differences between Tahiti and Darwin, Australia. The Southern Annular Mode (SAM) is the station-based index from Marshall[28], which tracks the zonally-averaged meridional pressure difference between 40°S and 65°S. We obtained the Amundsen Sea Low (ASL) relative central pressure from Hosking et al.[57]. Our ‘Best Fit’ index is the best fit (in the least-squares sense) of ONI plus ASL 12-month running integrals (lagged by 4–6 months) to the AS-averaged ice-shelf height record. The running integral (or running sum) on both indices ensures consistency between the climate records.

Atmospheric variables

Surface air pressure, wind fields, precipitation and temperature are from the ERA-Interim global atmospheric reanalysis[38] for 1979–2017 by the European Centre for Medium-Range Weather Forecasts (ECMWF). These products are provided on a grid with spacing of ~80 km. We converted the monthly averaged records to anomalies with respect to the climatological average (1979–2017) for each month. The ERA-Interim model estimates of ENSO precipitation variability agree well with ice-core records and in-situ radar from the catchments of Thwaites and Pine Island glaciers[24]. It also agrees well on a regional scale with satellite gravimetry observations from the GRACE mission[22].

Ice-shelf surface mass balance and firn density modeling

Time series of surface mass balance (SMB) over the ice shelves are from a simulation based on output of a regional atmospheric climate model (RACMO2.3[41]), which is forced by the ERA-Interim reanalysis, has a horizontal grid spacing of 27 km, and accounts for processes such as surface meltwater retention due to refreezing, evaporation, wind drift, and sublimation.

The modeled SMB, surface temperatures, and 10-m wind speeds have been used to force a transient run of a one-dimensional time dependent Firn Densification Model (IMAU-FDM) that models the vertical profile of snow and firn density at each location[58]. The product was evaluated by comparing the temporal mean of firn density against 750 in-situ measurements continent-wide (r = 0.88[40]), and comparisons with radar-derived accumulation[24].

Correlation analysis

We identified correlations between ice-shelf height changes and climate forcing by removing the long-term trend from each height record using locally weighted regression (LOWESS)[34], with a smoothing span equal to 1/3 of the record (~7.7 years, a compromise between retaining decadal fluctuations in the trend and excluding the shorter ENSO time scales), and then computing 12-month running means (hereafter referred to as ‘height anomalies’, δh(t)). For the wind components and climate indices we computed 12-month running integrals (the sum of the preceding 12 months). We correlated δh(t) with the time-integrals of forcing anomalies represented by ERA-interim wind fields and climate indices. This approach is analogous to correlating dh/dt with time-dependent anomalies and indices, but avoids the high-frequency noise amplification due to explicit differentiation of δh(t). We tested correlations with and without lags between ice-shelf height anomalies and the climate records. We constructed confidence intervals for each correlation coefficient, taking into account temporal autocorrelation effects, using the nonparametric stationary bootstrap[35] with an average block length proportional to the maximum estimated persistence time (τ) of each pair of time series being analyzed[36,59]: where ρ is the autocorrelation coefficient and dt is the time step.

Spectral analysis

We identified common oscillatory behavior in the AS sector by decomposing the set of 77 height-anomaly time series from the AS (Fig. 1) into orthogonal components, using multivariate Singular Spectrum Analysis (SSA)[44,45,60,61]. In SSA, pairs of eigenvectors (EOFs) with similar periodicity represent a phase-and-amplitude modulated temporal oscillation (e.g. EOFs 1 and 2 in Extended Data Fig. 5). Our spectral decomposition by SSA also included statistical tests to discriminate between potential oscillations and “white” (independent and identically distributed) noise through Monte Carlo simulations[46]. Free software for implementation of SSA, among other spectral techniques, is provided by the SSA-MTM Toolkit at http://www.atmos.ucla.edu/tcd/ssa.

We removed the long-term trend from each ice-shelf time series (without smoothing), and normalized (to unit variance) each record prior to applying SSA[44,62]. Reported results refer to an SSA window length of 6 years (window lengths up to 9 years provided similar results). The SSA identified two significant modes of variability in ice-shelf height: 4–5-year time scales, and the annual cycle (Extended Data Fig 5).

Similarity index

We examined the spatial pattern of ENSO influence on the Antarctic Pacific sector ice shelves, by evaluating the similarity of each ice-shelf-averaged height anomaly record to the ONI. We define the similarity index (analogous to correlation) as: where x is the ice-shelf height anomaly record and y is the ONI, both normalized to unit variance. For graphical clarity, we normalized the SI values to a range of 0 to 1.

Equivalent height and mass calculations

For mass calculations, we estimated changes in ice-shelf height due to: variations in sea-level atmospheric pressure (ΔhIBE) from ERA-Interim pressure fields, assuming that surface-height change is −0.01 m for +1 hPa change in pressure[39]; variations in surface mass-balance (ΔhSMB) from the RACMO2.3 model output[41] (Fig. 4) and assumed surface densities (see next paragraph); and variations in snowfall from ERA-Interim precipitation directly (Extended Data Fig. 3), which is natively provided in the product as meters of liquid-water equivalent per day. For each model, we converted surface mass-balance and precipitation values (anomalies with respect to the 1979–2015 and 1979–2017 means, respectively) to buoyancy-compensated snow/firn height, ΔhSMB, using surface snow/firn density values randomly sampled from a Normal distribution with mean and sigma equal to 495 and 25 kg/m3, respectively (see ‘Uncertainty assessment’ below), and ocean water density of 1028 kg/m3.

We chose the range of surface density values based on the limited information available from these ice shelves, for which direct surface measurements are scarce. Results from a semi-empirical FDM suggest that surface densities along the WAIS coast are higher than in the interior; some of the highest values in Antarctica are found along the AS coast (Fig. 7 in Ligtenberg et al.)[58]. These high densities are due to a combination of relatively high temperatures, strong winds and elevated precipitation rates. Extended Data Fig. 4 shows density profiles from 0 to 1 m depth over Pine Island Glacier as represented by the FDM. Modeled density ranges from ~480 kg/m3 at zero depth to ~510 kg/m3 around 1-m depth. Our range of likely surface (firn) densities (495 ± 25 kg/m3) for the AS ice shelves, extends ± 10 kg/m3 beyond the model-derived range.

Further justification for our choice of densities comes from the fact that radar altimeters are relatively insensitive to fluctuations in low-density surface mass balance. Part of the radar signal penetrates through fresh snow, and changes in surface density (due, for example, to seasonal variations in radiation and air temperature) can affect the radar backscatter and derived height measurement[33]. We found no correlation between changes in height and changes in backscatter at interannual time scales (12-month running means) for the AS ice shelves, suggesting that interannual estimates of surface height are not affected by changes in surface properties. This is not the case at seasonal time scales where correlations were significant[33]. The layer being tracked by the RA as the “surface” must be sufficiently compacted to limit radar penetration and backscatter-height correlated changes. This result implies that height changes driven by surface accumulation may be higher than our RA measurements show, which would lead to a more important role of SMB relative to BMB in the height variability.

Phase lag between height changes and climate drivers

Height changes lag the climate drivers by about 4–6 months, based on correlation estimates for a range of lags (with maximum correlations for 4 to 6 months). We suggest three potential causes of this phase lag. First, the oceanic ENSO signal in BMB will be affected by the lagged and filtered dynamic response of the coastal ocean to the atmospheric forcing, caused by advection time scales over the continental shelf, ocean mixing, and response to anomalies in surface-buoyancy forcing as sea-ice volume changes. A typical advection time scale is given by T = L/U, where L is a characteristic length scale for the continental shelf (~500 km) and U is a typical mean advective velocity of CDW (~5 cm/s)[63], giving T ~ 4 months. Second, the RA measurement of ice-shelf height is some integration of backscatter within the surface snow and firn layer, with a time lag being induced by radar signal sensitivity to firn compaction (the development of a subsurface strong reflector) that takes place over some time interval after the precipitation signal. Third, the peak signals of interannual climate modes including ENSO and SAM tend to be short-lived (roughly one year for large El Niño events (Fig. 1)), so that the effect of each event on an ice shelf depends on its phasing relative to the strong annual cycles of atmospheric, oceanic and sea-ice conditions. We accounted for the measured phase delays by calculating the mass budget with ΔhOBS lagged by 6 months relative to the ΔhSMB in Eq. (1).

Estimates of mass change between EN1997 and LN1999 for individual ice shelves are presented in Extended Data Table 1. The estimates for basal mass change are sensitive to the choice of surface density. Our estimated change in melt rate between El Niño and La Niña for Getz Ice Shelf is between 3.1 and 3.5 m/yr (Extended Data Table 1). These values are consistent with the magnitude of observed variability in Getz’s melt rate, with interannual changes on the order of 3 m/yr reported by Jacobs et al.[15].

Uncertainty assessment

We assessed the sensitivity of our results to measurement uncertainty using a Monte Carlo (MC) approach. Through an iterative process, we scaled each spatial field used in Eq. (1) by a percentage drawn from a Normal distribution with sigma equal to the field’s prescribed error. We also draw surface density values from a Normal distribution with mean and sigma equal to 495 and 25, respectively, as justified above. Our reported values with uncertainties are then the mean and standard deviation from 2000 MC realizations. This approach has the advantage of directly quantifying the impact on the results of prescribed errors in the input fields, propagating these errors according to the relationship between the variables.

On interannual timescales, the main components of the uncertainty in ΔhOBS are: the statistical error (e.g. precision of the instrument, poor geophysical corrections, improper range estimation from the return radar waveform); surface compaction rate (changes in firn-air content); and neglecting the dynamic term (ΔhDYN) in Eq. (1).

Since the spatial field used for ΔhOBS in Eq. (1) is the result of averaging tens-to-hundreds of crossover values (height differences) per ~30×30 km grid cell[9,33], and we further average height-change estimates over 12-month time intervals, the statistical error on each grid cell is small (<10%)[33].

Ku-band (~13.5 GHz) radar altimeters mostly penetrate the top layer of fresh snow. Radar measurements are, therefore, less sensitive to compaction rates than laser-derived measurements, which track the air-firn interface. Using laser-altimeter observations, Moholdt et al.[64] estimated compaction rates across Ross and Filchner-Ronne ice shelves on the order of −9 ± 10 mm/yr and 23 ± 10 mm/yr, respectively. If we assume that all of the compaction rate (an upper bound for uncertainty due to compaction) translates to the radar-measured height change, the largest of these values is less than 8% of the EN1997–LN1999 height change.

Accurate estimation of ΔhDYN in Eq. (1), as applied to the transition from EN1997 to LN1999, would require annual velocity changes for the 1997–1999 period over the entire AS ice-shelf area; such data do not exist. We can, however, estimate the order of magnitude of the height change that would result if all the measured ice-flux increase in the AS[65] went into dynamic thickening (i.e., an extreme upper bound), using: where ΔQ/Δt is change in flux, A is ice-shelf area, ρi and ρw are the densities of solid ice (917 kg/m3) and sea water (1028 kg/m3), respectively, and δt is the time interval in question (1.3 years). The only estimate of ice-flux anomalies (Q) around the EN1997–LN1999 period comes from Mouginot et al.[65], which provides Q for the glaciers that feed Pine Island (PIG), Thwaites, Crosson, and Dotson ice shelves at 1996.1 and 2001.0. Using values for PIG (+8 Gt over 4.9 years, and ice-shelf area of ~6000 km2), we obtain ΔhDYN ≈ 0.042 m, which is ~13.5% of the EN1997–LN1999 measured height change if we assume that all ice shelves in the AS sector have accelerated coherently at the same rate as PIG (which is clearly not the case[65]). Using PIG, Thwaites, Crosson and Dotson combined (+10 Gt over 4.9 years, and area of ~18,700 km2), we obtain ΔhDYN ≈ 0.017 m, which is ~5% of the EN1997–LN1999 measured height change.

While variability in precipitation is well represented in reanalysis products and surface-mass-balance models, its magnitude has substantial uncertainty[24,41], especially at the low elevations of ice shelves where surface mass fluxes are much larger than at higher altitudes of the ice-sheet interior. For the AS Embayment, previous studies[66,67] have attributed a 15% uncertainty (based on field measurements[68]) to the SMB product. This corresponds to about 21% ( ) of the EN1997–LN1999 SMB-induced height change.

The contribution of IBE to total height change in Eq. (1) is small (~0.05 m). Bracegirdle and Marshall[69] estimated differences between Era-Interim and observations for annual mean sea-level pressure around Antarctica over the period 1979–2008. They found a standard deviation of 0.31 hPa. This corresponds roughly to 9% ( ) of the EN1997–LN1999 IBE-induced height change.

To ensure that we captured all these uncertainties in our error assessment, we assumed large (upper bound) errors, and derived the (one-sigma) ranges used in our MC calculation as: ΔhOBS ± 30%, ΔhSMB ± 25%, ΔhIBE ± 15%, and ρfirm ≈ 495 ± 25 kg/m3.

We note that if accumulation values were higher, or if surface density was higher than the values we used based on a range of modeled densities, then precipitation-driven change in mass would be more important relative to melt-driven change than we report here. If we assumed a wider range of surface densities (for the “compacted” firn layer tracked by the RA), our estimated range for ENSO-driven change in basal melt would be wider. However, for all reasonable ranges that we have tested (range of firn densities up to 420–550 kg/m3), the fundamental conclusions of this study remain the same.

Complexity of oceanic and atmospheric forcing interaction

The changes in ice-shelf height and mass at interannual time scales arise from the combination of direct atmospheric terms, predominantly precipitation, and oceanic forcing affecting basal melt rates. These terms can be closely coupled through the ocean circulation’s response to wind stress and buoyancy forcing, both of which are also correlated with precipitation. However, ocean response may be delayed and filtered through advection (e.g. the time taken for CDW at the continental shelf edge to reach an ice shelf) and mixing. While our results specifically for ENSO indicate a clear and dynamically consistent relationship between SMB and BMB variability, this relationship will not necessarily hold for oscillations associated with other climate phenomena such as the Southern Annular Mode (SAM) and the Zonal Wave Three mode (ZW3). Our results from a relatively short (compared to typical ice-sheet and climate timescales of variability) 23-year observational window, only apply to time periods in which ENSO dominates.

Data availability

The ice-shelf height-related data generated and analyzed during the current study are available in the PANGAEA repository [https://doi.pangaea.de/10.1594/PANGAEA.882376]. The climate data that support the findings of this study are available in the ESR repository [ftp://ftp.esr.org/pub/datasets/ERA_Int]

Extended Data Figure 7

Change in pressure and precipitation Antarctica-wide

Circum-Antarctic change in sea level pressure (a) and precipitation (b) between the El Niño of 1997– 1998 and La Niña of 1999–2000. These large-scale patterns show that ENSO-driven anomalies are localized: pressure (the Amundsen Sea Low) changing over the Amundsen-Bellingshausen Sea, and precipitation changing substantially over the Amundsen Sea ice shelves. Fields were derived from the ERA-Interim reanalysis [38].