New insights into the decadal variability in glacier volume of a tropical ice cap, Antisana (0°29′ S, 78°09′ W), explained by the morpho-topographic and climatic context

We present a comprehensive study of the evolution of the glaciers on the Antisana ice cap (tropical Andes) over the period 1956–2016. Based on geodetic observations of aerial photographs and high-resolution satellite images, we explore the effects of morpho-topographic and climate variables on glacier volumes. Contrasting behaviour was observed over the whole period, with two periods of strong mass loss, 1956–1964 (−0.72 m w.e. yr⁻¹) and 1979–1997 (−0.82 m w.e. yr⁻¹), and two periods with slight mass loss, 1965–1978 (0.10 m w.e. yr⁻¹) and 1998–2016 (−0.26 m w.e. yr⁻¹). There was a 42 % reduction in the total surface area of the ice cap. Individually, glacier responses were modulated by morpho-topographic variables (e.g. maximum and median altitude and surface area), particularly in the case of the small tongues located at low elevations (Glacier 1, 5 and 16) which have been undergoing accelerated disintegration since the 1990s and will likely disappear in the coming years. Moreover, thanks to the availability of aerial data, a surging event was detected on the Antisana Glacier 8 (G8) in the 2009–2011 period; such an event is extremely rare in this region and deserves a dedicated study. Despite the effect of the complex topography, glaciers have reacted in agreement with changes in climate forcing, with a stepwise transition towards warmer and alternating wet–dry conditions since the mid-1970s. Long-term decadal variability is consistent with the warm–cold conditions observed in the Pacific Ocean represented by the Southern Oscillation index.

1 Introduction

The tropical Andes have experienced continuous atmospheric warming over the last 4 decades (Chimborazo et al., 2022). Although the pattern is not identical from one site to another (Aguilar-Lome et al., 2019; Pepin et al., 2015; Yarleque et al., 2018), it confirms the marked glacier shrinkage observed in this part of the Andes (Dussaillant et al., 2019; Masiokas et al., 2020). Furthermore, climate projections warn that warming will continue during the 21st century (Bradley et al., 2006; Urrutia and Vuille, 2009), thereby increasing the rate of glacier melt and leading to major shrinkage and even to the complete disappearance of many glaciers in the region in the coming decades (Rabatel et al., 2018; Yarleque et al., 2018).

Despite the atmospheric warming, the small Antisana 15α Glacier (0.3 km² in 2016) in Ecuador (Fig. 1), which has been surveyed since the mid-1990s, showed only a slight mass loss during the 1995–2012 period (−0.12 m w.e. yr⁻¹) (Basantes-Serrano et al., 2016). Antisana 15α is part of an ice cap, and thanks to its very high elevation (up to 5700 m a.s.l.), the glacier still has a large accumulation area (roughly 60 % of its total surface area), which is not the case for the vast majority of the small glaciers in this region (Rabatel et al., 2013a). These peculiar morpho-topographic characteristics associated with the almost year-long humid conditions in the inner tropics (which contribute a considerable amount of solid precipitation in the upper reaches of the glacier) are the main drivers of the limited mass loss of the Antisana 15α Glacier (Basantes-Serrano et al., 2016). However, despite being the longest glaciological time series in the inner tropics, the representativeness of the changes observed on the Antisana 15α Glacier needs to be assessed at the scale of the other small glaciers (e.g. Nicholson et al., 2013). To this end, high-resolution stereoscopic imagery is useful to evaluate the representativeness of the glacier mass changes at the scale of a mountain range (e.g. Rabatel et al., 2006; Soruco et al., 2009b).

In addition, field measurements taken on Glacier 15 (G15), Antisana, showed that the key periods of the year that explain more than 90 % of the variance of the mass balance are centred on the equinoxes, i.e. between March and May (MAM) and between September and November (SON) (Basantes-Serrano et al., 2016; Francou et al., 2004). Indeed, these quarters are those with maximum potential solar radiation and higher rates of precipitation caused by the transit of the Intertropical Convergence Zone (ITCZ) (Vuille et al., 2000). Previous studies also confirmed a strong relationship between the variability in the mass balance and the occurrence and intensity of El Niño–Southern Oscillation (ENSO) events at an annual time step (Favier et al., 2004; Francou et al., 2004; Rabatel et al., 2013a). Consequently, we expected a cumulative influence of seasonal climate conditions in the MAM and SON quarters on the long-term variability in the decadal mass balance.

The present study is the first comprehensive assessment of the multidecadal variability in the glacier mass fluctuations quantified on 17 glaciers covering an equatorial ice cap using the geodetic method. We explore the role of topographic features in determining the responses of individual glaciers and the impact that the change in surface area caused by the geometric adjustment of small glaciers may have on the mass balance of the ice cap as a whole. We then discuss the changes in ice masses resulting from the climatic conditions that have prevailed in the region since the mid-1950s. The paper is structured as follows: Sect. 2 describes the study site and the key climate features of the inner tropics. Section 3 provides a concise summary of remote sensing and climate data and data processing. Section 4 provides a detailed description of the changes observed in the glaciers in five sub-periods since the mid-1950s and the topographic–climatic drivers that may explain the glaciers' responses together with the potential influence of the large-scale ENSO episodes on the decadal variability in the tropical ice cap.

Our study provides a baseline to model the future effects of glacier retreat on the water cycle and mountain ecosystems in the inner tropics.

Figure 1 (a) The Antisana ice cap with its 17 glaciers. The blue polygon locates Glacier 15 where in situ observations have been available since 1995. The grey line represents the glacier boundaries in 1956, and the black lines represent the boundaries in 2016. Data are missing on the southeast side of the ice cap due to cloud cover (area shown in white). A false-colour (RGB 324) pan-sharpened composite of Pléiades 0.5 m ortho-images used as the background shows the distribution of the vegetation in the surrounding glaciers (© CNES 2016, Distribution Airbus D&S). The bar charts in (b) and (c) show precipitation seasonality throughout the year in the 1950–2018 period measured at two weather stations, M003 (M188) located on the western (northern) side of the volcano (see Fig. 2). Panels (d) and (e) are the seasonal mean precipitation maps for December–February (DJF) and June–August (JJA) derived from ERA5 reanalysis data for the 1950–2018 period (shading scale). The solid black line shows the average position of the ITCZ estimated from the ERA5 precipitation rates and the South American Monsoon System (SAMS, blue arrows) during the austral summer (DJF). The red star locates Antisana volcano.

2 Study site and climate settings

Antisana is an active stratovolcano capped by a ∼ 13 km² ice cap in 2016, located close to the Equator (0^∘29^′ S; ∼5700 m a.s.l.), around 40 km to the east of the city of Quito. Little information is available concerning the volcano's activity in the last 400 years; no recent morphological changes in the Earth's crust, thermal activity or local ice reductions due to hot streams have been identified that would suggest recent volcanic activity (Hall et al., 2017). The ice cap has an overall radial structure composed of 17 glacierised catchments identified through photogrammetric restitutions following the boundary of the water basins from the summit of the volcano to the valleys. Some of the glaciers have a local name, but most are named by numbers starting from the north and going clockwise. Hereafter, we use the numerical nomenclature for easy identification of the glaciers; for example, Antisana Glacier 15 is denoted G15 (Fig. 1a).

The climate features in the tropical Andes are mainly driven by the oscillatory position of the ITCZ and the presence of the South American Monsoon System (SAMS), illustrated in Fig. 1d and e, and the deep atmospheric convection originating from the very wet Amazon rainforest which creates precipitation regimes in this region which differ even at short distances (Fig. 1b and c). The Andes Cordillera controls the spatial patterns of the wet flows driven by the trade winds blowing from the Amazon basin (Espinoza et al., 2020; Garreaud, 2009). In fact, the site is characterised by a mixed seasonal rainfall cycle. On the western side of the volcano, a bimodal cycle is observed with a first maximum from March to May and a second maximum in November (Fig. 1b), whereas the eastern side shows unimodal seasonality with a maximum lasting from July to August (Fig. 1c). The climate patterns may occur simultaneously when the Amazon seasonality penetrates the Andes Cordillera rather deeply through the wide valleys due to the atmospheric circulation (Tobar and Wyseure, 2018).

Local precipitation occurs all year round and ranges from ∼1300 mm yr⁻¹ on the western slope to ∼3000 mm yr⁻¹ on the eastern slope of the volcano (Fig. 1b and c), whereas air temperatures display more regional behaviour and remain almost constant throughout the year. Furthermore, in situ meteorological observations recorded at elevations ranging from 1900 to 5000 m a.s.l. suggest a positive vertical precipitation gradient, with annual precipitation rates that can reach a maximum of ∼6000 mm yr⁻¹ (Emck, 2007; Garreaud, 2009) and thus have a significant effect on the accumulation rates over the course of the year (Basantes-Serrano et al., 2016). An increase (decrease) in air temperature causes rainfall (snowfall), thereby enhancing ablation (accumulation) processes at the surface of the glacier (Favier et al., 2004; Francou et al., 2004). In addition, the response of the glacier to climate variation is modified during ENSO events, with very positive (negative) mass balances during La Niña (El Niño) events (Francou et al., 2004).

Morpho-topographic and climate interactions play an important role not only in the behaviour of glaciers but also in the definition of the biogeography of the site (Fig. 1). Below 4750 m a.s.l., glacier retreats have left behind several morainic deposits on the western slope of the volcano, which become new ecosystems for certain plants and other specialised organisms capable of developing in particularly challenging environmental conditions (Cauvy-Fraunié and Dangles, 2019; Cuesta et al., 2019). Below 4650 m a.s.l., small páramo vegetation dominates the landscape. On the eastern side, glaciers extend down to 4400 m a.s.l., and the areas abandoned by the glacier are immediately occupied by a dense rainforest and shrub vegetation.

3 Data and methods

3.1 Processing of aerial and satellite data

Photogrammetry is a powerful tool to assess the glacier-wide mass balance of several glaciers at the scale of a mountain range. In the present study, we took advantage of the availability of high-resolution aerial and satellite stereoscopic data from different sources acquired on different dates. The stereoscopic data enabled us to map every glacierised basin on the Antisana ice cap at high spatial resolution (Table 1) and to quantify changes in the mass balance of the glaciers at a decadal time step since 1956.

Table 1 High-resolution remote sensing data used in this study.

* Geodetic reference used to evaluate the xyz consistency of the elevation dataset.

On the one hand, we used 90 raw aerial photographs for the Antisana volcano taken on five different dates in 1956, 1965, 1979, 1997 and 2009. The aerial surveys were conducted by the Instituto Geográfico Militar (IGM) of Ecuador. The aerial films were scanned at 14 µm resolution using an Intergraph PhotoScan TD system. All the calibration reports of each sensor were available; this information is essential to reconstruct the geometry of the sensor at the moment of the aerial acquisition. The bundle block adjustment was conducted based on a geodetic network of 21 ground control points (GCPs) collected using the differential global positioning system (DGPS) method in a static mode with an uncertainty of less than ∼0.20 m (see Fig. 1b in Basantes-Serrano et al., 2016). For each aerial survey, very dense 3D point clouds representing between 90 000 and 185 000 points were generated automatically on the non-glacierised areas by applying an image correlation technique which matches similar features in the overlap area between two images (Hirschmuller, 2005). Over the glaciers, where the snow cover, shadows and cloud cover prevented effective application of the correlation methods, 42 500 3D points per aerial survey were measured randomly by photogrammetric restitution, giving a total of 212 500 points and avoiding areas with low contrast and saturated areas. These 3D point clouds were subsequently used to compute differences in elevation in each study period. All the photogrammetric tasks were performed with IMAGINE Photogrammetry (formerly the Leica Photogrammetry Suite (LPS) (Hexagon Geospatial^®)).

On the other hand, a pair of Pléiades images of the same area taken from different perspectives enabled the geodetic mass balance survey to be extended to 2016. These images were acquired by the sun-synchronous Pléiades-1A platform with a viewing angle of 12.69 and 8.76^∘ in the across-track direction (ω) and of −8.66 and 10.16^∘ in the along-track direction (ϕ) with respect to the nadir. To check the ability of the stereoscopic pair to capture the mountain topography, we computed the ratio of the distance of two successive positions of the satellite (B) to its altitude above ground H, resulting in 0.33, which is optimum for processing numerical models in mountainous areas (Perko et al., 2019). The relative orientation of the Pléiades images was performed using the rational polynomial coefficients (RPCs) provided by the ancillary information. Finally, stereo-matching algorithms enabled the generation of a 3D point cloud by means of pixel-to-pixel correspondences between the images. The 3D mapping procedure for the stereoscopic pair was performed using the Ames Stereo Pipeline (ASP), free open-source software developed by the National Aeronautics and Space Administration (NASA; Beyer et al., 2018). The Pléiades images were acquired in the framework of the ISIS programme of the French Space Agency (CNES) in the framework of the Pléiades Glacier Observatory initiative (Berthier et al., 2014).

3.2 Geometric adjustment of the geodetic data

To assess spatial consistencies in the geodetic information, surface elevations in stable rocky areas in the forelands of the ice cap were compared following the iterative workflow proposed by Nuth and Kääb (2011). Co-registration is recommended before identifying changes in glacier elevation to ensure precise xyz alignment between the geodetic data. A conical structure like the Antisana volcano is perfect to correct elevation bias because, as long as there are sufficient examples of differences in elevation (hereafter dh samples) over stable rocky areas, it covers all terrain aspects. The elevation differences (dh) of the two misaligned digital elevation models (DEMs) are a function of the slope of the terrain (α) and the aspect (ω) of the dh samples as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}h}{\tan \left(\mathit{\alpha }\right)}}=a\cdot \cos \left(b-\mathit{\omega }\right)+c,\end{array}$$

where a is the magnitude of the misalignment, b is the direction of the misalignment and c is the mean elevation bias. Then, the Δxyz adjustment factors are computed as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle }{\mathrm{\Delta }}_x=a\cdot \sin \left(b\right),\text{(3)}& {\displaystyle }{\mathrm{\Delta }}_y=a\cdot \cos \left(b\right),\text{(4)}& {\displaystyle }{\mathrm{\Delta }}_z={\displaystyle \frac{c}{\tan \left(\stackrel{\mathrm{‾}}{\mathit{\alpha }}\right)}}.\end{array}$$

As the 2009 geodetic survey presents the best photogrammetric adjustment (Table 1), the DEM generated from these data was used as a reference to evaluate the bias in the geodetic data of the other acquisitions. Because high-spatial-resolution data were used in this study (Table 1), an elevation difference layer of 10×10 m cell size was computed for each geodetic comparison using the universal kriging technique based on an exponential semi-variogram model. After co-registration, analysis revealed that the geodetic datasets were slightly unaligned with respect to the reference, with a mean elevation difference ranging from −1.8 to 3.2 m on stable terrain. Consequently, we decided to perform 3D geometric adjustment to minimise the difference in the geodetic data. The results are presented in Table 2 and in the Supplement. In all cases, there was a reduction of around 8 % in the standard deviation and the mean difference in elevation was improved, and a normal distribution of the dh on the glacier forelands is better shaped with an average standard error of 0.03 m after 3D geometric adjustment (Fig. S1 in the Supplement).

Table 2 Differences in average elevation and its standard deviation over glacier-free terrain before (${\stackrel{\mathrm{‾}}{\mathrm{d}h}}^{\prime }$, σ^′) and after (${\stackrel{\mathrm{‾}}{\mathrm{d}h}}^{\prime \prime }$, σ^′′) the 3D geometric adjustment, computed for the total number of elevation differences (dh samples) off-glacier for the five study periods. The geometric solution on the three axes (Δ_x, Δ_y, Δ_z) and the average standard error after adjustment (SE) are given. All the statistics are expressed in metres (m).

The co-registration results showed a similar displacement vector of around 8 m in all five datasets. For the aerial-based geodetic data, this was attributed to the number and distribution of the GCPs used to solve the stereo orientation. For the Pléiades-based data, this is probably explained by the fact that the triangulation and geometric adjustment of Pléiades stereo-images were based on the RPC model without including GCPs. All the elevation data can thus float up to 10 m above or below the ground and have to be adjusted horizontally and vertically (Berthier et al., 2014).

3.3 Computing the geodetic mass balance

The geodetic method captures multidecadal changes in glacier mass (Cogley et al., 2011) and also makes it possible to document the spatial variability in changes in elevation at the surface of the glaciers. On the basis of this variability, we spatially optimised the 3D measurements at the surface of glaciers to compute changes in elevation using the geo-statistical framework proposed by Basantes-Serrano et al. (2018), which relies upon the spatial variability in the elevation change to densify a sampling network to optimise the quantification of the surface elevation change. The number of dh samples per glacier depends on their size; a minimum of 1000 sampling points per square kilometre was thus measured to estimate the glacier mass balance. With this approach, instead of trying to cover the entire glacier surface with 3D coordinates, we chose sites where significant elevation changes are expected, avoiding places with steeper slopes where large elevation outliers are probable (Berthier et al., 2014).

No measurements were made on about 10 % of the glacier surface comprising areas with low contrast, including saturated areas, or with a slope greater than 45^∘. dh outliers were removed considering a threshold of three normalised median absolute deviations (NMADs), as done in other studies (e.g. Braun et al., 2019; Brun et al., 2017). The procedure for removing outliers was applied to each individual glacier, and about 5 % of the differences in elevation outside the 3NMAD threshold were excluded. Based on the remaining dh samples, an exponential variogram model was fitted to generate an elevation difference grid of 10×10 m cell size using the universal kriging technique (Basantes-Serrano et al., 2018). Because all gridded surfaces have the same spatial resolution, we did not expect any curvature error in the topography of glaciers that could affect the estimation of the geodetic mass balance (Gardelle et al., 2012). Surface areas for each geodetic survey were manually digitalised in stereo mode following the boundaries of the glacierised catchments.

Finally, the glacier-wide mass balance expressed in metres water equivalent (m w.e.) for each glacier was computed by

$$\begin{array}{}\text{(5)}& B_{\mathrm{g}}=(\stackrel{\mathrm{‾}}{\mathit{\rho }}\cdot r^{\mathrm{2}}\cdot {\sum }_{i=\mathrm{1}}^p\mathrm{\Delta }h\left(x_i\right))/\stackrel{\mathrm{‾}}{S},\end{array}$$

where $\stackrel{\mathrm{‾}}{\mathit{\rho }}$ is the average density value, r is the pixel size, Δh is the change in glacier surface elevation at each location x_i, p is the number of pixels covering the glacier at its maximum extent, and $\stackrel{\mathrm{‾}}{S}$ is the glacier surface area averaged over the period between the date of the first aerial survey and date of the second aerial survey.

The mean annual rate of glacier-wide mass balance ${\dot{B}}_{\mathrm{g}.t}$ considering the number of years N for each study period is then calculated as

$$\begin{array}{}\text{(6)}& {\dot{B}}_{\mathrm{g}.t}={\displaystyle \frac{B_{\mathrm{g}}}{N}}.\end{array}$$

3.4 Uncertainty analysis

To determine the total uncertainty in the geodetic mass balance, we assessed the contribution of the following individual random and systematic sources:

• First, the uncertainty related to the DEMs used to compute the difference in elevations following the method proposed by Rolstad et al. (2009). Thus, we estimated the spatial auto-correlation of the grid of dh based on a semi-variogram analysis. For the Antisana ice cap, we found that the difference in elevation in non-glacierised areas is spatially uncorrelated beyond a horizonal distance (a₁) of about 1250 m in all geodetic surveys; this distance was computed from the variogram model. Because the averaging of the glacier areas in the first survey S_T is greater than the effective correlated area $S_{\mathrm{cor}}=\mathit{\pi }\cdot a_{\mathrm{1}}^{\mathrm{2}}$ in all five periods, the uncertainty in spatially averaged differences in elevation in metres is defined by

  $$\begin{array}{}\text{(7)}& {\mathit{\sigma }}_{\mathrm{d}h.\mathrm{g}}=\pm {\mathit{\sigma }}_{\mathrm{d}h}\cdot \sqrt{{\displaystyle \frac{S_{\mathrm{cor}}}{\mathrm{5}\cdot S_T}}},\end{array}$$

  where σ_dh is the standard deviation of dh over stable terrain in the vicinity of the ice cap.

• Second, regarding the uncertainty related to the density assumption, we analyse two extreme scenarios. First, we consider an average density recommended by Huss (2013) of $\stackrel{\mathrm{‾}}{\mathit{\rho }}=\mathrm{850}$ kg m⁻³ with a plausible uncertainty range of ${\mathit{\sigma }}_{\mathit{\rho }}=\pm \mathrm{60}$ kg m⁻³. This value is appropriate for a wide range of conditions and when no information on firn pack changes is available (Huss, 2013; Zemp et al., 2013). However, moderate mass gains occurred in the second study period for which the conventional density assumption may not be true. Taking advantage of firn compaction data in two shallow cores (mean depth ∼14 m) extracted from the summit of the Antisana volcano in February 1996 and November 1999, respectively (Calero et al., 2022; Williams et al., 2002), we propose a second scenario with an average density value of $\stackrel{\mathrm{‾}}{\mathit{\rho }}=\mathrm{564}\pm \mathrm{64}$ kg m⁻³, indicating that the mass gain or loss mostly comprised firn.

• Third, missing glacier elevation differences, i.e. where no elevation measurements were available, were attributed to (i) saturated areas because of the coarse grain size (low sensitivity) of the film used during the aerial surveys or the presence of bright snow patches, most of which were located in the lower flat areas; (ii) cloud cover which affects some parts of the glaciers; and (iii) sites with slopes greater than 45^∘. These three characteristics are mostly found in the accumulation zone of the Antisana ice cap. Recently McNabb et al. (2019) explored the effect of DEM voids on the geodetic mass balance and evaluated several interpolation methods to fill gaps in elevation data. These authors concluded that interpolation methods work quite well when data gaps are limited to small patches or when elevation changes are regularly distributed over the glaciers, as is the case in ablation zones. However, this is not necessarily true for the accumulation zones of Antisana glaciers, which present a noisy spatial pattern of either negative or positive surface elevation changes, linked to the interaction between ice dynamics and the harsh topography of the volcanic cone (Basantes-Serrano et al., 2016). In addition, it is worth mentioning that the dh coverage for all periods are evenly distributed over the glacier surface, which reduces the likelihood of inducing some spatial biases in the quantification of glacier elevation changes (Fig. S2 in Supplement). Thus, rather than filling dh voids using an interpolation method, which can introduce biases into the geodetic mass balance, we decided to estimate the uncertainty expressed in metres of water equivalent (m w.e.) by the following expression:

  $$\begin{array}{}\text{(8)}& {\mathit{\sigma }}_{S\mathrm{void}}={\mathit{\sigma }}_{\mathrm{d}h.\mathrm{g}}\cdot \left(\stackrel{\mathrm{‾}}{\mathit{\rho }}\right)\cdot \left({\displaystyle \frac{S_{\mathrm{g}.\mathrm{void}}}{S_{\mathrm{g}}}}\right),\end{array}$$

  where S_g.void is the portion of individual glacier with missing dh values and S_g is the total surface area of the corresponding glacier recorded during the first aerial survey.

• Fourth, uncertainty in surface area determination of glaciers σ_Sg considers a buffer zone of 1 pixel surrounding the glacier boundary and is computed by following the same approach as used to determine the uncertainty when no elevation measurements are available (see Eq. 8) but replacing S_g.void for S_g.b, which is the surface area of the buffer zone.

• Finally, to match the arbitrary timescale imposed by geodetic surveys and the hydrological timescale, it is necessary to perform temporal homogenisation that accounts for the ablation and accumulation processes which occur during the time lag (Thibert et al., 2008). Indeed, the dates of the aerial surveys used in this study (Table 1) do not match the hydrological year in the inner tropics (i.e. from January to December). However, no meteorological and glaciological observations at an annual time step are available in this region before the mid-1990s (Basantes-Serrano et al., 2016); therefore no time adjustment is possible, and we kept the original dates for the mass balance estimations. This is called the floating-date reference. To evaluate the systematic error due to the survey difference σ_t.ref, we assume constant monthly mass balance rates at the glacier surface based on the geodetic mass balance. Then, the monthly mass balance is multiplied by the number of months to match the hydrological year. The uncertainty due to the time reference is evaluated as the residual between the annual mass balance at the floating date and the annual mass balance at the fixed date.

Thus, the overall uncertainty in the geodetic mass balance of each individual glacier expressed in metres of water equivalent (m w.e.) was estimated using the propagation error approach of uncorrelated variables as follows:

$$\begin{array}{}\text{(9)}& {\mathit{\sigma }}_{B_{\mathrm{g}}}=\sqrt{{\mathit{\sigma }}_{\mathrm{d}h.\mathrm{g}}^{\mathrm{2}}\cdot \left(\stackrel{\mathrm{‾}}{\mathit{\rho }}\right)+{\mathit{\sigma }}_{\mathit{\rho }}^{\mathrm{2}}+{\mathit{\sigma }}_{S\mathrm{g}}^{\mathrm{2}}+{\mathit{\sigma }}_{S\mathrm{void}}^{\mathrm{2}}+{\mathit{\sigma }}_{t.\mathrm{ref}}^{\mathrm{2}}}.\end{array}$$

3.5 Climate variability and local drivers of glacier response

Meteorological observations are relatively recent in Ecuador and are very scarce above 4800 m a.s.l. Climate reanalysis provides an opportunity to obtain key information to reconstruct climate conditions in a specific region and to help understand the impacts of climate changes in the past when direct observations are scarce (Hersbach et al., 2020). To analyse the glacier response to regional climate variability, we considered monthly interpolated air-temperature and precipitation data at different pressure levels from the fifth generation of ERA5 reanalysis provided by the European Centre for Medium-range Weather Forecasts (ECMWF), which has a 0.25^∘ (∼30 km) grid resolution. ERA5 reanalysis replaces the former version of ERA-Interim reanalysis. As the 2009 geodetic it well represents the long-term climatology of this region (Fig. S3 in the Supplement); moreover it covers several geopotential levels and the entire study period.

Reanalysis data are prone to systematic biases due to the limited spatial resolution and simplified physics used to represent climate, particularly in mountain areas where orographic factors drive climate conditions. In practice, temperature (precipitation) measurements made near glaciers are somewhat higher (lower) than reanalysis data. To adjust the ERA5 dataset to observations, we compared the reanalysis data close to the grid cell to the station location and the in situ meteorological data at a monthly time step for the nine weather stations with the longest series that recorded climate variability from the western and the eastern sides of the cordillera between 1950 and 2018 (Fig. 2 and Table 3). This observational network is located below ∼4500 m a.s.l. and is maintained, homogenised and quality controlled by the Instituto Nacional de Meteorología e Hidrología (INAMHI).

This comparison showed that ERA5 is able to capture the western climate signal observed at M003 weather station (Izobamba), located 50 km from the volcano. The average annual rainfall in M003 is 1446 mm yr⁻¹ (±15 %), not far from the annual average recorded by one of the rain gauges installed on the western slope of Antisana, close to the terminus of G15, i.e. an average of 1235 mm yr⁻¹ between 1995 and 2018. On the eastern slope, ERA5 reanalysis was unable to capture the seasonality of precipitation of any of the weather stations located in the Amazon basin. However, it will be recalled that the estimated precipitation on the Amazonian side can reach 3000 mm yr⁻¹. On the other hand, the temperature recorded by the M188 weather station (Papallacta) was satisfactorily captured by ERA5 reanalysis.

Figure 2 Map of long-term weather stations located in the region. The city of Quito and the ice cap are represented by the hatched and solid grey polygons, respectively. The inset table lists the meteorological stations and the variables used in this study; the percentage of missing values at a monthly time step are reported. Cartographic information was provided by the Instituto Geográfico Militar of Ecuador (IGM).

To reduce the bias between the ERA5 data and reference observations, we retain M003 data and we applied the following approaches: (i) for temperature calibration we applied a bias-correction method which uses the mean differences and variability in the data series (Hawkins et al., 2013); and (ii) for precipitation calibration, we applied a quantile mapping approach which uses a statistical transformation function to adjust the distribution of ERA5 data to match the distribution of observations (Thrasher et al., 2012). To define the function, we applied a non-parametric transformation termed QUANT in the R package qmap (Gudmundsson et al., 2012). Once biases were reduced, a simple statistical metric enabled us to confirm that the mean difference between ERA5 data and in situ observations was improved by 90 % and the root mean square error (RMSE) reduced by 60 %, showing that the calibration methods performed quite well except for the RMSE of M188 station, which showed no improvement (Table 3).

Table 3 Temperature (T, ^∘C) and precipitation (P, mm) datasets used in this study. The mean difference and the RMSE between observations and reanalysis before bias correction and quantile mapping (bBC/QM) and after correction (aBC/QM).

3.6 Glacier clustering based on morpho-topographic features

To understand the relation between the thinning of each individual glacier and their morphometric characteristics over the entire period, we divided the glaciers into groups based on similarities in the elevation profile of annual changes in elevation. This requires a precise delineation of the boundaries of each of the glacial basins, which is difficult to achieve with medium-spatial-resolution imagery (e.g. Landsat images) of a conical structure like a volcano. However, we take advantage of the availability of high-resolution aerial imagery. Consequently, we applied the k-means classification technique (Wu, 2012). First, each glacier was divided into regular 50 m elevation ranges; we then computed the mean annual change in surface elevation for each elevation range; these values were then randomly divided into k groups, and the average for each group was computed. Next, the algorithm assigned each estimation to the closest mean value based on the Euclidean distance. The new mean value of each group was then updated and the process repeated iteratively until the sum-of-squares distance of the groups was minimised. Thus, it is assumed that values in the same group underwent similar changes in their surface elevation. To define the number of clusters, we performed k-means classification for different k groups, starting with 10 groups and continuing down to 1 group. We then compared the sum of squares of the groups, and the breaking point in the sum of squares is an indication of the optimal number of clusters. As a result, we divided the glaciers into two groups which match the exposure to the humid fluxes.

4 Results and discussion

4.1 Decadal changes in glaciers on the Antisana ice cap

Figure 3 gives an overview of the differences in elevation of the glaciers in the five sub-periods between 1956 and 2016. It should be noted that the five sub-periods are not representative of particular climate conditions but are based on the availability of remote sensing data. In all the periods, the strongest and most homogeneous changes in surface elevation occurred at lower elevations, then faded and became heterogeneous towards higher elevations. Such noisy patterns in the upper reaches of the ice cap (>5100 m a.s.l.) can be attributed to the rugged topography, where seracs, crevasses and steep slopes dominate the glacial landscape (e.g. average slope between 30 and 40^∘), combined with the very likely positive precipitation gradient (Basantes-Serrano et al., 2016), as discussed in more detail hereafter.

Figure 3 Spatial variability in changes in annual surface elevation of all glaciers on the Antisana ice cap in each study period. A decrease (increase) in the glacier surface elevation is shown by the red (blue) tone scale. The glacier boundaries in the first aerial survey are indicated by the grey lines, and areas for which no data are available are represented by the grey areas. The rate of elevation change and the surface area of the ice cap on the first date are given on the bottom left of each map (except for the last map illustrating the entire period for which the surface area measured on the last date is given as it is not mentioned elsewhere). All the maps have the same scale.

The average elevation change for the entire period (1956–2016) was $-\mathrm{30.01}\pm \mathrm{0.68}$ m with a minimum value of $-\mathrm{6.58}\pm \mathrm{0.70}$ m for G9; G1a, G5 and G16 had maximum values of $-\mathrm{54.83}\pm \mathrm{0.67}$ m, $-\mathrm{49.50}\pm \mathrm{0.68}$ and $-\mathrm{58.88}\pm \mathrm{0.66}$ m, respectively. A total reduction of 42 % in the surface area of the ice cap was found, from 23.93 km² in 1956 to 13.84 km² in 2016. In 2016, the average altitude of the glacier termini ranged from 4714 to 4794 m a.s.l. on the western slope and from 4334 to 4575 m a.s.l. on the eastern slope. Values were missing for some parts of the glaciers; the 2010–2016 period was the most affected with 26 % of the surface area of the ice cap not covered by measurements, representing an uncertainty of ±0.20 m.

Table 4 lists the glacier-wide mass balance for all the glaciers for each sub-period. Overall, during the entire study period, glaciers underwent an average cumulative mass loss of $-\mathrm{25.30}\pm \mathrm{0.58}$ m w.e. and an annual mass loss of $-\mathrm{0.41}\pm \mathrm{0.01}$ m w.e. yr⁻¹. G1a, G5 and G16 suffered the strongest mass losses with an average mass balance of $-\mathrm{0.74}\pm \mathrm{0.01}$ m w.e. yr⁻¹. On the contrary, G8, G9 and G10 showed moderate losses with an average annual mass balance of $-\mathrm{0.18}\pm \mathrm{0.01}$ m w.e. yr⁻¹. The annual mass balances of the remaining glaciers clustered between −0.29 and −0.57 m w.e. yr⁻¹. The cumulative mass balance computed from the five sub-periods is not exactly the same as the value given in the full period. The mean difference is equal to 0.05 m w.e., which is within the uncertainty in the mass balance. The biggest absolute difference can be seen in G1b and G4 with a mean value of 5.89 m w.e., followed by G12, G13 and G17 with a mean value of 3.25 m w.e.

Table 4 Glacier-wide cumulative mass balance in metres water equivalent (m w.e.) observed during the last 6 decades for the 17 glaciers on the Antisana ice cap. When mass losses occur the density conversion value is $\stackrel{\mathrm{‾}}{\mathit{\rho }}=\mathrm{850}\pm \mathrm{60}$ kg m⁻³; when mass gains occur the density conversion value is $\stackrel{\mathrm{‾}}{\mathit{\rho }}=\mathrm{564}\pm \mathrm{64}$ kg m⁻³. The * symbol indicates periods when more than 40 % of information concerning changes in the elevation of the surface area of the glacier was missing. Σ_period is the sum of glacier mass balance derived from the five sub-periods.

Now to look at the different sub-periods in more detail, Antisana glaciers experienced very negative to slightly negative or even positive mass balances, leading to a similar change in the surface area of all the glaciers in the ice cap (Figs. 4 and S4). Surprisingly, a very negative trend was observed in the 1956–1964 period, with an average annual mass balance of $-\mathrm{0.72}\pm \mathrm{0.05}$ m w.e. yr⁻¹. During this period, the northeastern glaciers (i.e. G17, G1a, G1b, G23), which represent about 40 % of the total surface area of the ice cap, had the most negative mass loss rates with an average of $-\mathrm{1.16}\pm \mathrm{0.05}$ m w.e. yr⁻¹, while the other glaciers lost mass at a 2-fold-lower negative rate ($-\mathrm{0.58}\pm \mathrm{0.05}$ m w.e. yr⁻¹). Despite the very negative rates observed in the first period, the reduction in glacier surface area was only 8 % and remained limited until the end of the 1970s.

In contrast to the 1956–1964 period, the 1965–1978 period presented a moderate gain in ice mass by all the glaciers, with an average positive annual mass balance of 0.10±0.04 m w.e. yr⁻¹, which helped offset part of the mass loss that had occurred in the previous period. Interestingly, during this period, only G5 showed a mass loss of the same order of magnitude as the gains by the other glaciers. After the late 1970s, the mass balance again became negative and a long period of ice mass loss started with an annual rate of $-\mathrm{0.82}\pm \mathrm{0.01}$ m w.e. yr⁻¹. The glaciers with the most negative mass balances were G1a, G1b, G5 and G16. This overall behaviour is in line with the results of previous studies highlighting a pronounced shrinkage of tropical glaciers since the late 1970s (Masiokas et al., 2020; Rabatel et al., 2006, 2013b; Soruco et al., 2009a). At the end of the 1990s, Antisana glaciers entered a period of limited mass loss, with an average annual mass balance of $-\mathrm{0.24}\pm \mathrm{0.03}$ m w.e. yr⁻¹, which has lasted for almost 20 years up to the present. The limited mass loss reported by Basantes-Serrano et al. (2016) for G15 was also confirmed for the other glaciers. By 2016, we observed the almost complete disintegration of G16 with a 95 % reduction in its surface area since 1956.

Figure 4 Annual mass balance in m w.e. yr⁻¹ over the different periods.

4.2 Sensitivity of the geodetic mass balance to the density assumption

In most of the geodetic studies, when there is no information available about changes in firn pack it is strongly recommended to use a conservative density value such as the one proposed by Huss (2013), especially in periods of mass loss. However, in our glaciers, the second period (1965–1978) is characterised by mass gain, and a density value close to the density of ice could lead to an overestimation of the mass balance. Assuming a density of 850 kg m⁻³ in both the accumulation and the ablation areas for the 1965–1978 period, the mass balance increases to 0.06 m w.e. yr⁻¹, which is within the uncertainty in the mass balance. In addition, during the 1998–2009 period, seven glaciers in the Antisana ice cap are close to equilibrium with a slightly positive or negative mass balance no matter what density scenario is assumed. Given the small difference between both assumptions, we decide to apply an average density value of 850 kg m⁻³ when mass losses prevail, and when positive conditions are present we use an average density of 564 kg m⁻³ according to observational data on the summit of the Antisana ice cap (Calero et al., 2022; Williams et al., 2002).

Figure 5 Antisana hypsometry in 1956 (red bars) and 2016 (cyan bars) and the elevation change rates for the study periods as a function of the altitude of each group of glaciers: (a, b) Group I composed of glaciers exposed to Amazon fluxes and (c, d) Group II composed of glaciers exposed to Pacific fluxes. The $\mathrm{d}h/\mathrm{d}t$ values are averaged over a 50 m elevation range.

Figure 6 (a) G8 shows an unusual pattern of elevation changes for the period 2009–2011 (also visible for the 2010–2016 period in Fig. 3) with a decrease in elevation (red) above 4800 m a.s.l. and an increase in elevation (blue) below 4800 m a.s.l., according to the aerial images acquired in September 2009 (b) and January 2011 (c). The dashed red lines represent the glacier outlines in September 2009, and the dashed black lines represent the glacier outlines in November 2011. The 0.5 m orthophoto (2011) used as the background was provided by the Ministerio de Agricultura y Ganadería (MAG) in the framework of the SIGTIERRAS project.

4.3 Morphometric drivers of changes in glaciers and effects on their interpretation

• Group I (Fig. 5a), with a B_g of −0.37 m w.e. yr⁻¹, includes very steep glaciers with slopes of more than 35^∘. The glacier terminus is located between 4100 m a.s.l. (in 1956) and 4334 m a.s.l. (in 2016). These glaciers are exposed to the east and south of the ice cap, facing the Amazon basin. Their ablation and accumulation zones are of the same width. The hypsometric changes show that the highest thinning rates occurred below 4900 m a.s.l., which is the mean altitude of the 0 ^∘C isotherm.

  This group includes G5, which showed the highest thinning rate ($-\mathrm{0.81}\pm \mathrm{0.04}$ m yr⁻¹) of the entire period. Also, during the 2010–2016 period, an anomalous increase in the surface elevation was detected in the tongue of G8. Thanks to the availability of the aerial photos taken in November 2011 by a high-resolution DMCii 143 sensor and after visual inspection of the orthophotos in 2009 and 2011, we were able to identify a surge-type event (Meier and Post, 196) that led to an advance of about 500 m of the glacier terminus with simultaneous thickening of the glacier tongue of about 12 m (Fig. 6). Indeed, above 4800 m a.s.l. G8 has an almost flat area (α>12^∘) which can be considered a reservoir area for the storage of ice masses if the glacier receives a significant amount of solid precipitation. In the case of surge events, after some time when the storage capacity is exceeded, the material flows down into the valley to the receiving zone that then advances rapidly. In the case of G8, this mass was transferred by a “channel” with a slope of up to 22^∘. Therefore, a large transfer of ice masses will appear as a drop in elevation in the departure zone (Fig. 6a, arrow a) and as a rise in elevation in the lower reaches of the glacier (Fig. 6a, arrow b). To our knowledge, surge events have never previously been reported in an inner-tropics glacier. Surging processes may be partially controlled by climate factors as well as by morphology, subglacial drainage patterns, glacier hypsometry and the magnitude of previous surge events In the present case, it could be hypothesised that sub-surface heating enhancing basal melt might be part of the triggers of this surge event, but no volcanic activity has been evidenced over the past 4 centuries.

• Group II (Fig. 5b), with a B_g of −0.49 m w.e. yr⁻¹, groups glaciers whose termini are located between 4550 m a.s.l. (in 1956) and 4700 m a.s.l. (in 2016). These fan-shaped glaciers have similar morphology with gentler slopes (5^∘ $<\mathit{\alpha }<\mathrm{15}$^∘) in the wider lower parts and steeper topography (α>35^∘) in the upper narrower zones. The change in surface area of these glaciers occurred between 4800 and 5000 m a.s.l., whereas upslope the change in the surface elevation of these glaciers was very limited. In this group, the average thinning of G1a and G16 ($-\mathrm{0.91}\pm \mathrm{0.04}$ m yr⁻¹) was strong over the entire period.

To understand the sensitivity of the glaciers to the morpho-topographic features in more detail, we performed a regression analysis between the average glacier-wide annual mass balance for the 1956–2016 period and morpho-topographic factors (Table 5). The balances for the entire study period were moderate but significantly correlated with the following variables: (i) maximum altitude; (ii) median altitude, a variable considered a proxy for the balanced-budget equilibrium-line altitude (ELA) (Braithwaite and Raper, 2009; Rabatel et al., 2013b); (iii) mean slope; (iv) the cosine of the mean exposure; and (v) the potential of incoming solar radiation under clear-sky conditions (Fu and Rich, 2003; Rich et al., 1994). The morpho-topographic characteristics of the glaciers explained 30 %–40 % of the observed long-term variation in the mass balance of the glaciers. However, when we considered the glaciers based on their group, the relationship between the mass balance and the maximum and median altitude, as well as solar radiation, was largely improved.

We observed that the slope played an important role in defining the rate of ice losses of the glaciers in Group I. Glaciers with a mean slope of less than 30^∘ and that extended above 5300 m a.s.l. presented a moderate imbalance. In fact, a relatively small ice loss is expected in steep slope sectors (Venkatesh et al., 2012) as long as the ice masses are present at altitudes where the conditions favour continued accumulation (see Table S1 in the Supplement). Also, around 60 % of the variance of the mass balance can be explained by solar radiation. When all the glaciers were included, we found a significant correlation between the mass balance and changes in surface area, and when the very unbalanced glaciers (G1a, G5 and G16) were removed, we found a lower but nevertheless significant correlation between mass balance and changes in surface area. The higher explained variance when the very unbalanced glaciers were included may be due to the morphometrical characteristics of these glaciers, whose maximum altitude is below 5300 m a.s.l., leaving a wider ablation area than the accumulation area and thus leading to a bigger reduction in surface area.

Table 5 Correlation coefficient R² for the regression analysis between the conventional mass balance and the morpho-topographic features and solar radiation for the period 1956–2016. Considering the mass balance B_m is for all the glaciers in the ice cap, and ^′B_m is excluding the outlier glaciers (G1a, G5 and G16). The level of significance is indicated with p value <0.01 (***), p value <0.05 (**) and p value <0.10 (*).

To test the sensitivity of the glacier mass balance to variations in the glacier surface area, we computed the glacier-wide mass balance for each sub-period by keeping the glacier surface area unchanged (hereafter glacier-wide fixed mass balance B_g,f).

This exercise was inspired by the work on the reference mass balance at an annual time step conducted by Elsberg et al. (2001), Harrison et al. (2005) and Huss et al. (2012), which was focused on the effect of climate in the glacier response without taking into account the change in geometry due to flow dynamics. We thus replaced the term $\stackrel{\mathrm{‾}}{S}$ in Eq. (5) with S_o, the surface area of the glacier at the date of the first geodetic survey in the study sub-period concerned. Therefore, one would expect the mass balance of the entire ice cap to be underestimated if the mean variation in the surface area of the small glaciers due to their imbalance with climate is not taken into consideration.

The cumulative mass balance of the whole ice cap over the entire period resulted in a discrepancy of 3.6 m w.e. between ∑B_g and $\sum B_{\mathrm{g},\mathrm{f}}$, which means that B_g,f is 0.11 m w.e. yr⁻¹ less negative than B_g (−0.41 m w.e. yr⁻¹). Considering the glaciers according to their group, the difference between B_g and B_g,f was 0.09 m w.e. yr⁻¹ for glaciers in Group I, which falls within the uncertainty in the mass balance quantification (Fig. 7a). This was not the case for glaciers in Group II where the difference between B_g and B_g,f was 0.15 m w.e. yr⁻¹ (Fig. 7c). At the scale of the sub-periods, the biggest differences between the two balances occurred in the 1979–1997 period with an overestimation of 0.11 m w.e. yr⁻¹ for Group I and of 0.20 m w.e. yr⁻¹ for Group II (Fig. 7b and d). For the other sub-periods, the difference is within the uncertainty in the average mass balance.

When G5 (Group I) and G1a and G16 (Group II) were excluded, the discrepancy between B_g and B_g,f was reduced (Fig. 7e to h), showing that (i) up to 30 % of the overestimation in the mass balance becomes apparent if the mean variation in the surface area of the ice cap is disregarded and (ii) most of this discrepancy may be due to a marked imbalance between the geometry of the small glaciers and the climate.

Figure 7 Cumulative mass balance and centred mass balance of the glaciers in Group I (including G5 in panels a and b and excluding G5 in panels e and f) and of the glaciers in Group II (including G1a and G16 in panels c and d and excluding G1a and G16 in panels g and h). The black line represents the conventional balance, and the red line represents the fixed mass balance. The grey-shaded areas represent the standard deviation of the mass balance.

Thanks to the availability of recent geodetic observations provided by satellites, several studies have been undertaken to estimate the glacier mass loss at a regional scale in recent decades (e.g. Braun et al., 2019; Brun et al., 2017; Dussaillant et al., 2019). Because no multitemporal inventory is available to quantify changes in the glacier surface area in the majority of the glacierised regions in the world, many regional studies (i.e. Braun et al., 2019; Brun et al., 2017; Dussaillant et al., 2019) assume the glacier surface area remains constant. As underlined by Dussaillant et al. (2019) at a regional or mountain range scale, the impacts of this assumption are limited, but at a glacier scale and particularly in periods with significant glacier retreat, not accounting for changes in the surface area of the glacier may result in overestimation of ice mass loss as is shown in this study (see Sect. S1 in the Supplement). This inaccurate estimation was as high as 30 % in the case of the smaller glaciers of the Antisana ice cap studied here (Fig. S4 in the Supplement). This fact is extremely important, for example, at a local scale, when assessing the contribution of glacier melt to the catchments that supply water to local communities.

4.4 Glacier response to climate variability

To analyse the decadal variability in glacier volume change, we included B_g but not glaciers with a strong geometric effect, and we computed the mean centred mass balances of the glaciers for each sub-period by subtracting the average mass balance for the 1956–2016 period. This resulted in centred mass balance values of −0.30 m w.e. yr⁻¹ ($\mathit{\sigma }=\pm \mathrm{0.31}$), 0.52 m w.e. yr⁻¹ ($\mathit{\sigma }=\pm \mathrm{0.10}$), −0.40 m w.e. yr⁻¹ ($\mathit{\sigma }=\pm \mathrm{0.43}$), 0.17 m w.e. yr⁻¹ ($\mathit{\sigma }=\pm \mathrm{0.25}$) and 0.14 m w.e. yr⁻¹ ($\mathit{\sigma }=\pm \mathrm{0.17}$) for 1956–1964, 1965–1978, 1979–1997, 1998–2009 and 2010–2016, respectively. Periods with negative mean centred mass balance values (i.e. negative anomalies) showed more dispersion between glaciers than periods with positive anomalies. However, the overall picture shows a homogeneous response of glaciers to climate variability. Although the periods for which the mass balance variability was reconstructed based on the available images using the geodetic method are constrained by the dates of the surveys, the mass balance reflects the decadal climate conditions and/or changes.

To link the regional climate to the decadal variability in the mass balance, we compared the centred geodetic mass balance with climate anomalies computed from ERA5 reanalysis data, i.e. anomalies in temperature, and in specific humidity as a proxy for precipitation anomalies (Fig. 9) (Ruiz-Hernández et al., 2021). To be able to see the footprint that ENSO events may have left on the decadal mass balances, we included in the comparison the Southern Oscillation index (SOI) averaged for each sub-period (Ropelewski and Jones, 1987). SOI is a standardised index based on changes in ocean temperatures across the eastern tropical Pacific during ENSO events. Negative values correspond to warm conditions during El Niño, while positive values correspond to cold temperatures during La Niña episodes. To understand the effect of seasonal variability in the long-term glacier response, we focused our analysis on the key quarters of the year, MAM and SON (see Sect. S2 in the Supplement). Anomalies were determined relative to the current normal climate period 1981–2010. For instance, a year is considered to be anomalous if temperature anomalies are greater than ±0.5^∘C and precipitation anomalies are less than ±10 %. Likewise, as a proxy for the energy available for melting, we included the estimated average anomaly of downward solar radiation (srad) per period taken from the ERA5 reanalysis as a proxy for shortwave radiation at the Earth's surface. It is worth mentioning that this is a broad approximation based on the average solar radiation over a grid box on the ice cap.

Figure 8 Climate anomalies relative to the current normal climate period 1981–2010 in (a) temperature and (b) specific humidity at pressure levels over the Antisana ice cap. (c, d) The same as in (a) and (b) but for the MAM quarter. (e, f) The same for the SON quarter. Pale to dark red represents warm and dry, whereas pale to dark blue represents cold and humid. The dashed line shows the maximum altitude of the Antisana ice cap. (g) The centred mass balance for the Antisana ice cap (black line) and the (1σ) standard deviation (grey-shaded area) for each sub-period; the orange (green) line represents SOI fluctuation averaged over the MAM (SON) quarter. (h) The centred mass balance compared with the 11-month moving average of the downward solar radiation anomaly over the Antisana ice cap (red line).

As is well known, glaciers in this region are particularly sensitive to humidity and as a consequence to precipitation and clouds (Mölg et al., 2009; Prinz et al., 2016; Sicart et al., 2005). The main patterns that emerge from Fig. 8 are described below:

• In the 1956–1964 period, a notable loss of ice mass was observed which can be attributed to the dry conditions that prevailed in the 1950s and were further accentuated in the 1960s. Dry conditions were even more marked in the SON quarter. Regarding temperature, this period was characterised by a negative anomaly ($\mathrm{\Delta }T<-\mathrm{0.2}$ ^∘C compared to the mean calculated over the period 1981–2010). Thus, while dry conditions may have limited accumulation in the upper reaches of the ice cap, a concomitant reduction in cloudiness may have been responsible for an increase in shortwave radiation (Fig. 8h), thereby enhancing melt due to the decrease in albedo. Sporadic snowfall in the MAM quarter may not have been sufficient to reduce the steep ablation rates.

• In the 1965–1978 period, the glacier mass balance was stable or even positive. These conditions are consistent with negative temperature anomalies and a mostly positive humidity anomaly, which likely triggered snowfall all year round, thereby increasing glacier surface albedo, reducing ablation and enhancing accumulation processes over the ice cap. Later, in a colder and more humid context, the probable increase in cloudiness could have reduced the energy available for melt (negative solar radiation anomaly, Fig. 8h), and snowfall events would continue to enhance accumulation processes until the mid-1970s. In the late 1970s, slightly warm and moderately dry conditions emerged which started the long negative period for the ice cap described below.

• In the 1979–1997 period, strong ice mass losses predominated. Warm conditions prevailed and were slightly more marked during the key quarters ($\mathrm{\Delta }T>+\mathrm{0.25}$ ^∘C). These conditions were combined with a long, slightly wet period, which was however interrupted by some dry years. Particular features were as follows. (i) A warm and humid context during MAM in the 1990s provided favourable conditions for rainfall on the tongues of the glaciers, thereby increasing melt rates. A very strong El Niño occurred in 1983 characterised by a positive temperature anomaly combined with very humid conditions favouring rainfall on the lower reaches of the glaciers, which would increase ablation rates via the shortwave radiation budget due to the low albedo. (ii) After 1985, warm conditions combined with the surplus humidity recorded in the SON quarter enhanced melt rates, which were further increased during the strong El Niño in 1987. The observed slight increase in incoming shortwave radiation may have helped maintain high melt rates.

• In the period following the late 1990s–early 2000s, slightly negative or even positive mass balances were documented. In humid conditions, presumably more continuous cloudiness over the volcano helped reduce shortwave radiation at the glacier surface. On the other hand, precipitation and slightly colder episodes could have maintained the snow cover long enough to protect the glaciers from the energy available for melting. One strong cold La Niña event occurred in 1999–2000 along with two warm El Niño events, one in 1998 and the other in 2015–2016. After 2015, conditions suddenly changed to become moderately warm and humid over the glaciers, increasing rainfall rather than snowfall and reducing the albedo and hence increasing ablation rates.

• The consistency between the decadal variability in the mass balance, the variability in the shortwave radiation available for melting, and the large-scale fluctuations in the sea surface temperatures between the western and eastern tropical Pacific represented by the SOI index is clearly visible throughout the study period (Fig. 8g).

5 Conclusion

This paper has provided a detailed view of the changes to glaciers in the Antisana ice cap, Ecuador, i.e. in the inner tropics, since the middle of the 20th century. Despite the contrasted decadal variability in the mass balance, an overall negative trend has prevailed over the last 7 decades with an annual mass loss rate of $-\mathrm{0.41}\pm \mathrm{0.04}$ m w.e. yr⁻¹. This study revealed that glacier mass losses were remarkably high ($-\mathrm{0.72}\pm \mathrm{0.06}$ m w.e. yr⁻¹) in the 1956–1964 period, whereas a moderate gain in ice mass was detected in the 1965–1978 period (0.10±0.05 m w.e. yr⁻¹) which partially offset glacier shrinkage in the preceding period. After the late 1970s, ice mass loss continued until the late 1990s ($-\mathrm{0.82}\pm \mathrm{0.04}$ m w.e. yr⁻¹). However, since the early 2000s, thinning of the glaciers of the Antisana ice cap has slowed down, as revealed by slightly negative mass balances ($-\mathrm{0.26}\pm \mathrm{0.05}$ m w.e. yr⁻¹). The overall change in the mass of the ice cap between 1956 and 2016 reduced its surface area by 42 %. G12 and G13, located on the western side, roughly followed the same general pattern observed in each sub-period, confirming that these glaciers are fairly representative of the decadal glacier response and could thus be used as benchmark glaciers for future studies.

Overall, 70 % of the thinning occurred below the average elevation of the 0 ^∘C isotherm (around 4900 m a.s.l.). We distinguished two types of glacier whose response to climate variability is determined by their morpho-topographic characteristics. The spatial pattern of the changes in glacier surface elevation is more marked on the eastern slopes of the ice cap where the glaciers are more directly influenced by the interplay between the flows of humidity from Amazonia and the rugged and heterogeneous ice cap morphology. Interestingly, a surge-type event was detected on G8, underlining the importance of including the effect of ice-flow dynamics as well as climate variations when studying glacier response. To our knowledge, no similar event has been reported in the tropics to date; thus more research is needed before being able to conclude on the internal (ice-flow dynamics) or external (climate, sub-surface heating due to volcanic activity) factors that triggered such an event.

Regression analysis between mass balance rates over the whole period and the morpho-topographic variables suggests that about 40 % of the mass balance spatial variability can be explained by the topography. However, glaciers located on the western side reacted more obviously to topographic features. We found that glaciers whose termini reached 4100 m a.s.l. in 1956 retreated further than glaciers whose termini were located at 4700 m a.s.l. Furthermore, small glaciers with a maximum elevation below 5300 m a.s.l. have almost disappeared (∼86 % loss of surface area since 1956), as reported for other glaciers in the region, and the remaining portions of ice will probably disappear within the coming decade. Additionally, our study shows that estimations of geodetic mass balance are prone to error when variations in glacier surface areas are neglected, particularly at the glacier scale and in the case of small glaciers when the loss is high. This can indeed lead to a misinterpretation of the glacier response to climate and may have a major impact on hydrological modelling used to estimate the contribution of meltwater in catchments that contain small glaciers.

Overall, the behaviour of glaciers on the Antisana ice cap is consistent with the regional climate signal, which has shown a stepwise transition towards warmer and alternating wet–dry conditions since the mid-1970s. Sub-periods with extreme mass balance anomalies have more scattered precipitation values than years with moderate or even slight mass balance anomalies.

Our results provide a starting point for future research on the physical processes that drive these changes in the complex topographical–climatic conditions found in the tropical Andes. Due to the importance of this site as a water resource for the population of the city of Quito, we also recommend analysis of the contribution of meltwater to the hydrological budget of the catchments, which is crucial for effective management of water resources under climate change.