Investigating the past, present and future responses of Shallap and Zongo Glaciers, Tropical Andes, to the El Niño Southern Oscillation

Abstract

Tropical Andean glaciers are highly sensitive to climate change and are impacted by the El Niño Southern Oscillation (ENSO). However, glaciological data are scarce, meaning that there are substantial knowledge gaps in the response of Andean glaciers to future anthropogenic and ENSO forcing and these are crucial to address, as glaciers represent a key water source for downstream populations and ecosystems. Here we integrated data from glaciological field studies, remote sensing, statistical analysis and glacier modelling to analyse the response of two Andean glaciers (Zongo and Shallap) to ENSO and their potential sensitivity to a range of climate forcing scenarios. Both glaciers retreated and experienced increasingly negative mass balance between the 1990s and the 2010s and responded strongly and rapidly to contemporary ENSO forcing, although this relationship evolved over time. Sensitivity experiments demonstrate that Shallap and Zongo are highly sensitive to ENSO forcing scenarios and the combination of ENSO and climate warming can cause rapid ice loss under the most extreme scenarios. Results also demonstrate the strong sensitivity of both glaciers to changes in the equilibrium line altitude, whereby rapid ice loss occurred when melt extended into present-day accumulation areas.

1. Introduction

Mountain glaciers represent important fresh water stores, releasing meltwater downstream for drinking water, agriculture, hydropower and ecosystems, yet most glaciers globally are shrinking in response to climate warming (Clarke and others, 2015; Zemp and others, 2019; Immerzeel and others, 2012). Tropical mountain glaciers are particularly sensitive to climate warming as their high altitudes and low latitudes expose them to strong radiative forcing (Chevallier and others, 2011), and because ablation occurs year-round (Veettil and Kamp, 2017). Glaciers of the Tropical Andes have retreated throughout the 20th and 21st centuries (Sicart and others, 2011), as average temperatures have increased by 0.26°C per decade (Rabatel and others, 2013). This is significant because tropical mountain glaciers are essential for sustaining base flows during the dry season: e.g. in Bolivia, one-third of the population live in mountain communities, where seasonal rainfall is supplemented by glacier melt (Rangecroft and others, 2015; Perry and others, 2017; Zekollari and others, 2019). In the near future, enhanced glacier melting may initially increase runoff, but as the glaciers continue to shrink, this water supply will reduce and will threaten water security in the Andes (Owen and others, 2009). At the same time, downstream populations and infrastructure are likely to expand, leading to increased demand and water stress (Rangecroft and others, 2015). It is therefore critical to understand how Andean Mountain glaciers will respond to climate warming and to predict their longevity as a crucial resource for mountain communities and ecosystems. Despite this, data on Andean glaciers are limited, partly due to their remote, high-altitude locations (Rabatel and others, 2013), and few studies have attempted to model their near-future behaviour (Réveillet and others, 2015; Yarleque and others, 2018).

Andean glacier retreat has accelerated since the 1970s (Schauwecker and others, 2014; Braun and others, 2019) but these changes have been spatially heterogeneous: low altitude glaciers (<5400 m) retreated faster than those at higher altitudes and are forecast to disappear during the 21st century (Fraser, 2012). Another source of variability in the response of Andean glaciers to climate change is their location in the inner, versus the outer, tropics. Inner-tropical glacier mass balances are controlled by changes in temperature, which influences the position of the equilibrium line altitude (ELA) and snow-line altitude (SLA) (Favier and others, 2004; Schauwecker and others, 2014). The outer tropics, by contrast, have relatively constant temperatures, but experience pronounced seasonality in precipitation, with minima occurring in June–August and maxima in January–March (Hofer and others, 2010). Consequently, outer-tropical ablation occurs throughout the year, but accumulation is limited to the wet season, which typically occurs between November and March and reaches its maximum in January–March (Rabatel and others, 2013). Mass balance in the outer tropics is therefore determined by the distribution and phase of precipitation, cloudiness and timing of the wet season onset (Favier and others, 2004).

Outer tropical glaciers can be further sub-divided into the wet and dry outer tropics, which have differing surface energy balance (SEB) and mass balance. In the wet outer tropics, the key control on SEB and mass balance is the onset of wet season precipitation and its impact on glacier albedo. Precipitation in the form of snow strongly reflects incoming solar shortwave radiation and therefore reduces melt rates; even a small delay in the arrival of seasonal precipitation can enhance glacier melt markedly (Wagnon and others, 1999; Sicart and others, 2011, 2008). By contrast, in the dry outer tropics, increased cloud cover and temperature are the key variables controlling the SEB and mass balance. Cloud cover impacts the amount of longwave radiation reaching the glacier surface, while temperature impacts the location of the ELA and SLA (Gurgiser and others, 2013; Maussion and others, 2015).

A further important influence on glacier mass balance in the Tropical Andes is the El Niño-Southern Oscillation (ENSO) (Sulca and others, 2018; Manzanas and Gutiérrez, 2019). ENSO has a periodicity of 2–7 years and influences climate in the Tropical Pacific and globally (Dietze and others, 2008). In its negative phase, El Niño, the tropical Andes experiences above-average temperatures and below-average precipitation, with the inverse occurring during La Niña (Dietze and others, 2008). The impact of ENSO on glacier mass balance varies between the wet and dry outer tropics (Favier and others, 2004; Vuille and others, 2008b). For example, at Zongo Glacier (Bolivia) in the wet outer tropics, strong El Niño conditions in 1991–1992 caused the peak wet season to be half its average duration (Francou and others, 1995). The delayed wet season onset and reduced duration of maximum precipitation intensity resulted in bare ice being exposed for longer, resulting in reduced albedo and enhanced ablation (Francou and others, 1995). In contrast, at Shallap Glacier (Peru) in the dry outer tropics, ENSOs affect temperature, which in turn controls the SLA and the phase of precipitation (Maussion and others, 2015; Schauwecker and others, 2017).

Despite these overall patterns, individual glaciers may respond differently to ENSO due to variable glacier area, elevation, thickness, hypsometry and topographic setting (Kaser, 1999; Grosshauser, 2012; Roe and Baker, 2014; Brahmbhatt and others, 2017; Mir and others, 2017; Christian and others, 2018). Furthermore, this relationship varies across the altitudinal range of individual glaciers, with ENSO having a greater effect at lower altitudes and in the ablation zone (Vuille and others, 2008b; Maussion and others, 2015; Silverio and Jaquet, 2017). This may have important implications for the future response of tropical Andean glaciers to ENSO: as glaciers retreat into higher altitudes, ENSO's influence may reduce. However, this may be offset by recent increases in ENSO intensity, enabling the ENSO signal to be transmitted to higher altitudes (Silverio and Jaquet, 2017).

Alongside spatial variability, temporal variability in the ENSO–glacier mass-balance relationship is also observed, with the relationship appearing to breakdown during certain periods (Levado-Casimiro and others, 2013; Schauwecker and others, 2014). For example, the strong El Niño conditions of 1982/1983 produced very little change in Andean glacier mass balance, and in 1993/1994, neutral ENSO conditions coincided with strong regional mass gain (Vuille and others, 2008a). Variable intensities and lengths of ENSO periods, the distribution of the ENSO peak in a given hydrological year and the gradual overriding of ENSO's climatic signal by climate warming are all potential causes for temporal variability in the ENSO–mass-balance relationship (Vuille and others, 2008b; Veettil and others, 2017a).

Increases in the intensity and frequency of ENSO since the 1970s have been linked to climate warming (Silverio and Jaquet, 2017; Veettil and others, 2017b). Pacific sea surface temperatures (SSTs) are forecast to rise by ~1.5°C during the 21st century, which is likely to increase ocean-atmospheric volatility and lead to ENSO conditions developing more regularly (Steinhoff and others, 2015). This may favour the dominance of El Niño events versus La Niña (Koutavas and Joanides, 2012; Cai and others, 2014; Veettil and others, 2017b). This could have significant impacts on glacier mass balance, as previous high-magnitude El Niño events (e.g. 2015–2016) were followed by several years of sustained glacier retreat in the Peruvian Andes (Seehaus and others, 2019). To date, there have been few numerical modelling studies of Andean glaciers, meaning that their future response to climate warming, and how this may interact with ENSO cycles, is highly uncertain.

To address these uncertainties, we assess the sensitivity of Shallap Glacier (dry outer tropics) and Zongo Glacier (wet outer tropics) to past, present and future ENSO conditions. Firstly, we use remote-sensing techniques to construct a geodetic mass-balance time series for Shallap Glacier from 1985 to 2020 and for Zongo Glacier, we utilise a glaciological mass-balance (GMB) time series from 1991 to 2017, produced by the GLACIO-CLIM Observatory. We then use statistical analysis to determine past and present responses of both glaciers to ENSO. The finite element ice flow model, Ùa (Gudmundsson, 2020), is then used to assess the sensitivity of Shallap and Zongo Glaciers to a range of different ENSO patterns that may develop under climate warming.

The respective locations of the study glaciers in the dry and wet outer-tropics produce substantial differences in their precipitation receipts (Sicart and others, 2011; Juřicová and Fratianni, 2018). Our study thereby offers an opportunity to examine whether glacier responses to ENSO varies across the outer tropics. Previous remote-sensing studies have linked terminus changes at both glaciers to ENSO climatic variability (Veettil and others, 2017a; Maussion and others, 2015). At Zongo, Réveillet and others (2015) conducted full-Stokes modelling to assess glacier evolution under different climate projections. However, previous work has not examined the temporal evolution of the ENSO–mass-balance relationship at either glacier, or at other outer-tropical glaciers more broadly. Furthermore, the volumetric response of glaciers in the region to ENSO and its potential future characteristics under a warmer climate has not been assessed. We aim to fill this knowledge gap, by modelling the response of Shallap and Zongo to current and future ENSO scenarios. Furthermore, we test the feasibility of deriving ice thickness and velocity for numerical modelling when field-data are unavailable.

2. Methodology

2.1 Study site

We focus on Shallap and Zongo Glaciers, as they are the only two glaciers in the outer-tropical Andes with the datasets required for numerical modelling (Fig. 1). Shallap is in the southern Cordillera Blanca, Peru, with an area of ~7 km² (Vignon and others, 2003) and spans an altitudinal range of 4700–5700 m a.s.l. (Maussion and others, 2015). Zongo is located in the Bolivian Cordillera Real, in the Huayna Potosi Massif (Wagnon and others, 1999). At ~2 km², Zongo is smaller than Shallap and spans an altitudinal range of 4900–6000 m.a.sl (Wagnon and others, 1999).

Fig. 1. Location of Shallap (b) and Zongo (c) in outer-tropical South America (a). Background imagery is sourced from Planet Labs. The outline of Shallap is shown in blue and Zongo in red.

Shallap and Zongo both experience a wet (September–March) and dry (May–August) season, with accumulation occurring primarily in the wet season and ablation occurring year-round (Sicart and others, 2008). However, their precipitation receipts differ, due to their location in the dry (Shallap) and wet (Zongo) outer-tropics (Sicart and others, 2011; Juřicová and Fratianni, 2018).

2.2 Data sources

The datasets used in the remote-sensing and modelling components of our study are summarised in Table 1.

Table 1. Summary of the datasets used for both Shallap and Zongo

Several datasets were additionally used to investigate the impact of the ENSO and climatic variability on both Shallap and Zongo. These datasets are summarised in Table 2 and discussed below in section 2.3.4.

Table 2. Summary of climate and weather/climate datasets used in this study

EIA, ENSO impact analysis; CCIA, climate change impact analysis; SA, statistical analysis

2.3 Remote sensing and statistical analysis

2.3.1 Changes in glacier length and area at Shallap Glacier

Changes in the length of Shallap's terminus position and overall area were calculated from Landsat imagery, between 1984 and 2020. We selected one image per year from the dry season (between June and August) to reduce the potential for snow patches to be misidentified as glacier ice (Martins and others, 2018). Area was then delineated manually from each Landsat image. Glacier terminus position was also digitised for each year, using the box method to account for uneven retreat (Lea and others, 2014).

There are two main error sources when mapping glacier area and frontal position. The first source is inaccurate geo-referencing of Landsat images. To account for this, static features (i.e. mountain-ridges) in each image were compared with the most recent image in our time series (2020) and all images aligned at the pixel resolution. Secondly, glacier margins can be obscured by snowfall, topographic shadowing and debris cover (Arnaud and others, 2001) and the digitised margin is subject to the interpretation of the person completing the digitising. To quantify this error, the area of Shallap and the terminus length, using the box method, were digitised five times from the 2020 image. We then calculated the standard errors from the five measurements of area and terminus position, yielding an error of ±2.02% (±141.674 m²) and ±0.220% (2.8 m) respectively.

2.3.2. Surface elevation change, ice thickness and geodetic mass balance at Shallap Glacier

We used ASTER DEMs to derive a time series of changes in Shallap's surface elevation and area at different altitudes from 2001 to 2020. While satellite images were available for every year from 2001 to 2020, DEMs covering the entire spatial extent of Shallap were only available for certain years. Therefore, elevation was determined using the nearest available DEM image. Contours were then generated at 200 m intervals across Shallap's surface for each DEM, which were then overlayed with satellite images with the nearest temporal match, enabling annual glaciated area to be calculated for each elevation band (ranging from <4800 to >5800 m) and for the entire glacier.

Surface mass balance for Shallap was only available for 2007–2008 (produced by Gurgiser and others, 2013) and no multi-decadal GMB time series exist for Shallap at the time of writing. Therefore, geodetic techniques were applied. Ideally, ASTER DEMs would be used to determine ice thickness changes and then GMB. However, ASTER DEMs covering Shallap are only available from 2001 to 2021 and the data quality is highly variable. Consequently, to derive a multi-decadal mass-balance time series, volume-area scaling was applied. Volume-area scaling uses a power–law relationship (Eqn (1)) to estimate volume changes that result from area changes. In Eqn (1), V is glacier volume, A area, c a dimensionless coefficient that quantifies volume changes resulting from area change and Y a second coefficient that dictates the degree to which volume scales area (Adhikari and Marshall, 2012).

(1)$$V = c \times \Delta A\Upsilon $$

Initially, values for c and Y of 0.0285 and 1.375, respectively, were used to calculate thickness/volume changes, as they have been previously suggested for temperate valley glaciers (Adhikari and Marshall, 2012). Thickness changes were also calculated from ASTER DEMs from 2001 to 2020 (where available) and were used to compare geodetic thickness changes for the initial values of c and Y. We then experimented with different values of c and Y until we obtained values that produced thickness changes in greatest agreement with those calculated from DEMs. This determined that c = 0.0205 and Y = 1.4 were most appropriate for Shallap.

To further calibrate our thickness changes, our best geodetic thickness changes and the difference between geodetic and DEM-calculated thickness changes across a central flowline were plotted to determine a cost-function that quantified the difference between the two thickness calculations (Fig. 2). The flowline was determined as the central point of the glacier from satellite imagery. A linear trend was applied to prevent overfitting and due to initial analysis suggesting a linear trend provided a good approximation to the data. Geodetic thickness changes were then calibrated using the equation from Figure 2 and used to calculate GMB (MB _GEO; Eqn (2)) from thickness change, Δh, ice density, ρ (assumed as 917 kg m⁻³) and time, Δt.

(2)$$MB_{GEO} = \displaystyle{{\Delta h\rho } \over {\Delta t}} = {\rm \;}\displaystyle{{\Delta V\rho } \over {\Delta A \times \Delta t}}$$

Fig. 2. Relationship between Geodetic calculated thickness changes and the difference between DEM and geodetic thickness measurements. The derived equation was used to calibrate geodetic thickness changes, used to calculate GMB.

With a GMB calculated for each elevation band, a specific mass balance (b _n) for Shallap is determined using Eqn (3) (Wang and others, 2014), where A represents the entire glacier area, b _i the elevation band GMB and s _i the elevation band area.

(3)$$b_n = {\rm \;}\displaystyle{1 \over A} \sum s_{i} b_{i} $$

Errors in MB _GEO and b _n are proportional to errors in the estimate of volume derived from Eqn (1), which is estimated by Eqn (4) (Huh and others, 2017) and the theory of propagation of uncertainty (section 2.5.2).

(4)$$\varepsilon _V = \sqrt {{\left({\mathop \sum \displaystyle{{\partial V} \over {\partial A}}\sigma_V} \right)}^2 + {( \sigma _A) }^2} $$

Percentage errors were calculated for total glacier volume for each year in the time series, with the average percentage error across the time series representing the standard error in volume, b _n and MB _GEO. Thus, error in a given volume, b _n and MB _GEO measurement is determined as ±8.31%

2.3.3 Zongo mass balance and length fluctuations

Zongo's mass balance, recorded at 100 m elevation intervals, and terminus length fluctuations are provided by the GLACIO-CLIM observatory (https://glacioclim.osug.fr) from 1991 to 2017. Mass balance was calculated from ablation stakes and snow pits, with annual length fluctuations determined by GPS monitoring of the terminus. Errors in terminus length are related to the error in GPS positioning, which was ±7 m (Réveillet and others, 2015). Errors in mass balance are discussed in section 2.4.5. Area change measurements were not conducted for Zongo because geodetic techniques were not required to determine mass balance. Instead, glaciological studies since 1991 have resulted in multi-decadal mass-balance measurements distributed across Zongo's entire surface.

2.3.4 Climate data

The ENSO was quantified through the Southern Oscillation Index (SOI) (https://www.ncdc.noaa.gov), which measures the sea-level pressure difference between Tahiti and Darwin, Australia. Pressure differences are a fundamental control on the distribution of warm ocean waters, and so the SOI captures the variations in Pacific SSTs, the driver of ENSO climatic variability (Melice and Servain, 2003). The PDO (Pacific Decadal Oscillation), downloaded from https://www.psl.noaa.gov/Timeseries/PDO/, was also analysed to determine whether its variability impacts the relationship between ENSO and glacier mass balance and length fluctuation. The PDO is another ocean-atmospheric climate phenomenon that impacts the Pacific Ocean and regularly interacts with the ENSO to alter Tropical Pacific climate (Li and others, 2020). Therefore, PDO forcing may amplify/obscure the ENSO signal and its impact on glacier dynamics.

The GLACIO-CLIM observatory collects hourly temperature and precipitation data from a meteorological station on Zongo's terminal moraine. However, data are only available from 2003 and are therefore unsuitable for analysing the ENSO–mass-balance relationship over Zongo's entire data collection period (1991–2017). Like GLACIO-CLIM, AClnn Weather Stations (www.chacaltaya.edu.bo/) collected meteorological data for Shallap through sensors located at 4790 m on the glacier's surface. However, data were only collected from 2010 and equipment failure results in low measurement accuracy and inconsistent data. Due to the poor quality and short records of data at the glaciers themselves, we obtained average annual temperature and number of rain days for La Paz, Bolivia (from 1991 to 2017), and annual average temperature and total annual precipitation for Huaraz, Peru (from 1985 to 2020) (https://www.climate-data.org). La Paz and Huaraz are located ~30 and ~20 km from Zongo and Shallap, respectively. Consequently, they may not be directly representative of conditions at each glacier but do provide the best-available overview of general climatic conditions at Zongo and Shallap. Due to the differences in altitude, the climate data were not used for statistical analysis, and only used to infer causes for changes in the relationship between glacier dynamics and the ENSO.

2.3.5 Statistical analysis

MATLAB functions were used to analyse the response time of the mass balance (total and each elevation band) and length fluctuations of Shallap and Zongo to the ENSO. Firstly, Change Points (Matlab 2021a) was used to identify statistically significant changes in the mean of the SOI, mass balance and length fluctuation time series and hence whether significant changes in the SOI were followed by changes in glacier mass balance/terminus fluctuation. Next, changes in the lag time between the ENSO changes and mass balance/terminus fluctuation were visualised with local climate data and the PDO to determine if changes in the lag time between mass balance and ENSO change points coincided with certain weather and/or PDO conditions, using MATLAB function Find Delays (Matlab 2021b). Finally, MATLAB function cross-correlation (Matlab 2021c) was used to identify the strongest average lag between the SOI and glacier mass balance/length fluctuation.

2.4 Numerical modelling

2.4.1 Model description and datasets

Úa is a finite element ice flow model, written in MATLAB, and was used to assess the sensitivity of Shallap and Zongo to the ENSO (Gudmundsson, 2020). Úa uses a vertically integrated formulation of the ice dynamic equations to solve ice flow for ice sheets, ice caps and mountain glaciers (Gudmundsson, 2020), using the Shallow Ice Stream Approximations (Hill and others, 2020). It has an adaptive mesh, meaning that key areas can be resolved at a higher resolution. Úa has been widely applied to assess the response of glaciers to climate change (e.g. Hill and others, 2020; De Rydt and others, 2021). Úa requires surface mass balance, surface velocity, ice thickness, surface elevation and basal topography as initial inputs. For both Zongo and Shallap, a number of these datasets were not available as pre-existing products and therefore had to be derived specifically for this work. The datasets used, and the approaches used to derive individual datasets, are described below.

2.4.2 Surface elevation and basal topography

Surface elevation data for Shallap were sourced from the ASTER DEM at 30 m resolution, downloaded from https://earthdata.nasa.gov. Error in surface elevation was ±9.5 to 10 m (Tachikawa and others, 2011). Zongo's surface elevation was computed from stereo-pairs of aerial images from 1997 and 2006 by Réveillet and others (2015). The data have a spatial resolution of 25 m and a date of 2006. The data were provided as a series of discrete points, which we converted into a surface elevation raster covering the entirety of Zongo's surface, using linear interpolation. Error in surface elevation was estimated as the error in GPS positioning used to geo-reference the stereo-images, which is ±7 m (Réveillet and others, 2015).

Basal topography for Shallap was calculated by subtracting the calculated ice thickness from the surface elevation. As a result, errors in basal topography are proportional to those in ice thickness (see section 2.4.5) and surface elevation (2.4.2). Zongo's basal topography (as of 2006) was provided by Réveillet and others (2015). Basal topography was calculated using ice thickness measurements, determined from ice-penetrating radar, at accessible reaches of the glacier (Réveillet and others, 2015). This produced point measurements of basal topography that were converted into a raster covering the entire glaciated area using linear interpolation. The error in basal topography is ±5 m (Réveillet and others, 2015).

2.4.3 Ice velocity

Millan and others (2019) derived ice velocity for the entire Cordillera Blanca using feature tracking and high-resolution satellite imagery for 2017–2018; the data cover the entire spatial extent of Shallap at a spatial resolution of ~50 m. Error in ice velocity was estimated using the mean of the standard deviation in the  x and y component of velocity (Millan and others, 2019); this is calculated as ±6.1 m a⁻¹.

Zongo Glacier is not covered by either the Millan and others (2019) data or ITSLive (https://its-live.jpl.nasa.gov/). Furthermore, directly measured ice velocities, from stake surveys, are available only for the lowest 300 m of Zongo glacier (https://glacioclim.osug.fr). Consequently, we generated ice velocities covering the entirety of Zongo using Im-GRAFT (Grinsted, 2021) within MATLAB, which calculates ice velocities from the displacement of prominent ice surface features (i.e. crevasses) between two satellite images. We used Planet Labs imagery (https://www.planet.com/) at 3 m resolution to generate the ice velocities, as this was the highest spatial resolution imagery available. Cloud cover and infrequent imaging over Zongo resulted in only three calculations of ice velocity being possible: 12/04/2020–11/12/2020, 17/08/2018–07/05/2019 and 16/08/2010–28/07/2011. However, results that produced velocity measurements over a substantial percentage of Zongo's surface (i.e. >80% coverage) were only obtained from the 17/08/2018 and 07/05/2019 image pair. We then used the velocities measured directly from ablation stake displacements in the lower terminus for 2017–2018 to calibrate the Im-GRAFT velocities, as these were closest in date to the Planet imagery. A calibration factor was then calculated as the ablation stake velocity minus our Im-GRAFT velocity measurement for each stake measurement and applied across our ImGraft-derived velocities, to ensure the best fit to measured data. Calibrated velocities fit the expected spatial distribution across the glacier surface, i.e. maximum velocities around the ELA, lower velocities up-ice and higher velocities down-ice (Phillips and others, 2013).

Errors in our ice velocity were calculated from the discrepancy between our calibrated values and ablation stake displacement measurements. This indicated errors in velocity of up to ±11.6 m a⁻¹, with a mean error of ±1.5 m a⁻¹, which transcribes to an average percentage error of ±14.0%. Consequently, errors in our ice velocity are high. However, with no ice velocity data in the Bolivian Andes and poor spatial and temporal coverage by high-resolution satellite imagery, our approach best utilises the limited resources in this data-sparse region.

2.4.4 Ice thickness

Ice thickness for Shallap was inverted for using the vertical velocity profile (Eqn (5)) (Hooke, 2019). U _s, ice surface velocity, and U _b, basal ice velocity, were determined using measured ice velocity data for Shallap provided by Millan and others (2019). Φ, the surface slope, was determined from surface elevation data and g, acceleration due to gravity, and ρ, ice density, defined as 9.81 m s⁻¹ and 900 kg m⁻³, respectively. Additional parameters in Eqn (5) such as S _f (a coefficient that represents the glaciers geometry as a ratio of its length and width), A (the ice flow rate factor) and n (an additional coefficient) were unknown for Shallap and had to estimated. We therefore generated multiple different thickness estimates for Shallap by redefining the estimates for A, n and S _f.

(5)$$h = \root {n + 1} \of {\displaystyle{{U_{\rm s}-U_{\rm b}} \over {( 2A/n + 1) {( S_f\rho g\Phi ) }^n}}} $$

The ice thicknesses derived from Eqn (5) were then used to invert for ice velocities in Úa, with the ice thickness distribution that returned ice velocities most concurrent with the measured velocities from Millan and others (2019) selected as the ice thickness distribution that would be used in the numerical modelling experiments. When solving for ice velocities in Úa, basal slipperiness, C (using Eqn (6)) and A needed to be estimated as if these variables were left as free for Úa to solve, the optimisation methods applied by Úa would lead to all ice thickness converging well with measured ice velocities. A prior estimate for C, using Eqn (6), was made, where U _b is basal velocity (assumed as 25% of surface velocity) and t _b basal shear stress (Ng and others, 2018), calculated using Eqn (7).

(6)$$C = \displaystyle{{U_{\rm b}} \over {t_{\rm b}}}-1$$

(7)$$t_{\rm b} = f\rho g{\rm sin}\Phi $$

Errors in thickness are proportional to the mean difference in measured and calculated ice velocities, which was determined as ±1.13 m a⁻¹. This corresponded to an average percentage error of ±3.1% in Shallap's ice thickness.

Zongo's ice thickness was calculated using QGIS 3.0 raster calculator to subtract the basal topography from the surface elevation. Error in ice thickness was estimated from the error in GPS positioning and error in the ‘radar-pick’ (Réveillet and others, 2015). Due to the propagation of uncertainty, mean errors in ice thickness were ±8.6 m.

2.4.5 Mass-balance data

Shallap's surface mass balance was acquired from Gurgiser and others (2013) and was determined from a glacier mass-balance model, with model outputs calibrated against ablation stake measurements. The modelled mass balance had a spatial resolution of ~1 m and produced mass-balance measurements every 2 months from 2007 to 2008, which were averaged to generate annual mass balance for the entire glacier and at elevation bands of 50 m. Error in mass balance was calculated by Gurgiser and others (2013) using Eqn (8), where b ₂₀₀₇ is the 2007 mass balance, b ₂₀₀₈ the 2008 mass balance and n the number of in-field measurements (Gurgiser and others, 2013). With a difference in mass balance of 0.83 m w.e. a⁻¹ and 20 in-field measurements, the error in mass balance is ±0.034 m w.e a⁻¹.

(8)$$SE = {\rm \;}\displaystyle{{b_{2007}-b_{2008}} \over n}$$

The GLACIO-CLIM Observatory provides annual mass balance for Zongo at 100 m elevation bands (https://glacioclim.osug.fr). Calculation methods vary across the altitudinal range of Zongo, with mass balance at 4900–5300 m calculated from ablation stakes, 5300–5500 m calculated from linear interpolation and >5500 m calculated from snow pits. Using surface elevation data, we interpolated mass balance across Zongo's entire surface. SMB data from 2006 were used, as it corresponds to the data collection period for basal topography, ice thickness and surface elevation. Errors in mass balance are unknown and not calculated by the GLACIO-CLIM Observatory. Therefore, errors are estimated using the same procedure as Gurgiser and others (2013) (i.e. Eqn (9)). Errors in mass balance were thereby found to be ±0.011 m w.e a⁻¹.

(9)$$SE = \displaystyle{{b_{2005}-b_{2006}} \over n}$$

2.4.6 Mesh generation

For both Shallap and Zongo glaciers, a computational mesh was defined and was locally refined to produce smaller element sizes around areas of high effective strain-rate gradients (where effective strain rates were iteratively derived from velocities in Úa). Different meshes with element spacings of 5, 10, 25, 50 and 100 m were produced. A convergence study was then conducted for Shallap and Zongo to ascertain the optimal balance between run time and resolution, and to ensure that the results were unaffected by mesh resolution. A 25 m mesh was used for both Shallap and Zongo glaciers because this value produced minimal change in output results in an adequate timeframe.

Due to dataset availability, modelling was not undertaken for the entirety of Shallap's surface. Specifically, we needed to invert for ice thickness, due to the lack of direct measurements, and the assumptions used in our calculation of ice thickness did not hold across Shallap, due to its large width (Farinotti and others, 2017). Additionally, a nunatak partitions Shallap into two flow units, and so we focused on Shallap's faster moving portion. By only modelling one of Shallap's flow streams, we must consider mass conservation and ice flow from adjacent parts of the glacier. To do this, we used a fixed boundary condition in Úa where ice velocity was fixed as the average for the flow stream (6.75 m a⁻¹). Modelling was undertaken for the entirety of Zongo glacier because of its simpler geometry and availability of required datasets. As such a mesh was defined for the entire glacier surface and natural boundary conditions were applied at its margins (i.e. ice is allowed to flow in and out).

2.4.7 Inversions

Inversions were performed to calculate distributions of A (the rate factor in Glen's Flow Law) and C (basal slipperiness) and utilise initial estimates of 2 × 10⁻⁸ for A and 0.001 for C across all nodes of the mesh. For Shallap, estimates for C were defined using Eqn (6), as used in the derivation of ice thickness. Inversions for A and C generate a cost function (Eqn (10)) that enables identification of the best estimates for A and C, which are when the difference between measured and calculated velocities, I, and regularisation term, R, are at their minima.

(10)$$J( {q( p ) , \;p} ) = I( {q( P ) } ) + R( p ) $$

Mean discrepancies between measured and calculated velocities were calculated as 0.687 and 0.804 m a⁻¹ for Shallap and Zongo, respectively.

2.5 Sensitivity experiments

2.5.1 Model scenarios

Our sensitivity experiments assessed responses of Shallap and Zongo to a range of potential future ENSO forcing scenarios (summarised in Fig. 3). Each transient experiment was run forward for 50 years, with an initial time step of 0.1 years. Adaptive time stepping was applied, so that the time step can automatically be decreased if the model is not converging or increased where the model is converging well.

Fig. 3. Different mass-balance scenarios for Shallap (a) and Zongo (b) in our modelling experiments. In experiment (v), the mass balance for experiment (iii) Super El Nino ELA Rise is used for the stronger mass El Nino mass balance.

The model scenarios were:

1. (i) Control run – neutral mass balance for 50 years; specifically, the 2006 mass balance for Zongo and 2007–2008 for Shallap. This is the control run, to which results of experiments (ii), (iii), (iv) and (v) were compared.

2. (ii) High-magnitude El Niño and La Niña events characteristic of the 1984–2020 period for Shallap, and 1991–2017 for Zongo, as identified using statistical analysis. For Zongo, the high-magnitude El Niño mass balance was defined as the average of mass balances from 1994, 1997 and 2015, and La Niña mass balance as the average from 1996, 2000, 2007 and 2014. Shallap's mass-balance data (Gurgiser and others, 2013) only produce mass balance for two years (2007–2008). Therefore, the GMB time series was analysed to determine strong El Niño and La Niña mass-balance years, identified as 2010 and 1998 for the El Niño, and 2011 and 2009 for La Niña. This produced an average mass balance for each event, which can be compared with the time series average to determine the subsequent change in mass balance that the El Niño and La Niña instigate. This is then applied to the 2007–2008 mass balance to approximate strong El Niño and La Niña mass-balance distributions.

3. (iii) Two ‘Super’ El Niño events, where mass loss is doubled. These are:

   1. (1) Niño-1: Accumulation is halved and ablation doubled, but the ELA remains at the same altitude as the regular El Niño mass balance.

   2. (2) Niño-2: Accumulation is halved, ablation is doubled and the ELA rises 50% up-ice compared to the regular El Niño mass balance.

   This enables us to test the sensitivity of Zongo and Shallap to changes in net mass balance and ELA rises.

4. (iv) Weaker El Niño at higher altitudes. The exact altitude that El Niño forcing may become weakly transmitted is unknown, so the entire altitudinal range of each glacier is examined at 100 m intervals. El Niño forcing (experiment (ii)) was then applied to the area below each 100 m altitudinal threshold and a ‘neutral’ mass balance (experiment (i)) applied above it.

5. (v) El Niño-dominated ENSO, with stronger and longer El Niño. In the future, El Niño events may become stronger and/or longer than La Nina events, leading to an El Niño-dominated ENSO cycle. These scenarios are conducted under a control mass balance between ENSO events and a simulated mass balance under climate change which we approximated by extrapolating the historic cumulative mass balance. We explore three ENSO scenarios:

   1. (3) More regular El Niño events than La Niña.

   2. (4) Stronger El Niño events than La Niña.

   3. (5) Stronger and longer El Niño than La Niña.

The phasing of ENSO events for experiment (v) are summarised in Table 3, where one ENSO cycle is displayed. Cycles are then repeated until 50-years have been reached. If a scenario cycle length is not an integer factor of 50, then the maximum number of whole cycles is used, with the required number of neutral mass-balance years called until 50-years of mass balances have been generated.

1. (vi) (v) Longer and stronger El Nino with climate change. Experiment (v) was repeated, but the mass balance called between ENSO events had climate warming added. To do this, we calculated the trend in measured cumulative mass-balance data over time and extrapolated it forward by 50 years (Fig. 4). While this is a simplification, the high R ² values (>0.95) and low RMSE values provide confidence, and the approach enables us to test sensitivity when detailed climate forecasts are not yet available for Zongo and Shallap. Cumulative mass balance for a given year in our forward runs is thus the neutral balance plus the extrapolated trend in cumulative mass balance. ENSO events are then applied as an additional increase or decrease in each modelled year's mass balance and were calculated as deviations from the control mass balance using our measured data.

Table 3. Phasing of ENSO cycles used in our sensitivity experiments

Blue refers to La Nina mass balances, black neutral mass balances, light red regular El Nino mass balances, the medium red super El Nino 1 mass balances and dark red super El Nino 2 mass balances.

Fig. 4. Modelled and measured cumulative mass balances for (a) Shallap and (b) Zongo. Changes in cumulative mass balance are used for estimating future mass balance under climatic change.

2.5.2 Modelling experiment errors

Standard error, σ, in modelled results will arise through the propagation of uncertainty, ɛ, in each of the input variables used in Úa (Eqn (11)).

(11)$$\sigma = \sqrt {{\left({\displaystyle{{\varepsilon_h} \over h}} \right)}^2 + {\left({\displaystyle{{\varepsilon_v} \over v}} \right)}^2 + {\left({\displaystyle{{\varepsilon_b} \over b}} \right)}^2 + {\left({\displaystyle{{\varepsilon_s} \over s}} \right)}^2 + {\left({\displaystyle{{\varepsilon_{as}} \over {as}}} \right)}^2} $$

Applying Eqn (11) and the propagation of uncertainty assumes errors in each variable are uncorrelated, a common simplification for complex datasets (Batstone, 2013). Furthermore, while the standard deviation is usually used for each variable in Eqn (11), standard deviations are unknown for some datasets. Therefore, percentage errors are used instead, another acceptable alternative for datasets where methods of derivation/calculation do not permit calculation of the standard deviation (Batstone, 2013). Using Eqn (11) and the assumptions outlined above, errors in modelled results are calculated as ±15.6 and ±17.5% for a given output result for Shallap and Zongo, respectively.

3. Results

3.1 Remote sensing

3.1.1 Changes of Shallap and Zongo Glaciers

Shallap and Zongo experienced strong decreases in cumulative mass balance and marked terminus retreat during the study period (Fig. 5). Shallap's cumulative mass balance decreased by 1.4 × 10⁴ ± 1162 mm w.e from 1986 to 2020, corresponding to an average annual specific mass-balance decrease of 404.4 ± 34 mm w.e a⁻¹ and an average thinning rate of 2.71 ± 0.22 m a⁻¹. Zongo displayed a similar cumulative mass-balance loss of 1.1 × 10⁴ ± 650 mm w.e from 1991 to 2017; however, annual average specific mass balance loss was higher, at 524.3 ± 41 mm w.e a⁻¹. Both glaciers underwent net terminus retreat: Shallap retreated by 870 ± 2 m (−25 ± 0.01 m a⁻¹) and Zongo by 322 ± 7 m (−12 ± 7 m a⁻¹), over their respective study periods. Statistical analysis determined that frontal position changes on Shallap lag mass-balance changes by 1 year (p-value<0.05). For Zongo, no statistically significant lag time between mass balance and length changes were detected (p-value>0.05).

Fig. 5. Cumulative mass balance (cmb) and cumulative terminus retreat (dl) for Shallap and Zongo over their respective study period.

3.1.2 Mass-balance response to the ENSO

The mass balances of both Shallap and Zongo responded strongly to the ENSO. Statistical analysis shows that ENSO results in a change in glacier mass balance after 1 year for Shallap, and within the same year for Zongo (p-value<0.05). Shallap and Zongo's mass balances consistently corresponded to the ENSO's variability, not just to high-magnitude events. Change points in both glaciers' mass-balance time series are concurrent with SOI change points and, despite ENSO maximum events occurring every 2–7 years, SOI and glacier mass-balance change points occur with intervals of 0–2 years for both glaciers.

The lag time between observed ENSO forcing and observed mass-balance response changed during both Shallap and Zongo's data collection periods. At Zongo, mass-balance changes occurred in the same year as ENSO changes (i.e. zero lag) from 1991 to 2003, but after 2003, mass-balance changes occurred 1 year after ENSO changes (i.e. a 1-year lag; Fig. 6a). This change in lag was coincident with average annual temperature in La Paz rising by ~1°C between 2003 and 2017 (Fig. 6a). At Shallap, the lag between ENSO forcing and mass-balance response varied through the study period: zero-lag response times were observed from 1988 to 1991, followed by a 1-year lag time for 1991–1999 (Fig. 6b) and then a period of decoupling between 2005 and 2010, which is discussed below. One-year lag times (from 1991 to 1999) occurred during a period of relatively cool temperatures and high precipitation receipts (Fig. 6b) despite pronounced ENSOs and a strong El Niño event from 1997 to 1999 (Fig. 6d). Furthermore, ENSO peaks occurred year-round and were not confined to the wet or dry season. We assessed the potential for the PDO counteracting the ENSO signal and found that during the period of 1-year lag times, the SOI and PDO were weakly correlated (r = −0.314; R ² = 0.094; p-value <0.05), suggesting PDO forcing is counteracting the ENSO and having a stronger influence on regional climate and Shallap's mass balance. Thus, Shallap's slower mass-balance response times to the ENSO occur when the PDO is dominating local climate/weather.

Fig. 6. Local climate conditions, SOI and mass-balance change points for Zongo (a) and Shallap (b). Displayed for reference are glacier mass balance and SOI timeseries, covering the data collection periods for Zongo (c) and Shallap (d).

For Zongo, we identified no periods of decoupling (periods where mass-balance change points precede SOI change points), with the ENSO continually influencing its mass balance between 1991 and 2017. In contrast, Shallap's mass balance decoupled from ENSO forcing between 2005 and 2010 (Fig. 6b). During this period, changes in mass balance preceded changes in SOI. This period of decoupling was concurrent with high average temperatures and low annual rainfall, which are usually associated with El Niño, but during this period three positive (La Niña) peaks occurred (Fig. 6d). Furthermore, from 2005 to 2010, mass balance was out of phase with the SOI compared to previous time periods, with more negative mass-balance fluctuations occurring during positive (La Niña) SOI values and less negative mass-balance fluctuations during negative (El Niño) SOI values (Fig. 6d).

During the period of decoupling, the SOI and PDO were more negatively correlated (r = −0.499; R ² = 0.283; p-value<0.05), suggesting positive PDO phases coincided with negative ENSO phases and vice versa. However, the SOI and PDO demonstrated even stronger negative correlation post-2010 (r = −0.567; R ² = 0.315; p-value<0.05), but mass balance responded rapidly to the ENSO. Therefore, our data suggest that PDO forcing alone may not have been responsible for the 2005–2010 period of decoupling. Instead, the timing of the ENSO peak may have forced the period of decoupling. The ENSO peak intensity during Shallap's period of decoupling occurs consistently during the wet season. From 2005 to 2010, six ENSO peaks occur, with three being in January and one in February, April and September. This timing of the peak ENSO intensity is not consistent with the longer record: between 1986 and 2020, ~40% of ENSO peaks occurred during the dry season and 60% during the wet season.

3.2 Numerical modelling

3.2.1 Future glacier responses to the ENSO

The results of all modelling experiments are summarised in Table 4. Specific details of each experiment and their implications are discussed in detail in the following sections.

Table 4. Summary of key sensitivity experiments results from modelling in this study

Results are expressed as a change when the original glacier (t = 0), compared when the glacier at the end of the model run (t = 50).

Control run – neutral mass balance (i). For Shallap, we define 2007–2008 as the neutral mass balance, to which other experiments will be compared. After 50 years of neutral mass balance, Shallap retreated by 670 m and ice volume reduced by 46% (Fig. 7a). Most of the glacier thinned, with mean elevation loss of −89.6 m, but central areas of the glacier thickened by an average of 26.1 m.

Fig. 7. Shallap ice thickness changes (dh) from the initial ice thickness (t = 0) at the end of each experiment's model run (t = 50).

Zongo's response to its neutral mass perturbation (2006) was terminus retreat of 58 m and minor area loss at high elevations (>5900 m a.s.l.; Fig. 8a). However, average elevation change across the glacier was +10.1 m and ice volume increased by ~13%. Thickening occurred at most elevations, with some thinning at the terminus, high altitudes and at the glacier margins (Fig. 8a). Overall, Zongo had a positive mass balance under its baseline (2006) conditions.

Fig. 8. Zongo ice changes (dh) from the initial ice thickness (t = 0) at each of the experiment's model run (t = 50).

Present-day high-magnitude ENSO response (ii). Zongo and Shallap displayed very different responses to the application of present-day, high-magnitude El Niño and La Niña mass balances, which were applied over the 50-year modelling period. Over 50 years, Shallap's terminus retreated in response to both El Niño (~1000 m) and La Niña (~170–450 m) events (Figs 7b, 7c). During the La Niña run, Shallap's terminus retreat was spatially variable, with the southern section of the terminus retreating less (~170 m), than the north (~450 m). Shallap also thinned and lost volume during both the El Niño and La Niña runs. For the El Niño run, average surface elevation change was −106.42 m (~80 m more than during the control run) and 55% of its volume was lost (9% more than during the control run), compared to −73.00 m and 37% for average elevation and volume change, respectively, for the La Niña run. Isolated areas of thickening and mass gain were detected between 5000 and 5200 m, where the ice is initially thickest (Fig. 7c).

Zongo's response to high-magnitude ENSO mass-balance perturbations was much less marked than at Shallap. In response to the El Niño run, the terminus retreated ~407 m (457 m more than during the control run) but average ice thickness (+0.602 m) and volume (~1%) increased (9.5 m and 12% less than during the control run). Maximum thickening occurred around the ELA (up to 282.1 m) and surface thickening occurred across the entire accumulation zone and proximal to the terminus (Fig. 8c). The La Niña run triggered an 18.5% increase in ice volume and average ice thickening of +14.9 m (5.5% and 4.8 m more than during the control run).

Super El Niño (iii). Two super El Niño scenarios were applied. First, in Niño-1 the magnitude of accumulation halves and ablation doubles but the ELA remains at the same altitude as in the regular El Niño used in experiment (ii). Second, Niño-2 is an El Niño event where the magnitude of accumulation halves, the ablation doubles and the ELA rises up 50% of the glacier's altitudinal extent: this represented a rise of 400 and 300 m up-ice for Shallap and Zongo, respectively. By simulating these two possible future El Niño behaviours, we can assess the relative sensitivity of Shallap and Zongo to accumulation/ablation changes and to an ELA rise.

Shallap's response was similar between experiments (ii) and (iii): its terminus retreated by ~100 m more in experiment (iii) and lost 4% more volume, giving a total loss of ~59% (Fig. 7d). Zongo, on the other hand, showed a far stronger response in experiment (iii). Terminus retreat was more than twice experiment (ii), at 920 m (Fig. 8d), and volume reduced by 26%, compared to 1% volume gain in experiment (ii).

Both Zongo and Shallap showed greater retreat and mass loss in response to Niño-2 compared to Niño-1. Shallap's terminus retreated 1613 m (~1100 m in Niño-1), average surface thinning was 150 m and volume reduced by 76% (~59% in Niño-1). Thus, Shallap is very sensitive to changes in the ELA and less sensitive to changes in net mass balance. Specifically, the central area (Fig. 7d) of comparatively thick ice between 5100 and 5300 m is crucial to its response, so that raising the ELA to fully encompass this area causes substantial ice loss, whereas simply increasing ablation and decreasing accumulation had only a minor impact. Zongo displayed a more complex response to the Niño-2 simulation. The west side of the terminus retreated far less than the east, 1145 and 1500 m (920 m in Niño-1), respectively (Fig. 7e). In addition, volume loss and mean surface thinning were 48% (26% in Niño-1) and 38.9 m across the entire glacier, reflecting a twofold increase from the Niño-1 response. Thus, Zongo is sensitive to both changes in ELA and mass balance.

Varied maximum altitude of ENSO forcing (iv). It has been suggested that climate warming may weaken the impact of ENSO at higher altitudes (Rabatel and others, 2013), but the altitude at which this may occur is uncertain. We therefore assess the sensitivity of both glaciers to changes in the altitude at which ENSO forcing occurs.

Shallap's volume loss increased by >0.02 km³ when the maximum altitude of ENSO forcing was changed from 5300 to 5400 m (Fig. 9a), suggesting that this altitude is critical for determining glacier response to the ENSO. Furthermore, Shallap retreated by ~150 m when forcing was moved from 5400 and 5500 m (Fig. 7a).

Fig. 9. Volume change (a) and retreat (b) for Shallap and Zongo in response to changes in the maximum altitude of ENSO forcing. Mass balance and area distribution for Shallap (c) and Zongo (d) are also shown. The black line in each figure refers to the glaciers ELA.

At Zongo, increasing the maximum elevation of ENSO forcing led to greater terminus retreat (Fig. 9b). Zongo gained mass (Fig. 9b) for all scenarios, except when ENSO was applied to all but the highest 100 m of Zongo's altitudinal extent, i.e. areas above ~5800 m (Fig. 9b). This suggests Zongo is sensitive to ENSO forcing reaching its high-altitude accumulation zone.

Stronger and longer El Niiño (v.1). It has been suggested that stronger and longer El Niño events, relative to La Niña, may form in the future as climate warms (Cai and others, 2014). Therefore, we simulated ENSO phases where (a) El Niño events last twice as long as La Niña; (b) El Niño events are twice as strong as La Niña; and (c) El Niño events twice as strong as La Niña and lasting twice as long. We do this under climate change and neutral climatic conditions, to assess how changing El Nino cycles interact with climate change.

Neither Shallap nor Zongo displayed sensitivity to longer El Niño events when the strength of El Niño events remained unchanged, and a neutral mass balance was called between ENSO cycles (Figs 7g, 8g). Zongo's volume increased by 12%, with no terminus retreat, meaning that volume change was very similar to that experienced during the control run (+13%) and was greater than for the El Niño run (+1%). A similar response occurred at Shallap: during the longer El Niño simulation; volume loss was 4% less than during the control run and the terminus retreated ~10 m less (Fig. 7g). At both Shallap and Zongo, therefore, mass gains are balanced by mass losses, as long as La Niña mass balances are of the same magnitude as El Niño, even if El Niño events persist for longer.

Both glaciers were sensitive to the stronger and the stronger and longer El Niño simulation. When El Niño phases were stronger, the terminus of Zongo retreated by ~1000 m and volume reduced by 41% (Fig. 8e). Furthermore, when El Niño events were stronger and longer, a further 3.5% ice-volume loss (48.5% in total) occurred, and the terminus retreated by a further ~450 m. Thus, the response of Zongo to the ENSO was dominated by the strength of El Niño events, rather than the length.

In response to the stronger El Niño, Shallap's volume decreased by ~0.2 km³ (46.8%), only 0.8% more than during the control run (Fig. 7e), and the terminus retreated an additional 160 m (170 more than during the control run). Thus, stronger El Niño events primarily resulted in terminus retreat, as opposed to overall ice-volume loss. Shallap underwent additional retreat and ice volume loss when El Niño events were both stronger and longer (Fig. 7h). The terminus retreated by an additional 170 m (1000 m in total and 330 m more than in the control run) and a further 3% of the glacier volume was lost (49.5% in total and 3% more than during the control run) compared to the stronger run. Thus, like Zongo, Shallap is also more sensitive to increases in the strength of El Niño events than the length.

Stronger and longer El Nino under climate change (v.2). We now recreate the stronger and longer El Niño simulation from section 3.2.1.5, but use a simulated mass balance under climate change, instead of the control mass balance, between ENSO cycles. When the El Niño events were stronger, with a rise in the ELA, Shallap disappeared entirely by the end of the 50-year simulation. Hence, Shallap is most sensitive to the ENSO when El Niño phases become stronger relative to La Niña. The length of El Niño events, even if they are stronger, had no additional impact on the speed of glacier recession (Fig. 10). In addition, Shallap recorded an additional 0.75 km³ of volume loss when its ELA was increased in response to the El Niño (as opposed to when just accumulation and ablation were changed and the ELA remained fixed), implying that Shallap is sensitive to ELA rises and less sensitive to changes in the magnitude of ablation and accumulation (Fig. 10).

Fig. 10. Total ice volume changes at a given model run time for Shallap under the ENSO conditions outlined in experiment (iv.2). Run 1 is with no ENSO forcing, run 2 is with regular ENSO forcing, run 3 is with longer regular El Niños relative to La Nina, run 4 is when El Nino phases are the same length as La Nina but stronger and with a fixed ELA, run 5 is when El Nino phases are longer than La Nina, stronger and with a fixed ELA, run 6 is when El Niños are the same length as La Nina but stronger and with a rise of the ELA and run 7 is when El Niños are longer than La Niñas, stronger and with a rise of the ELA. Errors are not displayed to assist with the clarity of the figure.

Like Shallap, it was also possible for Zongo to disappear entirely when the strength of the El Niño events increased relative to La Niña (Fig. 11). Zongo disappeared after 25 years when El Niño strengths increased with a rise in the ELA, and after 35 years when the strength of El Niño increased but the ELA remained the same (Fig. 11). Furthermore, the stronger El Niño event simulations (runs 4–7) resulted in mass loss that is more rapid and extreme than that of runs 1–3, suggesting that while Zongo is more sensitive to rises of the ELA and ablation occurring in the steady-state accumulation zone, it is also sensitive to changes in the amount of accumulation and ablation occurring in each respective zone. Moreover, similar to Shallap, Zongo was most sensitive to changes in the strength of El Niño, rather than length, with the results of experiments 2 and 3, 4 and 5, and 6 and 7 being the same, despite variable changes to the length of their uniform ENSO forcing (Fig. 11).

Fig. 11. Total ice volume changes at a given model run time for Zongo under the ENSO conditions outlined in experiment (iv.2). Run 1 is with no ENSO forcing, run 2 is with regular ENSO forcing, run 3 is with longer regular El Ninos relative to La Nina, run 4 is when El Nino phases are the same length as La Nina but stronger and with a fixed ELA, run 5 is when El Nino phases are longer than La Nina, stronger and with a fixed ELA, run 6 is when El Ninos are the same length as La Nina but stronger and with a rise of the ELA and run 7 is when El Ninos are longer than La Ninas, stronger and with a rise of the ELA. Errors are not displayed to assist with the clarity of the figure.

In runs 1–7, Zongo initially retreated rapidly, with retreat levelling off and reducing during the 50-year model run time (Fig. 11). By contrast, while the rate of volume loss does decrease over time for Shallap, it does not level off (Fig. 10).

In addition, the relative proportion of the glacier that is melted in response to run 1 (no ENSO forcing) is significantly higher for Shallap than for Zongo. Approximately 60% of Zongo was lost after the 50-year modelling period, substantially less than the 89% of glacier loss observed for Shallap (Figs 10, 11). However, the increase in volume loss was higher for Zongo than Shallap as increasingly stronger El Niño phases occurred throughout runs 2–7, with total glacier recession occurring faster for Zongo. As a result, Shallap is more sensitive to climatic warming and the subsequent gradual decreases in net mass balance. Zongo is more sensitive to increasingly stronger, and more negative, El Niño mass-balance perturbations.

4. Discussion

4.1 Observed mass-balance and terminus position changes at Shallap and Zongo

Both Shallap and Zongo exhibited negative mass balance during the study period. Shallap lost mass at a rate of 404.4 ± 34 mm w.e a⁻¹, from 1986 to 2020 and never exhibited positive mass balance (Fig. 5). As with other glaciers across the Tropical Andes, changes at Shallap are likely due to observed climate warming of ~0.26°C decade⁻¹ (Mark and Seltzer, 2005). Our data give an average thinning rate of 2.71 ± 0.22 m a⁻¹, which is lower than the average for the Cordillera Blanca at between −3.1 and −7.7 m a⁻¹ from the late 20th century to present (Hastenrath and Ames, 1995; Casassa and others, 2007; Emmer and others, 2015; Wigmore and Mark, 2017). We attribute Shallap's lower thinning rate to its comparatively large ice volume, despite its comparatively low maximum altitude of 5700 m.a.s.l (Fig. 1). Like Shallap, Zongo's cumulative mass balance became increasingly negative through the study period (Fig. 5), reducing at a rate of 524.3 ± 41 mm w.e a⁻¹. This is substantially less than the Bolivian Andes average of 736–1380 mm w.e a⁻¹ (Soruco and others, 2009) and suggests that Zongo is less sensitive to climate change than other Bolivian glaciers. This most likely reflects its relatively large size and high maximum altitude (~6100 m), which has enabled it to retain a substantial accumulation zone (Ramirez and others, 2001).

Zongo's rate of cumulative mass loss was ~120 mm w.e a⁻¹ greater than Shallap, suggesting Zongo is more sensitive to climate warming. Furthermore, at Zongo, we observed positive mass-balance values during the study period, which were seen at all elevations above 5000 m, and its ELA can rise/fall by hundreds of metres annually (Fig. 5). The higher sensitivity of Zongo, compared to Shallap, likely reflects the location of Zongo in the wet outer-tropics (Wagnon and others, 2001): here, the glacier SEB is strongly influenced by the timing of the wet season and subsequent changes to albedo, cloud cover and radiation receipts (Wagnon and others, 2001). Moreover, with mass balance dictated by precipitation, wet climate glaciers are also sensitive to temperature and its impact on the SLA (Wagnon and others, 2001).

Both Shallap and Zongo underwent significant net retreat during the study period, at rates of −25 ± 0.01 and −12 ± 7 m a⁻¹, respectively (Fig. 5). Similar retreat has been observed across the Andes and is likely due to the lack of thick, insulating debris coverage on the glaciers, their low latitudinal setting and/or high ice surface gradients (Benn and others, 2012; Shukla and Qadir, 2016). Shallap's more rapid retreat likely reflects its lower minimum altitude (~4700 m) compared to Zongo (~4900 m), as glacier retreat in the Tropical Andes is strongly altitude-dependant (Veettil and Kamp, 2017). The advance/retreat of Shallap's terminus lagged mass-balance perturbations by 1 year (section 3.1.2). This is quicker than the 3–5 years suggested for highly responsive glaciers in New Zealand's Southern Alps (Purdie and others, 2014) and the wider Peruvian Cordillera Blanca (Vuille and others 2008a). Unlike Shallap, no statistically significant lag time was found between Zongo's mass balance and terminus advance/retreat (section 3.1.2). This is unexpected, as Zongo is smaller, thinner and located in a wetter climate, which should lead to a more rapid response (Christian and others, 2018). We propose this behaviour is due to the sensitivity of Zongo to climate forcing. Its high sensitivity means that terminus advance/retreat never reaches equilibrium with mass balance, as Zongo is continuously responding to climate change over a range of timescales (Gabbud and others, 2016). To further investigate this, higher spatial and temporal resolution mass balance and terminus advance/retreat measurements are needed and processes such as sublimation require monitoring.

4.2 Past mass-balance response to the ENSO

The glacier-wide mass balances of both Shallap and Zongo responded rapidly to the ENSO, with Shallap's mass balance responding 1 year after the ENSO and Zongo's in the same year (Fig. 6). These rapid response times agree with the 1–3-month responses suggested for Ecuadorian glaciers (Veettil, 2012). The slightly shorter response time of Zongo, compared to Shallap, may reflect its smaller size and thickness (Christian and others, 2018). Furthermore, Zongo is located in the wet outer-tropics, which means that up to 70% of its mass balance is determined by wet season precipitation receipts and the wet season onset, which are strongly influenced by ENSO variability (Soruco and others, 2009; Sagredo and Lowell, 2012). Finally, our differing approaches for calculating mass balance may be responsible for the observed difference: Zongo's mass balance was derived from in situ glaciological and climatological measurements, whereas Shallap's mass-balance time series was derived from volume-area scaling. This utilises area changes to determine volume, and then mass balance, but area changes themselves will take time to respond to mass-balance change (Wang and others, 2014). Thus, our geodetic mass balance for Shallap likely represents the maximum response time.

At both Shallap and Zongo, mass-balance response times to the ENSO changed during our study period (Figs 6a, 6b). At Shallap, there was a 1-year response time between 1992 and 1999, which was far quicker than the response time of up to 8 years suggested previously by Veettil and Simões (2019). This may reflect the regional-scale approach of the Veettil and Simões (2019) study compared to our glacier-specific focus: Shallap likely has local micro-climates and specific glaciological characteristics (i.e. altitude, aspect, area) that mediate response times compared to the regional average variability (Akhtar and others, 2008; Immerzeel and others, 2012). Between 2005 and 2010, Shallap's mass balance temporarily decoupled from ENSO forcing, i.e. changes in mass balance preceded changes in ENSO (Fig. 6b). With the available data, it is difficult to determine the cause of this decoupling, but we propose five potential explanations:

1. (i) From 2005 to 2010, Shallap's mass balance was out of phase with the ENSO (i.e. El Niño events cause less negative mass balance) and the PDO and ENSO are also in antiphase. Thus, the PDO may have been dominating local climate, causing the ENSO–mass-balance relationship to breakdown. However, after 2010, the PDO and ENSO remained in antiphase and Shallap's mass balance responded to ENSO with no lag, suggesting that the PDO being out-of-phase with the ENSO was not the only cause of decoupling.

2. (ii) During the period of decoupling, five of the six ENSO peaks occurred in January–April, during the regional wet season. Glaciers in the dry outer-tropics are more sensitive to the ENSO if the peak occurs during the dry season, as their mass balance is more reliant on temperature given that precipitation is comparatively low all year round (Francou and others, 1995). Thus, if ENSO peaks occur in the wet season, the impacts on temperature are lower, meaning that the impact of the ENSO on mass balance is more limited.

3. (iii) The period of decoupling coincides with warmer temperatures in Huaraz (Fig. 6d), which may have overridden the ENSO signal (Vuille and others, 2008a; Juřicová and Fratianni, 2018) and limited its impact on glacier mass balance.

4. (iv) Other climate change/variability that has not been considered in this study, such as the ITCZ, may be influencing local climate.

5. (v) The mass balance of Shallap is responding to the ENSO but with a longer response time, but this is not evident in our data.

To determine which of these mechanisms drives the variability in the ENSO–mass-balance relationship in the future, we recommend the installation of continuous climate monitoring stations, in parallel with mass-balance measurements. Whatever the underlying mechanism, our results support the idea that ENSO impacts mass balance in the Tropical Andes in most, but not all years (Schauwecker and others, 2014), as opposed to a general breakdown in the ENSO–mass-balance relationship, as climate warming overrides the ENSO signal (Rabatel and others, 2013).

Unlike Shallap, Zongo displayed a continuous response to the ENSO from 1991 to 2017, but after 2003, mass-balance response times increased from 0 to 1 year (Fig. 6a). The timing of the ENSO peak had been suggested as a possible mechanism for changing the response time of Zongo's mass balance to the ENSO (Soruco and others, 2009) but we did not detect this. Instead, we suggest that temperature rises of ~1°C in La Paz since 2003 may have reduced the impact of ENSO on Zongo's mass balance and increased response time to 1 year. Thus, our study glaciers display different changes to the ENSO–mass-balance relationship over time: at Shallap, ENSO continues to impact mass balance in most years, but the relationship undergoes temporary breakdowns, whereas at Zongo, climatic warming has gradually overridden the ENSO signal. The differing response is likely due to differing rates of climate warming across the Tropical Andes (Rabatel and others, 2013), the variable responsiveness of glaciers in the wet and dry outer-tropics to climate variability (Vuille and others, 2018), and/or the differing elevations, areas and ice thicknesses of Shallap and Zongo.

4.3 Modelled future ENSO responses

4.3.1 Control run – neutral mass balance (i)

In response to the 2007–2008 (control) mass balance, Shallap shrunk substantially after 50 years, with the terminus retreating ~670 m and ice volume reducing by 46% (Fig. 7a). In contrast, Zongo displayed limited change in response to its 2006 (control) mass balance, with the terminus retreating ~58 m while ice volume increased by 13% (Fig. 8a). We attribute this difference to Shallap's lower minimum altitude (~4700 m) compared to Zongo (~4900 m), which agrees with the strong contemporary altitude-dependence of glacier retreat in the Tropical Andes (Veettil and Kamp, 2017). Furthermore, Shallap is located in the dry outer tropics and has a more negative mass balance (−0.7 m w.e a⁻¹), than Zongo (−0.3 m w.e a⁻¹), meaning that its accumulation is limited and so it is more difficult to replenish losses experienced at its terminus.

4.3.2 Present-day high-magnitude ENSO response (ii)

Shallap and Zongo had differing responses to mass-balance perturbations representative of high-magnitude ENSO events observed during the satellite era (section 3.2.1.2). Shallap underwent a 55% reduction in ice volume during its present-day El Niño simulation and a 37% reduction during the La Niña simulation (Fig. 7b). These volume losses are significantly higher than the 19% volume gain of Zongo during the La Niña simulation and 1% gain during the El Niño (Figs 8b, 8c). We suggest this is due to the greater average altitude and altitudinal range of Zongo, which increases its resistance to climate warming, by reducing sensitivity to rises in the ELA (McGrath and others, 2017). At Shallap, climate warming has significantly disrupted the mass-balance relationship with the ENSO. Temperature increases have lowered the ELA so far that even the ELA rises experienced during La Niña can no longer contribute mass gains.

4.3.3 Response to high-magnitude El Niño events

To analyse the sensitivity to Shallap and Zongo to different future ENSO conditions, two high-magnitude El Niño events were simulated: Niño-1 and Niño-2. Using exact climate projections of El Niño events would be ideal, but debate exists over the evolution of the ENSO under climate warming (Luo and others, 2018) and poor meteorological data for the Tropical Andes hinder climate model forecasting (Larsen and others, 2011). Therefore, we conducted a series of sensitivity experiments, based on potential future scenarios.

In response to Niño-1, Shallap lost an additional 4% of its volume and retreated by an extra 10% compared to experiment (ii) (Fig. 7d). This suggests that Shallap is relatively insensitive to major changes in accumulation and ablation. In contrast, Zongo's terminus retreated twice as much as in experiment (ii) and its volume decreased by 26%, compared to the 1% gain for experiment (ii) (Fig. 8d). This suggests that Zongo is far more sensitive to changes in mass balance. We attribute the different responses of Shallap and Zongo to Niño-1 to the impact of ENSO forcing on the ELA, along with glacier hypsometry and ice thickness. At Shallap, there is an area of comparatively thick ice located between 5000 and 5200 m, which is ~50 m above its equilibrium state ELA (Gurgiser and others, 2013). Thus, if the ELA remains below this altitude, Shallap retains a substantial ice reservoir and is comparatively insensitive to Niño-1. This is supported by results from Niño-2, when the ELA was raised to 5200 m, and 76% of the glacier's volume was lost at the end of the 50-year run period (Fig. 7e). Thus, the impact of ENSO forcing on the ELA appears critical to Shallap's response and demonstrates more broadly how relatively small ELA changes can instigate significant glacial loss, given the right glacier geometry and hypsometry (McGrath and others, 2020) (Fig. 12).

Fig. 12. Ice thicknesses distribution for Shallap (a) and Zongo (b). Shown are contour elevations of 200 m, the location of the ELA and the location of each glaciers retrospective ice-bulge (an area of localised ice thickness increase).

At Zongo, we observe a similar sensitivity to ENSO-induced changes in ELA in the Niño-2 simulation (Fig. 8e): during Niño-2, Zongo lost 48% of its volume, which was double the Niño-1 scenario (Fig. 8d). As at Shallap, Zongo has its thickest ice at 5100–5200 m, suggesting it is sensitive to ENSO forcing that causes the ELA to rise into this area. Zongo's hypsometry may also explain its comparative insensitivity to present-day high-magnitude ENSO events: the ELA did not notably increase in experiment (ii), meaning that most of the accumulation zone was not impacted by warmer temperatures under El Niño conditions. Furthermore, it experiences high accumulation rates, due to its location in the wet outer tropics, meaning it can more easily offset the enhanced ablation during El Niño events. Overall, our results demonstrate that glacier hypsometry and ice thickness are crucial in determining the response of tropical glaciers to ENSO forcing (Florentine and others, 2020). However, the ice thickness distribution of most Andean glaciers is poorly constrained (Andreassen and others, 2015), which limits our capacity to estimate their response to both ENSO forcing and climate warming and highlights the urgent need for reliable and comprehensive ice thickness measurements.

4.3.4 Varied maximum altitude of ENSO forcing (iv)

The maximum altitude of ENSO forcing influenced ice loss from Shallap and Zongo significantly during the 50-year run (section 3.1.1.4), which we attribute to both glaciers having ‘top-heavy’ hypsometries (i.e. having comparatively large accumulation areas vs ablation areas; Shukla and others, 2020). At Zongo, 5900 m was a critical altitude: mass loss increased substantially once ENSO forcing reached 5900 m (Fig. 9b). Similarly, at Shallap, applying ENSO forcing to 5300–5500 m was a key determinant of the glacier's response to ENSO. Shallap has a comparatively low ELA, meaning it has a relatively small ablation zone and large accumulation zone, but receives limited precipitation inputs due to its dry climate (Gurgiser and others, 2013). Most of Shallap's accumulation zone is between 5300 and 5500 m (Fig. 9c) and so this area is particularly sensitive to ENSO forcing. Thus, if El Niño conditions are applied to this zone (i.e. drier and warmer conditions and hence more negative mass balance) then the glacier undergoes rapid ice loss. As such our results from both glaciers suggest that the impact of ENSO forcing in the accumulation zone is greater than in the ablation zone, due to hypsometry and ice thickness distributions. This has important implications for their potential response to future warming. As climate in the outer-tropical Andes warms, we expect that ELAs will gradually move up-glacier (Carturan and others, 2020) and so ENSO forcing will reach higher altitudes. During El Niño events, ablation will occur at higher altitudes on Zongo and Shallap and encompass areas of thicker ice, which our results show is crucial for determining their response to ENSO. Over time, these reserves may be depleted by repeated El Niño events, meaning that both glaciers become more sensitive to both ENSO forcing and general climate warming (Fig. 13), which may hasten their eventual disappearance.

Fig. 13. Conceptual model of how climate warming and changes to the ENSO will impact the ELA of Shallap and Zongo, leading to rapid deglaciation. Grey lines are the steady-state ELA, red lines El Nino event ELAs and blue lines La Nina event ELAs.

4.3.5 Stronger and longer El Niño

Both Shallap and Zongo were relatively insensitive to longer/more frequent El Niño events in the future but were sensitive to stronger El Niño events (Figs 7, 8). This was the case for both neutral mass balances between El Niño events and for a simulated mass balance under climate change (Figs 7, 8). Thus, we suggest that an increased frequency of the El Niño phase under climate change, as has been forecast by many climate models (Timmermann and others, 1999; Williams and Patricola, 2018; Yang and others, 2021), will have limited influence on glacier evolution in the future if the strength of El Niño events does not increase. While the glacier–climate interactions that cause this are difficult to establish, we suggest this may be due to: (i) the climate warming signal overriding any volume changes associated with longer El Niño events; and/or (ii) the fact that we only analysed El Niño events occurring twice as regularly as La Niña or lasting twice as long, and even more frequent or longer perturbations may be needed to instigate a measurable glacier response.

Under stronger El Niño events, it was possible for Zongo to completely disappear after 50 years and for most of Shallap to be lost (Figs 10, 11), meaning that potential increases in the strength of El Niño events have significant consequences for glacier longevity. Furthermore, we use only extrapolated historic climate warming between El Niño events, meaning that deglaciation may be even more rapid if rates of climate change are higher in the future. Our results correspond well with those of Réveillet and others (2015) who used modelled climate forecasts to determine Zongo's future ice loss: both studies produced a reduction in the volume of ~50% by 50 years into the study.

Zongo showed greater sensitivity to increased El Niño strengths than Shallap, while Shallap appeared more sensitive to longer-term climate change. We suggest that this is because Shallap has a lower altitude (~4700 m compared to Zongo's ~4900 m) and is in a drier climate (dry outer tropics vs Zongo in the wet outer tropics climate zone) (Vuille and others, 2018), making it more difficult for it to compensate for increased ablation under a warming climate. At Zongo, the altitudinal range is greater, meaning that it can maintain an accumulation area, even under temperature rises of 3–4°C (Réveillet and others, 2015). However, the additional rise in the ELA under stronger El Niño events causes ablation in areas usually in the accumulation area. This causes ice loss from critical altitudes, which contain thick ice and/or a large portion of the accumulation zone by area. Thus, repeated strong El Niño events could hasten glacier loss, due to their impact on these critical ice reservoirs. Thus, we suggest that changes in the intensity of the ENSO is critical for predicting the future evolution of Zongo and Shallap, but the interaction between climate change and ENSO intensity is currently poorly constrained by climate models (Guilyardi, 2006; Da Rocha and others, 2014). This limits our capacity to forecast future ice losses from Shallap and Zongo, and from tropical Andean glaciers, where similar controls are likely to operate. Furthermore, our work highlights glacier hypsometry and thickness as key factors in the response of glaciers to both ENSO and climate warming, making it crucial to take into account both local climate and local glacier properties in forecasts of near future ice loss from the tropical Andes. Additionally, it should be noted that there are significant uncertainties in modelling the future evolution of ENSO and its characteristics at the high altitudes at which our study glaciers are found, which is further compounded by the lack of meteorological station data in these areas. Furthermore, it is uncertain how ENSO forcing will translate to SMB, highlighting the need for further modelling work on the impact of changes in ENSO on climate in the Andes.

5. Conclusion

Overall, our results demonstrate that Shallap and Zongo glaciers are highly sensitive to both anthropogenic climate change and ENSO forcing. Over the satellite era, both Shallap and Zongo experienced substantial retreat and increasingly negative cumulative mass balances. During this period, both glaciers displayed strong, rapid (0–2 years) responses to ENSO variability and did not just respond to high-magnitude ENSO events. However, the ENSO and mass-balance relationship showed temporal and spatial variability: at Zongo, we observed reduced sensitivity to the ENSO during periods of climate warming, while at Shallap, a temporary breakdown in the ENSO and mass-balance relationship occurred between 2005 and 2010. Our numerical modelling experiments indicate that Shallap and Zongo are both highly sensitive to forcing from anthropogenic climate change and ENSO and that the future evolution of ENSO in the region may be a critical control on ice loss from these glaciers. Furthermore, our results suggest that initial thickness distribution, hypsometry and altitudinal range exert a first-order control on response of the study glaciers to ENSO and climate warming. Specifically, increases in the ELA due to climate warming and/or El Niño can cause rapid ice loss, as they initiate melt in the large, high-altitude accumulation areas of our study glaciers. There is debate that the future impact of ENSO may diminish at the highest altitudes (i.e. above 5000 m) as climate warming begins to override the ENSO signal. If this occurs, then ice loss related to El Niño events may reduce, as these weakened events would have less impact on the ELA and hence on melt rates at high elevations. However, if the ENSO continues to impact temperature and precipitation regimes at high altitudes, then Zongo and Shallap will continue to remain vulnerable to enhanced and more extensive melt during El Niño events. Furthermore, our data indicate that under the most extreme scenarios of longer and stronger ENSO cycles, with climate change, the study glaciers could disappear entirely within 30 years. Taken together, our results suggest that tropical Andean glaciers are highly sensitive to both anthropogenic forcing and ENSO cycles and that their response to forcing is strongly modulated by glacier hypsometry. Thus, there is an urgent need for more extensive glaciological data in the region, particularly on ice thickness and mass balance, which is currently very sparse. These data are essential for accurately forecasting ice loss from tropical Andean glaciers during the 21st century and, in turn, the impacts of this ice loss on the downstream populations and ecosystems that rely on Andean glaciers for their water supplies.