Next Article in Journal
Predicting Water Availability in Water Bodies under the Influence of Precipitation and Water Management Actions Using VAR/VECM/LSTM
Previous Article in Journal
Effects of Climate Conditions on TP Outsourcing in Lake Kinneret (Israel)
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Distributed Energy Balance Flux Modelling of Mass Balances in the Artesonraju Glacier and Discharge in the Basin of Artesoncocha, Cordillera Blanca, Peru

by
María Fernanda Lozano Gacha
* and
Manfred Koch
Department of Geohydraulics and Engineering Hydrology, University of Kassel, 34125 Kassel, Germany
*
Author to whom correspondence should be addressed.
Submission received: 22 July 2021 / Revised: 7 September 2021 / Accepted: 7 September 2021 / Published: 17 September 2021

Abstract

:
A distributed energy balance model (DEBAM) is applied to estimate the mass balance of the Artesonraju glacier in the Cordillera Blanca (CB), Peru, and to simulate the ensuing discharge into its respective basin, Artesoncocha. The energy balance model calibrations show that, by using seasonal albedos, reasonable results for mass balances and discharge can be obtained, as witnessed by annually aggregated Nash Sutcliffe coefficients (E) of 0.60–0.87 for discharge and of 0.58–0.71 for mass measurements carried out in the period 2004–2007. Mass losses between −1.42 and −0.45 m.w.e. are calculated for that period. The elevation line altitudes (ELAs), which lie between 5009 and 5050 m.a.s.l., are also well simulated, compared to those measured by the Unidad Glaciologica de Recursos Hídricos del Perú (UGRH). It is demonstrated that the net radiation which drives the energy balance and melting processes is mainly affected by the amount of reflected shortwave radiation from the different surfaces. Moreover, the longwave radiation sinks between 63 and 73% of solar radiative energy in the dry season. Further sensitivity studies indicate that the assumed threshold temperature T0 is crucial in mass balance simulations, as it determines the extension of areas with different albedos. An optimal T0 between 2.6 and 3.8 °C is deduced from these simulations.

1. Introduction

Mass balance changes of tropical glaciers have produced a retreat of the respective glaciated areas over the last 40 years [1,2,3,4]. This glacier retreat results from alterations of physical processes between the atmosphere and the glaciers that are strongly affected by anticipated changes in various climatic factors. Modifications of glacier dynamics are indicators of the global warming trend in tropical regions and generate synergic impacts on the water system as a whole [3].
The case of the Cordillera Blanca (CB) is crucial, as CB is the largest mountain range with tropical glaciers on Earth. Glacier waters from the CB contribute to the total annual discharge in the basin of the Santa River in Peru. Some studies indicate that between 16 and 20% [5] or even up to 30–45% [6] of the discharge in each glaciated catchment of the Cordillera Blanca comes from meltwater. As a matter of fact, this meltwater is an essential water reserve for many downstream users along the Santa River feeding, in particular, the large Chavimochic agricultural irrigation project near the coastal city of Trujillo, which is considered one of highest-revenue producers in the whole country of Peru. Thus, it is of no surprise that adverse changes in the CB glaciers may have harsh consequences for the economic well-being of Peru, all of which stirs ongoing interest in appropriate glacier studies in the CB.
Furthermore, the importance of Andean tropical glaciers resides in the continuous water provision during the hydrological year and in the regulation of the water discharge in each season. Thus, the tropical glaciers act as a water storage system that supplies water during the whole year. This storage is vital, especially in the dry season, when precipitation reduces drastically. Therefore, changes in seasonal patterns of glacier discharge can impact water availability, namely in the dry season, and the occurrence of flooding in the rainy season [3,7,8]. Moreover, changes in the glacier dynamics pose severe risks to downward communities, due to destabilization of the rocks and glacier surfaces. Avalanches from slides of ice, snowpack, and rock can cause disasters such as the one that occurred in Chamoli, Indian Himalaya [9].
Consequently, the accelerated loss of mass of tropical glaciers has created a need to determine the climatic variables that drive these changes. Some studies carried out on the Andean Tropics explain to what extent climate variables, such as temperature, precipitation, humidity, and solar radiation, influence the physical dynamics behind mass balance [10,11,12,13,14,15,16,17,18,19,20,21]. Similar studies have focused on determining the effects of climate phenomena, such as ENSO, on tropical glaciers [13,22,23] or those of climate change on mountain hydrology [24,25].
The generalized retreat of glaciers across the Andes suggests a common pattern of causality [26]. However, research studies have led to different conclusions about which climatic variables are the dominant driving factors of the energetic fluxes in tropical glaciers. For the studied cases and regions, the different influencing factors that explain the glacier response to climate variations seem to depend on the specific locations within the tropics [11]. In addition, differences in the driving factors of the energy balance, specifically of the ablation and the accumulation zones, were found [14].
Research on energy balance fluxes on Andean tropical glaciers is relatively recent. The most studied glaciers regarding their energy dynamics are Antisana in Ecuador [10,27]; Zongo [18,19,20,21,28], Charquini [29], and Chacaltaya [12] in Bolivia; and Shallap [14,22] and Artesonraju in Peru [15,30].
For most studies, there is the common conclusion that the net radiation, namely, its contribution from net shortwave radiation, is the primary source of energy for melting in the ablation zone, and here, the albedo as the major influencing parameter [3,10,11,14,22]. There are differences, though, concerning which parameters affect the albedo the most. In most northern tropic glaciers, the so-called inner tropic glaciers, the temperature plays a more decisive role in the glacier dynamics. The former regulates the liquid and solid phases of precipitation and, consequently, the snowfall areas. However, these statements are valid only in white glaciers, while debris-covered glaciers show different energy flux dynamics [3].
This influence is not so explicit in the glaciers located in the outer tropic with a more subtropical climate, such as the Zongo glacier in Bolivia. As a matter of fact, the annual melting is directly related to the distribution of the precipitation. Solid precipitation extends over the whole glacier, reducing melting on the glacier surface [11,19], and temperature is poorly correlated with seasonal mass balances [3]. However, it cannot completely be ruled out either, that the temperature is not a determining variable in outer tropics glaciers, as it appears to be the case for the Shallap glacier in Peru. Here, energy simulations suggest that the threshold temperature can also have an important influence on albedo, by controlling the solid–liquid phase of precipitation [14].
For southern tropical glaciers, some studies state that the effects of temperature change on the mass balances are produced via a changing humidity pattern [20]. The latter is a climatic factor that influences sublimation.
Sublimation constitutes a further relevant physical process in certain subtropical and tropical environments, e.g., Kilimanjaro’s glaciers, where sublimation dominates ablation [31]. However, other studies contradict this conclusion about the essentiality of sublimation as a driver of the energy budget on Andean tropical glaciers, for instance, Zongo glacier, where the role of sublimation is not as important as that of the longwave fluxes [19]. Moreover, other study results for tropical glaciers, namely Shallap, indicate that the latent flux and sensible heat offset each other [14,22].
Liquid precipitation rarely occurs on the Zongo glacier with an altitude range of 5000–6000 m.a.s.l. Therefore, for such high-altitude glaciers, it is the amount of snowfall and not the temperature that determines the snowline [11]. This was also observed by Francou [12], who attributes seasonal variations in the Chacaltaya Glacier in Bolivia to precipitation and humidity changes. However, the author found a correlation between reanalyzed temperatures on interannual timescales and the mass balance of Chacaltaya, which he explained by the interrelation of temperature with humidity [12].
Despite the different and, sometimes, contradictory results from the aforementioned studies, regarding seasonal climate drivers, especially in the glaciers of the outer tropics, retreat and mass losses are common problems affecting all tropical glaciers. Increments in annual temperature along the Andes show, evidently, a direct negative correlation with mass balances [3]. Nevertheless, because of the diverse results from the various glacier studies, one must put forward the idea that temperature increases due to global climate change may induce different mechanisms that can act directly or indirectly with local climatic factors that, in turn, eventually impact mass balances.
Based on the above statements and the manifold of conclusions of the numerous glacier studies, it is clear that further investigation is required to examine the spatial and temporal dynamics of energy fluxes and the role of other climatic variables in tropical glaciers. Thus, this theme constitutes an ample study field, as the state of research shows that many aspects are not yet clear, such as the impact of patterns of seasonal climate drivers on the behavior of glaciers located in the inner and outer tropics [3,11]. The present research attempts to contribute to a better understanding of the physics of tropical glaciers and the effects of global climate change on their development through modeling recent mass balance changes in the Artesonraju glacier using a novel distributed energy balance model (DEBAM). The latter class of glacier models, although—because of a need for more variant climate/energy data—more complex and intricate to use than the more commonly employed simple temperature index models, appears to better mimic mass balances in glaciers, because of a better representation of the numerous physical phenomena controlling a glacier’s accumulation and ablation processes.

2. Materials and Methods

2.1. Study Area and Climate Characteristics

The glacier studied in this study is the Artesonraju, which is located in the basin of Artesoncocha, part of the Santa River basin in the Cordillera Blanca (CB). As mentioned, this river is the major water resource of the coastal region around the city of Trujillo, south of which, having carved its way for several hundred km through the western Andes chain, the river flows into the Pacific. The glaciated area of the CB is 527.62 km2. As a matter of fact, the CB is the mountain chain with the largest extension of glaciers in the tropics. Located at 8°12′–10°01′ S and 77°40′ W in the district of Ancash, Peru, this mountain chain comprises a total of 775 glaciers [2]. Among these, more than 200 have a maximum elevation above 5000 m a.s.l, with 30 even above 6000 m a.s.l. [1].
In the Artesoncocha basin, there are three snow-covered mountain peaks: Artesonraju (north), Paron (northeast), and Piramide (south) (see Figure 1). The Artesoncocha lake is a body of water of glacial origin located at the foot of the Artesonraju mountain. This lake is fed by the water of a recently formed small lake, called Laguna Nueva Artesoncocha Alta (situated on the tongue of the Artesonraju glacier) and some of the waters coming from the Caraz snow-covered areas. Artesoncocha lake has an area of 6.77 km2 at an elevation of 4288 m.a.s.l. [32]. The snout of the Artesonraju glacier is located at 8°57′07″ S, 77°37′56″ W. The glacier has an area of 3.5 km2 and a defined tongue. The altitude of the glacier reaches from between 4685 and 4720 m.a.s.l. up to 5979 m.a.s.l. The glacier is predominantly oriented towards the southwest and has an average slope of 54% [32].
The climate in the CB is mainly influenced by the easterly winds from the Amazon basin, the westerly winds that come from the Pacific which, in turn, are affected by the Humboldt Ocean current, local orographic effects, and last but not least, by the El Niño Southern Oscillation (ENSO), wherefore the multi-year periodicity of the latter leads to corresponding—though dampened—climate alterations in the mountain chain. The Intertropical Convergence Zone (ITCZ) produces tropical characteristics and some subtropical climate features in the dry season. Thus, only two main seasons, wet and dry, occur in the CB. The dry season is characterized by lower humidity and less precipitation than the wet season. However, the wet season can be characterized by two phases, the first one, the beginning of the rainy season, September and October (SO), when the precipitation starts with lower intensity, and the second one, from November to March (NDJFM), which comprises the core of the rainy season. The subsequent dry season lasts from April to August (AMJJA).
The seasonality of the Artesonraju glacier manifests itself in slight changes in temperature over the months of the year. In fact, tropical glaciers are known for their low variation of daily temperature during the year. Nevertheless, minor changes are noticeable in the records of Artesonraju. In particular, the minimum and maximum monthly mean temperatures have differences of 1.1 ± 0.2 °C, according to the data in the gauge station of Artesonraju at (4811 m.a.s.l.), whereas the precipitation follows the general yearly variation mentioned above. As such, the seasonal pattern of this glacier fits the one observed for the outer tropics quite well, i.e., subtropical influence in the dry season and tropical behavior in the wet season [33]. The mean wind speed in Artesonraju is larger in January and February and between August and November, and decreases between May and July, which correspond to the driest months [34].

2.2. Data

The climate time series of daily temperature, precipitation, humidity, wind, and energy fluxes of incoming shortwave radiation, shortwave reflected radiation, incoming longwave radiation, and outgoing longwave radiation were used for the timespan between 24 August 2004 and 28 August 2007. Moreover, mass balance- and discharge measurements were available from the Peruvian Autoridad Nacional del Agua (ANA) for the same period. This data is jointly managed with the Institute de Recherche pour le Dévelopement (IRD). In addition, the World Glacier Monitoring Service (WGMS) provided records of mass balance estimations [35]. The location of the gauge stations is shown in Figure 2, and details of the data characteristics are listed in Table 1.
The mass balance records available between 9 September 2003 and 18 September 2008 (i.e., for each hydrological year in that time period) are based on field measurements from the net of stakes and pits operated by UGRH and IRD. These measurements were taken in lapses varying between 20 to 68 days, with the number of stakes ranging from 29 in 2003 to 34 in 2008. Furthermore, there were two pits in the accumulation zone, where measurements were made once every year.
Energy fluxes of shortwave-incoming radiation, shortwave-reflected radiation, longwave radiation incoming from the atmosphere and longwave radiation at the surface, as well as the net radiation, are available every 30 min from 2004 to 2008. However, data gaps occurred in year 2006 in the period 9–19 March, 4–19 April, 25–31 May, 24 June 24 to 19 July, and 30 August to 1 November, and in the year 2007 in the periods 28 January to 13 March, 16–24 March, and 29 March 29 to 28 August. These gaps were filled in with the help of various interpolation and autocorrelation methods, specifically developed in the ®-environment by the authors, which take into account the seasonality of the corresponding time series and/or employ data from nearest neighbor stations. For further details, we refer the reader to [34].

2.3. Methods

2.3.1. Surface Distributed Energy Balance Model (DEBAM)

DEBAM General Modeling Approach

The distributed energy balance model (DEBAM) developed by Hock [36] allows to calculate the specific mass balance and discharge by first estimating the amount of melting derived from the interaction of the energy fluxes at the glacier surface. The energy gained is used for the fusion of the ice or the snow to water when the surface temperature is 0 °C or above. Another considered state of water change by the model is sublimation. Sublimation is always calculated estimated when there is a loss of latent heat from the surface. In the case of gains of latent heat at the glacier surface, condensation and refreezing/deposition are obtained.
The glacier calculation of the mass balance considers the accumulation, deposition, fusion, sublimation, and condensation in each grid point. It determines the specific feeding and loss of the Artesonraju glacier (3.5 km2), which belongs to the Artesoncocha basin (7 km2). The calculation of the discharge takes into account not only the Artesonraju glacier area but also the melting of snow zones from other nearby snowcapped mountains as runoff of some areas devoid of snow that contribute to the inflow to the Artesoncocha Lake (see Figure 1 and Figure 2).
First of all, the focus will be on the computation of each component of the energy balance equation. We shall explain how each flux is calculated and which parameters are considered. Afterwards, we shall elucidate how this energy is related to melting and how it is converted in water equivalent for the reckoning of discharge and the specific mass in the glacier.

Energy Balance Equation

The DEBAM model described as follows is based on the energy balance equation:
Q M = Q N e t + Q H S + Q H L + Q P + Q G ,
where Q M is the energy flux available for melting, Q N e t is the net radiation, Q H S is the flux of sensible heat, Q H L is the flux of latent heat, Q P is the heat flux from precipitation, and Q G is the subsurface energy flux. Units of all energy fluxes (power) are Wm−2. Downward fluxes towards the surface are considered positive (gain) and those upwards from the surface are negative (loss).
These flux variables will be discussed in detail in the following sub-sections with additional information provided in Appendix A.

Net Radiation

The net radiation Q N e t is the sum of net shortwave radiation and net longwave radiation fluxes
      Q N e t = ( S w i n c + S w r e f ) + ( L w i n c + L w s u f c )
where   S w i n c is the incoming shortwave radiation from the atmosphere, S w r e f is the shortwave radiation reflected by the surface, L w i n c is the longwave radiation incoming from the atmosphere, and L w s u f c is the longwave radiation emitted by the surface.
Further details about the estimation of each term of Equation (2) are listed in Appendix A.

Albedo

The reflected shortwave radiation S w r e f is defined to   S w i n c over the albedo α of the surface by:
  S w r e f = ( 1 α )   S w i n c
The albedo term α in Equation (3) turns out to be a very sensitive parameter in the DEBAM model. For the current simulations, the corresponding albedos are allocated for each glacier surface type that may be snow, ice, or firn, wherefore αsnow ranges between 0.75 and 0.98 for fresh dry snow, between 0.70 and 0.85 for old clean dry snow, and between 0.46 and 0.70 for old clean wet snow, αfirn ranges between 0.3 and 0.65, and αice between 0.06 and 0.46 [37].

Turbulent Fluxes

The bulk aerodynamic method is used for the estimation of turbulent fluxes (Tf) of latent ( Q H L ) and sensible heat ( Q H S ) . This method takes into account the differences in temperature, wind speed, and vapor pressure between the surface and a certain height in the atmosphere.
QHS is calculated as:
Q H S = c p k 2 ρ o P o P u 2 ( T 2 T S ) l n ( z / z o w ) l n ( z z o T )
where k is the Karman constant (=0.41), P is the atmospheric pressure (Pa), ρ o is the air density at P o (=1.29 kgm−3), T s is the surface temperature (°K), z o w and z o T are the roughness lengths for the logarithmic profiles of wind and temperature and z is the instrument height above the surface (2 m). Site values of roughness lengths of temperature and water vapor range between 1 and 10 mm for the Artesonraju glacier [30].
Q H L is calculated as:
Q H L = 0.623 L v / s k 2 ρ o P o u 2 ( e 2 e o ) l n ( z / z o w ) l n ( z / z o e ) .
where   L v / s is the latent heat of evaporation or sublimation, depending on which phenomenon is occurring. With a negative latent heat flux, sublimation occurs, whereas with a positive flux, condensation will occur for positive surface temperatures and deposition for negative ones. e is the vapor pressure, and z o e is the roughness length for a logarithmic profile of water vapor.

Ground Heat Flux and Rain Energy

The subsurface conductive heat transfer is minor if the temperature of the ice or snow is at its melting point [37]. For temperate glaciers (such as Artesonraju), which are characterized by a temperature at or near their pressure-melting point throughout the ice mass [38], the energy ground flux is usually neglected.
The rain energy Q R is calculated as:
Q R = ρ w c w R ( T r T s ) ,
where ρ w is the density of water, c w is the specific heat of water (J kg−1 K−1), R is the rainfall rate (msec−1), T r is the temperature of rain assumed to be equal to air temperature (K), and T s is the surface temperature (K).

2.3.2. Mass Balance Estimations

Accumulation

The glacier accumulation is calculated from the snowfall. The state of precipitation is calculated using the so-called threshold temperature snow/rain T0. In the temperature interval of T0 − 1 K and T0 + 1 K, the percentage of snow or rain is linearly interpolated.

Ablation

The ablation is due to melting and sublimation. The melting M, in water equivalent melt (m/day) is calculated as:
M = Q M ρ L f ,
where Q M is the energy flux available for melting, ρ is the water density, and L f is the latent heat of fusion (=334,000 J kg−1).
If the latent heat flux Q L is negative, no matter what the surface temperature is, sublimation S is assumed and computed as:
S = Q L ρ L s ,
where L s is the latent heat of sublimation (=2.849 × 106 J kg−1). If the latent heat flux is positive and the surface temperature is 0/<0, condensation /deposition is expected. The deposition is present when the surface temperature is <0.

Equilibrium Line Altitude (ELA)

The changes of the surface type of the grid quadrants around the equilibrium line altitude (ELA) are very important, as they eventually define how the glacier behaves. These changes are calculated from the initial surface maps at the beginning of the hydrological year and the initial conditions of the gains and losses of the surface at the start day of the simulation period. The snow area is determined by the solid precipitation according to threshold temperature rain/snow precipitation. The ice area, located below the ELA, is the zone with bare ice after the snow cover has melted away, and the firn area above the ELA is the one free of snow.

2.3.3. Discharge Estimations

The DEBAM glacier model uses a classical routing discharge model [39], which is based on the idea that different surfaces release water at a rate proportional to the storage amount which is then routed across the glacier surface grid in a time-delayed manner. More specifically, the sum of melt and rainfall is converted into discharge using a linear reservoir approach, which relies on the non-stationary water budget equation:
dS/dt = R(t) − Q(t),
where S is storage, expressed as S = k ∗ Q (linear reservoir, with k, the storage coefficient), R and Q are inflow and outflow into an area, respectively, and t is time. A separate linear reservoir model is set up for snow, ice, and firn areas, and specified by different storage coefficients k.
The computational, integrated version of Equation (9) is then [36]:
Q ( t 2 ) = i = 1 3 Q i ( t 1 ) e Δ t k i + R ( t 2 ) R i ( t 2 ) e Δ t k i
The cyclic (yearly period) input of melting and precipitation generates a corresponding dampened cyclic variation of runoff, with peak values delayed relative to the maximum input time [37]. After the input signal has ceased, the runoff declines exponentially with time, depending on the size of k, whose exact value for each surface type is determined in the model as part of the discharge calibration process.

2.3.4. Model Fit Evaluation

Finally, the quality of the DEBAM model in simulating the observed discharge (and other directly measured glacier mass balance parameters) in calibration and validation periods is evaluated by the commonly used Nash Sutcliffe coefficient (E) [40], which is defined as:
E = 1 i = 1 n ( O i P i ) 2 i = 1 n ( O i O ¯ ) 2
where Oi is observed and Pi predicted (simulated) data value. The theoretically possible range of E lies between 1.0 (perfect fit) and −∞. Acceptable ranges for a reasonable fit are usually 0.5 < E < 1 [41]. In contrast, an E < 0 is considered so bad that even the mean value of the observed time series would be a better predictor than that of the model. It should be noted that other measures of a model fit, such as the coefficient of determination R2 (with definition essentially similar to E) and/or the PBIAS have been proposed for evaluating general hydrological model fits, but their advantage over E is not always clear and appears to be more a matter of taste [41]. Another accuracy parameter that has been used here is the root mean square error (RMSE), for which, in contrast to the other parameters mentioned, no standard reference quality value is defined, as RMSE depends on the unit and the possible measurement error of the observable fitted.

2.4. DEBAM Input Data and Its Pre-Processing

The DEBAM model requires a large number of climate input variables and other parameters and/or physical constants, some of which need to be fine-tuned during the calibration of the model. Table 2 summarizes these input variables with a short description of how these are estimated, whether they are invariant or not across the model grid and how, for the latter case, the ensuing grid interpolation is carried out. The data can essentially be categorized in four groups: climate data, radiation data, hydrological data and geographic grid data. Target calibration data of the glacier model are measured glacier mass balances and discharge at the Artosoncocha Lake outlet. Important calibration parameters whose sensitivity to the model’s output will be discussed in detail in the results sections are the albedo, the reservoir storage constant and the rain/snow threshold temperature.

3. Results and Discussion

The accuracy of the DEBAM model in predicting the glacier’s changes is quantified with the Nash Sutcliffe coefficient (E) (Equation (11)) and the RMSE of the simulated versus observed discharge. Moreover, field readings of stakes and pits were compared with the cumulated mass loss/gain simulations at the corresponding locations. Annual mass balances and the ELA positions are also compared to local estimations included in WGMS (2015) [35].
In Section 3.1, Section 3.2 and Section 3.3, we first present the physically most palatable modeling results in terms of the glacier’s meltwater discharge as well as its mass-balance changes over the simulation period 2004–2007, and in Section 3.4 further details on the model optimization including a parameter sensitivity/uncertainty study. The latter is further elaborated on in Section 4.

3.1. Simulated Discharge

Although the discharge is the final outcome of the whole chain of model driving parameters (energy, meteorology, hydrology), it will be discussed first, as it is the ultimate measured parameter on which all other parameters of the model as described in the following sub-sections are calibrated. Figure 3 shows the performance of the calibrated model for the three hydrological years in terms of the Nash Sutcliffe coefficient (E) and the RMSE for the fitted discharge. Based on these values, the adjustment goodness is marked in the last column of Table 3 using the general quality classification for E of Cabrera (2009) [42]. The lowest E and maximum RMSE are obtained in the year 2005–2006, but the adjustment quality is still considered as good.
The observed and simulated discharge hydrographs for the various years analyzed (2004–2007) are illustrated in Figure 3. One may note from the different panels that the model cannot represent well the most prominent discharge peaks, especially, in the month December 2005–2006 and, to a lesser extent, some peaks of October and April of 2004–2005 and 2006–2007. For the hydrological years 2004–2005, 2006–2007, and the months of ASOND of the year 2007, the major discrepancy occurs in October. Significant bias in the estimation of some minima occurs mainly in the months of January and February 2005–2006, 2006–2007, and 2004–2005. These differences may be a consequence of disregarding the rapid fluctuations of albedo that could be produced by quick changes in precipitation [43,44]. However, the characteristic trends of the season are captured by the seasonal albedos, i.e., despite the use of seasonal albedos in the simulation, the performance of the model to simulate the discharge appears to be good enough. In fact, albedo parameterizations of the present DEBAM model [45] require an hourly resolution of the input data, but not all climate variables for the CB are available with that resolution. Albedo simulations may be particularly important in the rainy season when the reflection of shortwave radiation varies more strongly over the days.
Nonetheless, the use of seasonal albedos for each surface appears to be warranted for calculating the typical average discharge in an acceptable way and, consequently, the amount of water produced by the glacier. Moreover, the simulation does not estimate subsurface energy transfer and subsurface flows, due to a lack of data in that regard. For instance, the subglacial flow of melting through channels could be enhanced in the core of the wet season, when there is more water in the snow [46]. This uncertainty of the amount of subglacial flows could be at the origin of the observed underestimations of peak discharges in Figure 3.

3.2. Stake and Pits Measurements, Mass Balance Profiles, Annual Mass Balance, and ELAs

The mass balance is calculated only for the Artesonraju glacier, excluding some snow-covered areas of the Piramide, Paron, and Caraz peaks. These excluding areas also contribute to the discharge but do not feed specifically the Artesonraju glacier. The calculation area for the mass balance is 3.5 km2, with a DTM of 30 m resolution. This area corresponds approximately to the one used by the UGRH in its calculations [35]. The simulated annual mass balances are compared with the ones measured UGRH.
Figure 4a shows the cumulative mass at stakes/pits locations retrieved from the current simulation and Figure 4b the mass balance profiles for each of the three hydrological years. The estimated E and RMSE for year 2004–2005 are 0.71 and 1.59 m, respectively; 0.58 and 2.36 m for 2005–2006; and 0.69 and 1.9 m for year 2006–2007. For all three years, the model performance in simulating the cumulated mass balance is good. The error is quite similar for the three years.
One may also notice from the different panels of Figure 4 that for year 2004–2005, the energy balance model overestimates the amount of mass loss in the area of ablation. The contrary occurs in year 2006–2007, whereas in year 2005–2006, overestimations and underestimations of mass measurements occur alike. These outliers are most likely caused by the use of seasonal albedos in the model that cannot capture all changes that happen to the snowline position during the day. Moreover, the simulated mass in the accumulation area is very close to the measurements in the pits.
The ELA numbers listed in Table 4 indicate a good agreement between the UGRH reported and simulated ELAs, with minimal differences between the two of only up to 35 m. This agreement is essential because the ELA separates the accumulation zone where snow and firn prevail from the ablation zone, majorly formed by ice. The 2004–2007 yearly displacement of the simulated ELA reveal the same pattern as the one of the yearly snow line identified by Rabatel et al. (2012) [47] in Artesonraju (see Table 4). Meanwhile, the UGRH- ELA estimations show opposite displacement for the years 2004–2005 and 2005–2006.
The simulated mass balance profiles (see Figure 4b) unveil vertical mass balance gradients of ~4 m.w.e /100 m in the ablation area, but which decrease to ~0.5 m.w.e /100 m in the accumulation area. The estimations of the accumulation area ratio (AAR), defined as the ratio of the area of the accumulation zone to the glacier area, are also presented in Table 4 and indicate variations ranging between 59% and 65% of the total glacier area. In fact, changes in the glacier zone areas have a significant impact on the estimations of the mass balances. However, glaciers with high slopes are less sensitive to AAR changes, as is the case for the Artesonraju glacier.
The simulated, mostly negative annual balances agree rather well with those estimated by UGRH (see Table 4). Some minor differences occur in the year 2005–2006 when the simulated losses are lower than those estimated by UGRH. According to our model results, the largest yearly mass losses occur in year 2004–2005, whereas the UGRH-estimated losses are the largest in 2005–2006. However, these higher simulated losses for the first model year 2004–2005 are congruent to findings in the Zongo glacier, where the second-largest specific mass balances losses (after 1998–1997) were also measured in year 2004–2005 (−1.69 m.w.e.) [48]. Moreover, the location of the snow lines [47] shows steady increments of accumulated snow between 2004–2007. This trend indicates an increase in the accumulation area and, therefore, improvements in the glacier’s mass balance conditions towards the subsequent year 2007, all of which is in agreement with the current DEBAM simulations.

3.3. Seasonal Variations of Mass Balances and Association with Temperature and Precipitation

Figure 5 presents in a more detailed manner the seasonal variations of the simulated mean daily glacier mass balance, together with the mean daily temperatures and the mean daily precipitation. Generally, the steady mass losses obtained for all seasons in the year 2004–2005 and the beginning of the year 2005–2006 correlate clearly with seasonal increments of temperature and declines in precipitation for the periods NDJFM4-5, AMJJA4-5, and SO5-6, though this correlation is somewhat less stringent for the remaining seasonal periods between the years 2005 and 2007.
As the climate of the CB is strongly affected by ENSO tele-connection effects, the above findings can also be interpreted in terms of cyclic ENSO phenomena. Thus, a weak ENSO (El Niño) caused slight temperature increments and precipitation reductions in the CB in 2004–2005 [49]. For example, Figure 5 identifies the highest daily mass losses for season SO6-7, during which the precipitation was also significantly reduced. Indeed, one can observe that this season SO6-7 had the lowest precipitation from the beginning of the rainy seasons of all SO –periods in the 2004–2007 time span. Moreover, in this SO6-7 period, filled data of shortwave radiation incoming and longwave radiation was used. Ensuing small bias in this interpolated data produced increments in the calibrated threshold temperature T0 (see Section 3.4.2 later), which could lead to an overestimation of the mass balance loss in this SO6-7 period. Lagged cross-correlation analysis of ENSO- (El Niño SST) signals of Pacific Region 3.4 with climate time series in the CB by the authors [34] show a lag of one to four months between ENSO anomalies and their effects on the glacier areas of CB. Therefore, the El Niño of 2004–2005 probably extended its teleconnection until SO5-6 in the glacier, with lower precipitation and slightly higher temperature in that period.
In contrast, during the year 2005–2006, a Niña event occurred during October to April in Niño 3.4 Region [49], which can have a lagged effect in the glacier with a delay of up to 7 months. This Niña caused slight and steady increments of the precipitation that allowed the glacier to reduce losses in the first 9 months of year 2005–2006. Again, between August and January of the year 2006–2007, a weak Niño occurred in the Nino 3.4 region, which may have had a slight effect on the season AMJJA6-7. All these findings support the general idea that mass balances in the CB glaciers are significantly affected by El Niño/La Niña events [3,13,22,23,50], wherefore, in particular, the here observed and simulated seasonal mass balance of the Artesonraju glacier may be explained in part by intermittent occurrences of such ENSO phenomena.
The above results indicate that in the simulated 2004–2008 time period, the glacier is far from reaching equilibrium. During most of the seasons in that time span, the glacier experienced losses of mass, though with alternating negative and positive trends, triggered by ENSO phenomena. As the changes in the Artesonraju glacier’s simulated mass balance also have indirect effects on the observed and simulated discharge at the outlet of its catchment (Artesoncocha lake), further analyses of the possible cross-correlation of the measured discharge with various climatic factors resulted in a correlation coefficient of r = 0.35 of the former with temperature at the climate station of Artesonraju glacier, of r = 0.43 at climate station Artesoncocha, and of r = 0.24 with Rh (relative humidity), of r = 0.26 with Lwinc, and of r = 0.28 with precipitation. Although these correlation coefficients may not appear too significant, they are basically consistent with the findings above regarding the effects of different seasonal climate variables (temperature and precipitation) on the glacier’s mass balance.

3.4. Model Calibration, Optimization and Sensitivity Analysis

The optimization of the DEBAM model was made by fine-tuning pairs of different input parameters (see Table 2). Initial model simulations indicated that the assumed albedo value α of the different glacier surface types (snow, ice, firn), the threshold temperature rain/snow T0, had a very strong impact on the model results for mass balance and discharge—and for the latter, also directly the storage constants k for the three surface types of the routing reservoir model—as quantified by the ensuing E and RMSE-values (see Table 3). For better manageability and illustrative purposes, these three most sensitive parameters of the model were mostly tested in pairwise combinations, and the corresponding results are presented and discussed in the subsequent paragraphs.

3.4.1. Albedo

Different combinations of albedos for each of the glacial surfaces were tested, but by respecting the physical ranges of their possible values. The model sensitivity to the seasonal albedo selection was also investigated under different threshold temperature T0 conditions.
Figure 6 shows the response surfaces of the Nash Sutcliffe parameter E of the simulated discharge for different combinations of snow and ice albedos for each season of the three-year simulation period. Generally, one can notice for all three years similar tendencies of the E response surfaces for the same season.
More specifically, for the SO seasonal period, high calibrated values of the snow albedo are observed, varying from 0.90–0.95, whereas that of ice ranges between 0.35–0.5. This seasonal high albedo can be explained by the beginning of the wet season in which fresh snow is typical and clean ice prevails in the ablation zone.
For the subsequent wet season period NDJFM, the snow albedo decreased to a range of 0.80–0.88. The same holds for the ice albedo, which declined to 0.25–0.3 for years 2004–2006, but increased again to 0.45–0.5 for 2006–2007. In this wet season, more significant fluctuations of precipitation and a slightly higher temperature usually occur, so that the glacier melting is higher. The ensuing meltwater reduces the snow albedo directly and indirectly increases the grain size of the ice/snow surface which, in turn, reduces albedo [43]. The above albedo fluctuations are consistent with typical values from fresh snow and wet snow [37].
Lastly, for the AMJJA dry season period, the snow albedo ranges between 0.82 and 0.95. During that timespan of low precipitation, the water content is reduced, maintaining the snow accumulated in the previous season. Therefore, the snow albedo tends to increase compared with the previous season, which is also concordant with what was observed by Juen [15]. In fact, in the present case, two of the simulated seasons between 2005 and 2007 are characterized by (lower) albedos of old dry snow.

3.4.2. Threshold Temperature Rain/Snow T0

The threshold temperature rain/snow (T0) determines the surface or size of the areas of snow, ice, firn, and debris which, in turn, by their different physical features, influence the overall energy budget of the glacier. Generally, T0 is variable and depends on geographical and seasonal atmospheric conditions. The latter influence the state of precipitation through different mechanisms. For example, the thickness of the air layer that falling water droplets (leading to precipitation) must go through, the soil heat, the cooling influence on humidity, and the salt content of the particles affecting freezing point are some factors that influence the type of precipitation eventually reaching the soil [51,52,53]. A study with over 1000 climate stations in the United States of America showed that rain never occurs beneath a temperature of 0.8 °C and that snow is never observed, when the temperature exceeds 6.18 °C [54]. Observational data in the northern hemisphere indicate that the warmest threshold temperatures occur in certain mountain regions above 4000 m.a.s.l., e.g., Rocky Mountain 3.8 °C and Tibetan Plateau 4.5 °C [55]. Another study, specifically in the Cordillera Blanca of interest here, estimated a rain/snow threshold temperature range of 1.1–2.5 °C in the Shallap glacier [14].
Based on these studies, it is of no surprise that T0 turns out also to be a very sensitive parameter in the present model simulations, as illustrated by Figure 7, which shows the seasonal T0 sensitivity for the three hydrological years in terms of its effect on E. Figure 7 reveals that the best seasonal rain/snow threshold temperature T0 ranges between 2.6 °C and 3.8 °C, which are higher values than those found above in glacier Shallap [14]. These disagreements may be due to uncertainties in the spatial variation of other climate variables, such as precipitation, wind, relative humidity, vapor pressure, and daily snowfall variability. Moreover, any deviation of the assumed initial conditions from the real limits of the different surface types can be compensated by increasing the threshold temperatures.
Eventually, ENSO anomalies can also produce changes in T0. For the CB study region, an ongoing El Niño event usually leads to strong dryness—as is the case, for instance, for the one occurring during NDJFM and AMJJA in year 2004–2005, when higher seasonal daily mean temperatures and lower daily precipitation occurred (Figure 5). In fact, snowfall events at a lower relative humidity are more likely to fall as snow at higher T0, i.e., 4.5 °C (40–50%), 3.7 °C (50–60%), 2.8 °C (60–70%), 2.2 °C (70–80%), 1.4 °C (80–90%), 0.7 °C (90–100%) [55]. As for the glacier Artesonraju, the seasonality relative humidity shows reductions, especially in the dry season, to 62–75%, compared to 74–80% in the wet season, all of which is in tendency with the present results (Figure 7). In general, as the relative humidity in the Artesonraju glacier turns out to be slightly lower than that in the southern glaciers of the CB [34], the higher T0 modeled here are understandable. Notwithstanding, T0 could be overestimated in the model calibration, compensating for any uncertainty or inaccuracy in the input data, as mentioned above, i.e., bias in the delimitation of the different surface areas, but also uncertain precipitation losses and gradients.

3.4.3. Discharge Storage Constants

The optimal discharge storage constants k of the reservoir routing model (Equations (9) and (10)) obtained as part of the sensitivity study for the three glacier surfaces and the different seasons are listed in Table 5. Similarly to the two previously discussed sensitivity parameters, the spread of the storage constants across the three years simulated is relatively low for the same season. However, differences arise for the individual seasons. Thus, for the first season, SO, values of kfirn of 700–800 h, for the second (rainy) period, NDJFM, of 500–800 h, and for the third (dry) season, AMJJA, of 1100 h are obtained. The storage constants of firn are consistent with the cycle of firn formation in the hydrological year, i.e., maximum values at its end in the dry season (AMJJA), when a larger residence time prevails, owing to more snow accumulation (and conversion to firn) from the previous season. The lower kfirn—values obtained for season NDJFM can be explained by the fact that more melting water penetrates into the firn’s interstices. In the first (2004) SO season, on the other hand, the melting is still low, and the firn formed in the preceding season still prevails, which leads to a slightly larger kfirn-value.
For ksnow, the seasonal ranges are 390 h, 390–900 h, and 150–200 h for seasons SO, NDJFM and AMJJA, respectively; i.e., the ksnow-values for the two rainy periods SO and NDJFM are much larger than that of the dry season AMJJA. This is due to the fact that there is more snow during the two rainy seasons and which is maintained for a longer time period. Indeed, for tropical glaciers, the total snow accumulation of a year occurs in these two rainy seasons.
Finally, for kice, values ranging between 30 and 50 h for NDJFM (wet season) and of 250–300 h for AMJJA (dry season) are encountered. This is congruent with large melting and ensuing discharge in the wet season, resulting, consequently, in a lower residence time for ice, but lower melting and discharge in the dry season, i.e., a higher residence. However, for the season SO of 2006–2007, an extremely high value of kice = 500 h is found that may not be realistic and is most likely due to some gaps in the input data during that period which, though filled in, as discussed earlier, may lead to some modelling errors.
Table 5 shows further that the highest storage constants are obtained in the dry season (AMJJA), mainly for firn and ice surfaces. This is logical, as during that season, the glacier receives less energy input (southern winter) and, thus, less melting takes place, ergo, less discharge at the glacier’s outlet. Interestingly, the situation is somewhat opposite for know, where, for the dry season, the k values are the lowest. Indeed, as the solid precipitation is also reduced in that dry season, the sporadic snowfalls can melt faster, leading to a shorter snow storage time.
In conclusion of this section, it can be stated with some confidence that the optimally calibrated discharge storage constants represent the seasonality of accumulation and ablation of snow, ice, and firn, typical of the outer tropical glaciers, satisfactorily.

3.5. Contribution of Seasonal Energy Fluxes on the Melting Energy Balance

The mean seasonal contributions of the different radiation energy flux terms to the overall energy balance of the glacier (Equation (1) and following) for the whole 2004–2007 simulation period are presented in Figure 8a, which indicates that the net radiation (QNet), which eventually drives the total mass losses on the glacier’s surface, ranges between 15 and 73 Wm−2, depending on the season. In particular, QNet decreases in all dry seasons of the simulated years. This seasonal variation of QNet is consistent with the typical seasonal cycle of ablation of the outer tropic glaciers.
Figure 8b shows in more detail that lower shortwave balance (Swbal) and longwave balance (Lwbal) in the dry and cold seasons leads to a large decline of QNet, which, obviously, reduces glacial melting.
The net shortwave radiation balance (Swbal) has the most significant weight within the net radiation term (Equation (2)). Although there is more precipitation in the wet seasons (Figure 5), the calibrated albedos of the snow surfaces have similar seasonal values in the core wet season and the dry season (see Section 3.4.1). In the wet seasons, not only does more snow precipitation fall, but the latter also has a larger water content. This occurs during a time period of maximum net radiation and slightly increased temperatures (Figure 5). All this could reduce the snow albedo values for the wet season a bit. Nevertheless, the mean shortwave radiation reflected (Swref) (Equation (3)) of the whole glacier is comparatively larger in the two wet seasons than in the dry season. This can be explained by the fact that, during the dry season, the overall snow surfaces of the glacier shrink, reducing its average total reflection (lower Swref).
Swbal is generally higher in the rainy (51–86 Wm−2) than in the dry (43–63 Wm−2) season. Moreover, Swinc decreases slightly after the rainy seasons during the simulation years 2004–2007. This reduction in the subsequent dry seasons is to be expected, as with then lower precipitation and less cloudiness, diffuse incoming radiation diminishes.
As for the net longwave radiation (Lwbal), one can notice from Figure 8b that it is also reduced for the dry seasons. For the latter, the diagrams show, furthermore, that Lwbal makes up about 63%, 73%, and 66% of the Swbal energy received in years 2004–2005, 2005–2006, and 2006–2007, respectively. This overall dry-season reduction of Lwinc, necessarily, causes an even more significant drop of Lwbal during that time period, when there are more clear sky days.
It should be noted, finally, that these findings regarding the seasonal variations of longwave radiation balance in the Artesonraju agree well with those obtained for the Zongo glacier [19]. Further comparisons with other studies will be made in the subsequent section.
The seasonal mean of Tf (comprised of the sensitive heat QHS as well as the latent heat QHL, see Equations (4) and (5)) increases a little in the dry seasons of 2004–2005 and 2005–2006, compared to the wet seasons (Figure 8a). These seasonal increases (mainly of QHS) may be the result of temperature rises in the named seasons due to the then-occurring El Niño phenomenon.

3.6. Comparison of Simulated Energy Balances with Those of Other Studies in Tropical Glaciers

Averaging the seasonal flux terms of Figure 8b over the year, the corresponding values, as listed in Table 6, are obtained. Our estimations of the energy fluxes Swinc, Lwinc and Lwsurf are quite similar to other studies carried out in the Andean tropical glaciers. Energy fluxes in the glaciers of Zongo, Antisana, Shallap, and Artesoranju for different years [3,10,11,12,16,17,21,22,56] show that the annual averages of Swinc range between 206 and 239 Wm−2, Lwinc between 258 and 310 Wm−2, and Lwsurf between −277 and −311 Wm−2. Our annual averages for these variables are within the mentioned ranges.
On the other hand, the mean annual reflected shortwave radiation (Swref) range of 151–162 Wm−2 tends to be a little above the average range estimated in the mentioned studies of tropical glaciers (−86–137 Wm−2). Nonetheless, this discrepancy may be partly understood by considering that (a) Table 6 shows annual averages of shortwave reflected in the whole area which, keeping in mind that the AAR is between 59 to 65% (Table 4), the albedo in the accumulation area, which is predominantly high, raises the average, whereas some of the other studies presented site values; (b) different local climatic conditions could occur in the tropical glaciers, for instance, in the CB the influence of charged air masses that ascend from the eastwards located Amazon basin may also be felt in the Artesonraju glacier (indeed, there is a low correlation of the precipitation within some subbasins in the CB); (c) the use of seasonal albedos may influence the results obtained for Swref, since daily fluctuations are not captured; (d) the uncertainty in the assumed precipitation gradients is high, as there are no gauge stations in high elevations of the glacier. Notwithstanding, the simulated Swref values at the climate station are similar to the measurements there. The turbulent fluxes (Tf) found here (Table 6) are between the ranges found for other tropical glaciers, in which QHS is between −2 and 21 Wm−2 and QHL between −27 and −4 Wm−2. Generally, the turbulent fluxes Tf are comparatively less than the net radiation QNet, which is due to the fact that, over the whole year (but not for separate seasons), the mean daily sensible (QHS) and latent heat (QHL) fluxes counterbalance themselves pretty much, here, with their sum (=turbulent fluxes, Tf) hovering between +1 and 0 Wm−2.

3.7. Vertical Profiles of Seasonal Energy Balance Flux Terms across the Glacier’s Elevation

In the present sub-section, vertical variations (profiles) of the various simulated seasonal energy balance flux terms in the across the glacier elevation range are presented. To that avail, maps of the simulated mean seasonal distributed energy fluxes, with a horizontal resolution of 30 m × 30 m, are divided into elevation bands, each of which covering a range of 50–60 m of elevation, resulting in a total of 24 bands across the whole glacier altitude range. A mean energy flux is then calculated from the aggregated grid point values in one elevation band. Vertical profiles of the different seasonal energy flux terms obtained in this manner are shown in Figure 9, Figure 10 and Figure 11.
The vertical profiles of QNet (Equation (2)) Tf (Equations (4) and (5)) and the ultimate energy remaining for melting QM (Equation (1)) are illustrated for each season of the 2004–2007 simulation period in Figure 9. The three panels depict the results for the three seasons of the hydrological year, but for each year individually, allowing to better follow the seasonal time evolution of the energy flux terms across the analysis period.
One can notice from Figure 9 that QNet is much more prominent in magnitude and more variable (ranging between −1 and 177 Wm−2) in the lower-lying ablation zone than in the higher-lying accumulation zone (ranging between −15 and 43 Wm−2). This pattern is congruent with the typical minimization of mass losses in the accumulation zone and maximum mass losses in the ablation zone. Moreover, low QNet radiation occurs in the accumulation zone occasionally at the beginning of the rainy season (SO) and predominantly in the dry season (AMJJA).
The peculiar vertical and seasonal variations of QNet across the glacier’s elevation are a direct consequence of the underlying definition of QNet as the sum of the shortwave (Swbal) and longwave (Lwbal) radiation balance (Equation (2)). The vertical profiles of these two radiation fluxes are depicted in Figure 10. From this figure one can notice that high values of Swbal (between 30 and 194 Wm−2) occur in the ablation zone. In contrast, in the accumulation zone, these values are low (between 7 and 79 Wm−2) due to the increased albedo here (Equation (3)), which is characteristic of this zone (snow). On the other hand, Lwbal values differ less between accumulation (between −46 and 17 Wm−2) and ablation (between −66 and −12 Wm−2) areas, and this holds for all season. As these Lwbal values are predominantly negative across the whole glacier, particularly in the ablation area, where surface temperatures of 0 °C prevail, it becomes clear that the longwave radiation balance is not responsible for any glacier melting. This means, last but not least, that QNet in the ablation area is controlled by the shortwave radiation balance (Swbal).
The magnitude of the Tf flux terms in the two zones is less variable than that of the QNet. However, these fluxes tend to be higher (ranging between −16 and 30 Wm−2) in the ablation zone than in the accumulation zone (ranging between −40 and 15 Wm−2). In the latter, Tf is predominantly negative for both rainy seasons (SO, NDJFM), but also shows more significant variations (ranging between −40 and 15 Wm−2) and which are reduced again in the dry season (AMJJA) (between −2 and 15 Wm−2). This reduction suggests the strong influence of wind on Tf, rather than relative humidity and temperature. Wind increases in the wet season, whereas relative humidity and temperature decrease in the dry season [34]. Furthermore, as Tf is not high in the ablation zone (between −16 and 30 Wm−2), compared to the QNet radiation term (between −1 and 177 Wm−2), the former flux term does not drive the melting in this zone.
However, the situation is different for the accumulation zone, where the Tf profiles suggest that they may act as a significant sink of energy, especially in the beginning and the core of the wet season (see panels SO and NDJFM of Figure 9 discussed above); i.e., Tf offsets the effects of QNet in that zone. This energy sink represents the glacier mechanism to minimize the melting in the accumulation zone in the wet season (SO, NDJFM), in spite of higher radiative energy input during that period. Indeed, in the dry season (AMJJA), QNet is reduced due to lower back-emission of longwave atmospheric radiation (Lwinc) and slightly less incoming short wave solar radiation (Swinc).
A more detailed picture of how the turbulent flux Tf is divided across the glacier’s elevation into its two components, QHS (sensible heat) and QHL (latent heat), is presented in Figure 11 The sensible heat QHS varies between 2 and 29 Wm−2 in the ablation zone and between −47 and 15 Wm−2 in the accumulation zone; i.e., it is mainly negative there, particularly during the wet season. On the other hand, QHL is positive in SO, resulting in a loss of energy that enables deposition processes, but negative in NDJFM, generating energy for sublimation which requires eight times more energy than melting. Due to this, melting in the accumulation zone is prevented, although higher QNet is available in these two wet seasons. Additionally, as already mentioned for Figure 9, Tf is lower in the dry season which, by virtue of Figure 11, can now be explained by the fact that the two heat flux terms QHS and QHL counteract each other more or less over the elevation range of the glacier. Nevertheless, slightly higher values of Tf are obtained in the dry season of the two El Nino years 2004–2005 and 2005–2006, due to small increments of QHS.
In conclusion, the vertical profiles of the various energy balance equations terms presented in this sub-section allow to clearly distinguish the ablation (ice-covered) and accumulation (snow-covered) zones across the elevation of the glacier and their development over the various seasons of the analyzed 2004–2007 time period.

3.8. DEBAM Glacier Simulation Uncertainties and Their Possible Sources

The current DEBAM glacier simulations, in conjunction with the limited quantity and quality of the input and calibration data, are naturally fraught with some uncertainty and/or errors. One particular source of the latter may be the scarcity of mass-balance- and radiation measurements in the accumulation zone of the glacier. In the high mountains of the Cordillera Blanca, these glacier zones are usually steep and abrupt, making access for measurements difficult, let alone dangerous [33]. As the net and the instrumentation of the measurements are also costly, the latter generally cover just some parts of the glaciers [3]. Therefore, initial conditions for use in the numerical model, especially those describing accumulations, are still difficult to estimate over the whole area.
Another source of uncertainty is the spatial variation of the various climatic variables, such as precipitation, wind speed, and relative humidity. Precipitation was assumed with losses ranging between 20 and 30%. The wind speed in the DEBAM is assumed as constant across the whole glacier modeling area, which has direct consequences on the proper estimation of the turbulent fluxes (latent and sensible heat), so that these are, necessarily, fraught with certain inaccuracies. The latter will then propagate into the estimation of the net energy flux available for melting. Thus, exploring different parameterizations to account for the spatial wind variability may be an essential task for a better understanding of the role of sublimation, condensation, deposition, melting, and other variables, all of which depend strongly on the wind speed. Moreover, the direct interference of strong wind blowing leads to a reduction of the net snow accumulation at high altitudes, although it may produce accumulation in topographic depressions (so-called snowdrift) [3], both of which are out of the scope of the present model simulations.
The sporadic gaps in the measured energy flux time series filled with various statistical interpolation methods [34] are also sources of uncertainty. The possible errors generated in this way could affect the performance of the DEBAM model, as is reflected by the (relatively low) Nash Sutcliffe coefficients of the present simulations for the partly filled-in data of season SO of year 2006–2007.
Unknown daily variations of the albedo are another source of uncertainty. Indeed, variations of albedo are expected particularly in the wet season, due to the fact that strongly fluctuating precipitation amounts and moisture contents affect the permanence of the snow cover on the glacier’s surface. As a matter of fact, snow albedo falls exponentially after a few snowfall days [43,57]. However, such albedo changes are not considered in the current model.
The unknown amounts of energy transfers in the subsurface and the subglacial flows (see Equation (1)) are another source of uncertainty in the understanding of Artesonraju glacier dynamics and could contribute to an underestimation of discharge, especially in the wet season, when more water from melted snow flows downhill towards the outlet of Artesoncocha lake.
Finally, rock avalanches that produce significant changes of debris mass in high, abrupt mountains may affect the thermal environment on or near the glacier and so influence the dynamics of the latter. Such mass movements caused by massive slope failures (favored by rising temperatures, and so exacerbated by climate change) have been observed in high mountain glaciated areas of the alpine region [58].

4. Conclusions

By means of a novel distributed energy balance model (DEBAM), the present glacier modelling study sought to deepen the understanding of the various physical processes with their leading parameters that govern mass balance in tropical glaciers—here, the Artesonraju glacier in the CB, Peru. Thus, the results of the present modeling exercise contribute to the yet existing and occasionally contradictory conclusions of many other previous studies in this field, regarding the dominant drivers of the seasonal dynamics of tropical glaciers, in general.
Based on the scarcely available empirical input data, a time period of only three years (2004–2007) was modelled. In spite of this rather short calibration period, good-to-satisfactory results are obtained, with the Nash Sutcliffe coefficient E ranging between 0.60 and 0.87 for measured discharge at the Artesoncocha lake outlet (after fine-tuning the storage constants in the reservoir routing model in a seasonal manner) and of 0.58–0.71 for glacier stakes’ mass balance measurements. Moreover, annual mass balance estimations were compared with previous official calculations by UGRH [35]. The simulated annual mass balances show mass losses between −1.42 and −0.87 m.w.e., which turns out to be concordant with the official estimations. The annual and seasonal balances strongly depend on climatic variations associated with El Niño and La Niña phenomena, which occurred during the time period simulated.
A more detailed analysis of the various components of the energy balance equations indicates that the net radiation (sum of short- and longwave radiation balances) governs the energy balance and the melting processes, confirming results of some current studies [3,10,11,14,18,19,22]. More specifically, it is the reflected shortwave radiation which is the dominant component here. As the latter is governed by the albedo of the different glacier type surfaces (snow, ice, firn), further sensitivity studies with regard to that parameter were carried out. Although the daily short-term variations of albedo could not be simulated, the use of seasonal albedos for three different seasons (less rainy, rainy, dry) allowed to satisfactorily perform the modelling tasks. As a conclusion here, the seasonal albedos appear to reasonably account for the temporary characteristics of the surface conditions over the season under question.
On the other hand, longwave radiation still plays a role in determining seasonal glacier conditions, for instance, by reducing the radiation at the surface in the dry season when there are more clear sky days. This result is confirms statements of Sicart [19].
Vertical profiles of various seasonal energy balance flux terms across the glacier’s elevation reveal particularly interesting details on the factors controlling the temporal dynamics in the glacier’s accumulation and ablation zones. For example, the temporal analysis of the turbulent fluxes (sum of latent and sensible heat) reveals that their lowest values occur predominantly in the dry season. However, the turbulent fluxes are not mainly responsible for the glacial processes in the ablation zone, as they are comparatively lower than the net radiation.
In the wet season(s), in contrast, and for the accumulation zone, the turbulent fluxes could offset the net radiation there. This means that the turbulent fluxes then act as an energy sink, helping to maintain the glacier accumulation, despite an increase in net radiation in this (slightly warmer) period. Increments/decrements in wind speed in the wet/dry season are the most important cause of the major/minor transfers of turbulent fluxes.
In spite of these palatable results, it should be noted that the current DEBAM simulations do not consider the subsurface energy fluxes and flows, which could affect the estimations of turbulent fluxes and so increase the subsurface flows, especially in the wet season.
The role of the snow/rain threshold temperature T0 is investigated in a further sensitivity study and turns out to be a decisive parameter in the glacial modelling of (not only inner tropic glaciers, such as Shallap [14]). In fact, the outer tropic glacier Artesonraju studied here extends over the elevation band at which T0 occurs, with the latter determining so the snow area and, eventually (over the albedo), the magnitude of the reflected radiation. Thus, T0 controls the portion of accumulation and ablations areas. This is relevant because the physical surface characteristics regulate, in turn, the response to specific radiation fluxes, such as shortwave radiation reflected, longwave radiation emitted, and the turbulent fluxes. Other investigations of T0 in mountain regions show, furthermore, that it depends on other meteorological variables, such as humidity, pressure, and wind [55], so that all of these could then indirectly influence the spatial average of the various energy fluxes. In conclusion, the threshold temperature shows some variability that depends on geographical conditions or climate (ENSO) events, as is also found here during the 2004–2007 simulation period. More investigations in that regard are suggested to better determine the physical behavior/effects of threshold temperature in outer tropical glaciers.

Author Contributions

Conceptualization, M.F.L.G. and M.K.; methodology, M.F.L.G. and M.K.; validation, M.F.L.G.; formal analysis, M.F.L.G. and M.K.; investigation, M.F.L.G.; resources, M.K.; writing—original draft preparation, M.F.L.G.; writing—review and editing, M.K. and M.F.L.G.; visualization, M.F.L.G. and M.K.; supervision, M.K.; project administration, M.F.L.G. and M.K.; funding acquisition, M.K. All authors have read and agreed to the published version of the manuscript.

Funding

This research partly received external funding from the DAAD for collecting data in the study region.

Data Availability Statement

Datasets from climate data are available at Autoridad Nacional del Agua del Perú (ANA) web page (https://repositorio.ana.gob.pe/handle/20.500.12543/1497) (accessed on 12 September 2014), the Servicio Nacional de Meteorología e Hidrología del Perú (SENAMHI) (https://www.senamhi.gob.pe/?p=descarga-datos-hidrometeorologicos) (accessed on 9 September 2021), and the World Glacier Monitoring Service (WGMS) (WGMS, 2015) for the data provided for this research. Model is available at (https://github.com/regine/meltmodel) (accessed on 9 September 2021).

Acknowledgments

Thanks to the Autoridad Nacional del Agua del Perú (ANA), the Servicio Nacional de Meteorología e Hidrología del Peru (SENAMHI) and the World Glacier Monitoring Service for the data provided for this research.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Appendix A.1. Estimation of the Terms of Net Radiation

Shortwave radiation:
The shortwave radiation is the result of interaction between the incoming shortwave radiation and the reflected shortwave radiation taken from Hock [36,43,45].
The incoming shortwave radiation   S w i n c , also called global radiation G, is the sum of the direct radiation I and the diffuse radiation D.
The reflected shortwave radiation S w r e f is due to the albedo of the surface. Thus, the net shortwave radiation is calculated as:
S w N = ( 1 α )   S w i n c = ( 1 α )   ( I + D )
where:
I is the direct radiation
D = D s + D t is the diffuse radiation, as the sum of the portions that come from the sky
D s and from the adjacent terrain D t .
Direct radiation I:
The calculation of the potential radiation for clear sky conditions (I) in high slopy areas for estimating diffuse radiation uses a calculation of radiation from sky and from adjacent terrain [36,43,45].
The calculation of the potential radiation for clear sky conditions ( I ) , Equation (A2), takes into account the following aspects: scattering, reflection and absorption of the solar beam by the atmosphere’s transmissivity, inverse relation between solar radiation and atmospheric pressure, variation of the path length with the sun’s altitude by means of local zenith angle.
I = I 0   ( r m r ) 2   Ψ ( p P 0 c o s Z ) c o s Z
where:
I 0   is the solar constant (=1368 Wm−2)
r   is the distance between the earth and the sun (m refers to the mean)
Ψ is the transmissivity
p is the atmospheric pressure which depends on elevation h through the barometric formula
p 0 is the standard atmospheric pressure
Z   is the zenith angle
The transmissivities Ψ used in the current simulations are taken from the study of Baigorria et al. (2004) on seasonal transmissivities in Peru and vary between 0.5 in the warm/rainy season and 0.6 in the cold/dry season.
The zenith angles Z are calculated according to (Oke, 1987) cited by Hock (2012) [45] as:
c o s Z = s i n Φ s i n δ + c o s Φ c o s δ c o s h
δ = 23.4 c o s [ 360   ° ( t j + 10 ) / 365 ]
h = 15 ( 12 t )
where:
Φ is the latitude
δ is the solar declination
h is the hour angle
t j is the Julian date
t is the local apparent solar time
In addition, topographic conditions for the calculation of the direct radiation I are taken into account by the model by incorporating the slope and the exposition of the plane in each grid. Thus, the radiation on the slope of a plane is calculated as:
I s l o p e = I c o s θ / c o s Z
c o s θ = c o s   β c o s Z + s i n β s i n Z c o s ( s l o p e )  
c o s =   s i n δ c o s Φ + c o s δ s i n Φ cos h / s i n Z
where:
θ is the angle of incidence between the normal to the slope and the solar beam
β is the slope angle
is the solar azimuth angle
s l o p e is the slope azimuth angle
The calculations of shading, the exposition correction factor and potential direct radiation are done for subintervals. The correction factor ( c o s θ / c o s Z ) is calculated as a weighted mean of all correction factors for all the subintervals within the time step.
Total diffusive radiation D:
The total diffusive radiation D is:
D = D s + D t
and is the sum of the diffusive radiation D s that comes from the sky and D t of that comes from the terrain. D s is computed as:
D s = D 0 S f
where:
D 0 is the diffuse radiation from an unobstructed sky,
Sf is the sky view factor indicating the portion of the visible sky of an element.
These parameters are calculated for each grid cell, using the Relief Visualization Toolbox (RVT) [59].
The diffusive radiation from the surrounding terrain D t is computed as:
D t = α m S w i n c ( 1 S f )
where αm is the mean albedo of the surrounding terrain, and the other variables are as defined. Since D t depends on the local sky view factor of a slope element, it is variable across the grid.
There has been much research and debate on how to estimate or compute D and/or Do (e.g., [36]). When measurements of the global radiation S w i n c at a climate station are available, Hock [36] proposed an empirically corroborated relation of the diffuse radiation D to the global radiation or ( S w i n c ):
D S w i n c = { 0.15 : S w i n c I T o A   0.8   0.929 + 1.134 S w i n c I T o A + 5.111 ( S w i n c I T o A ) 2 + 3.106 ( S w i n c I T o A ) 3 : 0.15 < S w i n c I T o A < 0.8 1 : S w i n c I T o A 0.15
where I T o A is the radiation at the top of the atmosphere (which is equal to the solar constant affected by the instantaneous and the mean solar earth distance and the cosine of the zenith).
Equation (A12) shows clearly how the ratio D/ S w i n c increases with the decrease of S w i n c due to cloud cover and vice versa. Using Equation (A12), the total D for measurements can be computed and, employing D = Ds + Dt, with Dt known by Equation (A9), Ds = D − Dt, i.e., Do = Ds/Sf. With this computed Do at the climate station (which it is invariant for the whole area), D for an arbitrary grid element is estimated by the sum of Equations (A10) and (A11) and using the appropriate sky view factor Sf in these equations.

Appendix A.2. Longwave Net Radiation

The net longwave radiation ( L w N ) is calculated as:
L w N = L w l i n c   + L w s u r f ,
with longwave incoming radiation L w l i n c   calculated as:
L w i n c = L w l i n c s + L w l i n c t
describing the sum of the longwave radiation that comes from the sky, L w l i n c s , and from the adjacent terrain, L w l i n c t .
L w l i n c s from the sky is estimated by:
L w l i n c s = L 0 S f
where:
L 0 is the sky irradiance from an unobstructed sky,
S f is the sky view factor.
L 0 is calculated at the climate station through L w l i n c x   = L w l i n c   (measured) −   L w l i n c t (A14) and then by L 0 = L w l i n c x     / S f (A15). Lo is then used invariant for the whole grids. S f which is different in each grid, shows the portion of the visible sky.
The longwave radiation that comes from the terrain, L w l i n c t , is based on a parameterization Plüss and Ohmura [60], and discussed further by Hock [36], which considers the part of the sky which is obstructed, the air temperature, and the temperature from the emitting surface. More specifically, L w l i n c t is computed as:
L w l i n c t   =   ( 1 S f ) π s r ( L b + a T a + b T s ) ,
where:
L b is the emitted radiance of a black body at 0 °C; 0   ° C   ( = 100.2   Wm 2 sr 1 ) ,
a and b are constants (a = 0.77 Wm−2sr−1°C−1 and b = 0.54 Wm−2sr−1°C−1).
T a and T s are the temperatures of the atmosphere and surface, respectively.
Longwave outgoing radiation emitted by the surface L w s u r f from the surface calculated for each grid which is computed as:
L w s u r f = ε σ T s 4
where:
ε is emissivity of the surface,
σ is the Stefan Boltzmann constant (=5.67 × 10−8 Wm−2K−4), T s is the surface temperature (K).
At infrared wavelengths, snow and ice act as near-perfect radiators, with typical emissivity’s ε between 0.94 and 0.99 [37], and most of the cases are taken approximately to 1.
The parameterization of the surface temperature T s assumes melting conditions at a temperature of 0 °C. Moreover, when the energy balance results in a negative value, the surface temperature is lowered, until the energy balance reaches equilibrium.

References

  1. Ames, A.; Francou, B. Cordillera Blanca: Glaciares En La Historia. Bull. Inst. Fr. des Etudes Andin. 1995, 24, 37–64. [Google Scholar]
  2. Autoridad Nacional del Agua. Inventario de Glaciares Del Perú. Inventario. Nacional de Glaciares y Lagunas; Autoridad Nacional del Agua: Huaraz, Peru, 2014; Available online: http://www.ana.gob.pe/media/981508/glaciares.pdf (accessed on 9 September 2021).
  3. Francou, B.; Rabatel, A.; Soruco, A.; Sicart, J.; Silvestre, E.; Ginot, P.; Cáceres, B.; Condom, T.; Villacís, M.; Ceballos, J.L.; et al. Glaciares de Los Andes Tropicales: Víctimas Del Cambio Climático; Comunidad Andina (CAN): Quito, Ecuador, 2013; p. 100. [Google Scholar]
  4. Rabatel, A.; Francou, B.; Soruco, A.; Gomez, J.; Cáceres, B.; Ceballos, J.L.; Basantes, R.; Vuille, M.; Sicart, J.E.; Huggel, C.; et al. Current State of Glaciers in the Tropical Andes: A Multi-Century Perspective on Glacier Evolution and Climate Change. Cryosphere 2013, 7, 81–102. [Google Scholar] [CrossRef] [Green Version]
  5. Schaner, N.; Voisin, N.; Nijssen, B.; Lettenmaier, D.P. The Contribution of Glacier Melt to Streamflow. Environ. Res. Lett. 2012, 7, 034029. [Google Scholar] [CrossRef]
  6. Mark, B.G.; Seltzer, G.O. Tropical Glacier Meltwater Contribution to Stream Discharge: A Case Study in the Cordillera Blanca, Peru. J. Glaciol. 2003, 49, 271–281. [Google Scholar] [CrossRef] [Green Version]
  7. Villacis, M.; Cadier, E.; Pouyaud, B. Relaciones Hidrológicas Entre El Glaciar y Los Páramos de Los Andes Tropicales Del Ecuador: Su Papel En La Disponibilidad de Recursos Hídricos. In Proceedings of the IV Simposio Internacional Sobre Cambios Globales, La Paz, Bolivia, 30 September–1 October 2010; pp. 1–8. [Google Scholar] [CrossRef]
  8. Schoolmeester, T.; Johansen, K.S.; Alfthan, B.; Baker, E.; Hesping, M.; Verbist, K. Atlas de Glaciares y Aguas Andinos. El Impacto Del Retroceso de Los Glaciares Sobre Los Recursos Hídricos; UNESCO: Paris, France; GRID-Arendal: Arendal, Norway, 2018. [Google Scholar]
  9. Shugar, D.H.; Jacquemart, M.; Shean, D.; Bhushan, S.; Upadhyay, K.; Sattar, A.; Schwanghart, W.; McBride, S.; de Vries, M.V.W.; Mergili, M.; et al. A Massive Rock and Ice Avalanche Caused the 2021 Disaster at Chamoli, Indian Himalaya. Science 2021, 373, 300–306. [Google Scholar] [CrossRef]
  10. Favier, V.; Wagnon, P.; Chazarin, J.P.; Maisincho, L.; Coudrain, A. One-Year Measurements of Surface Heat Budget on the Ablation Zone of Antizana Glacier 15, Ecuadorian Andes. J. Geophys. Res. Atmos. 2004, 109, 1–15. [Google Scholar] [CrossRef] [Green Version]
  11. Favier, V.; Wagnon, P.; Ribstein, P. Glaciers of the Outer and Inner Tropics: A Different Behaviour but a Common Response to Climatic Forcing. Geophys. Res. Lett. 2004, 31, 1–5. [Google Scholar] [CrossRef] [Green Version]
  12. Francou, B. Tropical Climate Change Recorded by a Glacier in the Central Andes during the Last Decades of the Twentieth Century: Chacaltaya, Bolivia, 16°S. J. Geophys. Res. 2003, 108, 4154. [Google Scholar] [CrossRef]
  13. Francou, B.; Vuille, M.; Favier, V.; Cáceres, B. New Evidence for an ENSO Impact on Low-Latitude Glaciers: Antizana 15, Andes of Ecuador, 0°28’S. J. Geophys. Res. Atmos. 2004, 109, 1–17. [Google Scholar] [CrossRef]
  14. Gurgiser, W.; Marzeion, B.; Nicholson, L.; Ortner, M.; Kaser, G. Modeling Energy and Mass Balance of Shallap Glacier, Peru. Cryosphere 2013, 7, 1787–1802. [Google Scholar] [CrossRef] [Green Version]
  15. Juen, I. Glacier Mass Balance and Runoff in the Tropical Cordillera Blanca, Perú. Ph.D. Thesis, University of Innsbruck, Innsbruck, Austria, 2006; pp. 3–195. [Google Scholar] [CrossRef]
  16. Kaser, G.; Georges, C. Changes of the Equilibrium Line Altitude in the Tropical Cordillera Blanca (Perú) between 1930 and 1950 and Their Spatial Variations. Ann. Glaciol. 1997, 24, 344–349. [Google Scholar] [CrossRef] [Green Version]
  17. Kaser, G. Glacier-Climate Interaction at Low-Latitudes. J. Glaciol. 1999, 25, 1–26. [Google Scholar] [CrossRef] [Green Version]
  18. Sicart, J.E.; Wagnon, P.; Ribstein, P. Atmospheric Controls of the Heat Balance of Zongo Glacier (16°S, Bolivia). J. Geophys. Res. D Atmos. 2005, 110, 1–17. [Google Scholar] [CrossRef] [Green Version]
  19. Sicart, J.E.; Hock, R.; Ribstein, P.; Litt, M.; Ramirez, E. Analysis of Seasonal Variations in Mass Balance and Meltwater Discharge of the Tropical Zongo Glacier by Application of a Distributed Energy Balance Model. J. Geophys. Res. Atmos. 2011, 116, 1–18. [Google Scholar] [CrossRef]
  20. Wagnon, P.; Ribstein, P.; Francou, B.; Pouyaud, B. Annual Cycle of Energy Balance of Zongo Glacier, Cordillera Rreal, Bolivia. J. Geophys. Res. 1999, 104, 3907–3923. [Google Scholar] [CrossRef]
  21. Wagnon, P.; Ribstein, P.; Francou, B.; Sicart, J.E. Anomalous Heat and Mass Budget of Glaciar Zongo, Bolivia, during the 1997/98 El Nino Year. J. Glaciol. 2001, 47, 21–28. [Google Scholar] [CrossRef] [Green Version]
  22. Maussion, F.; Gurgiser, W.; Großhauser, M.; Kaser, G.; Marzeion, B. ENSO Influence on Surface Energy and Mass Balance at Shallap Glacier, Cordillera Blanca, Peru. Cryosphere 2015, 9, 1663–1683. [Google Scholar] [CrossRef] [Green Version]
  23. Vuille, M.; Kaser, G.; Juen, I. Glacier Mass Balance Variability in the Cordillera Blanca, Peru and Its Relationship with Climate and the Large-Scale Circulation. Glob. Planet. Chang. 2008, 62, 14–28. [Google Scholar] [CrossRef] [Green Version]
  24. Vergara, W.; Deeb, A.; Leino, I.; Kitoh, A.; Escobar, M. Assessment of the Impacts of Climate Change on Mountain Hydrology: Development of a Methodology through a Case Study in the Andes of Peru; World Bank: Washington, DC, USA, 2011. [Google Scholar]
  25. Condom, T.; Escobar, M.; Purkey, D.; Pouget, J.C.; Suarez, W.; Ramos, C.; Apaestegui, J.; Tacsi, A.; Gomez, J. Simulating the Implications of Glaciers’ Retreat for Water Management: A Case Study in the Rio Santa Basin, Peru. Water Int. 2012, 37, 442–459. [Google Scholar] [CrossRef]
  26. Pierrehumbert, R. Tropical Glacier Retreat. Available online: https://www.realclimate.org/index.php/archives/2005/05/tropical-glacier-retreat/ (accessed on 9 September 2021).
  27. Wagnon, P.; Lafaysse, M.; Lejeune, Y.; Maisincho, L.; Rojas, M.; Chazarin, J.P. Understanding and Modeling the Physical Processes That Govern the Melting of Snow Cover in a Tropical Mountain Environment in Ecuador. J. Geophys. Res. Atmos. 2009, 114, 1–14. [Google Scholar] [CrossRef]
  28. Sicart, J.E.; Ribstein, P.; Wagnon, P.; Brunstein, D. Clear-Sky Albedo Measurements on a Sloping Glacier Surface: A Case Study in the Bolivian Andes. J. Geophys. Res. Atmos. 2001, 106, 31729–31737. [Google Scholar] [CrossRef]
  29. Lejeune, Y.; Wagnon, P.; Bouilloud, L.; Chevallier, P.; Etchevers, P.; Martin, E.; Sicart, J.E.; Habets, F. Melting of Snow Cover in a Tropical Mountain Environment in Bolivia: Processes and Modeling. J. Hydrometeorol. 2007, 8, 922–937. [Google Scholar] [CrossRef]
  30. Winkler, M.; Juen, I.; Mölg, T.; Wagnon, P.; Gómez, J.; Kaser, G. Measured and Modelled Sublimation on the Tropical Glaciar Artesonraju, Perú. Cryosphere 2009, 3, 21–30. [Google Scholar] [CrossRef] [Green Version]
  31. Mölg, T.; Cullen, N.J.; Hardy, D.R.; Kaser, G.; Klok, L. Mass Balance of a Slope Glacier on Kilimanjaro and Its Sensitivity to Climate. Int. J. Climatol. 2008, 28, 881–892. [Google Scholar] [CrossRef]
  32. Instituto Nacional de Investigación en Glaciares y Ecosistemas de Montaña INAIGEM. Reconocimiento de Peligros Naturales En La Laguna Nueva “Artesoncocha Alta”; Instituto Nacional de Investigación en Glaciares y Ecosistemas de Montaña (INAIGEM): Huaraz, Peru, 2016. [Google Scholar]
  33. Kaser, G.; Georges, C. On the Mass Balance of Low Latitude Glaciers with Particular Consideration of the Peruvian Cordillera Blanca. Geogr. Ann. Ser. A Phys. Geogr. 1999, 81, 643–651. [Google Scholar] [CrossRef]
  34. Lozano, M. Energy Budget- and Temperature Index-Based Calculations of Seasonal Mass Balances and Discharge in the Tropical Glacier Artesonraju and Artesoncocha Basin in the Cordillera Blanca, Peru. Ph.D. Thesis, Universität Kassel, Kassel, Germany, 2020. [Google Scholar] [CrossRef]
  35. Zemp, M.; Gärtner-Roer, I.; Nussbaumer, S.; Hüsler, F.; Machguth, H.; Mölg, N.; Paul, F.; Hoelzle, M. Global Glacier Change Bulletin No. 1 (2012–2013); Global Glacier Change Bulletin 1; World Glacier Monitoring Service (WGMS): Zurich, Switzerland, 2015. [Google Scholar] [CrossRef]
  36. Hock, R. Modelling of Glacier Melt and Discharge. Ph.D. Thesis, ETH Zurich, Zurich, Switzerland, 1998. [Google Scholar] [CrossRef]
  37. Cuffey, K.M.; Paterson, W.S.B. The Physics of Glaciers, 4th ed.; Academic Press: Cambridge, MA, USA, 2010. [Google Scholar]
  38. International Association of Classification Societies (IACS). Glossary of Glacier Mass Balance and Related Terms; UNESCO In-ternational Hydrological Programme (UNESCO IHP): Paris, France, 2011; p. 114. [Google Scholar]
  39. Baker, D.; Escher-Vetter, H.; Moser, H.; Oerter, H.; Reinwarth, O. A Glacier Discharge Model Based on Results from Field Studies of Energy Balance, Water Storage and Flow. IAHS Publ. 1982, 138, 103–112. [Google Scholar]
  40. Nash, J.E.; Sutcliffe, J.V. River Flow Forecasting through Conceptual Models Part I—A Discussion of Principles. J. Hydrol. 1970, 10, 282–290. [Google Scholar] [CrossRef]
  41. Krause, P.; Boyle, D.P.; Bäse, F. Comparison of Different Efficiency Criteria for Hydrological Model Assessment. Adv. Geosci. 2005, 5, 89–97. [Google Scholar] [CrossRef] [Green Version]
  42. Cabrera, J. Calibración de Modelos Hidrológicos. Módulo Introductorio de curso de Hidrología. Imefen. Uni. Edu. Pe 2009, 1, 1–7. [Google Scholar]
  43. Hock, R.; Holmgren, B. A Distributed Surface Energy-Balance Model for Complex Topography and Its Application to Storglaciären, Sweden. J. Glaciol. 2005, 51, 25–36. [Google Scholar] [CrossRef] [Green Version]
  44. Hock, R. Glacier Meteorology Energy Balance. Summer School in Glaciology Fairbanks/McCarthy 7–17 June 2010. Available online: https://glaciers.gi.alaska.edu/sites/default/files/mccarthy/Notes_Energybal_Hock.pdf (accessed on 9 September 2021).
  45. Hock, R.; Tijm-Reijmer, C. A Mass-Balance, Glacier Runoff and Multi-Layer Snow Model DEBAM and DETIM Distributed Energy Balance and Distributed Enhanced Temperature Index Model Users Manual. 2012. Available online: https://regine.github.io/meltmodel/ (accessed on 9 September 2021).
  46. Pellicciotti, F.; Carenzo, M.; Helbing, J.; Rimkus, S.; Burlando, P. On the Role of Subsurface Heat Conduction in Glacier Energy Balance Modeling. Ann. Glaciol. 2009, 50, 16–24. [Google Scholar] [CrossRef] [Green Version]
  47. Rabatel, A.; Bermejo, A.; Loarte, E.; Soruco, A.; Gomez, J.; Leonardini, G.; Vincent, C.; Sicart, J.E. Can the Snowline Be Used as an Indicator of the Equilibrium Line and Mass Balance for Glaciers in the Outer Tropics? J. Glaciol. 2012, 58, 1027–1036. [Google Scholar] [CrossRef] [Green Version]
  48. Soruco, Á.; Vincent, C.; Francou, B.; Rabatel, A. Comparación de métodos para estimar el balance de masa del Glaciar de Zongo, Bolivia (16° S, 68° O). Geoacta 2014, 39, 154–164. [Google Scholar]
  49. Null, J. El Niño and La Niña Years and Intensities. Available online: https://ggweather.com/enso/oni.htm (accessed on 12 November 2018).
  50. Arnaud, Y.; Muller, F.; Vuille, M.; Ribstein, P. El Niño-Southern Oscillation (ENSO) Influence on a Sajama Volcano Glacier (Bolivia) from 1963 to 1998 as Seen from Landsat Data and Aerial Photography. J. Geophys. Res. 2001, 106, 17773. [Google Scholar] [CrossRef]
  51. Dai, A. Temperature and Pressure Dependence of the Rain-Snow Phase Transition over Land and Ocean. Geophys. Res. Lett. 2008, 35, 1–7. [Google Scholar] [CrossRef]
  52. Feiccabrino, J.; Graff, W.; Lundberg, A.; Sandström, N.; Gustafsson, D. Meteorological Knowledge Useful for the Improvement of Snow Rain Separation in Surface Based Models. Hydrology 2015, 2, 266–288. [Google Scholar] [CrossRef] [Green Version]
  53. Ye, H.; Cohen, J.; Rawlins, M. Discrimination of Solid from Liquid Precipitation over Northern Eurasia Using Surface Atmospheric Conditions. J. Hydrometeorol. 2013, 14, 1345–1355. [Google Scholar] [CrossRef] [Green Version]
  54. Auer, A.H. The Rain versus Snow Threshold Temperatures. Weatherwise 1974, 27, 67. [Google Scholar] [CrossRef]
  55. Jennings, K.S.; Winchell, T.S.; Livneh, B.; Molotch, N.P. Spatial Variation of the Rain-Snow Temperature Threshold across the Northern Hemisphere. Nat. Commun. 2018, 9, 1–9. [Google Scholar] [CrossRef] [Green Version]
  56. Sicart, J.E.; Hock, R.; Six, D. Glacier Melt, Air Temperature, and Energy Balance in Different Climates: The Bolivian Tropics, the French Alps, and Northern Sweden. J. Geophys. Res. Atmos. 2008, 113, 1–11. [Google Scholar] [CrossRef] [Green Version]
  57. Sugathan, N.; Biju, V.; Renuka, G. Influence of Soil Moisture Content on Surface Albedo and Soil Thermal Parameters at a Tropical Station. J. Earth Syst. Sci. 2014, 123, 1115–1128. [Google Scholar] [CrossRef] [Green Version]
  58. Fischer, L.; Purves, R.S.; Huggel, C.; Noetzli, J.; Haeberli, W. On the Influence of Topographic, Geological and Cryospheric Factors on Rock Avalanches and Rockfalls in High-Mountain Areas. Nat. Hazards Earth Syst. Sci. 2012, 12, 241–254. [Google Scholar] [CrossRef] [Green Version]
  59. Zakšek, K.; Oštir, K.; Kokalj, Ž. Sky-View Factor as a Relief Visualization Technique. Remote Sens. 2011, 3, 398–415. [Google Scholar] [CrossRef] [Green Version]
  60. Plüss, C.; Ohmura, A.; Plüss, C.; Ohmura, A. Longwave Radiation on Snow-Covered Mountainous Surfaces. J. Appl. Meteorol. 1997, 36, 818–824. [Google Scholar] [CrossRef]
Figure 1. Locations of glacier Artesonraju, Artesoncocha basin limits, and Artesoncocha lake.
Figure 1. Locations of glacier Artesonraju, Artesoncocha basin limits, and Artesoncocha lake.
Climate 09 00143 g001
Figure 2. Locations of climate station (CS) used in the study. Satellite image of Artesonraju glacier taken from Google Earth 2019.
Figure 2. Locations of climate station (CS) used in the study. Satellite image of Artesonraju glacier taken from Google Earth 2019.
Climate 09 00143 g002
Figure 3. Simulated and measured yearly (a) (2004–2005), (b) (2005–2006), (c) (2006–2007) discharge at the outlet of Artesoncocha lake.
Figure 3. Simulated and measured yearly (a) (2004–2005), (b) (2005–2006), (c) (2006–2007) discharge at the outlet of Artesoncocha lake.
Climate 09 00143 g003
Figure 4. (a) Measurements in stakes/pits of mass losses/gains and simulated mass losses/gains in the same locations, (b) mass balance profiles for each hydrological year between 2004–2007, (c) distribution of areas in the elevation range.
Figure 4. (a) Measurements in stakes/pits of mass losses/gains and simulated mass losses/gains in the same locations, (b) mass balance profiles for each hydrological year between 2004–2007, (c) distribution of areas in the elevation range.
Climate 09 00143 g004
Figure 5. Seasonally averaged daily means of (a) simulated mass balance, (b) temperature, and (c) precipitation over the simulated 2004–2007 hydrological years.
Figure 5. Seasonally averaged daily means of (a) simulated mass balance, (b) temperature, and (c) precipitation over the simulated 2004–2007 hydrological years.
Climate 09 00143 g005
Figure 6. Seasonal snow and ice albedo sensitivities for the three hydrological years in terms of their E–response surfaces.
Figure 6. Seasonal snow and ice albedo sensitivities for the three hydrological years in terms of their E–response surfaces.
Climate 09 00143 g006
Figure 7. Seasonal threshold temperature rain/snow T0 sensitivity for the three hydrological years in terms of its effect on E. Red points show the points located above 3 °C.
Figure 7. Seasonal threshold temperature rain/snow T0 sensitivity for the three hydrological years in terms of its effect on E. Red points show the points located above 3 °C.
Climate 09 00143 g007
Figure 8. (a) Simulated mean seasonal net radiation and turbulent fluxes on the Artesonraju glacier; (b) simulated mean seasonal contributions of the various radiative flux terms.
Figure 8. (a) Simulated mean seasonal net radiation and turbulent fluxes on the Artesonraju glacier; (b) simulated mean seasonal contributions of the various radiative flux terms.
Climate 09 00143 g008
Figure 9. Vertical profiles of elevation-band averaged energy flux terms for melting, net radiation and turbulent fluxes for the three seasons September–October (SO), November–March (NDJFM), and April–August (AMJJA) for each of the hydrological years 2004–2007.
Figure 9. Vertical profiles of elevation-band averaged energy flux terms for melting, net radiation and turbulent fluxes for the three seasons September–October (SO), November–March (NDJFM), and April–August (AMJJA) for each of the hydrological years 2004–2007.
Climate 09 00143 g009
Figure 10. Similar to Figure 9, but for shortwave, longwave, and net radiation fluxes.
Figure 10. Similar to Figure 9, but for shortwave, longwave, and net radiation fluxes.
Climate 09 00143 g010
Figure 11. Similar to Figure 9, but for sensible and latent heat, and turbulent fluxes.
Figure 11. Similar to Figure 9, but for sensible and latent heat, and turbulent fluxes.
Climate 09 00143 g011
Table 1. Climate, discharge, and mass balance data used in the glacier model.
Table 1. Climate, discharge, and mass balance data used in the glacier model.
DataTime ResolutionLocationElevationAccuracy
TemperatureDaily8°58′11.7″ S
77°38′13.6″ W
4811 m.a.s.l.+/−0.2 °C
PrecipitationDaily8°58′11.6″ S
77°38′14″ W
4836 m.a.s.l.
Relative humidityDaily8°57′55″ S
77°38′12″ W
4838 m.a.s.l.
WindDaily8°57′55″ S
77°38′12″ W
4838 m.a.s.l.
Shortwave radiation incomingEvery 30 min. For simulation, daily averages are used.8°57′55″ S
77°38′12″ W
4838 m.a.s.l.+/−3%
Shortwave radiation reflectedEvery 30 min. For simulation, daily averages are used8°57′55″ S
77°38′12″ W
4838 m.a.s.l.+/−3%
Longwave radiation incomingEvery 30 min. For simulation, daily averages are used.8°57′55″ S
77°38′12″ W
4838 m.a.s.l.+/−3%
DischargeEvery 30 min. For simulation, daily averages are used.8°58′37″ S
77°38′40″ W
Artesoncocha Lake Outlet
4254 m.a.s.l.
Mass balance measurementsBetween 20 and 68 days.
Table 2. DEBAM input variables and methodology of their estimation across the glacier model grid.
Table 2. DEBAM input variables and methodology of their estimation across the glacier model grid.
Data TypeVariableMethod of Estimation and/or Spatial Interpolation
Climate
Data
TemperatureSpatially interpolated using the gradient of temperature 0.65 °C/100 m.
PrecipitationSpatially interpolated using precipitation changes with elevation between 20 and 35% × 100 m−1 up to 4900 m.a.s.l. precipitation change with elevation was between 5 and 10% above 4900 m.a.s.l.
Relative humiditySpatially constant.
Wind SpeedSpatially constant.
Radiation DataShortwave radiation incomingSpatially interpolated according to estimations of direct (I) and diffuse radiation (D).
Direct radiation: estimations and interpolation use, local data of atmospheric pressure, atmospheric transmissivity, geographical location, altitude, slopes, zenith angle, solar declination, azimuth angle and sky view factor (see Equations (A2)–(A8)).
Diffuse radiation: estimations consider measured shortwave incoming radiation at the climate station; the diffuse radiation is also derived using radiation at the top of the atmosphere. These estimations use an empirically corroborated equation of the diffuse radiation D to the shortwave incoming radiation (taken from Hock [36]). The interpolation also assumes a sky view factor for each grid cell (see Equations (A9)–(A12)).
Shortwave radiation reflectedSpatially interpolated using distributed Swinc. Seasonal albedos used as calibration parameters.
Longwave incoming radiationSpatially interpolated using Equations (A13)–(A17). Interpolation uses measured longwave radiation. For interpolation, the sky view factor is used.
Longwave outgoing radiationSpatially interpolated using Equation (A17), the subsurface temperature is derived by iteration closing the Energy Balance Equation (1)
Turbulent fluxesSpatially interpolated using Equations (4) and (5), using relative humidity, wind speed, and atmospheric pressure. Turbulent fluxes (sum of latent and sensible heat) are used to determine whether sublimation, deposition, or fusion occurs.
Hydrological DataDischargeArtesoncocha Lake outlet, used for calibration. Estimation using Equation (10).
Storage constants k used as calibration parameters.
Geographic Grid
Data
ElevationDTM30 m
Glacier and basin areaFrom DTM calculations using ARCSWAT
SlopeInterpolated from DTM using ArcGIS and used for direct radiation
ExpositionInterpolated from DTM using ArcGIS and used for direct radiation
Sky view factorInterpolated From DTM calculated using SVF program used to calculate direct radiation and longwave radiation.
(Snow) Water equivalent Spatially interpolated simulation from measurements with kriging for the first day of simulation. Snow water equivalent changes temporally with precipitation data and threshold rain/snow temperature. Used to allocate albedo and discharge constants according to surface type. Threshold rain/snow parameter T0 used as calibration parameter.
Firn, ice, and debris areasDelimited at the beginning of the simulation of hydrological year, according to the ELA reported by UGRH and WGMS. Used to decide location of firn area when snow has been removed. It also reduces melt over debris areas.
Table 3. Nash Suffcliffe coefficients (E) and root mean square error (RMSE) for discharge and mass balances with quality classification categories of Cabrera (2009) [42].
Table 3. Nash Suffcliffe coefficients (E) and root mean square error (RMSE) for discharge and mass balances with quality classification categories of Cabrera (2009) [42].
YearDischargeMass Balance
ERMSE
(m3s−1)
Adjustment
Quality
ERMSE
(m3s−1)
Adjustment Quality
2004–20050.870.12Very Good0.711.59Good
2005–20060.610.22Good0.582.36Good
2006–20070.600.16Good0.691.9Good
Table 4. Annual simulated and UGRH measured [35] mass balances and ELAs, observed snow line [47] and calculated accumulation-area ratio from simulated data.
Table 4. Annual simulated and UGRH measured [35] mass balances and ELAs, observed snow line [47] and calculated accumulation-area ratio from simulated data.
YearMass Balance (UGRH)
(m.w.e)
Mass Balance
(Simulated)
(m.w.e)
ELA (UGRH)
(m.a.s.l)
ELA (Simulated)
(m.a.s.l)
Snow Line
(Observed)
(m.a.s.l)
AAR
(%)
2004–2005−1.55−1.4250155050492559
2005–2006−1.67−1.3750505021489063
2006–2007−1.52−0.4549865009488065
Table 5. Optimally calibrated seasonal storage constants for the three kinds of glacier surfaces.
Table 5. Optimally calibrated seasonal storage constants for the three kinds of glacier surfaces.
YearSeasonkfirn
(h)
ksnow
(h)
kice
(h)
2004–2005SO700390100
NDJFM50039035
AMJJA1100150250
2005–2006SO70039050
NDJFM80090050
AMJJA1100150300
2006–2007SO800390500
NDJFM80045030
AMJJA1100200250
Table 6. Simulated annually averaged daily energy flux terms in the Artesonraju glacier.
Table 6. Simulated annually averaged daily energy flux terms in the Artesonraju glacier.
YearEnergy Flux Components (Wm−2)
SwincSwrefSwbalLwincLwsurfLwbalQNetQHSQHLTf
2004–2005219−15167278−301−23448−70
2005–2006222−15665281−304−23438−80
2006–2007216−16254285−305−2035−41
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Lozano Gacha, M.F.; Koch, M. Distributed Energy Balance Flux Modelling of Mass Balances in the Artesonraju Glacier and Discharge in the Basin of Artesoncocha, Cordillera Blanca, Peru. Climate 2021, 9, 143. https://0-doi-org.brum.beds.ac.uk/10.3390/cli9090143

AMA Style

Lozano Gacha MF, Koch M. Distributed Energy Balance Flux Modelling of Mass Balances in the Artesonraju Glacier and Discharge in the Basin of Artesoncocha, Cordillera Blanca, Peru. Climate. 2021; 9(9):143. https://0-doi-org.brum.beds.ac.uk/10.3390/cli9090143

Chicago/Turabian Style

Lozano Gacha, María Fernanda, and Manfred Koch. 2021. "Distributed Energy Balance Flux Modelling of Mass Balances in the Artesonraju Glacier and Discharge in the Basin of Artesoncocha, Cordillera Blanca, Peru" Climate 9, no. 9: 143. https://0-doi-org.brum.beds.ac.uk/10.3390/cli9090143

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop