Simulations of changes to Glaciar Zongo, Bolivia (16° S), over the 21st century using a 3-D full-Stokes model and CMIP5 climate projections

Abstract

Bolivian glaciers are an essential source of fresh water for the Altiplano, and any changes they may undergo in the near future due to ongoing climate change are of particular concern. Glaciar Zongo, Bolivia, located near the administrative capital La Paz, has been extensively monitored by the GLACIOCLIM observatory in the last two decades. Here we model the glacier dynamics using the 3-D full-Stokes model Elmer/Ice. The model was calibrated and validated over a recent period (1997–2010) using four independent datasets: available observations of surface velocities and surface mass balance were used for calibration, and changes in surface elevation and retreat of the glacier front were used for validation. Over the validation period, model outputs are in good agreement with observations (differences less than a small percentage). The future surface mass balance is assumed to depend on the equilibrium-line altitude (ELA) and temperature changes through the sensitivity of ELA to temperature. The model was then forced for the 21st century using temperature changes projected by nine Coupled Model Intercomparison Project phase 5 (CMIP5) models. Here we give results for three different representative concentration pathways (RCPs). The intermediate scenario RCP6.0 led to 69 ± 7% volume loss by 2100, while the two extreme scenarios, RCP2.6 and RCP8.5, led to 40 ± 7% and 89 ± 4% loss of volume, respectively.

1. Introduction

The Bolivian Cordillera Oriental contains more than 1800 glaciers, covering a total area of ~450 km² (Soruco, 2008). Water resources from these glaciers are extremely important for the Bolivian Altiplano, where the administrative capital, La Paz, is located. Runoff from glaciers is an essential source of fresh water for direct consumption by the inhabitants, as well as for irrigation water for agriculture and water for the generation of hydropower. Soruco (2008) and Soruco and others (2015) determined that glacier meltwater contributes ∼15% of the water resources of La Paz annually, and up to 27% in the dry season.

Several studies have shown that tropical glaciers are particularly sensitive to climate change because tropical climate conditions maintain the lower reaches of the glaciers under ablation conditions all year round, so small changes in temperature and precipitation can have a major impact on their mass balance (Francou and others, 2003; Vuille and others, 2008; Rabatel and others, 2013). All glaciers in the tropical Andes have retreated significantly since the Little Ice Age maximum (Rabatel and others, 2008; Jomelli and others, 2009). In a recent review of glaciological surveys in the tropical Andes, Rabatel and others (2013) showed that in recent decades glacier shrinkage has accelerated with changing climate conditions and particularly with increasing temperatures. In this context, the question of the future of glaciers in the Andes, particularly in the Bolivian cordilleras, is extremely important for scientists and society.

A wide range of future climate projections for the next century, under different forcing scenarios derived from different climate models, is available from the Coupled Model Intercomparison Project phase 5 (CMIP5) (http://www.cmip-pcmdi.llnl.gov). These climate model results provide a solid basis for studying the response of glaciers to ongoing climate change.

To quantify future changes in glacier volume, the surface mass balance needs to be estimated, and ice dynamics have to be accounted for by transferring mass from accumulation to ablation zones. Approaches used to estimate mass balances in previous studies range from the sensitivity of the mass balance to temperature (e.g. Vincent and others, 2014) to the spatial distribution of surface mass balance calculated using a model driven by climate variables (e.g. Huss and others, 2010), despite considerable climateforcing uncertainties. We preferred to use a simple approach to parameterize changes to the surface mass balance. In tropical conditions, sublimation is too significant to use a simple degree-day model (Sicart and others, 2008), and the use of a full surface energy-balance model remains questionable because some of the input parameters are not accurately reproduced by the climate models. Different approaches to ice dynamics have been used in previous studies on mountain glaciers. Empirical approaches do not aim to reproduce glacier dynamics and fluctuations with high precision or in great detail and can thus rely on parameterization of varying complexity (Oerlemans and others, 1998; Vincent and others, 2000, 2014; Sugiyama and others, 2007; Huss and others, 2010; Carlson and others, 2014). Physical approaches rely on the use of numerical ice-flow models to solve the equations describing the flow of ice. Numerous ice-flow models of varying degrees of complexity have been developed and applied to both real and synthetic situations (Bindschadler, 1982; Hubbard and others, 1998; Le Meur and Vincent, 2003; Schäfer and Le Meur, 2007; Gudmundsson, 1999; Le Meur and others, 2007; Jouvet and others, 2008). One drawback of these approaches is that they require a large number of in situ observations to calibrate and validate the model.

In the present study, we used a physical approach. Thanks to the increase in computational power, solving the full set of mechanical equations that govern the flow of ice in glaciers or ice sheets without approximations (i.e. the full-Stokes equations) does not limit potential applications. Here we use the 3-D full-Stokes ice-flow model Elmer/Ice (Gagliardini and others, 2013) for an application on a mountain glacier in the tropical Andes, i.e. Glaciar Zongo, Bolivia.

Glaciar Zongo is a part of the GLACIOCLIM observatory (http://www-lgge.ujf-grenoble.fr/ServiceObs/SiteWebAndes/index.htm) and has been monitored continuously for more than two decades. Available observations include topography (surface and bedrock), changes in surface elevation, the spatial and temporal distribution of surface mass balance and surface velocities.

All these observations are described in Section 2. The model set-up, climate forcing and surface mass-balance parameterization are detailed in Sections 3 and 4. Finally, the calibration/validation of the model over the period 1997–2010 and projections for the 21st century are described and discussed in Section 5.

2. Study Area and Glaciological Data

Glaciar Zongo is located in the Bolivian Cordillera Real, about 30 km from La Paz (Fig. 1). The glacier covers an area of 1.9 km² and flows on the southeastern side of the Huayna Potosi summit from 6000 to 4900 m a.s.l. An extensive network of glaciological, hydrological and meteorological measurements on Glaciar Zongo has been developed since 1991, making it the glacier with the longest measurement series in the tropics (Rabatel and others, 2013). The observations used in this work to calibrate and validate the glacier model include glacier topography (bedrock topography and surface elevation), changes in surface elevations, surface mass balances and surface velocities.

Fig. 1. Glaciar Zongo, Bolivia. The inset map shows the location in Bolivia. Glacier outlines in 1983, 1997 and 2006 are shown, as well as the glaciological measurements network and the 2012 IPR thickness profiles (adapted from Soruco and others, 2009).

Digital elevation models (DEMs) at a resolution of 25 m were computed from stereo-pairs of aerial photographs for 1997 and 2006 (Soruco and others, 2009), and allowed accurate estimation of changes in the surface elevation over the entire surface area of the glacier during the period 1997–2006. The limits of the glacier, and especially the position of the front, were calculated using the same aerial photographs and differential GPS measurements.

A surface DEM was used to calculate bedrock topography in ice-free areas. In the safely accessible parts of the glacier, the ice thickness was measured in 2012 using icepenetrating radar (IPR). Our IPR is a geophysical instrument specially designed by the Canadian company Blue System Integration Ltd in collaboration with glaciologists to measure the thickness of glacier ice (Mingo and Flowers, 2010). It comprises a pair of transmitting and receiving 10 MHz antennas that allow continuous acquisition, and is georeferenced with a GPS receiver. Several profiles were made on the glacier tongue, which revealed a maximum thickness of 120 m at the foot of the seracs zone. Two profiles made in the upper part of the glacier showed a maximum thickness of 70 m (Fig. 1). In the seracs zone and in the highly crevassed upper reaches of the glacier, ice thickness was estimated taking into account the surface slope of the glacier and checking for overall agreement with the closest radar measurements taken above and below the seracs zone.

The surface mass balance of Glaciar Zongo has been measured since September 1991 using the glaciological method (Rabatel and others, 2012). The annual mass balance was obtained from measurements made at 21 stakes in the ablation area and in three snow pits/cores located in the accumulation area (Fig. 1). Because of the inaccessibility of some areas, including the steep seracs zone in the middle of the glacier, which represents about one-third of the glacier area, measurements cannot be made over the entire glacier surface. For the period 1997–2006, the annual mass-balance time series obtained using the glaciological method was compared and calibrated against the mass balance computed from the surface DEMs using the geodetic method (Soruco and others, 2009).

Velocities were measured on the network of stakes located in the ablation zone of the glacier (below 5200 m a.s.l.). Mean annual surface velocities measured at each stake were computed from the position of the stakes that have been measured precisely every year since 1991 using geodetic techniques (i.e. theodolite or differential GPS). In the upper and middle parts of the glacier, the flow was estimated using the movements of crevasses and seracs, which are easy to identify on high-resolution satellite images. This method is less accurate than geodetic surveys (accuracy in meters rather than in centimeters) but is sufficient in our particular case, as in these inaccessible areas the surface velocities are between 40 and 50 m a^–1.

3. Ice Flow Model

Glaciar Zongo flow dynamics were modeled using the parallel finite-element model Elmer/Ice (Gagliardini and others, 2013). This model was originally designed to solve complex ice-flow problems for ice sheets and has since been applied to model mountain glacier flow in several studies (e.g. Zwinger and others, 2007; Jay-Allemand and others, 2011; Adhikari and Marshall, 2013; Zhao and others, 2014).

3.1. Mesh

Because of the small aspect ratio of the glacier, the three-dimensional (3-D) mesh was constructed in two steps. First we meshed a two-dimensional (2-D) footprint of our domain. In all the simulations, the horizontal position of our mesh nodes was fixed. As a consequence, in the upper part of the glacier, where fluctuations in the margins of the glacier are not expected, the model domain is given by the glacier margin measured in 1997. In the lower part of the glacier, the model domain was extended onto the glacier foreland to allow the glacier to advance freely. This footprint was meshed with unstructured triangular elements at a spatial resolution of ~25 m. This choice was a compromise between spatial resolution and computational cost. The 2-D mesh was then extruded vertically using six layers between the bedrock and the free surface. The vertical resolution of the bottom layer was set at approximately half the resolution of the top layer to better resolve vertical gradients near the bedrock. The resulting 3-D mesh contained 32 784 nodes.

The free surface elevation at the mesh nodes was obtained from a cubic spline interpolation (Haber and others, 2001) of the 1997 surface DEM. Owing to the lack of bedrock-elevation measurements for some parts of the glacier, bedrock elevation at the mesh nodes was obtained using natural-neighbor interpolation (Fan and others, 2005).

3.2. Field equations

The 3-D velocity field u =(u_x , u_y, u_z ) is the solution of the Stokes equations that expresses conservation of momentum (Eqn (1)) and the conservation of mass in an incompressible fluid (Eqn (2)):

(1)

(2)

where τ is the deviatoric stress tensor, P is the isotropic pressure, g = (0, 0, –9.81 m s^–2) is the gravity vector and ρ = 917 kg m^–3 is the density of the ice. We ignored variations in firn density in the accumulation area; this assumption is discussed in Section 5.3.

The glacier-free surface z _s (m) evolves with time following the kinematic equation

(3)

where b (m w.e.) is the surface mass balance. As the finite-element mesh cannot have a null thickness, a lower limit of 2 m above the bedrock elevation was applied to z _s in Eqn (3).

3.3. Ice-flow law

We used a viscous isotropic nonlinear flow law, Glen’s law (Glen, 1955), to link the deviatoric stress tensor τ to the strain-rate tensor D:

(4)

The ice effective viscosity μ in Eqn (4) is given by

(5)

where n is the Glen exponent, is the second invariant of the strain-rate tensor and A is a rheological parameter that depends on temperature T following an Arrhenius law:

(6)

where the activation energy Q is equal to 139 kJ mol^–1, the gas constant R is equal to 8.314 J K^–1 mol^–1 and the rate factor A ₀ is equal to 1.916 10³ Pa^–3 s^–1. Temperature measurements in an ice core drilled at 6300 m a.s.l. in the accumulation area of nearby Glaciar Illimani revealed a vertical temperature gradient in the ice and firn of ∼1°C (100 m)^–1 (Gilbert and others, 2010). For the sake of simplicity, given the limited thickness of Glaciar Zongo, we ignored the variation in temperature with depth and assumed that the temperature inside the glacier is equal to the surface temperature. The surface temperature was estimated from the mean air temperature measured at a meteorological station located at 5150 m a.s.l. on the glacier itself. We used a lapse rate of –0.55°C (100 m)^–1 (Sicart and others, 2011) to account for the altitude effect. The glacier temperature was limited to a maximum of 0°C, as the estimated mean air temperature can be positive on the lowest part of the glacier. This temperature distribution was kept constant in all the simulations; these assumptions are discussed in Section 5.3.

3.4. Boundary conditions

As the lower part of the glacier is temperate (Francou and others, 1995), a certain amount of sliding of the ice on its bed is to be expected. A linear friction law relating the basal shear stress τ_b to the basal velocity u _b was applied on the lower boundary:

(7)

The basal friction parameter β was adjusted by comparing modeled and observed surface velocities. Because of the scarcity of observations, we used a simple approach to quantify it by trying to reproduce the observed surface flow velocities. β was parameterized as a function of the bedrock elevation, z _b, only. For the highest part of the glacier (z > 5600 m a.s.l.), where not much basal sliding was expected, a constant value of 0.1 MPa was imposed for β. For the lower part, an initial guess, β _ini, was taken from simple approximations: according to the shallow-ice approximation (Greve and Blatter, 2009) the basal shear stress in Eqn (7) is computed as

(8)

where h is the thickness of the ice and α is the surface slope.

Following Eqn (7), we can then estimate β _ini as

(9)

where u _b is assumed to be the surface velocity measured at the stake. β _ini was set for the 12 elevation bands using the mean values obtained from Eqn (9) and then linearly interpolated at each mesh node. The mean differences between modeled and observed velocities at each stake from 1997 to 2006 are given for each elevation band in Table 1 (note that the accuracy of velocity measurements made using a differential GPS is a few centimeters). These differences ranged from 7 to 13 m a^–1, which was more than the measured velocities. To reduce these differences and obtain a better match between observed and modeled velocities, the friction parameter was manually adjusted using a trial-and-error method. The differences in surface flow velocities obtained with this adjusted friction parameter are listed in Table 1. These differences are lower than the values obtained with the initial distribution of the friction parameter (Fig. 2). The adjusted distribution of the friction parameter was kept constant in both current and future simulations. This assumption is discussed in Section 5.3.

Table 1. Differences between modeled and observed surface velocities (m a^–1) averaged for each elevation band, with the initial friction parameter distribution β _ini (Eqn (9)) and the adjusted friction parameter

Fig. 2. (a) Polynomial mass-balance function (dashed blue line) computed from 1997–2006 surface mass-balance measurements (the black dots show the average mass balance for the period 1997–2006 for each elevation band, and blue shading shows the standard deviation). (b) Mean surface velocities measured over the period 1997–2010 plotted against modeled velocities computed with βadjusted for the six elevation bands given in Table 1. Blue error bars represent spatial and temporal variability of velocity measurements. Red error bars correspond to spatial variability of simulated velocities.

For the lateral sides of the domain, two kinds of boundary conditions were applied: in the upper part, where the glacier extent is not expected to change, a null velocity was applied; in the lower part, a stress-free condition was applied as the model domain was extended compared with the margins of the glacier in 1997; this condition was also applied to non-glacierized areas and did not influence the results. The limit between the two parts is located at the average equilibrium-line altitude (ELA) in the period 1997– 2006, i.e. 5300 m a.s.l.

3.5. Relaxation

It is well known that changes in surface elevation show high and unphysical values at the beginning of prognostic simulations due to uncertainties in the initial model conditions (Seroussi and others, 2011). For that reason, we allowed a relaxation of the free surface before all prognostic simulations, in the same way that Gillet-Chaulet and others (2012) had done. The length of the relaxation was a compromise between leaving sufficient time to remove the largest transient effects but not too long, so that the free surface remains close to observations. We chose a relaxation period of 8 years. The surface mass balance b in Eqn (3) was kept constant over the relaxation period and was given by observations made in 1995. We chose these observations because the mass balance integrated over the whole glacier surface was nearly zero, resulting in only small changes in the glacier total volume during the relaxation period.

The glacier geometry at the end of the relaxation period was compared with the observed glacier geometry in 1997 (Fig. 3). The difference between the modeled and measured glacier areas was <1%. Changes in glacier thickness between the beginning and the end of the relaxation period were ~1 m in the lower part and <5 m in the upper part.

Fig. 3. Simulated thickness of Glaciar Zongo in 1997, 2006 and 2010. The glacier outline observed in 1997 is represented by a dotted line in each diagram. The front position observed in 2006 and 2010 is represented by bold lines.

4. Climate Forcing and Mass-Balance Parameterization

4.1. Climate forcing for the 21st century

Temperatures were taken from the CMIP5 results to force the relationship between temperature and the ELA used to quantify the surface annual mass-balance distribution over the 21st century. We selected nine global climate models and three representative concentration pathways (RCPs). Scenarios RCP2.6, RCP6.0 and RCP8.5 are based on radiative forcing of 2.6 W m^–2, 6.0 W m^–2 and 8.5 W m^–2, respectively. Scenario RCP2.6 is the most optimistic scenario in terms of greenhouse gas emissions, while scenario RCP8.5 is the most pessimistic. Scenario RCP6.0 is considered intermediate. For each of the nine climate models used, we considered the modeled annual average 2 m air temperature for the 21st century for the coordinates 15–18° S, 68–70° W, which cover the location of Glaciar Zongo. The modeled temperature data were converted into anomalies by subtracting the 1997–2006 average.

4.2. Mass-balance parameterization

For the calibration/validation period, the annual surface mass balance was derived from observations in 12 elevation bands from 4900 to 6000 m a.s.l.

For the future simulations, we proceeded in two steps:

• First, the vertical profile of surface mass balance computed from the available annual measurements averaged over the period 1997–2006 was fitted by a polynomial function to obtain a mass-balance function (Fig. 2a); this function gave a better representation of the observed data. Note that this function has been used for the validation period to test its efficiency compared with the use of direct surface annual mass-balance measurements, giving similar results.

• Second, we combined the ELA sensitivity to temperature with the temperature anomalies given by the climate models (Section 4.1) to obtain a projection of the annual ELA for the 21st century. Note that the sensitivity of Glaciar Zongo ELA to temperature was estimated by Lejeune (2009) to be 150 ± 30 m °C^–1. To link changes in ELA to changes in the spatial mass-balance distribution over the glacier surface, we assumed that the relationship based on the measurement period will continue over the 21st century. Then the mass-balance function was shifted according to the annual position of the ELA to obtain the annual spatial mass-balance distribution for all elevations of the glacier surface. Note that the changes in surface area and surface elevation of the glacier are taken into consideration for future simulations because the mass balance was computed for each point of the mesh at each time step and both the glacier surface area and the glacier surface elevation were updated at each time step.

5. Results and Discussion

5.1. Recent period (1997–2010)

Prognostic simulations were run for the period 1997–2010. Two independent datasets based on field measurements were used to constrain and validate the model. First, the model was forced and calibrated with available observations of topography, surface velocities and surface mass balance. The adjusted values of the basal friction parameter led to simulated surface velocities close to the field measurements (see Section 3.4). Then the model was validated over the recent period using different datasets from those used for calibration: the position of the front and the total volume loss.

Observed and modeled glacier outlines in 2006 and 2010 are compared in Figure 3. The retreat of the glacier front between 1997 and 2010 was well reproduced by the model. The difference between the simulated and measured glacierized areas was <1.5% in 2006 and <1% in 2010. The observed retreat of the front was 140 m between 1997 and 2006 and 50 m between 2006 and 2010. Results of the simulations were 125 m and 60 m, respectively. Modeled changes in surface elevation between 1997 and 2006 are shown in Figure 4 and compared with the measurements obtained by photogrammetry (Soruco and others, 2009).

Fig. 4. Changes in surface elevation between 1997 and 2006: (a) simulated by the model, (b) measured by photogrammetry (adapted from Soruco and others, 2009). (c) Difference between absolute model and observed surface elevation changes.

Measurements revealed a maximum shift from ~40 m for the glacier tongue down to ~15 m below the seracs zone. These model results are in good agreement in terms of both spatial distribution and magnitude. Indeed, the volume loss computed by photogrammetry was 1.10 × 10⁷ m³ while the simulated volume loss was 0.95 × 10⁷ m³. This 13% difference in volume loss is within the range of accuracy of the photogrammetric method.

5.2. Simulations of the 21st century

Simulations of the 21st century started in 2006. The state of the glacier in 2006 was projected by the simulations covering the calibration period (1997–2006), described in Section 5.1.

Volume projections for Glaciar Zongo with temperature anomalies obtained from the nine climate models for scenario RCP2.6 and for a sensitivity of the ELA of 150 m °C^–1 are shown in Figure 5. All model results showed similar trends, with a continuous decrease in glacier volume down to ~50% of its present value, until the period 2035– 60, followed by stabilization of the volume of the glacier until the end of the century. This stabilization reflects the temperature stabilization projected by the models with a time lag of about 10–15 years. The scatter of the volumes projected by the nine models is large and increases with time, with cumulative volume losses of 11–39% in 2030, 18–48% in 2050 and 32–59% in 2100.

Fig. 5. (a) Simulated Glaciar Zongo volumes from 1997 to 2100 using the RCP2.6 scenario. Results are given for the nine climate models used in this study and an ELA sensitivity of 150 m °C^–1. The horizontal dashed line indicates half the 2006 volume of Glaciar Zongo. (b) Temperature anomalies (compared with the 1997–2006 period) from 2006 to 2100 projected using the RCP2.6 scenario and the nine CMIP5 models.

Volume changes computed using temperature data from the three RCP scenarios are shown in Figure 6. The solid lines show the average of the nine-model forcing, and the shaded area shows the 1σ interval, where σ is the standard deviation of the nine model outputs. The three scenarios projected similar continuous mass loss until approximately 2035. As already mentioned, the change in volume projected by scenario RCP2.6 then remained stable for the rest of the century, whereas the two other scenarios projected a continuous decrease until the end of the century. As expected, RCP8.5 projected the highest rate of volume loss, with the glacier nearly disappearing in 2100. The volume losses projected by scenarios RCP2.6, RCP6.0 and RCP8.5 were respectively 27 ± 7%, 24 ± 7% and 29 ± 8% in 2030; 35 ± 8%, 34 ± 7% and 44 ± 7% in 2050; and 40 ± 7%, 69 ± 7% and 89 ± 4% in 2100.

Fig. 6. Simulated Glaciar Zongo volumes from 1997 to 2100 using scenarios RCP2.6 (purple), RCP6.0 (green) and RCP8.5 (orange). Solid lines show the average of the nine models, and the shaded area shows the ±1σ interval. The horizontal dashed line indicates half the current volume of Glaciar Zongo.

Lastly, we assessed the sensitivity of the volume changes computed by the model considering the sensitivity of the ELA to temperature changes. The average changes in volume obtained with the three RCP scenarios for ELA sensitivities to temperature of +120 m °C^–1, +150 m °C^–1 and +180 m °C^–1 are shown in Figure 7. As expected, the trends are similar and the volume loss increased with an increase in the sensitivity of the ELA. The average volume losses in 2100 obtained with the three values of the sensitivity were respectively 35%, 40%, 45% with RCP2.6; 60%, 69%, 77% with RCP6.0; and 82%, 89%, 94% with RCP8.5.

Fig. 7. Simulated Glaciar Zongo volumes from 1997 to 2100 using scenarios RCP2.6 (purple), RCP6.0 (green) and RCP8.5 (orange). Solid lines represent an ELA sensitivity of 150 m °C^–1. Dashed lines represent a sensitivity of 180 m °C^–1 and 120 m °C^–1. The horizontal dashed line indicates half the current volume of Glaciar Zongo.

5.3. Limitations of our study

While this study is among the first to use a full-Stokes model for a mountain glacier application, modeling ice dynamics still relies on three main assumptions: (i) basal sliding is parameterized using a linear friction law where the coefficient is a function of altitude only; (ii) the viscosity field was calibrated from current observations of the surface air temperature, and thermomechanical coupling is ignored; (iii) our glacier is composed only of ice, and variations in density are ignored.

Sensitivity tests for the first two assumptions were made for the validation period (1997–2010, Section 5.1). Variations of ±10% for the friction coefficient and the temperature led to variations in volume of <1% and <0.05%, respectively, in 2010.

In the case of a mountain glacier with no calving front and a surface mass balance parameterized as a function of altitude only, ice dynamics affect glacier volume only by altering the surface area of the glacier and free surface elevation. This explains the low sensitivity of the variations in volume to the ice dynamics parameters and justifies the use of pragmatic assumptions. Progress is underway to model these processes more accurately (e.g. Gilbert and others, 2014), and these processes will be incorporated in our simulations when the complexity of the coupling between ice dynamics and surface mass balance requires it.

Finally, uncertainties in the simulations of glacier changes for the coming decades mostly depend on climate changes and their impact on the surface mass-balance quantification. Our approach for surface mass-balance parameterization is simple but relies on direct surface mass-balance measurements. Using such a parameterization together with ELA shifts for future simulations may lead to uncertainties regarding interannual surface mass-balance variations. As shown in Figure 2, temporal mass-balance variability is higher in the lower reaches of the glacier than in the upper ones. Consequently, assuming a linear shift of the surface mass-balance profile leads to overestimated annual values for extreme positive/negative mass-balance years. Impacts of this approach are hardly quantifiable at secular scale, and in any case are assumed to be limited compared with the large scatter of temperature changes related to the different scenarios and climate models.

6. Conclusion

This paper describes the first application of a 3-D full-Stokes model on a mountain glacier in the tropical Andes to take into account the ice dynamics in ongoing and future glacier changes. Glaciar Zongo is a perfect candidate for such an application because the glacier has been continuously monitored since the early 1990s by the GLACIOCLIM observatory, therefore all the data needed to calibrate and validate the model are available. In simulating the future mass balance, our approach was based on the sensitivity of the ELA to temperature changes and assumed that the observed vertical profile of surface mass balance remains unchanged. Regarding glacier dynamics, some assumptions were made concerning basal sliding, the internal temperature of the ice and the density of the material as their impact is very limited in simulations of the glacier changes compared with the impact of changes in climate conditions and hence in the mass balance.

For the 1997–2010 validation period, the modeled surface area and changes in volume were in good agreement with observations. Consequently, the model was then used to simulate future changes in Glaciar Zongo during the 21st century using nine different climate models and three RCP scenarios. In the near future (2030), the variability between the different climate models was greater than the variability between the scenarios, and the volume loss ranged from 24 ± 7% to 29 ± 8%. The changes in volume computed with the three RCP scenarios differed significantly after 2050. Scenario RCP2.6 projected a stable volume in the second half of the century with a value corresponding to ~60% of its present volume, while the other two scenarios projected continuous and sustained volume loss over the whole century, with the near disappearance of the glacier projected in 2100 by the most pessimistic scenario, RCP8.5.

Future research could focus on (1) simulations of mass-balance changes using a complete surface energy-balance model; (2) considering a certain amount of snow/firn with different density values and rheological properties; (3) considering fracturing of the ice that would make it possible to represent specific zones of glaciers, like the seracs zone.