Slight mass loss revealed by reanalyzing glacier mass-balance observations on Glaciar Antisana 15α (inner tropics) during the 1995–2012 period

Abstract

In this paper, we reanalyze the glacier mass balance on Glaciar Antisana 15α over the 1995–2012 period. Annual glacier mass balances were quantified on the basis of monthly glaciological measurements using an adaptation of Lliboutry's statistical approach. The geodetic mass balance was computed between 1997 and 2009 giving a cumulative balance of −1.39 ± 1.97 m w.e. and a slightly negative adjusted annual glaciological mass balance (−0.12 ± 0.16 m w.e. a⁻¹). Despite a careful analysis of uncertainties, we found a large discrepancy between the cumulative glaciological and the geodetic mass balances over the common period, of 4.66 m w.e. This discrepancy can mainly be explained by underestimated net accumulation in the glacier upper reaches, which could be due to the peculiar climate conditions of the equatorial zone with year round accumulation, thereby preventing clear identification of annual layers. An increase of ~70% in measured rates of net accumulation would be needed to balance the glaciological and geodetic mass balances; a hypothesis confirmed by estimated ice flux in the vicinity of the ELA. Consequently, the vertical gradient of precipitation may be higher than previously estimated and the accumulation processes (including the role of frost deposition) need to be carefully analyzed.

1. INTRODUCTION

Several studies have highlighted the high sensitivity of tropical glaciers to climate change (e.g. Kaser and Osmaston, 2002; Bradley and others, 2009). These glaciers are a proxy to better understand climate variations at annual to decadal scales. A major glacier recession is documented for the tropical Andes since the late 1970s (Rabatel and others, 2013) linked with climate variability at the interannual scale, mainly forced by El Niño Southern Oscillation (ENSO), and, at the decadal scale, linked with the trend to increasing atmospheric temperatures (Francou and others, 2004; Vuille and others, 2008). To document and understand the interannual to decadal changes in these glaciers, there is a need for accurate long-time series of glacier mass balances measured on several reference glaciers. This is also a major consideration in the monitoring strategy of the Global Terrestrial Network for Glaciers (http://gtn-g.org) focused on: ‘extensive glacier mass-balance studies within major climatic zones for improved process understanding and calibration of numerical models’.

In the Ecuadorian Andes, a monitoring program began in June 1994 on Glaciar Antisana 15α (Francou and others, 2000). Today, it is a joint program of the GLACIOCLIM observatory and the Instituto Nacional de Meteorologia e Hidroligia (INAMHI, Ecuador). Mass-balance data are freely available through the GLACIOCLIM website (http://www-lgge.ujf-grenoble.fr/ServiceObs/index.htm) and the World Glacier Monitoring Service (WGMS, 2012). Up to now, data collected on Glaciar Antisana 15α have been used to: (1) study the physical processes governing ablation at the glacier surface (Favier and others, 2004), (2) understand the response of an equatorial glacier to regional climate variability, particularly when conditions are forced by ENSO variability (Francou and others, 2004); and (3) assess the contribution of glacier melt to river flow and water availability (Villacis, 2008).

However, specific annual mass balances on Glaciar Antisana 15α are very sensitive to systematic and random errors depending on the reference surface, the spatial distribution of the direct measurements (stakes/pits/cores), the spatial variability of the surface mass-balance components and the interpolation/extrapolation method used to cover areas for which no measurements are available. These uncertainties can lead to unreliable interpretation of the long-term relationship between climate and the glacier. The aim of this study was to reanalyze the glaciological mass-balance series obtained from in situ observations, to quantify the uncertainties, and to adjust the series on the basis of the geodetic mass balance obtained using aerial photogrammetry. To this end the analytical framework proposed by Zemp and others (2013) was applied; this approach is summarized in a generic scheme of six basic steps (see Fig. 1 in Zemp and others, 2013). The standard terminology and notations proposed by Cogley and others (2011) are used.

2. STUDY SITE AND CLIMATE SETTINGS

Antisana (0°28′S; 78°09′W) is an irregular cone-shaped stratovolcano located 40 km east of Quito (capital of Ecuador) in the Cordillera Oriental (Fig. 1a). The volcano is covered by 17 glaciers as listed by Hastenrath (1981). The present paper focuses on Glaciar Antisana 15α (5700–4850 m a.s.l., 0.28 km², 1.6 km long) located on the northwestern side of the volcano (Fig. 1b). Since the mid-1950s, the lower part of Glaciar Antisana 15 has been divided into two tongues 15α and 15β (Francou and others, 2000).

Fig. 1. (a) 2009 orthophoto-map of Glaciar Antisana 15, red and orange contours are the glacier extents in September 2009 and August 1997, respectively. The dashed black contour is the restituted area; 50-m interval contour lines are shown. The black dots show the ablation stakes in 2010, the red stars show the location of the accumulation pits, the blue dashed circle is the area where accumulation measurements and snow probing were made, the red line at ~5000 m a.s.l. shows the uppermost GPR cross section, the white points represent the four rain gauges, and the yellow triangle shows the location of the weather station. Coordinates are given in UTM zone 17S, WGS84. (b) Green triangles show the location of the 21 GCPs used to carry out the photogrammetric adjustment, and the red circles show the locations of the meteorological stations around the volcano. Glaciar Antisana 15 is shown in red. Inset (c) map of Ecuador with the location of Antisana Volcano (red dot) and Quito, capital of Ecuador (blue dot).

As described by Francou and others (2004), the seasonal cycle of precipitation is bimodal and there is no clear dry season. One precipitation maximum occurs from March to May followed by a minimum from July to October, and a second maximum occurs in November (Fig. 2). Annual precipitation measured on the glacier foreland does not exceed ~1300 mm a⁻¹. Cloud cover follows roughly the same pattern as precipitation, whereas air temperature remains fairly constant throughout the year and interannual variability is ~3°C (standard deviation (SD) over the last 50 years). The 0°C isotherm, the ‘freezing line’, shifts up and down across the ablation zone of the glacier (between 4800 and 5100 m a.s.l.) and a minor variation in air temperature can influence melt processes at the glacier surface by determining the precipitation phase (snow or rainfall), thereby affecting the surface albedo in the ablation zone; an increase of 1°C in temperature can move the snow/rain limit 150 m up the glacier (Francou and others, 2004; Rabatel and others, 2013).

Fig. 2. Climate conditions on the Glaciar Antisana 15α foreland. The gray bars represent the average monthly precipitation from a precipitation index (Pglacier Index) computed from the four rain gauges shown in Figure 1 (Manciati and others, 2014). The black dots and the line represent the average monthly temperatures recorded at the weather station on the glacier moraine at ~4850 m a.s.l. (Fig. 1) over the 1995–2009 period. Vertical bars represent the one SD.

3. DATA

3.1. Glaciological in situ observations

A total of 15 to 20 stakes in the ablation zone of the glacier have been measured monthly since June 1994. These stakes are located between the snout (currently ~4850 m a.s.l.) and the ELA0 (~5080 m a.s.l.), where the gently sloped glacier tongue provides safe conditions to manage the stake network. The stakes are located in the center and on the lateral sides of the glacier tongue to capture the spatial variability of the mass balance (Fig. 1a).

Access to the accumulation areas of the glaciers is challenging in this part of the Andes since all glaciers are located on volcanoes with rugged topography in their upper reaches. From 5100 m a.s.l. to the summit, Glaciar Antisana 15α is in the form of a narrow corridor with steep sides and is subject to frequent avalanches. These avalanches are caused by seracs falling from the large crevasse surrounding the dome-shaped summit. Two snow pits/cores were dug in this zone to estimate annual accumulation on a relatively gentle slope although containing some crevasses (between 5400 and 5650 m a.s.l., at the limit with Glaciar Antisana 15β). These accumulation measurements were completed by probing the snow in the vicinity (~50 m radius) (Fig. 1a). In each snow pit/core, snow stratigraphy and variations in density with depth were measured. Accumulation was measured at the end of the hydrological year, which is not well defined close to the equator, but was fixed in late December/early January in a period in which precipitation is generally low. It is worth noting that the lack of precipitation seasonality in the equatorial Andes makes it difficult to identify a well-marked crust, i.e. a layer whose density is close to the density of ice at the base of the current year's accumulation.

As part of monitoring, annual changes in glacier geometry (i.e. snout retreat, changes in surface area) and surface flow velocities at the level of the ablation stakes were quantified using geodetic observations obtained with the Differential GPS (DGPS).

In addition, ground penetrating RADAR (GPR) measurements were conducted in September 2000 (J Ramirez, unpublished data) to determine the thickness of the ice in four cross sections from the glacier snout up to ~5000 m a.s.l. Only the uppermost cross section located near the equilibrium line (Fig. 1a) was used to estimate the ice flux and validate our hypothesis concerning net accumulation in the upper reaches of the glacier.

3.2. Aerial photographs

Aerial photographs acquired in 1997 and 2009 over Antisana Volcano were collected. These data comply with the peculiarities required for glaciological applications: (1) a minimum snow cover in the ablation zone, (2) no cloud cover over the glacier and minimum cloud cover over the surrounding areas, (3) images of the required brightness and contrast, and (4) known optimal geometric characteristics (i.e. available calibration reports). The films were scanned at 14 μm resolution using an Intergraph PhotoScan TD system. The characteristics of the aerial photographs for each survey are listed in Table 1.

Table 1. Characteristics of the aerial photographs and topographic restitutions on the Glaciar Antisana 15α and its catchment

3.3. Meteorological data

Precipitation measured at five rain gauges located on the glacier foreland, known as P0 (4850 m a.s.l.), P2 (4780 m a.s.l.), P3 (4550 m a.s.l.), P4 (4450 m a.s.l.) and P5 (3930 m a.s.l.) were used. These datasets were analyzed and validated by Manciati and others (2014). We also used precipitation series from four meteorological stations named Papallacta, Cuyuja, Baeza and Cosanga (Fig. 1b). All these data series were used to interpret the mass-balance processes occurring on Glaciar Antisana 15α and their impacts on the specific annual mass balance.

4. METHODS

4.1. Construction of DEMs using aerial photogrammetry

Aerial photogrammetry is widely used in glaciology to estimate fluctuations in volume (e.g. Rabatel and others, 2006; Bauder and others, 2007; Thibert and others, 2008; Soruco and others, 2009). This method makes it possible to obtain changes in the surface elevation over a long period (~10 a), and has the advantage of covering the entire glacier surface.

All photogrammetric tasks were carried out using the software Orima-DT from the Leica Photogrammetry Suite (LPS) of the Intergraph Company^©. The aerial photographs were adjusted using an average of 25 tie points for each pair of stereo images. For each geodetic survey, a bundle block adjustment was performed on the basis of the selected aerial photos using a geodetic network of 21 ground control points (GCPs) around the volcano observed on stable rocky features by rapid static DGPS survey during several campaigns in 2009 and 2011 (Fig. 1b, the 3-D RMS error, (RMS)xyz of GCPs network is 0.20 m). Setting up GCPs on the eastern side of the volcano was challenging, as access is very difficult, but a good spatial distribution of control points was achieved. Which enabled us to obtain the optimal geometry to adjust the photogrammetric block without jeopardizing internal accuracy (Table 1).

The topographic mapping by photogrammetric restitution of the glacierized and surrounding areas and the editing of these features were made by applying a discrete data collection (isolated points). The collection of three-dimensional (3-D) features and terrain editing tasks on the pairs of stereo images were simplified by the fact that all the aerial surveys were made around noon, thus reducing shadows and enhancing visual contrast between ice, snow and off-glacier areas.

The DEMs of 10 m cell (pixel) size were created using the minimum curvature interpolation method based on the point cloud for each date.

To assess the uncertainties in the photogrammetric restitutions, we compared five topographic profiles located on the stable non-glacierized terrain on the western side of the volcano. Three topographic profiles were acquired by a rapid static DGPS survey and two were collected by photogrammetry for each survey (Fig. 3). The comparison of these profiles gave an average vertical RMS error, RMSz of ±1.1 m. In addition, the small dome-shaped summit of the volcano allowed us to evaluate changes in elevation in the upper part of the accumulation zone where changes in the thickness of the glacier are expected to be very small. Indeed, because the summit is exposed to strong eastern winds, part of the snow may be drifted on the leeward side to the volcano. We compared changes in elevation at 18 sites measured on the flat dome-shaped summit giving an RMSz of ±1.5 m. This quasi-flat area represents ~1% of the entire 15 km² glacierized area of the volcano.

Fig. 3. Difference in altitude on cross sections measured on non-glacierized terrain between photogrammetric measurements in 1997 (green dots) and 2009 (red dots) and DGPS surveys (a–c); and between photogrammetric measurements in 1997 and 2009 (d, e).

We also evaluated the uncertainty concerning the differences in elevation by comparing the two DEMs on stable terrain. The remaining uncertainty in the difference in elevation was calculated, giving an RMSΔz and a SD of ±2.1 m (Fig. 4). Note that negative elevation changes (red pixels) evidenced some dead-ice debris-covered patches in 1997 that decreased in volume leading to a change in the morphology of the catchment. A few remnants of these structures were observed during field campaigns and were confirmed by the changes that were visible on the orthophotos. These features were not taken into consideration in order to consider only stable terrain. The standard error in the difference in elevation was 0.04 m (5131 pixels) ensuring good agreement between 1997 and 2009 DEMs. Some authors (e.g. Nuth and Kaab, 2011; Zemp and others, 2013) recommend performing a 3-D coregistration of successive DEMs when systematic errors are identified related to: (1) the displacement in the xyz location of the successive DEMs; (2) an elevation-dependent bias; and (3) differences in the sensor acquisition geometry. In our case, the two DEMs were consistent: the difference in elevation on stable terrain was sufficiently accurate, and the same geodetic reference was used for each aerial survey during the aerotriangulation process. As a consequence, 3-D co-registration of the DEMs was not required.

Fig. 4. Map of changes in elevation over the non-glacierized terrain between [DEM1997 – DEM2009] for 5131 pixels (spatial resolution =10 m). The normal distribution, mean and RMS of the elevation changes are presented in the inset graph.

4.2. Quantification of the mass balance using the geodetic method

Algebraic subtraction between two successive DEMs gives the change in surface elevation for each pixel over the period concerned. Changes in surface elevation are then converted into changes in volume resulting in a hypothesis concerning the density of mass lost/gained over the average glacier surface area in the time span of the geodetic surveys. The cumulative mass balance is computed using the equation (Zemp and others, 2013).

(1) $$B\hbox{geod}_{[\rm{year}_2 - \rm{year}_1]} = \displaystyle{{{\bar \rho} \times {r^2} \times \mathop \sum \nolimits_{i = 1}^n \Delta {h_n}} \over {\dot {\bar S}}},$$

where ${\bar \rho} $ is the average density value of the volume change, r ² is the pixel size, Δh _n is the change in the glacier surface elevation for each pixel, n is the number of pixels covering the glacier at its maximum extent, and $\dot {\bar S}$ is the average glacier area over the period [year₂ − year₁].

4.3. Quantification of the mass balance using the glaciological method

The glaciological method was used to compute the specific annual mass balance (Cuffey and Paterson, 2010).

In previous studies (e.g. Francou and others, 2004; Manciati and others, 2014), the glacier was divided into irregular elevation ranges that varied from year to year and depended on the location of the measurement sites. In addition, the mass-balance quantification only included a limited number of ablation stakes selected as representative of a specific elevation range. The mass balance for each elevation range was weighted by its relative surface area. For the elevation ranges in which no direct measurements were available, e.g. between the uppermost ablation stake at ~5080 m a.s.l. (and up to ~5200 m a.s.l. in a few years) and the lowest snow pit/core (~5370–5470 m a.s.l.), a linear interpolation was made. The same consideration was assumed from the lowest to the highest snow pit/core (~5650 m a.s.l.). Between the highest snow pit/core and the summit, a constant net accumulation value equal to the value of the highest snow pit/core was assumed. Note that ~60% of the glacier area is not monitored by direct measurements and that most of this area is in the accumulation zone.

For the present study, all the raw field measurements were processed to compute the specific annual mass balance. This represented 3288 stake measurements in the ablation zone and 50 pits/cores in the accumulation zone over the study period from August 1994 to December 2012. The specific annual mass balance was computed by referring to a fixed date (January 1) and to the surface area of each 50 m elevation range for each year. For the computation, we used the following steps:

• Adjustment of reference surface area. For this adjustment, the annual fluctuations in glacier extent and changes in the hypsography were taken into consideration. For the lower reaches of the glacier (i.e. below 5100 m a.s.l.) the annual DGPS observations of the glacier outline revealed nonlinear temporal changes in this area of the glacier. Above 5100 m a.s.l. the topographic map made in 1997 and 2009 by photogrammetry restitution revealed no important change in glacier topography (location of the elevation contours and glacier boundary). Accordingly, the reference surface area was adjusted as follows: (1) above 5100 m a.s.l., the glacier was divided into 50 m elevation ranges using the 1997 topographic map, (2) below 5100 m a.s.l. we used the annual DGPS measurements of the glacier outline and the elevation lines from the two photogrammetric restitutions (1997 and 2009) to calculate the surface area of the 50 m elevation ranges using the formulation,

  (2) $${{\bar S}_{{\rm e}.t}} = {S_{{\rm ref}.t}} + {\rm \;} \displaystyle{{S^{\prime}} \over N} \times (t - {t_0}) + \displaystyle{{S^{\prime \prime}} \over N} \times (N - t + {t_0}),$$

  where ${{\bar S}_{{\rm e}.t}}$ , is the adjusted elevation range area for the elevation range and the corresponding year t. S _ref.t , t is the surface area that was common to each elevation range in 1997 and 2009. Time t ₀ is the first photogrammetric year, 1997. S′ and S″ are the surface areas between the upper/lower elevation contours of the elevation range area in 1997 and 2009. (N) is the number of years, i.e. 12.

  From 2009 to 2012, the glacier outline and the elevation lines resulting from the 2009 geodetic survey were not affected above 5100 m a.s.l. Below this elevation, the annual glacier tongue outlines obtained by DGPS surveys were included. In addition, to compensate for the lack of DGPS observations of the glacier outline in 1995 and 1996, we kept the same trend for these two years as that observed after 1996, leading to a change of 1% in the surface area.

• Interpolation/extrapolation method used to estimate the specific mass balance. A major source of uncertainty in the quantification of the glacier-wide mass balance is the lack of measurements in inaccessible areas due to crevasses, seracs and the risk of avalanches. In fact, the main problem involved in computing the specific annual mass balance is the interpolation/extrapolation technique used for the unmeasured areas. As shown by Soruco and others (2009), simple interpolations can lead to significant errors in the computation of the glacier mass balance. In order to better account for the spatial variability, we adapted the statistical approach developed by Lliboutry (1974) (e.g. Vincent and others, 2000; Thibert and others, 2008; Eckert and others, 2011). This approach assumes that the mass balance can be decomposed into a spatially independent temporal variation β _t and a spatial pattern independent of time α _j , thus:

  (3) $${b_{\,jt}} = {\alpha _j} + {\beta _t} + \varepsilon, $$
  where b _jt is the specific mass balance at the elevation range j during the year t, and ε is an error term.

  However, in this classical approach, the temporal term (β _t ) is assumed to remain constant over the entire surface area of the glacier, which is not the case for Glaciar Antisana 15α. Indeed, Figure 5 shows that temporal variability of the mass balance of this glacier decreases markedly with elevation. About 80% of the variance originates in elevation ranges below 5300 m, representing roughly one-third of the total surface area. The centered annual mass balances computed for distinct elevation ranges confirmed high variability of the mass balance in the ablation zone compared with the accumulation zone, which shows very low variability (Fig. 6), similar to that observed by Soruco and others (2009) on Glaciar Zongo in the outer tropics.

  Consequently β _t cannot be assumed to be independent of altitude. Note that the interannual variability of the surface mass balance above 5600 m a.s.l. remains within measurement uncertainty. To account for the dependence on elevation, the specific mass balance (b _jt ) measured at each range was normalized using an elevation factor [Z _max − Z _j  /Z _max − Z _min], before using Lliboutry's approach to compute the specific mass balance. In this elevation factor, Z _j is the altitude of each elevation range considered, Z _min is the altitude of the lowest elevation range of the glacier, and Z _max is the altitude at which temporal fluctuations can be disregarded. As a consequence, the term β _t represents the temporal fluctuation of the lowest parts of the glacier (below 5080 m a.s.l.). Lliboutry's approach was applied to the mean specific mass balance computed from the stakes located in each elevation range of the glacier for each year.

Fig. 5. Multi-annual average mass-balance profiles b(z) of Glaciar Antisana 15α. Average mass-balance curve interpolated linearly (black dashed line) and computed using Lliboutry's approach (red line).

Fig. 6. Centered specific mass balances for different elevation ranges over the 1995–2012 period.

4.4. Quantification of the ice flux using the kinematic method

The GPR measurements made it possible to identify the bedrock topography along the cross section. At the level of the uppermost cross section at ~5000 m a.s.l. (close to the ELA) the average ice thickness was 48 ± 4.3 m resulting in a surface area of the cross section of 1.09 ± 0.14 × 10⁴ m². At the same elevation, the annual surface flow velocity was estimated from the displacement of two stakes, giving an average velocity of 30 ± 4.9 m a⁻¹. Following the approach presented by Azam and others (2012), the ice flux through this cross section was estimated to be 2.8 ± 0.29 (×10⁵) m³ a⁻¹.

5. RESULTS

5.1. Geodetic mass balance

Glaciar Antisana 15α showed a geodetic mass balance (B _geod.t ) of −1.99 m w.e. between 1998 and 2009. Figure 7 shows the spatial distribution of the changes in the surface elevation over the glacier. Note that in the ablation zone, below 5050 m a.s.l. where the terrain is homogeneous and has an average slope of ~30%, the lowering of the glacier surface regularly decreased with elevation (maximum ice thinning reached ~20 m close to the glacier snout). On the other hand, above 5050 m a.s.l., a complex heterogeneous distribution of changes in surface elevation was observed. Indeed, the steepest accumulation zone (average slope, 50%) of the glacier comprises several seracs and crevasse sites, so that a large fall of seracs will appear as a lowering of the elevation of the surface in its departure zone (Fig. 7, arrow A) but as a rise in the elevation of the surface in the deposit zone (Fig. 7, arrow B). This fall and rise tends to compensate for another and are balanced by the ice flow at a longer timescale (e.g. several months to several years, depending on the velocity of the glacier flow). This phenomenon is related to glacier dynamics and does not depend on climate factors. Above and beyond this dynamic phenomenon, part of the heterogeneity in the distribution of changes in surface elevation in the accumulation area can also result from wind erosion/transportation/deposition of snow, a process that is extremely difficult to quantify.

Fig. 7. Spatial distribution of changes in surface elevation (in m) between 1997 and 2009 on Glaciar Antisana 15α, 50-m interval contour lines are shown. Surface elevation lowering (from pale to dark orange) are concentrated along the ablation zone, whereas surface elevation rising (from pale to dark blue) or even slight surface lowering can be seen in the upper reaches. Black arrows show the zones with unstable seracs (a) and possible deposit areas (b).

5.2. Specific mass balance computed using the glaciological method

Using the glaciological method and Lliboutry's approach for interpolation/extrapolation, the average specific mass balance of Glaciar Antisana 15α over the 1995–2012 period was −0.64 m w.e. a⁻¹, which corresponds to a cumulative mass balance of −11.57 m w.e. Table 2 presents the entire data series.

Table 2. Glaciological time series of Glaciar Antisana 15α from 1995 to 2012. Columns (A) and (B) show the specific mass balance obtained with the glaciological method and with Lliboutry's approach before (A) and after (B) the adjustment using the geodetic method. The last two columns list changes in surface area (C) and the cumulative glacier snout retreat (D)

Between 1998 and 2009, the cumulative glaciological mass balance (B _glac.t ) obtained using Lliboutry's approach was −6.05 m w.e., which is 5.69 m w.e. less negative than the result obtained by simple linear interpolation (i.e. −11.74 m w.e.), described in Section 4.3.

5.3. Assessment of uncertainties in the glacier mass-balance series

Before comparing the different estimated mass balances, systematic errors must be considered (Zemp and others, 2013). Table 3 summarizes the corrected surface mass balances (glaciological and geodetic methods) for the period 1998–2009 and the corresponding random errors.

Table 3. Cumulative mass balances computed from the glaciological (B _glac.t ) and geodetic (B _geod.t ) methods for Glaciar Antisana 15α. Errors (±σ) corresponding to each method are presented. Data are missing (–) when the errors could not be quantified due to lack of information

5.3.1. Correction of the geodetic mass balance and assessment of uncertainties

According to Zemp and others (2013) and Huss (2013), the geodetic mass balance is prone to uncertainties related to the density assumption, the difference in survey dates and internal ablation/accumulation, which thus need to be taken into consideration:

• On a glacier, density is not constant over time and space (Huss, 2013). Accordingly, if we use an unchanged density conversion (dc) value, such as glacier ice (900 kg m⁻³), we will likely overestimate the mass balance (Zemp and others, 2013) because the loss of mass in the ablation zone is primarily ice, whereas in the accumulation zone, either snow/firn or ice can be gained or lost (Zemp and others, 2013). To assess the sensitivity of the geodetic mass balance to the density assumption, three tests were performed. Test I: mass changes involve only glacier ice over the entire surface with density =900 kg m⁻³ according to Sorge's law, and assuming steady-state conditions (Bader, 1954). This assumption has been widely used by the scientific community in similar studies (e.g. Thibert and others, 2008; Soruco and others, 2009). Test II: same as Test I but using a density value of 850 kg m⁻³ to account for changes in snow/firm over the glacier surface (e.g. Sapiano and others, 1998; Huss, 2013). Test III: A zonal density assessment, considering the density values of 600 kg m⁻³ above the ELA and 900 kg m⁻³ below the ELA (e.g. Kääb and others, 2012; Vincent and others, 2013).

  The cumulative geodetic mass balances given by the three tests I, II and III, were −1.99, −1.88 and −2.20 m w.e., respectively. The maximum difference between the three tests was low (~0.31 m w.e.). Because the geodetic balance given by Test I is very close to the average of the three tests, we used a density value of 900 kg m⁻³.

• The geodetic mass balance has to be adjusted when the dates of the aerial surveys do not match the dates of the field measurements (Cox and March, 2004; Thibert and others, 2005); this is called the survey difference (sd). In Ecuador, the annual mass balance was computed from 1 January to 31 December. However, the dates of the aerial photographs taken in 1997 (3 August 1997) and 2009 (13 September 2013) did not correspond with the established hydrological year. Consequently, mass-balance processes occurring at the glacier surface between the date of the aerial survey and the end of the year had to be included. It is possible to use precipitation data and a degree-day model to assess the accumulation and ablation processes. However, because the use of a degree-day model for tropical glaciers has been called into question (Sicart and others, 2008), we used a simpler approach. We computed the average monthly mass balance over the period between the two aerial surveys by dividing the geodetic mass balance by the number of months. This approach assumes constant monthly ablation/accumulation rates at the glacier surface, which is acceptable in the case of an equatorial glacier (Francou and others, 2004). As a result, mass-balance correction values were 0.07 and 0.05 m w.e. for 1997 and 2009, respectively, leading to the total cumulative mass balance from 1998/01/01 to 2009/12/31.

• Internal ablation (A _i ) is the amount of water lost that is not linked with climate conditions but resulting from three sources:

  Geothermal heat fluxes at the glacier/bedrock interface. The Antisana Ice Cap is located over an active volcano. However, little information is available about the volcano's activity over the past 400 years (personal communication from ML Hall, 2014), and there is no evidence for thermal activity or a local decrease in ice due to hot streams on the glaciers and the surrounding terrains (personal communication from P Ramon, 2014). Because the geothermal heat fluxes are not measured on Antisana Volcano, we considered an average geothermal heat flux of ~73 ± 10 mW m⁻² (Davies and Davies, 2010). Using the specific latent heat of fusion, basal melting gives 0.0069 m w.e.  a⁻¹.

  Heat conversion of the potential energy loss from ice motion. This component was estimated using the approach proposed by March and Trabant (1997). Considering an average surface ice velocity of 45 m a⁻¹ (±40%), we assumed the same velocity over the surface of the entire glacier. Thus, the vertical component of the glacier surface displacement (ΔZ) is about ~32 m a⁻¹. Using an average ice thickness of ~50 m (unpublished GPR measurements performed in 2000), the ice density value, the acceleration of gravity and the specific latent heat of fusion, we estimated an annual melting rate of 0.04 m w.e. a⁻¹.

  Heat conversion of the potential energy loss from water flowing through and beneath the glacier. These parameters are difficult to quantify accurately due to the lack of data, and were consequently not taken into account in the quantification of the glacier-wide mass balance.

• Internal accumulation (C _i ) results from refreezing of percolated melt water or the freezing of capillary-trapped water in the cold snow or firn layers in the accumulation zone (Thibert and others, 2008; Cuffey and Paterson, 2010). We assumed internal accumulation to be limited because Glaciar Antisana 15α is a temperate glacier (Francou and others, 2004) for which accumulated snowfalls can be consider more important for computing the mass balance. At very high elevations (>5500 m a.s.l.) where cold firn may exist, melt water from the surface could refreeze in the firn. However, on one hand, melting is negligible at such elevations, and ablation, when it occurs, is mainly through sublimation, which, according to Favier and others (2004) is likely to be low. On the other hand, the presence of ice lenses in the snow pack, resulting from refreezing of the surface melt in the firn, was taken into consideration in the measurements of density in the snow pits and hence in the measurements of net accumulation.

Finally, the corrected geodetic mass balance (B _geod.corr.t ) of Glaciar Antisana 15α between 1998 and 2009 was −1.39 m w.e., corresponding to a slightly negative average annual mass balance of −0.12 m w.e.  a⁻¹.

Following Zemp and others (2013), the random error in the geodetic mass balance ( ${\sigma _{{B_{{\rm geod}.t}}}}$ ) combines the remaining elevation error between two DEMs (σ _Δ.DEM.t ) and the spatial autocorrelation error in the elevation difference ( ${\sigma _{{B_{{\rm geod}.{\rm autocorr}.t}}}}$ ). These values were estimated on the stable non-glacierized area following the approach proposed by Rolstad and others (2009).

We considered that the random error in the density assumption ( ${\sigma _{{B_{{\rm geod}.{\rm dc}.t}}}}$ ) is the maximum difference given by the density tests (see above). The random error related to the difference between surveys ( ${\sigma _{{B_{{\rm geod}.{\rm sd}.t}}}}$ ) accounts for the random error in the geodetic mass balance distributed equally over the study period and attributed to the months between the date of the aerial survey and the end of the hydrological year.

Regarding random errors associated with basal ablation ( ${\sigma _{{B_{{\rm int}.{\rm bas}.t}}}}$ ), Davies and Davies (2010) proposed a range of 5–10 mW m⁻² as an estimated error in the geothermal heat flux, which represents ±20% of the value concerned. In addition, random errors in internal ablation ( ${\sigma _{{B_{{\rm int}.{\rm abl}.t}}}}$ ) must include uncertainties on ice flow velocity (ranging from 25 to 65 m a⁻¹ for Glaciar Antisana 15α) and the uncertainties on ice thickness estimated from the GPR measurements performed on the glacier, which were estimated to be 5 m.

Finally, the overall random error of the corrected geodetic mass balance ( ${\sigma _{{B_{{\rm geod}.{\rm corr}.t}}}}$ ) was obtained using the expression (Zemp and others, 2013):

(4) $${\sigma _{{B_{{\rm geod}.{\rm corr}.t}}}} = \sqrt {\sigma _{{B_{{\rm geod}.t}}}^2 + \sigma _{{B_{{\rm geod}.{\rm dc}.t}}}^2 + \sigma _{{B_{{\rm geod}.{\rm sd}.t}}}^2 + \sigma _{{B_{{\rm int}.{\rm abl}.t}}}^2 + \sigma _{{B_{{\rm int}.{\rm bas}.t}}}^2}. $$

5.3.2. Assessment of uncertainties on glaciological mass balance

We considered three main sources of uncertainties:

1. (1) At the point location ( ${\sigma _{{B_{{\rm glac}.{\rm point}.t}}}}$ ) we used the classical propagation error law considering the uncertainties in ablation and accumulation measurements suggested by Gerbaux and others (2005), the total number of point observations in each elevation range of the glacier and the total number of years in the study period.

2. (2) The random error due to the spatial integration of the surface mass balance ( ${\sigma _{{B_{{\rm glac}.{\rm spatial}.t}}}}$ ) combines two components: (i) the contribution of each elevation range to the glacier-wide surface mass-balance quantification following the analysis of variance proposed by Soruco and others (2009) and (ii) the SD of the remaining difference between the observations and the results of Lliboutry's approach.

3. (3) Uncertainties associated with changes in the surface area were substantially reduced by considering the annual measurements of the outline of the glacier and the changes in hypsometry (Section 4.3). The random error in the surface reference ( ${\sigma _{{B_{{\rm glac}.{\rm ref}.t}}}}$ ) was computed from the random errors associated with the annual glacier outline measured with DGPS and with changes in elevation quantified by photogrammetry.

Because most tropical glaciers are temperate (Francou and others, 1995) superimposed ice is not significant (Cuffey and Paterson, 2010), and because the flux divergence is not relevant in the computation of the glacier-wide mass balance (Cuffey and Paterson, 2010), errors related to these components can be disregarded.

Despite the better agreement of the cumulative mass balance (computed from the glaciological method using Lliboutry's approach) with the geodetic mass balance, a significant discrepancy persisted. Most of this difference could be attributed to systematic errors related to accumulation measurements (snow height and snow density profile). Such systematic errors are discussed in detail in Section 6 below.

5.4. Comparison and adjustment of the glacier mass-balance series

Comparing the cumulative corrected geodetic (−1.39 m w.e.) and glaciological (−6.05 m w.e.) mass balances over the period 1998–2009 gave a difference of 4.66 m w.e. (Fig. 8). This difference is greater than the calculated uncertainties. Assuming a normal distribution of the discrepancy (Δ) the null hypothesis (Ho:Δ_{(Bglac.t−Bgeod.t)} = 0) can be stated. Considering a confidence level of 95 and 90%, Ho was accepted, but the β-probability of 68 and 56% respectively shows that some systematic differences remain between the two series and that there is a good chance of having accepted a false statement. As a consequence, according to the β-probability test, we considered that the glaciological mass balance needed to be adjusted.

Fig. 8. (a) Specific annual mass balance for Glaciar Antisana 15α adjusted with the geodetic mass balance. (b) Cumulative glaciological mass balance computed using a simple linear interpolation (red triangles and dashed line), using Lliboutry's approach (green triangles and dashed line), and after calibration (black dots and line) using the geodetic mass balance (gray triangles). Uncertainties associated with each method are shown (vertical bars and light gray-shaded area).

To this end, the average annual mass balance was computed for the period of time common to the two methods (1998–2009). Then the glaciological mass-balance time series was centered by removing the 1998–2009 average from each annual value. Lastly, the average annual geodetic mass balance was added to each year of the centered glaciological mass-balance time series to obtain the adjusted glaciological annual mass balance. Finally, the glaciological mass-balance time series was adjusted for the whole period, 1995–2012, assuming the same bias observed in the geodetic mass balance over the 1998–2009 period (Table 2).

Furthermore, due to the absence of seasonality in the Ecuadorian Andes, it was not necessary to distinguish between winter and summer mass balances. Consequently, the calibration of the seasonal mass balances is not applicable here.

6. DISCUSSION

6.1. From glaciological observations to estimation of the mass balance

It should be kept in mind that conducting glaciological observations at Antisana Ice Cap, just as on the other ice-covered volcanoes in Ecuador, is challenging due to the complex topography of the surface of the glaciers. As mentioned above, direct measurements are not made on 60% of the surface area of the glacier. This can be assumed to contribute to uncertainties in the computation of the specific mass balance of the Glaciar Antisana 15α and may consequently lead to unreliable quantification of other glaciological variables such as the ELA and AAR.

We showed that using an adaptation of Lliboutry's approach to compute the specific mass-balance on Glaciar Antisana 15α, the difference between the cumulative glaciological and geodetic mass balances for the period 1998–2009 was 4.66 m w.e. In fact, the interpolation method applied in previous studies for the estimation of the glaciological mass balance on this glacier was the source of 55% of the total discrepancy. Unlike the linear interpolation method, Lliboutry's approach allowed us to respect the profile of the surface mass-balance gradient with elevation. This profile was obtained considering all available measurements and we can assume this profile is similar from one year to another. The results for each year obtained with the two approaches are presented in supplementary materials. For instance, if the mass-balance measurements are missing between 5000 and 5350 m a.s.l., as was the case in several years (supplementary materials), Lliboutry's approach makes it possible to estimate the surface mass balance for this missing data at each elevation range without affecting the shape of the mass-balance profile with elevation.

A test was performed to validate the performance of the statistical approach, using the year 1998, for which numerous mass-balance measurements are available and have been carried out at almost every elevation range. First, the glacier-wide surface mass balance was computed using all the measurements with the glaciological method using the linear interpolation and Lliboutry's approach, giving respectively −1.43 and −1.39 m w.e. a⁻¹. Second, the measurements between 5050 and 5450 m a.s.l. were removed from the dataset, and the glacier-wide surface mass balance was recalculated using the two approaches, giving respectively −2.03 and −1.45 m w.e. a⁻¹. This test clearly showed that the glacier-wide surface mass balance depends, to a great extent, on the amount and distribution of available measurements in the case of the linear interpolation, whereas Llibroutry's approach is less sensitive and allows more reliable quantification even when measurements are missing for several elevation ranges.

6.2. Uncertainties in estimated accumulation

As previously mentioned, on Glaciar Antisana 15α, surface accumulation measurements are made in two snow pits/cores with additional snow probing in between (Fig. 1a). Due to the almost complete absence of seasonality in precipitation, obtaining accurate measurements is challenging. Indeed, the seasonal layer that distinguishes one annual cycle from the previous is almost invisible in most years and may lead to uncertainties when estimating accumulation. Accumulation may be underestimated, leading to a too negative specific glacier mass balance. We tested the sensitivity of the mass balance to the accumulation value by assuming the amount of the annual precipitation measured on the glacier foreland instead of the accumulation values measured on the glacier. This resulted in a cumulative specific mass balance of −6.0 m w.e. (−0.49 m w.e. a⁻¹), which did not reduce the discrepancy between the cumulative glaciological and geodetic mass balances. This is logical because of a possible rain gauge under-catch and an increase in precipitation with elevation (see below).

If we consider that the remaining difference between the cumulative geodetic and glaciological mass balances is related to the inaccurate quantification of accumulation resulting in underestimation of the accumulation amount, a 70% increase in annual net accumulation will be necessary to balance the cumulative geodetic and glaciological mass balances.

This hypothesis is in good agreement with the ice flux estimate made from the kinematic method (Section 4.4). Indeed, we also calculated the ice flux at the level of the same cross section using the average annual surface mass balance over the entire study period. Because over the study period the average elevation difference above the selected cross section is equal to 0.18 m (close to 0), we neglected the dynamics changes and we considered the ice flux to be similar to the surface mass balance in the accumulation zone. Finally, an average annual mass balance above 5000 m a.s.l. of 1.07 m w.e. a⁻¹ is needed to balance the two ice flux estimates, which is 60% higher than the mass balance measured in the upper part of the glacier.

Such an increase in the net accumulation implies that precipitation on the glacier is underestimated and is higher than the precipitation measured in the glacier foreland. Former studies conducted in Ecuador or in Bolivia have already discussed the problem of precipitation captured by the rain gauges located on the glacier forelands, mostly related to wind and relief effects. Indeed, Wagnon and others (2009) suggested that due to the strong easterly winds, the deficit in precipitation captured by the rain gauges on Antisana 15α foreland could be as high as 50%. On Glaciar Zongo, precipitation deficits ranging from 30% to 50% were mentioned by Soruco and others (2009).

Laraque and others (2007) showed the strong spatial variability in precipitation between the western and the eastern sides of the Cordillera Central in Ecuador. In fact, the precipitation regimes are very complex in this region: a bimodal regime with two maxima in March–May and November and a minimum in June–October have been documented in the western side (leeward) of Antisana Volcano (Francou and others, 2000). In contrast, a unimodal regime with a maximum in June has been observed on the eastern Amazonian side (windward) (Laraque and others, 2007; Espinoza and others, 2009). Antisana Volcano is thus located in a transition zone between the two different regimes. Precipitation data from meteorological stations are very scarce in the region. Table 4 presents annual precipitation recorded at the seven available locations in the close vicinity of the volcano (Fig. 1b). It will be noted that precipitation decreases strongly from east to west and with increasing elevation. However, based on data from the available stations, the two effects are almost impossible to distinguish and more data are needed to clearly understand the impact of each variable. On the other hand, the vertical gradient of precipitation appears to be a critical factor. Although a negative gradient is apparent in the data from the seven stations, the comparison of the measurements in the pluviometers located on the glacier foreland (~4700 m a.s.l.) with the estimated accumulation values in the upper reaches of the glacier required to balance the geodetic and glaciological balance, showed a positive vertical gradient of ~40 mm (100 m⁻¹), i.e. precipitation increases with elevation on the glacier.

Table 4. Annual precipitation measured at meteorological stations located close to the volcano (between 1973 and 2013 depending on the station concerned)

In addition, the amount and distribution of precipitation are strongly affected by local conditions such as topography and wind exposure. As an example, in June–September when the winds are strong (Favier and others, 2004) the easterlies transport convective clouds that cling to the mountain and lead to condensation and then heavy precipitation on the eastern slopes of the volcano (Pourrut, 1994; Laraque and others, 2007) due to the cooling effect of adiabatic expansion (Vuille and others, 2000). This moisture and related precipitation may reach the summit and the upper western slopes in the accumulation zone and result in low-intensity precipitation and/or frost deposition, whereas clear-sky conditions may prevail on the foreland of western glacier tongues. The impacts of such meteorological conditions on the accumulation are poorly documented and not yet well understood.

One way to increase our understanding of the accumulation processes would be to perform accumulation measurements several times a year using accumulation stakes in the few flat sites above 5200 m a.s.l. We suggest considering periods of the year when wind speed and moisture are closely related. A first measurement could be made in May at the end of period of low winds and important cloud cover and thus precipitation; a second could be made at the end of October to document the processes that occurred during the very windy period with potentially high frost deposition in the upper reaches of the volcano; and finally a third at the end of the hydrological year would document the processes that occurred during the second precipitation maximum.

7. CONCLUSION

Glaciological mass-balance measurements over the period 1995–2012 were reanalyzed using the geodetic mass balance on Glaciar Antisana 15α. This glacier is considered as a ‘benchmark glacier’ for the inner tropics in South America. Once the geodetic and glaciological mass balances were computed and the uncertainty assessment was performed, the comparison between the results given by the two methods showed that the surface mass balance was overestimated with the glaciological method. It is worth emphasizing that:

• The method used to extrapolate the glaciological data in the zones in which no measurements are available plays an important role in quantifying the specific mass balance. We adapted the statistical approach developed by Lliboutry (1974) to interpolate/extrapolate the data and quantify the specific annual mass balance. Indeed, applying Lliboutry's approach, the difference between the glaciological and geodetic mass balances has reduced by ~55%.

• Because of the absence of precipitation seasonality, a seasonal surface layer that makes it possible to distinguish between two hydrological years is rarely visible. As a result, measurements of accumulation may be highly inaccurate, which will have a strong impact on the specific mass balance. Assuming that the remaining difference between the geodetic and glaciological mass balance computed applying Lliboutry's approach is related to the annual accumulation measurements, an increase in 70% of the net accumulation will be needed to balance the two estimations. This hypothesis is supported by an ice flux estimate in the vicinity of the ELA from GPR and surface flow velocity measurements. Other authors have shown that precipitation measured at rain gauges located on the glacier foreland may be underestimated by ~50%. However, future research is required to increase our knowledge of accumulation processes, especially the role of the physical factors that control snowfall on the glacier (atmospheric temperatures, water vapor content and wind drift).

• The cumulative mass balance of Glaciar Antisana 15α over the period 1995–2012 revealed a trend of −0.25 m w.e. a⁻¹ with some balanced or even positive years (e.g. the mean annual mass balance between 1999 and 2000 was 1.14 m w.e. a⁻¹ resulting in an immediate ~3% increase in the surface area of the glacier). Over the period 1995–2012, the glacier mass loss is related to a 15% decrease in surface area and a 235 m retreat of the glacier snout.

• Our results show that Glaciar Antisana 15α was subject to a less pronounced mass loss during the past decade than the other glaciers in other mountain ranges, such as the Andes (e.g. −0.70 ± 0.10 m w.e. a⁻¹; Rabatel and others, 2013) or the Alps (e.g. −1.00 ± 0.04 m w.e.  a⁻¹; Huss, 2012).