Modeling energy and mass balance of Shallap Glacier , Peru

We calculated the distributed surface mass and energy balance of Shallap Glacier, Cordillera Blanca, Peru (9 S, 77 W, 4700–5700 m a.s.l., ∼ 7 km2), on hourly time steps for two years (September 2006–August 2008) using a process-based model and meteorological measurements as input. Model parameter combinations were optimized against 21 temporal readings of 20 stakes in the ablation zone of the glacier. Uncertainty caused by model input parameters and parameterization schemes was estimated using a leaveone out cross-validation scheme, which yields values of root mean square deviation (RMSD) of surface height change < 1 m (< 10 % of the measured amplitude) for all stakes. With the best parameter combination (smallest RMSD) applied, the modeled annual surface mass balance of the glacier was −0.32± 0.4 m w.e. (water equivalent) for September 2006–August 2007 and 0.51 ± 0.56 m w.e. for September 2007–August 2008. While the mass balance above 5000 m was similar in both years ( 10.33± 0.68 m w.e.) due to similar annual sums of solid precipitation, a difference of 1.97± 0.68 m w.e. was calculated for the lower parts of the glacier. This difference is associated with more frequent occurrence of higher snow line altitudes during the first year, which was mainly caused by a higher fraction of liquid precipitation due to higher mean air temperatures. As the net shortwave budget was found to be the main source for ablation throughout the year at Shallap Glacier, lower surface albedo especially caused by lower solid precipitation amounts explains most of the difference in modeled ablation and mass balance between the two years.


Introduction
In the semi-arid outer-tropical climate of the Cordillera Blanca range of the Peruvian Andes, glacier meltwater is cru-cial to sustaining regional runoff during the dry season from May to September (Juen et al., 2007;Kaser et al., 2010).Andean glaciers, and the hydrological stores they represent, have been shrinking substantially over recent decades (Rabatel et al., 2013;Vuille et al., 2008;Mark and Seltzer, 2003;Kaser and Georges, 1999;Ames, 1998).Investigations based on measured hydrographs (Baraer et al., 2012) and semiempirical models (Juen et al., 2007) suggest glacier retreat is changing the future runoff regime, resulting in seasonal decreases of river discharge of up to 30 % if the glaciers are lost entirely.Concurrently, growth in population, as well as hydropower, agriculture and tourism activities, increases the hydrological demand and vulnerability of local communities to changes in year-round water availability (Bury et al., 2010;Mark and Seltzer, 2003).Therefore, there is a critical need to develop reliable projections of the glacier contribution to runoff in the Cordillera Blanca in order to inform regional planning of hydrological resource usage.
A crucial first step in determining glacier runoff under changing climatic conditions is to investigate glacier mass balance and the processes through which climate conditions drive glacier response.Where sufficiently high quality input data are available, this is best done through application of process-based surface energy and mass balance models with high temporal and spatial resolution (e.g., Mölg et al., 2008Mölg et al., , 2009aMölg et al., , 2012;;MacDougall et al., 2011;Machguth et al., 2009;Hock, 2005;Klok and Oerlemans, 2002), which provide crucial additional information to empirical temperature index (e.g., Carenzo et al., 2009) or semi-empirical energy balance models (e.g., Kaser, 2001;Juen et al., 2007).
Mass and energy balance studies in the tropical Andes have been carried out within the inner tropics at Antisana Glacier, Ecuador (0 • 28 S), and in the outer tropics at Zongo Glacier, Bolivia (16 • S).At Antisana Glacier, point energy balance studies in the ablation zone show that net shortwave radiation is closely coupled to surface ablation throughout the year (Favier et al., 2004a) and variations in glacier-wide ablation are controlled by surface albedo and the position of the snow line, which are both governed by the occurrence and spatial distribution of solid precipitation (Francou et al., 2004).As air temperature controls the phase of precipitation, the sensitivity of the annual mass balance to both precipitation amount and temperature is high.Seasonal variations in energy fluxes are small and related to a slight seasonality in wind speeds.In the outer tropics at Zongo Glacier there is a clear seasonality in meteorological conditions dominated by changes in atmospheric humidity and resultant cloud cover conditions and atmospheric long wave emissions.While the surface energy balance is still dominated by the net shortwave radiation, low atmospheric humidity and atmospheric long wave emissions during the dry season result in high energy losses via the net long wave flux, which reduces the energy available for melting the glacier and consequently restricts ablation to sublimation (e.g., Wagnon et al., 1999;Favier et al., 2004a;Mölg and Hardy, 2004;Nicholson et al., 2012).In contrast to Antisana Glacier, the sensitivity of Zongo Glacier mass balance to air temperature was found to be low because air temperatures at this site are sufficiently low that precipitation over the whole glacier area was solid (Favier et al., 2004a).Instead, glacier mass balance here is most strongly controlled by atmospheric humidity, cloudiness and precipitation, especially at the onset of the humid season (Sicart et al., 2011).Kaser (2001) generally suggests that the mass balance in the inner tropics is most sensitive to variations in air temperature, whereas in the outer tropics sensitivity to atmospheric moisture is highest due to its multiple impacts on mass balance.
The Cordillera Blanca of Peru lies between the latitudes of Ecuador (inner Tropics) and Bolivia (outer Tropics).As no high spatial and temporal resolution mass and energy balance investigations have been carried out in this region so far, it is not clear how the seasonal and annual characteristics of glacier mass and energy balance relate to those identified in Ecuador and Bolivia.Furthermore, the atmospheric variables that exert the strongest control on glacier mass balance, which is crucial information for estimating the possible impact of climate change on the glaciers of the Cordillera Blanca, remain unknown.Thus, the aims of this study are to (1) use a process-based model to compute and compare the characteristics of glacier mass balance of Shallap Glacier (9 • S), Cordillera Blanca, Peru, for two hydrological years, and (2) investigate the energetic drivers of the surface mass balance.

Study area, measurements and methods
Shallap Glacier (9 • 20 S, −77 • 20 W) lies west of the divide of the Cordillera Blanca range of the Peruvian Andes (Fig. 1).The outer tropical climate is characterized by relatively con-stant annual temperature and a distinct variation in precipitation (e.g., Juen et al., 2007;Favier et al., 2004b;Kaser, 2001) with a wet season from October to April and a dry season from May to September (see Supplement, Fig. S1).On the basis of satellite data (SPOT XS Sensor) obtained in 2002 (Georges, 2004), the glacier at that time spanned 4700 to 5740 m a.s.l., with more than half of the 7.035 km 2 glacier surface area below 5200 m.In the absence of more recent information, the 50 m resolution digital elevation model (DEM) derived from these satellite images was used as the glacier surface for our mass balance modeling.

Automatic weather stations
In 2002, an automatic weather station was installed on the southern glacier moraine (4950 m a.s.l.) approximately 150 m above the glacier surface (hereafter AWS M ) (Fig. 1).Continuous hourly data of temperature, humidity, wind speed, global radiation and precipitation are available from January 2006 to November 2009 and from July 2011 to March 2012.Data from AWS M were used as input for the surface energy and mass balance model applied in this study (Sect.2.2).From July 2010 to September 2012 a second weather station (hereafter AWS G ) was operated in the ablation zone of the glacier (4796 m a.s.l.), recording hourly temperature, humidity, wind direction and speed and net radiation at 2 m above the glacier surface.No data are available from February 2012 to April 2012 due to a defect in power supply.Data from AWS G were used to improve the surface albedo module of the mass balance model (Sect.2.2) and to develop a function for transferring temperature from the moraine to the glacier surface.The temperature transfer function is based on a multiple regression considering temperature, relative humidity and wind speed measured at AWS M (Gurgiser et al. 2013).Figure S2 (Supplement), showing measured and transferred temperature, and information on the meteorological sensors used at both AWS are included in the supplement.

Stake measurements
Stake measurements (between 4758 and 4824 m a.s.l) in the ablation zone of Shallap Glacier were conducted by the Unidad de Glaciologia y Recursos Hidricos (UGRH) of the Peruvian Autoridad Nacional de Agua (ANA) from August 2005 to August 2008.At 20 stake locations (Fig. 1), relative surface height change was measured over 21 intervals spanning 14 to 64 days (dates are shown in Fig. 4).As neither snow depth nor snow density were recorded regularly, measurements cannot be readily transferred into mass units but as the stake measurement area was snow free on 29 August 2006, 21 August 2007, and 27 August 2008, annual mass balances can be derived between these dates, assuming a constant density of ice (900 kg m −3 ).The available stake measurements dictated the duration of our investigation period, and consequently in order to have two complete years of mass balance we define each mass balance year as 1 September-31 August, which is slightly offset from the traditional hydrological year for the region, defined as 1 October-30 September.
All available surface height measurements at all of the 20 stakes were used for model evaluation.Stake elevations were measured by the UGRH in 2006 and 2007.The average of these point elevations at each stake, which show a mean absolute difference of 10 m compared to the framing DEM grid points, were used for model optimization and evaluation runs (Sects.2.2.2 and 2.2.3).

Model design
Glacier surface energy and mass balance was investigated using a process-based model (Mölg et al., 2008(Mölg et al., , 2009b(Mölg et al., , 2012) ) that has been previously applied to Shallap Glacier (Gurgiser et al., 2013).Here we provide only a brief overview of the model as it is fully described in these prior publications.The model was run at hourly time steps from 1 June 2006 to 31 August 2008 (to include three months of spin-up which proved to be sufficient to produce reliable snowpack patterns) using inputs of temperature, humidity, wind speed, global radiation and all-phase precipitation measured at AWS M .Temperature is modified from AWS M to represent conditions above the glacier surface using a transfer function (Sect.2.1.1).The reference altitude of meteorological model input is 4800 m.Realistic value ranges for 18 free model parameters were obtained from the literature or field measurements at Shallap Glacier and those parameters to which the model showed high sensitivity were optimized simultaneously (Sect.2.2.2).
Within the model the mass balance is calculated as the difference between the sum of accumulation terms (solid precipitation, surface water deposition and refreezing of liquid water in the snowpack) and the sum of ablation terms (surface sublimation/evaporation and surface and subsurface melt).Surface melt and sublimation are derived from the surface energy balance (all terms in W m −2 ) formulated as where SWin is the incoming shortwave radiation, α the surface albedo, LWin and LWout are the incoming and outgoing long wave radiation, QS and QL are the turbulent sensible and latent heat fluxes, and QG is energy flux into the glacier body.F is the residual energy flux.If F is positive www.the-cryosphere.net/7/1787/2013/The Cryosphere, 7, 1787-1802, 2013 and surface temperature is 273.15K, the melt energy QM equals F .Otherwise F warms or cools a defined surface layer (Mölg et al., 2009a).
As described in detail in Gurgiser et al. (2013), measured SWin at AWS M is separated into direct and diffuse components, following Hock and Holmgren (2005).Then SWin(x, y) is calculated as a function of topographic shading, slope and aspect.LWin is calculated as a function of air temperature, water vapor and cloud cover, following a parameterization scheme developed by Sicart et al. (2010) for Zongo Glacier, whereas the parameter values were adapted for the local conditions; the skill of the parameterization was demonstrated to be similar at our site (Gurgiser et al., 2013).LWout is derived from modeled surface temperature and the Stefan Boltzmann law assuming emissivity to be unity.The turbulent heat fluxes QS and QL are calculated using the bulk method and published values of surface roughness lengths for snow and ice (Table A1).The effect of stable conditions on turbulent exchange is accounted for following Braithwaite (1995, Eq. 11).The vertical air temperature gradient is assumed to be −0.55 • (100 m) −1 , based on measurement at the Zongo Glacier (Sicart et al., 2011).In the absence of information for Shallap Glacier, we again followed Sicart et al. (2011) by assuming that precipitation does not vary with elevation, as discussed in Sect.4.3.The solid fraction of measured all-phase precipitation, which is derived from linearly interpolating between two threshold temperatures (Table A1), is converted into a snow height using an optimized value for fresh snow density and correcting the snow height for the surface slope of the glacier grid point.Within the model subsurface, the thermodynamic energy equation is solved numerically on a vertical grid with 14 layers (from 0.09 to 3 m total depth).Ice temperatures in the subsurface are unknown, but the temperature below the lowest layer is assumed to be 273 K as measurements at AWS G suggest that the surface temperature is close to the melting point throughout the year, and this is expected to be the best guess for most of the glacier area because of very high SWin throughout the year.QG is the sum of QPS, the fraction of net shortwave radiation that penetrates into ice or snow, and the conductive heat flux at the glacier surface, QC.The model accounts for compaction of snow and refreezing of liquid water, as described in Mölg et al. (2012).Due to its small magnitude compared to the other energy fluxes, the heat flux from precipitation is not considered in the selected model setup, nor is mass redistribution due to wind and gravity.
The surface albedo parameterization is based on Oerlemans and Knap (1998): where α is the surface albedo with subscripts "s" and "i" referring to snow and ice, respectively.d is the actual snow height and d * the albedo depth scale smoothing the transition between snow and ice albedo when d approaches d * .
The surface albedo of snow used within Eq. ( 2) is calculated as where subscripts "os" and "fs" refer to old snow and fresh snow respectively, t is the time elapsed since the last snow fall event and t * is the characteristic timescale (hereafter albedo timescale) over which albedo transitions exponentially from that of fresh to old snow (Oerlemans and Knap, 1998).Within this scheme, a snow fall event is defined as when snowfall accumulated over consecutive hours exceeds 1 cm.The findings of Gurgiser et al. (2013) suggest model skill at Shallap Glacier, where frequent light snowfall on aging snowpacks are common, could be improved by accounting for the albedo of underlying snow (Machguth et al., 2008).Therefore the albedo scheme was modified so that if the snow from a single event ablates away before reaching the albedo of old snow, α s is set back to the value preceding the snow fall event.
The albedo timescale, t * , required in the aging function was determined by applying a fixed value of d * (1 cm) and running the model for the location of AWS G , driven by the data from AWS M , and then comparing the model performance against measured surface albedo at AWS G over a 6 month period from 1 August 2012 to 31 January 2013.Highest correlations and lowest root mean square errors were found for t * between 1 and 4 days, which reflects the typically fast decay of snow albedo when ablation occurs throughout the measurement period (Mölg et al., 2008;Kaser, 2001).To account for a slower decrease of snow albedo in higher altitudes due to less surface melt, we applied a vertical gradient of the albedo timescale of 1 d (100 m) −1 (discussed in Sect.4.3).The skill of this albedo parameterization using a median t * value of 2.5 days (Fig. 2) is comparable to other schemes (Brock et al., 2000).

Values for free model parameters
The mass balance model (Sect.2.2.1) was optimized using a Monte Carlo approach (e.g., Mölg et al., 2012).Firstly, a single parameter sensitivity test was carried out to determine to which of the 18 free model parameters modeled surface mass balance shows greatest sensitivity.This was estimated by computing, for each of the 20 stake locations, a standard model run using median values for all parameters listed in Table A1, and then two further model runs using the upper/lower range limits (Table A1) while keeping all other parameters constant.The difference of the cumulative mass balance at the end of two years (September 2006to August 2008) between the standard run and the max/min perturbation runs for each point and parameter is an indication of model sensitivity to the parameter (Fig. 3).Although the influence of some parameters is interdependent, this single parameter sensitivity was used to identify the most crucial free model parameters within the considered value range as those that change the cumulative mass balance at the end of the time series by more than ±10 % (using the median of the values for the 20 stake locations).This procedure identified 7 key parameters for subsequent optimization (Fig. 3; Table A1): (1) the temperature range over which precipitation phase changes from 100 % solid precipitation to 100 % liquid precipitation; (2) the density of solid precipitation that affects snow height and surface albedo in the case of very thin snow layers that are within the influence of the albedo depth scale (see Oerlemans and Knap, 1998); (3) the fraction of net shortwave energy in the snow layer that influences the distribution of energy to the surface or subsurface; (4) the ice, (5) old snow, and (6) fresh snow albedo; and (7) the albedo timescale.
These 7 parameters were optimized by allowing them to vary randomly within their specified ranges over 1000 model runs in which the other parameters were held constant at their median values and the optimum parameter combination was that with the minimum root mean square deviation (RMSD) between measured and modeled surface height change for all stake points.Parameter values of this run, hereafter referred to as "best parameter combination", are given in Table A1.Following Machguth et al. (2008), the evolution of the standard deviation, and the standard deviation of the standard deviation, of modeled mass balance were used to evaluate the number of simulations required for the optimization and, as fluctuations in both variables approach zero well below 1000 simulations, we can assume that the upper and lower limits of the possible solution space are reached and the sensitivity test yields stable results.

Model uncertainty evaluation methods
Model parameter transferability in space and time was estimated with a leave-one-out cross validation approach (e.g., Wilks, 2011;Hofer et al., 2010) calculated on the 1000 model runs generated for the optimization procedure.To investigate model parameter uncertainty in space, we first calculated the root mean square deviation (RMSD) between measured and modeled surface height change (SHC) for each point as where r is the number of the respective model run (1000 in total), i is the stake location (20 in total), p is the measurement period (21 in total) and n p is the number of periods (21).Then we determined a RMSD for the stake location x by selecting the model run r that produced the smallest mean RMSD at all the other locations and calculating the RMSD of that run at the location x.Thus, this RMSD is obtained independently from measurements at location x and gives an estimate of the model error when the model, using the parameter www.the-cryosphere.net/7/1787/2013/The Cryosphere, 7, 1787-1802, 2013 values optimized at the stake locations, is applied to a location without measurements (Supplement Fig. S3a).This was repeated for all 20 stake locations.Even though we do not model the glacier mass balance beyond the stake measurement period (Sect.2.1.2),we have also evaluated model transferability in time.This was done in a manner analogous to the determination of transferability in space described above, only that measured time intervals rather than measured locations were left out.Consequently, the error we calculated for a period y is the one that can be expected if the model parameters optimized during the availability of measurements are applied at a time when no measurements are available (Supplement Fig. S3b).
For modeling the distributed mass balance of the whole glacier, we selected the run r where the mean RMSD of all stake locations is minimal (Sect.2.2.2).However, the stake measurement area only covers a small part of the glacier (Fig. 1), and we also apply our model outside of the evaluation zone for calculating the spatially distributed mass balance of the whole glacier.Therefore, we investigated whether or not the RMSD depends on the distance between the stake locations, i.e., whether there is a systematic spatial pattern in the error magnitudes.A correlation between this distance and the RMSD would indicate that we would have to be aware of growing model uncertainty outside of the evaluation zone.For completeness, we also investigated whether or not there is a correlation with the temporal distance, which would mean that our best parameter combination yields inferior results when applied before or after the measurement period.
Model error for the annual specific mass balance of the whole glacier for each year was calculated as where CMB is the cumulated mass balance at the end of the year, which is obtained by converting measurements of annual surface height change into mass units (see Sect. 2.1.2) and directly as an output from the mass balance model, respectively.n i is the number of stake locations.Uncertainty in the difference in mass balance between the two years was calculated as where n y is the number of considered years.
There are several other sources of uncertainty in our model results that could arise from uncertainties in spatial gradients of temperature, precipitation and albedo timescale (Sect.4.3), topographic features outside the stake area extracted from the DEM, or currently unconsidered processes such as redistribution of snow and ice and modified ablation due to debris cover.Further measurements and studies will be need to thoroughly investigate and quantify uncertainty arising from such sources or to implement further processes in the model.

Model parameter transferability in space and time
RMSD values obtained for the spatial model parameter transfer (Sect.2.2.3, Supplement Fig. S3a) are all < 1 m, with a mean value of 0.5 m.This is a markedly better performance than obtained in previous modeling (Gurgiser et al., 2013), indicating that the modifications to the albedo scheme described in Sect.2.2.1 improve model performance.RMSD values are < 10 % of the cumulated surface height change at the end of the 2 yr at all points.The time series of surface height change at each stake, modeled using the best parameter combination, optimized at all other stake locations, indicates that although the modeled amplitude is generally slightly smaller than the measured one, the variability and the total surface height change are captured quite well (Fig. 4).This is also the case for annual time steps (Fig. 5).In the first year the model tends to slightly underestimate the surface lowering at the lower stakes and overestimate it at the higher stakes.In the second year, when surface lowering was less, this tendency is not evident, but the correlation between measured and modeled surface height change is slightly worse.
Parameter transfer in time (Supplement Fig. S3b) causes a maximum error of 1.3 cm day −1 in period 9 and a mean value over all periods of 0.6 cm day −1 .As the distribution of RMSD at all periods is random and does not show trends over consecutive periods (Supplement Fig. S3b), the model performance does not favor any particular part of the time series.

Skill of best parameter combination for spatial distributed model run
The best model parameter combination (Sect.2.2.3 and Table A1) yields mean RMSD of 0.6 m (or 4 % of the surface height change amplitude) between modeled and measured surface height change (calculated with Eq. 4) at the 20 stake locations for all periods.The mean RMSD between measured annual mass balance and that modeled with best parameter combination (following Eq. 5) at the 20 stake locations at the end of each year is 0.40 m w.e.(water equivalent) (r = 0.97) for year 1 and 0.56 m w.e.(r = 0.81) for year 2, which is small with respect to the all-stake mean of measured amplitudes of −6.8 m w.e.(year 1) and −3.8 m w.e.(year 2).We found no significant correlation (p < 0.05, |r| < 0.13) between RMSD and location or time when the best parameter combination was applied; therefore the computed model uncertainty is not expected to increase when the model is employed outside the stake area, or before and after stake measurements.

Characteristics of surface mass balance
The spatially distributed surface mass balance was calculated for all glacier grid points using the best model parameter combination for the period September 2006 to August 2008.The specific mass balance was negative (−0.32 m w.e.) in the first year and positive (+0.51 m w.e.) in the second year (Table 1, Fig. 6a and b).Spatially averaged values of modeled equilibrium line altitudes (ELA) of 4985 m a.s.l. and 4953 m in years 1 and 2 respectively (Table 1), and accumulation area ratio (AAR) of 0.70 and 0.74 in year 1 and year 2 respectively are relatively similar in both years, despite the specific mass balance being markedly different (Table 1).Specific mass Table 1.Summary of glacier mass balance characteristics for the two hydrological years and difference between year 1 and year 2. The annual model uncertainties are calculated with Eq. ( 5) (yearly values) and Eq. ( 6) ( Y1-Y2) as described in Sect.2.2.3.ELA means equilibrium line altitude, AAR is the accumulation area ratio and MB is the mass balance.balance (Table 1, Fig. 6c) above the ELA was similar for both years ( 0.33 m w.e.), coinciding with similar annual sums ( 0.118 m w.e.) of precipitation (Table 2), but the differences below the ELA were large ( 1.97 m w.e.).These patterns are also reflected in the differences in the vertical mass balance gradients between the two years (Fig. 7).Generally, the gradients were similar in both years (Table 1 and Fig. 7) and showed some deviations around the ELA, reflecting differences in the persistence of snow cover in this zone between the two years (Fig. 8a).In the first year the vertical mass balance gradient tends to zero above 5400 m, while in the second year this occurs slightly lower at 5200 m, which reflects the effective limit of melt over the glacier surface within our model.The slight kink in the gradient below 4800 m is a result of the model structure (Sect.2.2.1) as we set a lower limit for the albedo timescale (1 d).
Precipitation was high during the wet season from October to April in both years (Fig. S1) with precipitation sums of 1.40 m w.e. and 1.50 m w.e. in year 1 and year 2, but small precipitation events also occurred in the dry season from May to September with total precipitation sums of 0.17 m w.e. and 0.18 m w.e.(Table 2).During the wet season the glacier was regularly entirely snow covered while the snow line altitude generally increases progressively during the dry season (Fig. 8a).Above 5000 m ablation rates were generally small.Below 5000 m ablation occurred during both seasons with varying magnitude, and no clear seasonality.Within both dry seasons, individual accumulation events caused the glacier to become snow covered, which coincides with dramatically reduced or no ablation for a time (Fig. 8c).
Interannual differences in both seasons contribute to the contrasting annual mass balances of the two years (Table 2).The first year experienced a period of intense melt (∼ 50 mm w.e.d −1 ) during February of the wet season that coincided with a higher snow line than was experienced during the wet season in the second year, and also in the dry season of the first year the snow line rose higher than in the second year.Similarly the frequency and duration of reduced ablation were greater in the second year than in the first year, with the contrast in ablation in the months of January-March being particularly marked.

Characteristics of surface energy balance
The drivers of the surface mass balance were analyzed on the basis of the surface energy fluxes modeled with the best parameter combination.As is the case with the interannual mass balance, daily means of the energy balance terms and surface albedo averaged over 10 m elevation intervals (Fig. 9) show greatest diurnal variability on the glacier tongue, so energy fluxes averaged over the glacier surface area below 5000 m are examined in more detail (Fig. 10; Table 2).
Extraterrestrial solar radiation is highest (Supplement Fig. S1) in the wet season, but increased cloud cover means there is little difference in incoming solar radiation between the two seasons.The highest daily mean net shortwave flux (> 200 W m −2 ) occurred during a clear sky interval during Table 2. Seasonal mean values of SWin (W m −2 ), the energy balance terms (all W m −2 ), precipitation (mm w.e.), temperature ( • C), snow line altitude SLA (m), wind speed ws (m s −1 ) and relative humidity RH (%) with differences from year to year for all area below 5000 m.WS is the wet season (October-April) and DS is the dry season (May-September).QM S is the melt energy in the subsurface.P T is all phase precipitation and P L is the fraction of liquid precipitation from P T .T P is the mean temperature during precipitation.the wet season of the first year (Figs.9a and S1), but these high values were restricted to a short time interval and to the lowest elevations of the glacier where albedo was also low-est (Fig. 9f).Instead, net shortwave radiation in the lower parts of the glacier is generally higher during the dry season (Table 2, Fig. 10) because of a combination of lower surface albedo and higher SWin (Table 2, Supplement Fig. S1).

Parameter
In general, net shortwave energy decreased quickly with altitude from the glacier terminus (Fig. 11) until just above the mean ELA (∼ 5000 m) where the glacier surface was always above the snow line (Fig. 8a) and thus modeled albedo (Fig. 9f) was consistently high.In contrast, the net long wave radiation budget varies little with altitude but shows a clear seasonality over the whole glacier surface.Mean seasonal values over the lower glacier range between −25 and −30 W m −2 during the wet season and are around −50 W m −2 during the dry season when maximum mean daily energy losses of up to −100 W m −2 can occur on the clearest days.In seasonal and spatial averages, energy loss via net long wave radiation offsets only between a third and a half of the energy gained via net shortwave radiation (Table 2).Seasonally averaged daily means of sensible heat flux averaged over the lower glacier are small during the wet season (∼ 10 W m −2 ) and enhanced in the dry season (∼ 20 W m −2 ), during which multiple short periods of enhanced wind speeds and temperature gradients between air and the surface occur.Some peaks in QS of almost 100 W m −2 were modeled above the snow line, but these high values were not relevant for melt as they occur only when net shortwave radiation is low due to the high albedo of the snow-covered surface.Also, latent and long wave net fluxes are strongly negative and acting to lower the surface temperature of the glacier, thus increasing the temperature gradient between air and surface ,which forces high sensible heat flux.Seasonally and spatially averaged latent heat flux over the lower glacier is very small during the wet season (∼ −5 W m −2 ) and increases during the dry season (∼ −25 W m −2 ) due to numerous stronger sublimation events (Figs. 9 and 10) during days with low moisture and enhanced wind speeds.The seasonally averaged difference in QL is driven by seasonal differences in both wind speed and nearsurface vapor pressure gradient, which are ∼ 1 m s −1 and ∼ 0.25-0.35hPa higher in the dry season compared to the wet season.Sublimation and evaporation accounted for 4 and 8 % of the wet season ablation in the first and second years, respectively, and 17 and 22 % in the dry season, whereas the variation in sublimation fraction is dominated by changes in the melting rather than the sublimation.The sum of the turbulent fluxes is minor compared to net all-wave radiation because positive QS is offset by negative QL.The spatially averaged sum of mean turbulent fluxes is positive in the wet season, but higher sublimation rates in the dry season more than offset the increased QS so that the summed turbulent fluxes become negative.
Ground heat flux was always negative and small averaged over the whole glacier (−3 W m −2 ).However, in the lowest parts of the glacier the maximal magnitudes of ground heat flux were similar to the turbulent fluxes (Fig. 9e).As a higher fraction of the shortwave radiation absorbed at the surface penetrates into ice than snow (Table A1), QPS is most negative where the glacier surface is snow free.Thus, the higher QG at lower elevation is driven by higher QPS where snow cover is less frequent and consequently lower albedo (Figs.8a and 9f) increases net shortwave flux.The energy lost from the surface through QG was a source for melt energy in the subsurface (Figs. 10 and 11).
Daily mean energy flux values averaged over the glacier area below 5000 m highlight the high variability of net shortwave energy and its control on total melt energy over the glacier ablation zone (Fig. 10).Thus, high melt rates in the ablation zone of Shallap Glacier can be seen to be contingent on high absorption of solar radiation.

Discussion
The magnitudes of the calculated specific mass balance are comparable to mass balance estimates obtained for the same period from measurements available at two other glaciers in the Cordillera Blanca (Artesonjaru, Yanamarey) and within the range of reconstructed annual mass balances for several catchments in the Cordillera Blanca (Rabatel et al., 2013).The modeled ELAs are close to those expected from conceptual approaches (Kaser and Osmaston, 2002), and the high values of the modeled AAR compared to those typical for mid-latitude glaciers (Porter, 1975) are also consistent with values expected for tropical glaciers (e.g., Kaser and Georges, 1997).The calculated vertical mass balance gradients (Table 1) for Shallap Glacier (Fig. 7) were high in the lower parts of the glacier (∼ 3 m w.e.(100 m) −1 ) in both years and comparable to the gradient of the lower Zongo Glacier (Sicart et al., 2007).Such high vertical mass balance gradients on the glacier tongue can be explained by frequent changes in snow cover in the ablation zone throughout the year (Sicart et al., 2007;Kaser and Osmaston, 2002).Above the equilibrium line the vertical mass balance gradients were very small (< 0.2 m w.e.(100 m) −1 ).This form of vertical mass balance and the occurrence of year-round ablation drives the high mass turnover that characterizes tropical glaciers (Kaser, 2001).Thus the modeled mass balance for Shallap Glacier captures all the expected features of a tropical glacier.

Topographic signature in mass and energy balance
The similarity of the AARs between a year with negative and positive mass balance at Shallap Glacier (Sect.3.2, Table 1) is strongly influenced by the local topography (Kaser and Osmaston, 2002).Above 5000 m the glacier slope increases (Fig. 1) and parts of the glacier above 5000 m are exposed to the southwest.This combination of high slope angle and aspect result in slightly less incident solar radiation (up to −10 W m −2 over the whole study period) in higher altitudes compared to the flatter lower ablation zone.Furthermore, the slower decrease of snow albedo with altitude as set in our model structure (Sects.2.2.1 and 4.3) reduces the amount of shortwave energy absorbed at higher elevations.
The steep surface slope above ∼ 5000 m (Fig. 1), and resultant rapid decrease in air temperature upglacier, means that the zone of transition between areas with exclusively solid and areas with solid and liquid precipitation phase is very narrow.Thus, the transition from ablation to accumulation area is abrupt at Shallap Glacier (Figs. 6 and 7).Finally, the steep slope just above the ELA means that any increase in ELA will have a small impact on AAR.Taken together, these factors imply that the sensitivity of the AARs to increasing temperatures is expected to be relatively low at Shallap Glacier.

Energetic drivers of ablation
At Shallap Glacier, very high ablation can occur in the wet season.For example, the highest modeled daily ablation (∼ 50 mm w.e.d −1 ) occurred at the terminus during February of the wet season of year 1.The high ablation rates were a result of low albedo values coinciding with very high val-ues of incoming solar radiation due to clear skies coinciding with peak solar radiation (Supplement Fig. S1).However, as both surface and subsurface melt tend to peak in dry season conditions, overall, seasonal mean melt energy was similar or higher in the dry season than in the wet season due to a more positive net radiation budget (Table 2).This is in contrast to findings at Zongo Glacier where melt rates are small or zero in the dry season due to strongly negative long wave net radiation, which offset shortwave energy gains (Sicart et al., 2011).Over the ablation zone at Shallap Glacier, even during periods of peak net long wave energy losses (∼ −90 W m −2 ) and strong net turbulent energy losses (∼ −40 W m −2 ) during the dry season, for example at the end of June in year 1 (Fig. 10), high net shortwave radiation (∼ 160 W m −2 ) means that surplus energy for melting remains (∼ 30 W m −2 ).This also holds when we consider the whole glacier surface.In general, melting is mostly inhibited by reduced net shortwave and only rarely by strong long wave emission and sublimation, which consume the available energy for ablation (e.g., end of June in the second year, Fig. 10).
Comparison of the two years (Table 2, Fig. 8c) shows that melt rates were higher in both the wet and dry seasons of year 1.The largest differences from season to season and year to year were calculated for the net shortwave budget, again indicating that annual differences in ablation were a result of differences in surface albedo and snow cover (Fig. 9f).This is confirmed by the contrasting seasonal sums of precipitation and their distribution in solid and liquid phase.Generally, precipitation was slightly lower in both seasons of the first year and a higher fraction fell as rain.These differences explain why the snow line was frequently higher in the wet season of the first year (Fig. 8a, Table 2), and hence, surface albedo was lower and the resultant impact on net shortwave caused higher ablation rates over the course of the season.
The smaller snowpack developed during the wet season of the first year was removed more rapidly during the subsequent dry season so that during the dry season the snow line rose to higher elevations in the first year than in the second year (Fig. 8), causing higher mean values of melt energy (Table 2, Fig. 9f).This highlights that the precipitation sums during the dry season were too small to form thick snowpacks and that precipitation, and also temperature, during the wet season impact the persistence of the snowpack (and higher albedo values) and subsequently affect the ablation rates during the dry season.Therefore the temperature sensitivity of annual mass balance and ablation are expected to be highest for the wet season temperature.
The similar evolution of the cumulative differences in snow line altitude and total melt energy in Fig. 12 again confirms that during our study period (i) the magnitudes of the net turbulent fluxes and net long wave radiation, which have similar or equal energy transfer properties for snow and ice, were small or compensated by additional shortwave net energy and (ii) that the seasonal and annual variability in melt energy was controlled by the shortwave net budget, thus snow cover, which was mainly a function of precipitation and temperature.So, in contrast to our findings for the ELA and AAR (Sect.4.1), annual melt rates and the specific mass balance on Shallap Glacier were sensitive to temperature, especially in the wet season when most of the annual precipitation occurs.

Uncertainties
The model uncertainty estimates presented in Sect.3.1 account for uncertainties in model parameterization schemes, parameter values and ranges (Table A1) and errors in the AWS measurements.The RMSD values calculated for the spatially distributed model run (Sect.3.1.2)were within the range of previously published studies (e.g., MacDougall and Flowers, 2011;Hock and Holmgren, 2005).Higher RMSD and lower correlation in annual surface height change in year 2 (Fig. 5) are likely to be results of more frequent snow fall and snow cover (Table 2).Observations indicate that snow cover ablates in a spatially irregular manner due to smallscale irregularities of the glacier surface, which results in higher scatter in the measured surface height change for year 2 compared to year 1 (Fig. 5) and, as the model cannot resolve small scale snow irregularities, lower model skill.
Possible errors in mass balance due to the applied vertical albedo timescale, temperature, or precipitation gradient are not captured in our evaluation as the available data extend over only a small elevation range.Thus, we only discuss potential impacts here.The dominance of melting processes for the aging of snow surfaces and related albedo decrease (e.g., Gardner and Sharp, 2010) requires accounting for the less frequent occurrence of melting and thus increasing aver-age snow albedo at higher altitudes (Brock et al., 2000).We here applied an albedo timescale gradient of 1 day (100 m) −1 (starting at the model input reference altitude of 4800 m), which yields values between 1 d (fixed lower limit) in the lowest parts of the glacier and 10.3 days at the highest elevations (5700 m) when applied to the timescale of 1.3 days (Table A1).As the accurate value of the gradient is uncertain, we calculated additional runs for gradient values of 0.5 day (100 m) −1 and 1.5 day (100 m) −1 .The difference in surface mass balance is −0.28 m w.e. for both years (lower bound) and +0.19 and +0.11 m w.e. for year 1 and year 2, respectively (upper bound).For all of the albedo timescale gradients evaluation, the timescale values over the elevation range of the glacier are within the range of values used in former studies (e.g., Sicart et al., 2011;Mölg et al., 2008;Oerlemans and Knap, 1998).Niedertscheider (1990) calculated a "general" precipitation gradient of 3.5 mm (100 m) −1 for the Cordillera Blanca based on 21 stations between 1380 and 4550 m a.s.l.which, in fact may also contain a horizontal distance effect.Our study site is well above the elevation range covered by these measurements.On Zongo glacier, Bolivia, no significant gradients in precipitation were found between 4750 and 5700 m (Sicart et al., 2007;Ribstein et al., 1995), which reflects the likely profile of no, or slightly negative, vertical precipitation gradients on high tropical mountains (Weischet, 1965).Therefore, we accept values of our precipitation gauge located at 4950 m (Fig. 1) as a best guess for the whole glacier.If precipitation does in fact increase with elevation at this site, it would yield more positive mass balance values, and we calculated that a precipitation increase of +5 % (100 m) −1 would increase mass balance by 0.3 m w.e. in both years.Mass balance is also affected by the separation of solid and liquid precipitation which was here done by linearly interpolating between two threshold values.However, possible errors were captured within our evaluation method as most variability in precipitation phase occurred at the lowest parts of the glacier where the stake measurements were conducted.
As the role of turbulent energy fluxes is minor (Sect.4.2) the vertical gradient of air temperature primarily impacts the mass and energy balance through its effect on incoming long wave radiation and the phase of precipitation, which affects accumulation, but also ablation by impacting the position of the snow line.Thus, variable lapse rates (Petersen and Pellicciotti, 2011) could have an impact on some glaciers, but for our study site the area of the glacier that was most sensitive to temperature (Sect.4.1) was very close to the altitude where temperature was measured (Fig. 1).If a temperature gradient of −0.65 • (100 m) −1 (e.g., Kaser, 2001) is applied (instead of −0.55 • (100 m) −1 ) mass balance would increase by 0.14 m w.e.(meters water equivalent) in year 1 and 0.08 m w.e. in year 2.
Wind speed is kept constant in space even though it is expected to vary due to topographic impacts on different scales.In the ablation zone, where ablation and mass balance would be most sensitive to changes in available melt energy due to variations in the turbulent fluxes, our model evaluation accounts for possible errors in wind speed.We generally expect the magnitude of the wind speeds there to be very similar to those used as model input (AWS M ) as the wind speeds averaged over 214 days of simultaneous records at AWS G and AWS M only varied by 0.1 m s −1 and the correlation for hourly values is 0.89.We also keep relative humidity constant in space as no data were available to calculate gradients.However, the deviation between the mean values measured at AWS G and AWS M are again low (1 %) and the correlation coefficient for the hourly values was high (0.95), which means that the values should be well suitable for the most sensitive lower parts of the glacier.
For future studies, additional validation data, e.g., proglacier runoff data or time series of snow line altitude, would be desirable to provide an independent model evaluation beyond the abilities of a cross validation approach, and particularly to test the robustness of the applied spatial gradients.

Conclusions
The albedo parameterization used in the surface energy and mass balance model has been improved from an earlier study (Gurgiser et al., 2013), and now better accounts for the snowfall characteristics of tropical glacier sites with frequent small snowfalls coinciding with continuous ablation.
Cross validation reveals that the RMSD of transferring parameterizations in space is < 10 %.This model evaluation benefits from the fact that the validation period included interannually contrasting mass balance conditions and stake data that come from the elevation range of the glacier, which is characterized by the steepest vertical mass balance gradients.This means that despite the small areal coverage of validation data, the stake locations can be expected to capture characteristics of the wider glacier area.
Results of the cross validation and error assessment of the model indicate that this modeling system is robust and performs well in comparison to previously published model results (e.g., MacDougall and Flowers, 2011;Hock and Holmgren, 2005).Furthermore, as there is no evidence of positive or negative trends between measurements and model solutions in space or time (Sect.3.1.2),both parameter transfer and model application outside the measurement area should yield stable internal model uncertainty for short timescales as long as the atmospheric conditions are comparable with the ones captured in this study period.For investigations on longer time periods, the characteristics of model input from AWS measurements or regional climate model output should be compared to the characteristics of the meteorological data used as model input for this evaluation (Supplement Fig. S1).If they are markedly different, the results for the parameter transfer in space and time presented here may no longer be representative.
Glacier-wide specific mass balance (Table 1 and Figs.6-8) was negative (−0.32 m w.e.) in the first year and positive (+0.51 m w.e.) in the second year, with the interannual difference in mass balance being almost exclusively produced in the portion of the glacier below 5000 m (Table 1, Fig. 6c).The calculated ELA, AAR and shape of the vertical mass balance profiles (Fig. 7) agree with the ranges obtained from conceptual approaches or measurements (Sicart et al., 2007;Kaser, 2001).The geometry of Shallap Glacier, with steep slopes above the present day mean ELA (approximately 5000 m) means that the sensitivity of the AAR to increasing temperature is low.
Seasonal means of all surface energy fluxes were strongest during the dry seasons .Averaged over the glacier area < 5000 m, net shortwave was the dominant energy source throughout the year, and net all-wave radiation remained positive as net long wave energy losses did not offset net shortwave gains in either season.Turbulent heat fluxes partially compensated each other, and the negative ground heat flux from the surface to the subsurface was mostly balanced by subsurface melt, so these fluxes did not offset the positive radiation flux and consequently ablation rates were high in all seasons.This is also the case at Antisana Glacier, but is in contrast to Zongo Glacier where net all-wave radiation can be strongly negative in the dry season and melting is restricted to the humid season (Favier et al., 2004b;Sicart et al., 2011).In addition, the model results demonstrate that the importance of sublimation during these two years is smaller than expected (e.g., Kaser, 2001).This study cannot verify whether this reflects longer-term changes in atmospheric moisture characteristics over the Cordillera Blanca.
The more negative mass balance rates for the first year are, on both annual and seasonal scales, related to lower total precipitation and higher fraction of liquid phase due to higher mean temperatures, especially during the wet season (Table 2).Thus, variations in snow cover and the altitude of snow line explain most of the difference in available melt energy in year 1 compared to year 2 (Fig. 12).As melting is the dominant ablation process throughout the year, these variations also account for the difference in annual ablation.Whereas the mass balance in the altitudes above 5000 m, where essentially all precipitation is solid, was primarily sensitive to annual precipitation sums, the mass balance of the lower glacier area was sensitive to both precipitation and temperature through its control on snow line altitude.The predominant role of the snow line altitude and sensitivity of mass balance to air temperature is a feature more similar to the mass balance regime of the inner tropical glaciers of Ecuador than the outer tropical glaciers of Bolivia (Favier et al., 2004a, b).Identifying this temperature and precipitation phase sensitivity is useful for estimating likely impacts of a changing climate on the glaciers in the Cordillera Blanca.

Fig. 1 .
Fig. 1.Overview of Shallap Glacier located on the west side of the Peruvian Andes.Study sites on Antisana and Zongo glaciers are indicated in the inset map.

Fig. 2 .
Fig. 2. Comparison of daily mean measured surface albedo with modeled surface albedo computed with optimized albedo model parameters.RMSD is the root mean square deviation, r 2 the coefficient of determination and n is the number of considered days.To exclude the effects of surface slope on albedo measurements, only hourly values with solar zenith angle above 60 • were considered for calculating the daily means.

Fig. 3 .
Fig.3.Sensitivity of modeled mass balance to the selected free model parameters set to their upper/lower limit compared to the results obtained with median values for all parameters at each of the 20 model points framing the stake locations.Points show the median values, the edge of the horizontal bars show the 25 % and 75 % percentiles, the whiskers extend to the extremes not considered as outliers and the outliers are plotted as small circles.

Fig. 4 .Fig. 5 .
Fig. 4. Results of the leave-one-out-cross-validation: measured (blue dots) and modeled (black solid lines) cumulative surface height change for the study period at each stake.Date as dd.mm is shown on the upper horizontal axis.Individual model parameter combinations for calculating surface height change for each point were derived using the best parameter combination evaluated as run with the smallest mean RMSD considering all other points.For details see Sect.2.2.3.

Fig. 6 .
Fig. 6.Annual mass balance for (a) September 2006 to August 2007 and (b) September 2007 to August 2008, showing the ELA as a black line.Given uncertainty values only reflect uncertainties in internal model parameter values (Sect.2.2.3).(c) Difference in annual mass balance between year 2 and year 1.Differences below (in sense of altitude, Fig. 1) the white line are significant above the 95 % confidence level.

Fig. 7 .
Fig. 7. Mean annual mass balance vs. altitude.Mass balance values are calculated as annual means for 10 m elevation increments.The grey areas indicate the estimated error ranges for each year (Sect.2.2.3).

Fig. 8 .
Fig. 8. Time series (September 2006-August 2008) of surface mass balance (a) with the dry and wet season duration highlighted by the grey and blue bars and mean snow line (snow depth ≥ 1 cm) indicated as black line.Time series of snow accumulation and total ablation in (b) and (c).White areas in (b) correspond to no accumulation.All values are calculated as daily means for 10 m steps in altitude.Note the different scales of the color bars.

Fig. 9 .
Fig. 9. Time series (September 2006-August 2008) of modeled energy balance terms described in Sect.2.2.1 (a-e) and of modeled surface albedo in (f).All values are calculated as daily means for 10 m steps in altitude.The grey and blue bars in (a) indicate the dry and wet season.Note the different value range in the color bar in (a) and (f).

Fig. 10 .
Fig. 10. 3 day running means of energy balance terms as spatial averages over the glacier area < 5000 m.The melt energy in the subsurface is shown as dotted line.Melt energies are multiplied by −1 to aid clarity.The grey and blue bars indicate the dry and wet season.

Fig. 11 .
Fig. 11.Time series of energy balance terms as mean values over the whole study period (September 2006-August 2008) vs. altitude.QM S is the melt energy in the subsurface.