1. Introduction
The Bolivian Cordillera Oriental contains more than 1800 glaciers, covering a total area of ~450 km2 (Reference SorucoSoruco, 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. Reference SorucoSoruco (2008) and Reference SorucoSoruco 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 (Reference Francou, Vuille, Wagnon, Mendoza and SicartFrancou and others, 2003; Reference VuilleVuille and others, 2008; Reference RabatelRabatel and others, 2013). All glaciers in the tropical Andes have retreated significantly since the Little Ice Age maximum (Reference Rabatel, Francou, Jomelli, Naveau and GrancherRabatel and others, 2008; Reference Jomelli, Rabatel, Brunstein, Hoffmann and FrancouJomelli and others, 2009). In a recent review of glaciological surveys in the tropical Andes, Reference RabatelRabatel 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. Reference Vincent, Harter, Gilbert, Berthier and SixVincent and others, 2014) to the spatial distribution of surface mass balance calculated using a model driven by climate variables (e.g. Reference Huss, Jouvet, Farinotti and BauderHuss 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 (Reference Sicart, Hock and SixSicart 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 (Reference OerlemansOerlemans and others, 1998; Reference Vincent, Vallon, Reynaud and Le MeurVincent and others, 2000, Reference Vincent, Harter, Gilbert, Berthier and Six2014; Reference Sugiyama, Bauder, Zahno and FunkSugiyama and others, 2007; Reference Huss, Jouvet, Farinotti and BauderHuss and others, 2010; Reference CarlsonCarlson 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 (Reference BindschadlerBindschadler, 1982; Reference Hubbard, Blatter, Nienow, Mair and HubbardHubbard and others, 1998; Reference Le Meur and VincentLe Meur and Vincent, 2003; Reference Schäfer and Le MeurSchäfer and Le Meur, 2007; Reference GudmundssonGudmundsson, 1999; Reference Le Meur, Gerbaux, Schäfer and VincentLe Meur and others, 2007; Reference Jouvet, Picasso, Rappaz and BlatterJouvet 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 (Reference GagliardiniGagliardini 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 km2 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 (Reference RabatelRabatel 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.
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 Reference SorucoSoruco 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 (Reference SorucoSoruco 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 (Reference Mingo and FlowersMingo 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 (Reference RabatelRabatel 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 (Reference SorucoSoruco 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 (Reference GagliardiniGagliardini 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. Reference Zwinger, Greve, Gagliardini, Shiraiwa and LylyZwinger and others, 2007; Reference Jay-Allemand, Gillet-Chaulet, Gagliardini and NodetJay-Allemand and others, 2011; Reference Adhikari and MarshallAdhikari and Marshall, 2013; Reference ZhaoZhao 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 (Reference Haber, Zeilfelder, Davydov and SeidelHaber 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 (Reference Fan, Efrat, Koltun, Krishnan, Venkatasubramanian, Demetrescu, Sedgewick and TamassiaFan and others, 2005).
3.2. Field equations
The 3-D velocity field u =(ux , uy, uz ) 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)):
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
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 (Reference GlenGlen, 1955), to link the deviatoric stress tensor τ to the strain-rate tensor D:
The ice effective viscosity μ in Eqn (4) is given by
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:
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 0 is equal to 1.916 103 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 (Reference Gilbert, Wagnon, Vincent, Ginot and FunkGilbert 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 (Reference Sicart, Hock, Ribstein, Litt and RamirezSicart 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 (Reference Francou, Ribstein, Saravia and TiriauFrancou 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:
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 (Reference Greve and BlatterGreve and Blatter, 2009) the basal shear stress in Eqn (7) is computed as
where h is the thickness of the ice and α is the surface slope.
Following Eqn (7), we can then estimate β ini as
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.
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
(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 (Reference SeroussiSeroussi and others, 2011). For that reason, we allowed a relaxation of the free surface before all prognostic simulations, in the same way that Reference Gillet-ChauletGillet-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.
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 Reference LejeuneLejeune (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 (Reference SorucoSoruco and others, 2009).
Changes in surface elevation between 1997 and 2006: (a) simulated by the model, (b) measured by photogrammetry (adapted from Reference SorucoSoruco 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 × 107 m3 while the simulated volume loss was 0.95 × 107 m3. 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.
(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.
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.
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. Reference Gilbert, Gagliardini, Vincent and WagnonGilbert 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.
Acknowledgements
This study was funded by the French Institut de Recherche pour le Développement (IRD) through the Andean part of the French glacier observatory service, GLACIOCLIM (http://www-lgge.ujf-grenoble.fr/ServiceObs/SiteWebAndes/index.htm), and was conducted in the framework of the International Joint Laboratory GREAT-ICE, a joint initiative of the IRD and universities and institutions in Bolivia, Peru, Ecuador and Colombia. Contributing authors from the Laboratoire de Glaciologie et Géophysique de l’Environnement (LGGE) acknowledge the contribution of LABEX OSUG@2020, ANR grant No. ANR-10-LABX-56. We are grateful to all those who took part in field campaigns for mass-balance measurements, to Benjamin Lehmann (IRD) for his help with the radar field campaign, to Gerhard Krinner (LGGE) for providing the CMIP5 data for our region of interest, and to Laurent Mingo (Blue Ice Integration Ltd) for assistance with the IPR (http://www.radar.bluesystem.ca/). This work was granted access to the high-performance computing resources of CINES under the allocation 2014-016066 made by GENCI. Finally, we thank G. AËalgeirsdóttir (Scientific Editor) and two anonymous referees for constructive comments that helped to improve the paper.
