2.1. Measurements

2.1.2. Automatic weather stations

Between 2002 and 2012 an automatic weather station (Fig. 1; AWS Moraine, hereafter AWS_M) was operated on the south moraine of Glaciar Shallap, 150 m south of the stake measurement zone at 4950 m a.s.l. All measured variables and sensors are listed in Table 1. Data are available from January 2006 until November 2009 without gaps. Hourly data from this station were used as input for our MB model as the measurements span the period of the available stake measurements that define the study period.

Since July 2010 a second AWS (Fig. 1; AWS Glacier, hereafter AWS_G) has been operating on the glacier surface at 4796 m within the stake reading zone (Section 2.1.2). Data from this station were used to (1) determine statistical transfer functions between the moraine station and the glacier station, (2) develop an incoming longwave radiation parameterization and (3) assign some model parameters pertaining to the glacier surface (Table 6 in the Appendix). A data gap in all records occurred between February and April 2012 when water entered the logger box and measurements ceased. Although the surface height was measured at this station with a sonic ranger (SR50), these data are not used here because, despite frequent maintenance, strong ablation over the whole measurement period resulted in both turning and tilting of the sensor.

Concurrent temperature, humidity and wind-speed data from both stations from July 2011 to February 2012 (214 days) were used to develop a local transfer function for air temperature from AWS_M (4950 m) to AWS_G (4796 m). Best results were obtained by the multiple regression

where is the calculated hourly air temperature at AWS_G, T _M, RH_M and ws_M are the air temperature, relative humidity and wind speed measured at AWS_M and the coefficients a ₀₋₃ are −1.825, 0.925, 0.022 and 0.191 respectively. Calibration of the regression coefficients was done using all available values. However, as we use this transfer function beyond our calibration period, we evaluated our results by applying a leave-one-out cross-validation (Reference Hofer, Mölg, Marzeion and KaserHofer and others, 2010; Reference WilksWilks, 2011), taking into account autocorrelation of the time series, to obtain an upper limit for the RMSD between measured and calculated values.

This was done by first calculating the hourly anomalies ΔT, ΔRH and Δws with respect to the mean daily cycle of the entire period. Then, for each hourly value, we calculated a multiple regression based on all other values, except those within the autocorrelation timescale (~1 week, |r| < 0.2). This leads to 4538 regression value estimates, and a reconstructed time series of the same length, where each value was derived independently of the measured value. The performance of this time series therefore provides a reference for the error magnitude of the transfer function that can also be expected outside the calibration period. The resulting RMSD between measured and calculated temperature anomalies is 0.76°C whereas the standard deviation of the measured anomalies is 1.2°C. These results indicate that the transfer function is robust, and for our study site has greater skill than constant linear temperature gradients in the range −0.4 to −0.9°C (100 m)⁻¹ which had RMSD > 1 °C. Vapor pressure (e) on the glacier was computed from RH_M using as an input. The RMSD between e calculated with T _M and is <0.2 hPa and the difference between the mean values is <0.1 hPa.

2.2. Energy- and mass-balance model

To calculate the relative surface height change (driven by climatic mass balance) for the two points on the glacier, we applied an energy-balance-based MB model. The model is described in detail elsewhere (Reference Mölg, Cullen, Hardy, Kaser and KlokMölg and others 2008, Reference Mölg, Cullen and Kaser2009a, Reference Mölg, Maussion, Yang and Scherer2012), so here we present only a concise overview. The model was run in point mode at hourly time-steps for the period 1 June 2006 to 27 August 2008, with a 3 month spin-up period in order to develop appropriate initial snow conditions. Values of temperature transferred to the glacier (Section 2.1.1), and hourly measurements of relative humidity, wind speed, solar radiation and all-phase precipitation at AWS_M serve as input. Values for the 24 parameters used in the model are given in Table 6 in the Appendix.

The model calculates accumulation as the sum of solid precipitation, surface water deposition and refreezing of liquid water in the snowpack. Solid precipitation is extracted from all phase precipitation by assuming a linear interpolation of the percentage of liquid and solid precipitation between an upper temperature threshold, above which all precipitation is liquid, and a lower temperature threshold, below which all precipitation is solid (Table 6). Total modeled ablation includes surface melt, sublimation and subsurface melt. Surface melt and sublimation are based on the surface energy balance (all terms in W m⁻²):

where S ↓ is the incoming shortwave radiation, α is the broadband surface albedo, L ↓ and L ↑ are the incoming and outgoing longwave radiation, Q _S and Q _L are the sensible and latent heat fluxes, Q _C is the conductive heat flux,Q _PS is the part of shortwave radiation that penetrates into the subsurface,Q _RPC is the heat flux from precipitation and F is the residual energy flux. If F is positive and surface temperature is 273.15 K, the melt energy Q _M equals F.

Absorbed solar radiation is the main energy source for ablation on tropical glaciers (e.g. Reference Mölg and HardyMölg and Hardy, 2004; Reference Sicart, Wagnon and RibsteinSicart and others, 2005). The original model computes S ↓ from theoretical considerations in conjunction with a cloud cover factor (Reference Mölg, Cullen, Hardy, Winkler and KaserMölg and others, 2009b), but to minimize the use of input parameterizations we chose measured S ↓ as model input. To extrapolate the horizontal radiation measurements at AWS_M to the stake locations, we first separated S ↓ into its direct and diffuse components following Reference Hock and HolmgrenHock and Holmgren (2005) and then calculated S ↓ (x, y) for the glacier gridpoints containing the two evaluation stakes accounting for topographic shading, slope and aspect. Reflected shortwave radiation from the surroundings is not considered because the sky-view factors for the stake measurement zone at Glaciar Shallap are high, and the southern slopes are usually exposed dark rock with infrequent snow cover.

The albedo module within the MB model is based on the scheme of Reference Oerlemans and KnapOerlemans and Knap (1998) which calculates surface albedo as a function of fresh snow albedo, firn albedo, ice albedo, a timescale and a depth scale. In our model, the daily threshold used to determine a snowfall event in the original model (Reference Oerlemans and KnapOerlemans and Knap, 1998) is replaced with a threshold applied to the sum of consecutive hours with snowfall (parameter P24 in Table 6). For example, if there are five consecutive hours with a snowfall rate of 0.3 cm h⁻¹, the sum is 1.5 cm and the snow albedo will be set to the fresh snow value as soon as the threshold value is exceeded.

L ↓ was calculated following the method of Reference Sicart, Hock, Ribstein and ChazarinSicart and others (2010) for Glaciar Zongo, Bolivia, using data from AWS_G to optimize the model for local conditions. The equation for all-sky longwave radiation L ₀ ↓ then reads

where c is a coefficient for clear-sky emissivity (1.24; see Fig. 2), T is the air temperature, τ _Atm is atmospheric solar transmissivity and σ is the Stefan–Boltzmann constant. L ↓ is then calculated by considering the sky view factor and the emission of the surrounding terrain (Reference Sicart, Hock, Ribstein, Litt and RamirezSicart and others, 2011, equation 11). A comparison with measured L ↓ based on 7989 hourly values shows that the fit yields a RMSD of 25 W m⁻² and a r ² of 0.59. For daily means RMSD is 16 W m⁻² and r ² is 0.81, comparable to the values obtained by Reference Sicart, Hock, Ribstein and ChazarinSicart and others (2010). L ↑ is calculated from surface temperature following the Stefan–Boltzmann law and assuming emissivity to be unity. Surface temperature change for every time-step is calculated from F which warms or cools a defined layer thickness due to energy storage change (Reference Mölg, Cullen and KaserMölg and others, 2009a).

The computation of the turbulent fluxes Q _S and Q _L is based on the bulk method, with the option of using a choice of two different stability correction functions (Table 6 in Appendix). Characteristic surface roughness lengths for snow and ice were adopted from published literature (Table 6). The subsurface module solves the thermodynamic energy equation numerically on a vertical grid with 14 layers at depths of 0.09, 0.18, 0.3, 0.4, 0.5, 0.6, 0.8, 1, 1.4, 1.8, 2.2, 2.5 and 3 m. At the bottom boundary we prescribe a constant temperature of 272 K because for the local climate we assume ice temperatures to be close to the pressure-melting point in the ablation zone throughout the year. The module considers Q _PS as a fraction (values in Table 6) of net shortwave radiation that penetrates the subsurface and is attenuated exponentially with depth (Reference Bintanja and Van den BroekeBintanja and Van den Broeke, 1995). Q _C is determined by the temperature difference between the surface and the first layer in the subsurface. As observed temperatures during precipitation are always close to 0°C, Q _RPC is not considered here.

As noted by Reference Mölg, Maussion, Yang and SchererMölg and others (2012), the model accounts for densification of the snowpack through compaction and refreezing of liquid water. This leads to a physically based modeled surface height change that we use here as the metric to be evaluated against the available measurements of surface height change at the stakes. The 24 model parameters (Table 6) can be varied within their physically meaningful ranges (e.g. Reference Mölg, Maussion, Yang and SchererMölg and others, 2012). Surface albedo constants (Reference Oerlemans and KnapOerlemans and Knap, 1998) were estimated from AWS_G measurements and lie within the range of previously reported values (e.g. Reference Klok and OerlemansKlok and Oerlemans, 2002; Reference Mölg, Cullen, Hardy, Kaser and KlokMölg and others, 2008; Reference Sicart, Hock, Ribstein, Litt and RamirezSicart and others, 2011). All other parameter ranges were chosen on the basis of previously published values and these, and their sources, are specified in Table 6.

2.3. Model evaluation method

The relative mass-balance sensitivity to each of the 24 parameters was found by computing (1) a standard model run for each stake position over the two consecutive hydrological years using median parameter values and (2) two further model runs for the upper and lower limit value of each parameter (Table 6 in Appendix) while keeping the other parameters constant at the median range values. We then computed the difference in cumulative calculated mass-balance values over the entire period between the two perturbed runs and the standard run. As all calculated model sensitivities are influenced by the chosen median value and range of the parameters, this simple test yields only estimates of the relative mass-balance sensitivity of our model to each parameter.

The parameters that resulted in a cumulative mass-balance change of more than ±5% in this single parameter sensitivity study were P3, P4, P5, P12, P15, P19, P20, P21, P22, P23 and P24 (Table 6; Fig. 5, further below). Within 1000 simulations for each of the two evaluation stake locations in each hydrological year, these 11 parameters were assigned randomly from a uniform distribution of their respective prescribed ranges (Table 6). Following Reference Mölg, Maussion, Yang and SchererMölg and others (2012), results of a model run, and its associated combination of parameters, are accepted if the RMSD between modeled and measured relative surface height at each stake reading is <10% of the measured cumulative amplitude and the deviation is <10% at the end of the time series. The best parameter combination (BPC) of the accepted runs is the one which results in the model output with the smallest RMSD between measured and modeled surface height. The model transferability in space and time is then tested by comparing RMSD and annual deviations between the results of different model runs and the measurements for each period.

Fig. 5. Difference in relative surface height change of model output for lower/upper limit of each model parameter value within the defined range (Table 6 in Appendix) subtracted from the standard run with median values for all parameters

3.1. Characteristics of measured surface height change

Annual surface point mass balance in ice equivalent can be determined directly from all stake measurements, as the surface at the stakes at the beginning and end of each hydrological year was snow-free. Cumulative surface mass balance at the points was more negative in the first year, with annual values of −9.15 and −4.75 m ice eq. at SH11 and SH22 respectively, than in the second year, for which cumulative surface mass balance was −5.65 and −2.7 m ice eq. at SH11 and SH22 respectively (Fig. 3). Less negative mass balance in the second year is accompanied by a 127 mm higher annual precipitation sum and a 0.6°C lower annual mean temperature at AWSM (Table 2). Mean vertical surface mass-balance gradients within the stake area were very high in both years (6.3 and 5.2 m ice eq. (100 m)⁻¹, respectively) although the terrain and surface features do not change markedly across this area. Similarly high gradients in surface mass balance within the ablation zone of other tropical glaciers in South America were attributed to strong vertical gradients in snow-cover frequency (Reference KaserKaser, 2001; Reference Sicart, Ribstein, Francou, Pouyaud and CondomSicart and others, 2007; Fig. 4). It is likely that the large gradient in surface mass balance within the stake measurement zone on Glaciar Shallap is also related to frequently inhomogeneous snow cover within the stake area during both years of our study period.

Mean daily relative surface height change between each measurement at the two reference points (Fig. 4) shows a strong variability for the first year, with two lowering maxima in measurement intervals 1 and 5 and two periods with reduced lowering in the wet season (intervals 3 and 4; intervals 7–9). The periods of reduced surface lowering coincide with the formation of several snowpacks that persisted for several days to weeks and which were also documented in the data of the UGRH (Section 2.1.2). The second year is different because surface lowering was strongly reduced, most likely due to frequent snow cover, between December and April (intervals 16–20) at both stakes. Strong surface lowering therefore only occurred in the dry season at the beginning and at the end of the second year (intervals 13 and 21). To support the argument about the general influences of snow-cover patterns, modeled snow depth is also shown in Figure 4. Although modeled snow depth is not always consistent with the measured surface height change, it provides an indication of the timing of snow-cover events that influenced the surface height change over the study period.

3.2. Energy-balance components and parameter sensitivity

The modeled energy-balance components (Table 3) show that at Glaciar Shallap net shortwave radiation is the greatest source of energy to the glacier surface, as has previously been well documented for other tropical glaciers (e.g. Reference Wagnon, Ribstein, Kaser and BertonWagnon and others, 1999; Reference Favier, Wagnon, Chazarin, Maisincho and CoudrainFavier and others, 2004a; Reference Francou, Vuille, Favier and CáceresFrancou and others, 2004; Reference Mölg and HardyMölg and Hardy, 2004; Reference Nicholson, Prinz, Mölg and KaserNicholson and others, 2012). Sensible heat flux is the next largest source of energy to the glacier but on average reaches a maximum of only 10% of the energy contribution of shortwave radiation. Only 30–40% of the net shortwave energy receipt is offset by the negative net longwave flux, and although a portion of the net shortwave penetrates the glacier surface and warms the subsurface, surface melting consumes the largest proportion of the net shortwave energy contribution at both stake locations. The main differences in the modeled mean energy-balance components at the two stake locations are in the reflected and penetrating short-wave components, which indicates that the snow-cover conditions at the two stake locations were different over the course of the study period, with SH22 having on average a higher surface reflectance than SH11.

As a consequence of the dominance of net shortwave in the surface energy balance, results of the sensitivity tests for SH11 and SH22 (Fig. 5) show that model output is most sensitive to parameters that are related to the shortwave energy budget. These are especially ice albedo (P19), fresh snow albedo (P20), old snow albedo (P21), albedo timescale (P22), albedo depth scale (P23) and the snow event threshold for albedo (P24). In addition, modeled mass balance is sensitive to parameters that affect the amount of solid precipitation (P3 and P4), which also impact net shortwave radiation, as well as parameters that affect the development of the snowpack and the penetration of shortwave radiation into the snow (P12 and P15). The importance of snowpack- and albedo-related parameters to modeled mass balance is in agreement with former studies of glacier energy balance both within and outside the tropics (e.g. Reference Favier, Wagnon and RibsteinFavier and others, 2004b; Reference Mölg and HardyMölg and Hardy, 2004; Reference Mölg, Cullen and KaserMölg and others, 2009a, Reference Mölg, Maussion, Yang and Scherer2012; Reference MacDougall and FlowersMacDougall and Flowers, 2011; Reference Sicart, Hock, Ribstein, Litt and RamirezSicart and others, 2011). Modeled mass balance also shows moderate sensitivity to the representative layer thickness (P5; Reference Klok and OerlemansKlok and Oerlemans, 2002) due to its impact on simulated surface temperature (Reference Mölg, Cullen, Hardy, Winkler and KaserMölg and others, 2009b).

By contrast, mass-balance sensitivity to parameters relevant for turbulent fluxes (P6–P11 in Table 6 in the Appendix) is generally low. Comparison of the data collected at AWS_M and AWS_G indicates that this is unlikely to be because the conditions at the glacier surface are not well represented by the input data (Section 2.1.2). It is possible that the sensitivity to the turbulent fluxes is reduced because the sensible and latent heat flux have different signs during the observation period (Table 3); however, even the direct comparison of the mean turbulent fluxes from the standard run (median parameter values) with the means from the sensitivity runs for parameters P6–P11 shows maximum absolute differences of <1 W m⁻². For 98% of the hourly values the absolute differences are <10 W m⁻². In this study, temperature and precipitation gradients have little impact on modeled results, but this is because our stake sites are separated by <60 m in altitude.

Model sensitivity to the selected parameters differs slightly between the two model points (Fig. 5). At SH11, snow cover is less frequent and usually thinner. Therefore model sensitivity to ice albedo (P19) is higher. In the case of the shallow snow cover, both parameters P12 (density of falling snow) and P23 (albedo depth scale) have an enhanced effect on the modeled surface height change. At SH22 the sensitivity to fresh snow albedo (P20) dominates because the surface is snow-covered for more of the study period at this location (Fig. 4).

3.3. Model performance at point scale

The results of our simulations described in Section 2.3 are plotted in Figure 6. For SH11, 306 runs passed the criteria for acceptable performance, while at SH22 only ten parameter combinations produced acceptable model output. The main reason for this difference in the number of accepted model runs is the smaller amplitude in cumulative surface mass balance measured at SH22. If we use the amplitude of SH11 to calculate the relative RMSD, 234 runs fulfill the criterion at SH22. The best model run for SH11 has a RMSD of 1.5% (0.21 m) relative to the measured surface height amplitude. The difference between the cumulative measured (−14.06 m) and modeled surface height change is −0.14 m at the end of the time series. The best model run for SH22 has a RMSD of 6.4% (0.48 m) between measured and modeled surface height, and the difference between measurement (−7.47 m) and model at the end of the series is +0.9 m. The parameter values for the best-fit parameter combinations (BPCs) resulting in model runs with the best fit to measured surface height change at both SH11 and SH22 (BPCSH11 and BPCSH22 respectively) are listed in Table 4.

Fig. 6. Comparison of cumulative relative surface height change from stake measurements and 1000 model runs for SH11 (a) and SH22 (b). Parameters to which modeled mass balance is most sensitive (Section 3.2) are assigned randomly from uniform distributions within defined ranges (Table 6 in Appendix). Accepted runs are those that meet two criteria: (1) RMSD between modeled and measured relative surface height at each stake <10% of the measured cumulative amplitude and (2) surface height difference <10% at the end of the model run. The best-fit model run is the one with the minimum RMSD between measured and modeled surface height, and the RMSD shown is for the best-fit run.

3.4. Parameter transferability in space

In Section 3.3 the model performance was evaluated and BPCs were determined for each of the two stake locations separately. Here we test the parameter transferability in space by applying the BPCs for SH11 to SH22 and vice versa. As the BPCs for the two points are different, the RMSD between model and measurements increases when the parameter values of the other stake location are used. The increase in RMSD compared to the results for the best runs is 4% (0.26 m) for SH11 and 10.3% (0.77 m) for SH22 over the entire period.

At location SH11 (Fig. 7a), when the best parameter combination for SH22 (BPCSH22) is used, the modeled surface evolution only deviates markedly from that modeled with BPCSH11 in the second year from interval 14, after which surface lowering is reduced until the middle of interval 18. This is due to the smaller albedo depth scale (P23) which keeps surface albedo higher for longer (e.g. Reference Mölg, Cullen, Hardy, Kaser and KlokMölg and others, 2008) and therefore reduces the net shortwave energy receipts (Fig. 7c). Subsequently, surface lowering is reduced in both runs until the end of interval 19 and increases again toward the end of the study period. Figure 7c shows that net shortwave radiation controls the differences in melt energy between the two runs. The additional errors between measurement and model as a result of applying the transferred parameter combination (BPCSH22) instead of the best one (BPCSH11) are 5 mm ice eq. for year 1 and 735 mm ice eq. for year 2.

Fig. 7. (a. b) Modeled cumulative relative surface height change from 29 August 2006 to 27 August 2008 with best parameter combinations (BPCs) for each location at SH11 (a) and SH22 (b). Measured values are shown as black dots. Bars at the upper borders of the plots show calculated surface albedo for both runs at hourly time-steps (upper bar is for SH11 run, lower bar for SH22 run). (c, d) Cumulative differences in daily means of energy-balance terms (SW_net and LW_net are the net shortwave and net longwave radiation fluxes; the other abbreviations are explained in Section 2.2) for SH11 (c) and SH22 (d). In (c) the values calculated with BPCSH11 are subtracted from BPCSH22, and vice versa in (d).ΔSW_in and ΔLW_in are zero for all runs. For easier readability, ΔQ _M is multiplied by −1 (since Q _M in the energy balance is negative by sign convention).

At location SH22 the model solutions diverge early in the first year during interval 3 (Fig. 7b). Albedo time series and modeled snow depth show that BPCSH11 cannot reproduce the snow cover that was documented by measurements and caused the strongly reduced surface lowering measured during intervals 4 and 5. As surface albedo remains low for SH11 settings during this period, the differences in melt energy and ablation between the two runs are large (Fig. 7d). However, from interval 5 on, both runs show a similar pattern of surface height evolution. For the second year, mass balance is overestimated by model runs with both BPCs but the pattern of surface height change produced by the runs is similar, so the differences between the runs are smaller than in the first year. As with SH11 above, the differences in net shortwave radiation and in the penetrating flux of shortwave radiation dominate the difference in melt energy (Fig. 7d). The additional errors due to the parameter transfer (BPCSH11 instead of BPCSH22) are 1326 mm ice eq. for year 1 and 542 mm ice eq. for year 2.

Overall, Figure 7c and d illustrate that the parameter transfer in space primarily alters the net shortwave energy fluxes, and thus different snow-cover patterns evolve (S ↓ is the same for all runs). As the sum of modeled solid precipitation differed by only 2%, the differences in snow-cover patterns are mainly driven by different rates of snow ablation and removal in the model, which are in turn a result of the different combinations of albedo constants. On short timescales (e.g. a few months) these differently modeled patterns can cause very different ablation rates, while on longer (annual) time series the differences can be compensated (e.g. Fig. 7b and d) as evidenced by a tendency for convergence in cumulative melt energy differences and modeled surface height evolution with both BPCs toward the end of the records. For periods with only sporadic snow cover (year 1 at SH11) the results for both parameter combinations are very similar.

3.5. Parameter transferability in time

In Section 3.3, we presented model performance with constant parameter combinations for both hydrological years. Here we target parameter transfer in time. Figure 8 and Table 5 show the best runs for each point and each year, and compare them to runs with parameter settings for the same stake but optimized for the other year.

Fig. 8. (a, b) Modeled cumulative relative surface height change with BPC for each year at SH11 (a) and SH22 (b). Measured values are shown as black dots. Bars at the upper borders of the plots show calculated surface albedo for both runs at hourly time-steps (upper bar is for SH11 run, lower bar for SH22 run). (c, d) Cumulative differences in daily means of energy-balance terms (SW_net and LW_net are the net shortwave and net longwave radiation fluxes; the other abbreviations are explained in Section 2.2) for SH11 (c) and SH22 (d). In (c) the values calculated with BPCSH11Year2 are subtracted from BPCSH11Year1; in (d) the values calculated with BPCSH22Year2 are subtracted from BPCSH22Year1.ΔSW_in and ΔLW_in are zero for all runs. For easier readability, ΔQ _M is multiplied by −1 (since Q _M in the energy balance is negative by sign convention).

At SH11, Figure 8a and c illustrate that during the first year the modeled energy fluxes are very similar and the model performs well in reproducing measurements with the settings optimized for either year. However, in year 2, the best parameter combination for year 1 (BPCSH11Year1) clearly underestimates surface lowering, so the RMSD becomes large. The different evolution can be explained by the higher surface albedo between intervals 14 and 17 in the second year, which is related to persistent snow cover. In contrast, the best parameter combination for year 2 (BPCSH11Year2) produces good results for both years. The additional error between measurement and model when using BPCSH11Year2 instead of BPCSH11Year1 for year 1 is therefore only 18 mm ice eq. For year 2 the additional error (BPCSH11Year1 instead of BPCSH11Year2) is 2069 mm ice eq.

At SH22, the BPCs per year at location SH22 result in much poorer model performance if they are transferred to the other year. This corresponds to a very different evolution of albedo (Fig. 8b), and thus high deviations in net shortwave radiation and melt energy (Fig. 8d) through time in the model output using the different BPCs for each year. The additional error between modeled and measured mass balance due to the parameter transfer in time (BPCSH22Year2 instead of BPCSH22Year1) is 3179 mm ice eq. for year 1 and 1657 mm ice eq. for year 2 (BPCSH22Year1 instead of BPCSH22Year2).

As was also found to be the case for parameter transfer in space, parameter transfer in time primarily affects the shortwave energy budget (Fig. 8c and d). As all other fluxes are much smaller (Section 3.2), the differences in net shortwave radiation control the differences in melt energy.

3.6. Impact of snow cover on model performance and transferability

Throughout the model experiments in the former sections, model results correspond to the measurements most closely in the first year at SH11, where snow cover was infrequent (73 out of 357 days with snow depth >1 cm in the first year; 212 out of 372 days in the second year) in all model runs. At SH22, frequent snow cover was experienced throughout both the study years (524 of 729 days with snow depth >1 cm) and model performance was poorer than at SH11. Thus, we can deduce that the model error, and especially that associated with transferability, increases considerably during periods when there are frequent snowfall events and when there is persistent snow cover over several days or weeks in the ablation zone.

This finding is consistent with the fact that modeling snow ablation, and its concomitant surface height evolution, is harder than modeling ice ablation (e.g. Reference Mölg, Cullen, Hardy, Winkler and KaserMölg and others 2009b, Reference Mölg, Maussion, Yang and Scherer2012) due to the more complex structure of snow and associated processes (e.g. large density variations, refreezing of meltwater, retention of liquid water), which pose challenges even in detailed snowpack models (Reference EtcheversEtchevers and others, 2004). In addition, uncertainties in precipitation measurement and in distinguishing between liquid and solid precipitation (e.g. Reference Wagnon, Lafaysse, Lejeune, Maisincho, Rojas and ChazarinWagnon and others, 2009) can be a large source of uncertainty in the model accumulation, but we have no adequate data to assess this uncertainty.

The increase of model error in the presence of more frequent snow cover in the ablation zone could also be related to deficiencies in parameterized albedo over snowpacks, associated with our model not capturing the full complexity of the underlying processes. Many published albedo models using simple parameterizations or process-resolving snowpack models are still subject to sizeable errors compared to observations (Reference Brock, Willis and SharpBrock and others, 2000; Reference EtcheversEtchevers and others, 2004), and in our model it is obvious that as soon as surface albedo shows high fluctuations, model performance decreases and parameter transfer can lead to large errors (Figs 7a and b and 8a and b).

An example is plotted in Figure 9, where calculated snowfall intensity, snow depth and surface albedo are shown for 15 days in May/June 2008 for BPCSH22 applied at SH22. Until 29 May, only minor precipitation events were measured and the snowpack thickness and albedo were decreasing steadily. In the evening of 30 May a snowfall event of 1.6 cm exceeded the albedo threshold value (P24 in Table 4), so model albedo is set to the fresh snow value (P20 in Table 4). A similar event occurred the following day. Then surface albedo decreased until the end of the considered period. However, even though the total amount of fresh snow has ablated in the model before 3 June (right red dot), surface albedo remains much higher than it was before the fresh snow event. That is, the albedo scheme overestimates surface albedo in case of fresh snow falling on older snow if the fresh snow melts away before its exponential ageing reaches the albedo value of the former snow surface. This effect is especially noticeable in the case of frequent light snowfall events on an existing snowpack, which is often the case at SH22 where >40% of all 336 calculated snowfall events (sum of consecutive hours with ongoing snowfall) result in snow cover only 1–2 cm thick. As frequent, light, snowfall events are typical for glaciers in the outer tropics, especially in the dry season, we surmise that a similar problem is likely to affect surface albedo calculated for Glaciar Zongo (Reference Sicart, Hock, Ribstein, Litt and RamirezSicart and others, 2011) using the parameterization of the US Army Corps of Engineers (1956). As albedo schemes are in many cases the largest source of errors (Reference KlokKlok, 2004), modifying the range of albedo input parameters or modifying the parameterization used in the MB model might lead to improved performance. For example, Reference Brock, Willis and SharpBrock and others (2000) developed a parameterization scheme for surface albedo of a mid-latitude glacier as a function of cumulative daily maximum temperatures greater than 0°C which performed better than previously published parameterizations. Another potentially helpful modification of the applied albedo scheme (Reference Oerlemans and KnapOerlemans and Knap, 1998) could be to account for the snow albedo before the event (Reference Machguth, Purves, Oerlemans, Hoelzle and PaulMachguth and others, 2008) or to keep track of the surface albedo prior to a snowfall event, and reset it to this value if the fresh snow layer vanishes before it has completed the model-defined ageing time. More complex albedo schemes are not easily applicable at Glaciar Shallap because they require extensive high-quality measurements for calibration, and further alternatives require accurate information on snow grain size (e.g. Reference Wiscombe and WarrenWiscombe and Warren, 1980; Reference Gardner and SharpGardner and Sharp, 2010), that are not available at this site.

Fig. 9. Calculated snowfall rate, snow depth and surface albedo for BPC for SH22 between 22 May and 5 June 2008. Red dots indicate model snow depth of 44.8 cm.