Degree-day melt models for paleoclimate reconstruction from tropical glaciers : calibration from mass balance and meteorological data of the Zongo glacier ( Bolivia , 16 ◦ S )

Introduction Conclusions References


Introduction
The physical mechanisms controlling the mass balance of mountain glaciers are primarily determined by the climatic conditions in the accumulation and ablation area (e.g.Ohmura et al., 1992).Additionally, because of their relative small size, alpine glaciers respond quickly (at the timescale of 10 0 to 10 1 yr, which is fast compared to the historical and geological timescales) to changes in main atmospheric variables, such as temperature or precipitation (e.g.Oerlemans, 2005).These two characteristics imply that mountain glaciers are very useful and precise proxies for paleoclimatic reconstructions (e.g.Blard et al., 2009;Flowers et al., 2007;Mark et al., 2005;Oerlemans, 2005).This is particularly true in high elevation areas where other climatic archives (such as pollen records) are scarce or inexistent.Because paleo-glaciers are by definition vanished, this climatic reconstruction must however be done through numerical models in which past glacial extents are considered as inputs and paleotemperature and paleoprecipitations as outputs.Accurate reconstruction of these main paleoclimate variables requires a proper modeling of the past glacier mass-balance (mass balance = accumulation -ablation).The computation of accumulation is generally quite straightforward (e.g.Hock, 2005;Lejeune et al., 2007;Vincent, 2002), so the main difficulty of paleoclimate studies is to model properly the ablation of ice and snow at the glacier surface.
In the recent decades, a large variety of melt models have been developed to link climatic variables with ablation ranging from (i) simple temperature-index models, also called positive-degree-days (PDD) models (e.g.Braithwaite, 1995;Braithwaite and Olesen, 1985;Hock, 1999Hock, , 2005;;Johannesson et al., 1995;Vincent, 2002) that are based on the empirical relationship existing between ablation and the sum of positive air temperatures, to (ii) sophisticated energy balance models (e.g.Anslow et al., 2008;Brun et al., 1989;Greuell and Konzelmann, 1994;Hock and Holmgren, 2005;Klok and Oerlemans, 2002;Molg et al., 2003) that compute melt from the energy flux budget at the glacier surface.
Given the complexity and the large variety of surface processes involved in the mass balance of tropical glaciers, the most accurate models are those based on an energy balance-budget.However, these models generally require a large number of input data, such as, i.e., wind velocity, relative humidity, cloud cover, all these parameters being generally difficult to know over geological timescales, even for the Quaternary period (last 2 million years).Thus, although energy balance models are physically more relevant than PDD models, this sort of model is generally difficult to use in paleoclimatic studies.This issue could however be overcome through careful sensitivity tests, but such approach is time-consuming and can yield significant uncertainties associated with the modeling results (e.g.Plummer and Phillips, 2003).Moreover, although PDD models assume significant simplifications of several physical processes, it has been demonstrated that they are able to yield robust results at the scale of a watershed, at monthly and annual timescales (Hock, 2005;WMO, 1986).PDD models also have the advantage to require only a limited number of input data that are generally widely available, such as monthly air temperature.The classical version of PDD models (Braithwaite and Olesen, 1985) assumes a proportional relationship between ablation, A (in mm water equivalent (mm w.e.) per unit of time) and the averaged positive temperature, T p (in is the degree-day melting factor (in mm w.e.• C unit −1 of time).However, this equation does not include the influence of solar radiation, and this simple model has thus a poor ability in reproducing day-night fluctuations, the influence of insolation variations over geological time scales, nor the influence of topography.This drawback can possibly hamper the model ability in reproducing hourly change in melting water flux, but it might also affect the accuracy of paleoclimatic reconstructions.This observation led several authors (e.g.Hock, 1999;Pellicciotti et al., 2005) to improve the classical PDD model by including the effect of potential direct solar-radiation in the key relationship linking air temperature and ablation.These improved PDD models were however calibrated from mass balance and meteorological data from mid-latitude glaciers (from Sweden and Switzerland, respectively) located at low elevation.The validity of these so-calibrated parameters is uncertain for tropical glaciers.Indeed, several studies (e.g.Sicart et al., 2005) have shown that, in the Tropics, ice melting is more sensitive to radiative budget, which implies that albedo (through snow precipitation) play a significant role on the ablation amount.Nonetheless, despite this particularity, ablation seems to be relatively well correlated with inter-annual temperature variations (e.g.Sicart et al., 2008;Francou et al., 2003), which suggests that PDD models could be used to reproduce annual mass balance of tropical glaciers.
The main goal of the present study is to take advantage of the high resolution mass balance and meteorological dataset obtained on the Zongo glacier (Bolivia, 16  (Braithwaite, 1995), (ii) the PDD model that includes short wave radiation to calculate ablation (Hock, 1999) and (iii) a PDD model in which short wave radiations are totally decoupled from temperatures (e.g.Pellicciotti et al., 2005).Calibration of the model parameters is performed by searching the best fit between the modeled and the measured mass balance for the 9 hydrological years (from 1 September to 31 August) between 1997 and 2006.Some of these models use different melt factors for ice and snow.This is an empirical way to include the effect of albedo, which is one of the main parameters controlling of the shortwave budget at a glacier surface (Sicart et al., 2005).
The motivation of this study is thus two-folds: first, to assess the temperature-index models capacity in reproducing the yearly mass-balance of a well-studied tropical glacier, and, second, to establish a robust calibration of these empirical models that can then be used for paleoclimatic quantitative reconstructions.Indeed, the parsimony of PDD models makes this approach very useful to interpret the past extents of ancient glaciers.This is of particular importance in certain high elevation tropical areas where ancient glaciers are sometimes the sole plaeoclimatic proxies.

The Zongo glacier
The Zongo glacier (16 • 15 S, 68 • 10 W) is located in the Huyana Potosi massif, Cordillera Real (Bolivia), in the East border separating the wet Amazonian Basin and the dry Altiplano plateau, about 30 km North of La Paz city (Fig. 1).The meteorological conditions in the Huayna Potosi massif are typical of outer tropical conditions (Kaser et al., 1996) and 5200 m) and to the East in the lower part (between 5200 and 4900 m) (Sicart et al., 2005(Sicart et al., , 2007;;Wagnon et al., 1999).In 2006, the Zongo glacier itself covered about 1.96 km 2 , for a total length of 3 km and a width of 0.75 km (Soruco et al., 2009) (Fig. 1).The Zongo has the characteristics of a temperate glacier (Francou et al., 1995).The hydrological year has been defined as starting on 1 September and ending on 31 August (Ribstein et al., 1995).During the wet season (summer, from September to March) solid accumulation and ice-snow melt occurs simultaneously.Contrarily, during the dry season (winter, from April to August) solid accumulation is sparse, and ice-snow melt is reduced but sublimation is significant (Favier et al., 2004).

Meteorological data
The meteorological data of the nine hydrological years studied in this article (  1).Such variability probably results from the seasonal fluctuations of the relative humidity and cloud covers, which are important factors controlling the warming of air column by long and short-wave radiation (Kageyama et al., 2005).Although this seasonal variability is small and that the annual lapse rate is not significantly different from the previously reported values (Klein et al., 1999), the following calculations are performed taking into account these variable monthly lapse rates.
Since September 2004, precipitation has been continuously measured at the ORE site at 5050 m (cumulative sum every ten minutes) using a Geonor T-200B sensor (opening area of 200 cm 2 ), which is the reference sensor to measure both solid and liquid precipitation thanks to its weighing device.The correction factor proposed by (Forland et al., 1996), depending on air temperature and wind velocity, was applied as in the works of (Lejeune et al., 2007) and (Wagnon et al., 2009) to take into account the difficulty of the gauge to properly collect solid precipitation in a windy environment (precipitation is systematically underestimated).For both 2004-2005 and 2005-2006 hydrological years, applying this correction factor resulted in a 40 % increase in the total amount of precipitation recorded by the Geonor sensor (Lejeune, 2009).Close to the ORE site, at 5080 m, a totalizer gauge (big cylindrical waterproof tank containing a thin layer of oil to avoid evaporation and of opening area of 2000 cm 2 ) has been manually measured every month since 1996.There is a good agreement between the corrected Geonor monthly precipitation and the monthly precipitation measured using the totalizer for both hydrological years 2004-2006 (R 2 = 0.94), but with a systematic under-estimation of 20 % for the totalizer.Consequently, for monthly precipitation before September 2004, we used the totalizer measurements corrected by a factor 1.2.Moreover, precipitation data recorded from the snow pits located at different elevations on the Zongo glacier (Fig. 1) do not indicate significant change in precipitation against elevation.The precipitation data measured on the moraine at 5050 m were thus used over the whole elevation range over which the Zongo mass balance is modeled (i.e. between 4900 and 6000 m).Mean annual precipitation range from 1044 to 1529 mm yr −1 for the 1997-1998 and 2000-2001 hydrological years, respectively (Table 1).All years are characterized by significant precipitation seasonality.Indeed, more than 70 % of precipitation on the Zongo occurs between November and March (Table 1).This seasonality arises from the location of the Cordillera Real: indeed, this mountain range belongs to the outertropics, in a zone of the Andes that is very sensitive to the South American monsoon (Garreaud et al., 2009).

Mass balance data
Mass balances of Zongo glacier are available from glaciological (since 1991 to 2006) hydrological (from 1975 to 2006) and geodetic methods (in 1956, 1963, 1975, 1983, 1997 and 2006) (Soruco et al., 2009).Annual mass balances of Zongo glacier have been systematically measured since 1991 (Francou et al., 1995).However, only the period from 1997 to 2006 have been corrected using the geodetic method by photogrammetry (Soruco et al., 2009).Therefore, for the present study, we have only used the mass balances data monitored over these 9 hydrological years.Only a short description of the glaciological data acquisition is given here, because the methods have already been thoroughly described in previous articles (Francou et al., 1995;Sicart et al., 2007;Soruco et al., 2009).
The traditional glaciological method (Paterson, 1994) applied to the Zongo glacier aims to estimate the loss or gain in mass over the whole glacier surface from field measurements.Annual mass balances have been obtained using in average 14 stakes of monthly measurements over the ablation zone and 3 snow pits and drilled cores for yearly measurements over the accumulation zone.Moreover, snow height measurements are also required in the ablation area, because snow accumulation can occur at any time in the year over the entire glacier.These measurements were performed assuming that snow and ice densities are 400 and 900 kg m −3 , respectively (Soruco et al., 2009).In the accumulation area, the snow density of the snow pits and drilled cores were measured at the end of the hydrological year.Finally, the specific glacier mass balances were computed using actual surface areas derived from photogrammetric measurements (Soruco et al., 2009).
The annual masse balance data were finally calculated by summing the monthly data for the 9 hydrological years between 1997 and 2006 (Table 2).
3 Description of the models

General considerations
The present study aims to evaluate the efficiency of 3 different types of models.Models 1 and 2 belong to the same category: both are derived from the classical PDD model (Braithwaite and Olesen, 1985).Model 3 is inspired from the enhanced versions developed by (Hock, 1999) and model 4 is a derivate of the model of (Pellicciotti et al., 2005).These 2 last types of models have the particularity to include the short wave solar radiations.The main difference between our study and these previous works is that the computation is here performed at a monthly timescale.Indeed, the Zongo mass balance data are only available at this temporal resolution.Moreover, the main motivation of the study is not to provide a detailed modeling of ablation at an hourly time scale, but rather to test a robust and simple approach to derive general paleoclimatic conditions from tropical glacier paleoELAs.Finally, (Sicart et al., 2008) has demonstrated that temperature-index models are not appropriate to simulate daily melting but that they can be relevant to describe the ablation of tropical glaciers at longer timescales (monthly or yearly).Although the model are run on a monthly time-step, all ablation parameters are reported here as daily values, to allow a direct comparison with the ablation parameters obtained from other calibration studies.For all models, the required input data and the calibrated parameters are described in Table 3.
All models share a common structure: monthly snow accumulation, and snow/ice ablation are computed each 50 m over the relevant elevation range (from 4950 to 6050 m for the Zongo).Annual accumulation and ablation are calculated by summing the monthly values.Finally, the mass balance elevation law, MB (mm w.e.yr −1 ) is computed at each elevation subtracting ablation (A in mm w.e.yr −1 ) from accumulation (S in mm w.e.yr −1 ):

Modeling of snow accumulation
The annual snow accumulation, S (mm w.e.yr −1 ), is computed in a similar way in all models.S is the sum of the monthly snowfalls: that temperatures have a normal distribution around the monthly mean temperature T mi , with a standard deviation σ i (Johannesson et al., 1995): where P i (mm w.e.month −1 ) is the monthly precipitation.Lejeune et al. (2007) has shown that, in the tropical Andes, 100 % of falls is snow below −1 • C, and that 100 % is rain above 3 • C. T s is thus supposed to be 1 • C in the following.

Modeling of ablation
Annual snow ablation, A (mm w.e.yr −1 ), is computed in a similar way in all models.A is the sum of the monthly ablation: A i (4)

Calculation of the monthly positive temperature T Pi
In each ablation model described below, the monthly positive temperature T Pi ( • C) is calculated following the positive degree-day approach, which assumes that temperature follows a normal distribution characterized by the monthly mean T mi ( • C) and the standard deviation σ i ( • C).The average positive temperature is thus computed as the integral of this distribution function on the interval (0,+∞) • C (Johannesson et al., 1995): This model is the classical PDD model (e.g.Braithwaite and Olesen, 1985).In this high simplicity model, only one melting factor is involved, MF (mm w.e. • C day −1 ).Monthly ablation is computed as: In this model ablation of snow and ice is calculated with different degree-day factors: -MF snow (mm w.e. • C −1 day −1 ) for snow (if the glacier surface is snow).
At each monthly time step, total ablation, A i , is the sum of snow ablation and ice ablation: Ice melting only occurs if the ablation of snow was sufficient to remove all the snow cover.In order to correctly quantify the fraction of time during which the ablation of snow occur each month, the model first calculates a parameter α that is defined as: As, for each month, the amount of total snow input is known (S i ), "α" is in practice computed as follow: -If a >1, the snow ablation is smaller tan accumulation, thus, in this case, only partial melting of the snow cover occurs.Monthly ablation A i is then: -If a <1, the total monthly ablation A i is sufficient to remove all the snow cover and start melting the glacier ice.In this case, A i is then calculated as:

Model 3:
This model is inspired from the one developed by (Hock, 1999).As in model 2, different melting factors are used for snow and ice.The particularity of this model is to incorporate the clear-sky short wave solar radiations.Similarly to the classical PDD model, ablation only occurs under positive temperature.In this model, snow ablation is calculated following: Ice ablation is obtained from: where MF (mm w.e. • C −1 day −1 ) is the temperature melt factor (MF is similar for snow and ice in model 3), a (mm w.e.d study, I i was calculated assuming clear-sky conditions, according to the equation (Hock, 1999): where I 0 is the solar constant (1368 W m −2 ), R i the mean monthly distance Sun-Earth, R m the mean Sun-Earth distance, ψ = 0.75 the atmospheric clear-sky transmissivity, P and P 0 are respectively the local and sea-level atmospheric pressure, Z i the mean monthly zenith angle and θ the angle of incidence between the normal to the grid slope and the solar beams.θ was calculated for each altitude step (between 4950 and 6050 m) by using a 30 m digital elevation model of the Zongo surface.This type of models (i.e.models 3 to 4) thus takes into account the glacier geometry to compute the actual amount of incident solar beams reaching the glacier surface.
In practice, the model follows the same algorithm as Model 2 to compute, for each month, the respective proportion of snow and ice melting: α = A snow /(A snow + A ice ).

Model 4:
This model is directly derived from the one proposed by (Pellicciotti et al., 2005).The slight difference with model 3 is that the ablation arising from solar radiation, I, is fully independent from temperatures.The particularity of the model presented here is to allow ablation, even when temperatures are negative (i.e.T p = 0 • C).
a is the solar radiation melting factor expressed in mm w.e.d −1 W −1 m 2 .I is the short wave incoming solar radiation calculated as for model 3 (Eq.14, Sect.3.3.4).For each hydrological year, the best-suited models parameters were determined by minimizing the efficiency criterion R 2 .This criterion is calculated as follow (Nash and Sutcliffe, 1970): where MB is the mass balance, the subscripts "obs" and "sim" denote "observed" and "simulated", respectively and i denotes the elevation step.
In practice, the parameters optimization was realized for all models and each hydrological years following the method described in (Hock, 1999).The R 2 values are mapped as a function of the considered ablation parameters value (Fig. 2).The retained best parameters are those yielding maximum R 2 values.Note that, for each of the 4 models, we also calculated the best R 2 for the whole 9 hydrological years (from 1997 to 2006) taken together.This approach permitted to obtain the most relevant ablation parameters for the whole dataset.

Results: calibrated parameters
All the calibrated parameters from the 9 hydrological years are displayed in Table 4.

Model 1:
The calibrated melting parameters MF range between 9 (2003-2004) and 15.5 mm w.e.• C −1 d −1 (1998-1999).The efficiency criterion R 2 ranges between 0.77 and 0.99, which underlines that, despite its simplicity, this model is quite efficient in reproducing the observed mass balance (Table 4, Figs. 3 and 4).The altitudinal gradient of the simulated mass balance is indeed in very good agreement with the observed one (Fig. 3).Besides, tropical glaciers are characterized by much stronger vertical mass balance in their ablation area than other glaciers (e.g.Kaser, 2001;Wagnon et al., 1999) and the fact that this gradient is well simulated gives confidence in the ability of the model to reproduce the glacier mass balance.Similarly, the modeled and observed ELA are generally closely matching, within the observation uncertainty (<50 m).The model optimization performed from the whole dataset yielded a MF of 11.9 ± 1.3 mm w.e.
• C −1 d −1 , with a R 2 of 0.92. .The efficiency criterion R 2 ranges from 0.82 to 0.99, which indicates that this model structure has also a good ability to fit the observed mass balance.However, the R 2 obtained from the whole dataset (0.93) is not significantly better than those of models 1 and 2. Thus the introduction of the parameter "a" did not result in a significant improvement of the model efficiency (Fig. 3).The "all-years" calibration yielded a value of 8.8 ± 1.0 mm w.e.
for "a", with a R 2 of 0.93.

Inter-model comparison
All models including the solar radiations (model 3 and 4) are characterized by a significant interannual variability of the "a" factors calibrated from each hydrological years (Table 4).Two main explanations could be invoked to explain such a large discrepancy: first, these low complexity models do not take directly into account the fluctuations of cloudiness, nor those of albedo.Only clear-sky radiations (input I in W m −2 ) are considered in Model 3 and 4.Under cloudy conditions, the true incoming short wave radiations may be significantly lower, and, so, the energy available for melting ice or snow.The observed interannual variability of these empirical "a" factors might thus result from unrecognized fluctuations in cloudiness or other atmospheric processes.Although actual cloudiness or albedo data are available for the 2003-2006 hydrological years (Sicart et al., 2008), this study did not test models including the true global short-wave radiations, nor the long wave budget.Indeed, the aim of this study is to keep as small as possible the number of input data.Moreover, past fluctuations of cloudiness and albedo are by definition very difficult to constrain.Second, the R 2 value has only a low sensitivity to the "a" parameter for certain hydrological years.Such a smooth definition of the best R 2 may lead to quite ambiguous results, implying a large uncertainty in the calibrated "a" factors.This observation is consistent with the conclusion that the introduction of the short-wave radiation does not improve significantly the model efficiency.
The optimized temperature-index factors MF are not affected by such a large dispersion.Indeed, the dispersion is still significant but all MF values belong to the same order of magnitude.In certain cases, as for Model 1, the standard deviation of the yearly MF factors is not higher than ∼15 % of the average value.A comparison of the MF values of "solar radiation models" (Model 3 and 4) with those of pure "temperatureindex" models indicates that the introduction of the solar radiation parameters "a" induces a reduction of the calibrated MF parameters values.Indeed, in Model 3, the melting term a × I is in the same order of magnitude as MF (Table 4, Fig. 5).For model 4 (A = MF × T + a × I), the ratio (MF × T )/(a × I) is highly correlated with the elevation (Fig. 5b).This is due to the fact that, at lower elevation, higher temperature implies that the term MF × T is significantly higher than a × I. On the other hand, above 6000 m, the ratio (MF × T )/(a × I) is lower than 1.This characteristic is a specificity of Model 4. It allows ablation even under negative temperatures.
Important observations can be obtained by comparing the efficiency criterion of each model (Table 4).In spite of its simplicity, Model 1 is characterized by a remarkably high R 2 values ("All years R 2 " is 0.92).This indicates a significant ability of this model in reproducing the observed mass balance (Fig. 3)."All years R 2 " of model 2, 3 and 4 are equal to 0.93, respectively.Consequently, the introduction of solar radiations or the use of different melting factors did induce only a slight improvement of the model performance (Figs. 3 and 4, Table 4).
Correlation coefficients have also been computed considering separately the positive and negative mass balances (R 2 pos and R 2 neg , respectively).These efficiency criterions allow a detailed diagnostic of the ability of each model to fit the data in the accumulation or in the ablation zones (Table 4).For all models R 2 neg (0.87) are significantly better than R 2 pos (0.48 to 0.50), which shows that these models are more efficient to model the mass balance in the ablation zone than in the accumulation zone.R 2 neg are similar for each of the 4 models, suggesting that all of them have a comparable ability in fitting the data in the ablation zone.However, R 2 pos are slightly different between each model.The highest value for R 2 pos is obtained with model 4 (R 2 pos = 0.50).This is probably because this model allows ablation in the accumulation zone.Indeed, Fig. 3 indicates that model 4 offers a much better fit between model and data above the ELA, while all of the 3 other models systematically overestimate the observed mass balance in the accumulation zone.Indeed, models 1 to 3 cannot yield ablation under negative temperatures.This observation is also confirmed by the slope of the regression line determined plotting the observed against the calibrated mass balance (Fig. 4): indeed models 1, 2 and 3 are characterized by a slope that is smaller than 1, while the MB simulated from model 4 define an almost perfect 1:1 line.
By comparing the efficiency of different type of models on a glacier in Switzerland, (Pellicciotti et al., 2005) obtained very similar results: models accounting separately the influence of temperature and short-wave radiations perform better.This characteristic may be potentially useful to take into account ablation mechanisms such as sublimation or snow erosion by wind.Although it is very difficult to parameterize these mechanisms, the possibility to take into account these processes in low-complexity models such as model 4 may represent a significant improvement of classical PDD models.
It is thus tempting to suggest that model 4 should be favored for next climatic studies based on tropical glaciers.We also tested model having different "a" parameters for snow and ice, but this modification did not induce any improvement in the model performance and thus these models are not described in the present study.Moreover, considering solar-radiations in the ablation calculation (model 3 and 4) can be very useful in certain situations.First, this type of model is able to take into account specific effects, such as the influence of the glacier orientation, slope or shading in entrenched valleys.The way the solar radiation is calculated in the ablation equations of model 3 and 4 indeed permits to incorporate geometrical effects that can potentially modulate the incoming short wave radiation.These models can thus be used as efficient tools to evaluate and quantify the influence of topographic effects.Second, these models have the ability to incorporate the past fluctuations of solar radiations resulting from orbital forcing.This aspect is important for paleoclimatic reconstructions.

Models 1 and 2
As previous calibration studies from tropical glacier are almost inexistent, it is tough to compare directly the temperature-index factors obtained here from the Zongo glacier with previously published melting index.The closest-to-equator calibration studies were located in the Himalaya (Kayastha et al., 2000;Singh and Kumar, 1996;Singh et al., 2000) between 27 • N and 35 • N latitude.In a review paper, (Hock, 2003) showed that these Himalayan glaciers yielded ablation factors ranging between 5.9 and 11.6 mm ) for snow and from 6.6 to 16.9 mm ) for ice respectively.
These values are not significantly different from those obtained here from model 1 (MF = 11.9 ).Some of the MF ice factors obtained here are higher than previously determined ice ablation factors.However, this observation should be

Model 3
This model is very close to the one used in our study dedicated to the last glaciation of the Hawaiian summits (Blard et al., 2007).In this initial study, model parameters were calibrated from the 2003-2005 hydrological years record from the Zongo glacier.This study yielded optimized values of 3.7 mm for MF ice and ∼4 × 10 −3 for "a" (Blard et al., 2007).The new MF parameters are comparable, but the new "a" parameters is twice larger than the one obtained from the initial calibration study of (Blard et al., 2007) (Table 4).This discrepancy probably results from the fact that the 2 studies were performed by using different calibration years: 9 hydrological years (1997 to 2006) are used in the present calibration study, while only 2 yr (2003 to 2005) were considered in (Blard et al., 2007).More important, the precipitation data used here were significantly revised (by ∼+20 %) to account for a bias in the precipitation gauge that was not considered in (Blard et al., 2007).The calibration of this type of model from glaciers in Sweden (Hock, 1999) and in Switzerland (Pellicciotti et al., 2005) also yielded ablation factors in the same order of magnitude than those obtained here.

Model 4
Only one study tested this model at an hourly timescale on a mountain glacier in Switzerland (Haut Glacier d'Arolla) (Pellicciotti et al., 2005).The MF factors calibrated in this study are significantly higher (about ×10) than those obtained by Pellicciotti et al. ( 2005), while the "a" factors are significantly lower.However, this comparison may be considered with caution since Pellicciotti et al. (2005) used measured net radiation and not clear-sky radiations.
The ablation processes of high latitude glaciers are quite different than those occurring in the Tropics (Hock, 2003;Sicart et al., 2008).These discrepancies in calibrated parameters may thus reflect differences in the involved physical ablation processes.This underlines that any paleoclimatic or projection studies should benefit from using parameters calibrated in similar conditions than those of the studied zone, as already suggested by Hock (2003).

Physical limitations of the temperature-index models
In spite of any improvement that could be made, these PDD models remain empirical models.They are not based on a real energy balance approach and they have not a good ability in catching the complex interactions between the glacier surface and the atmospheric variables (relative humidity, albedo, cloudiness, wind velocity) at a short timescale (Sicart et al., 2005(Sicart et al., , 2008)).However, the goal of the present study is to test and valid ablation models for paleoclimatic reconstruction, an exercise that requires keeping as low as possible the number of input data.The different models presented here should thus be used with caution, especially if the goal is to model the glacier mass balance under very different atmospheric conditions.For example, in very dry and low precipitation conditions, at high elevation (above 5500 m in the dry tropical Andes) sublimation becomes the main ablation process (Favier et al., 2004;Rupper and Roe, 2008;Wagnon et al., 2003).In the case of the Zongo glacier (precipitation > 1 m yr −1 ), the contribution of sublimation on total ablation is less than 10 % (Fig. 6 in Rupper and Roe, 2008).Sublimation however becomes the major ablation process under drier conditions.Indeed, the sublimation/melt ratio is higher than 1 when precipitation is lower than 0.5 m yr −1 (Rupper and Roe, 2008).In that particular case, classical PDD models (1 and 2) and models 3 will not be able to correctly capture the real amount of ablation.The structure of model 4 clearly represents a more realistic approach to calculate the mass balance of tropical glaciers fed by very low precipitations.2141 Ideally, the calibration of the "a" parameter of model 4 should be done using mass balance data from a glacier developed in a drier tropical zone.Unfortunately, there is a lack of mass balances records in such areas.

Sensitivity tests: precipitation-temperatures curves
We used the 4 calibrated models to determine the precipitation-temperature solution able to reproduce the present day-observed ELA at 5400 m (Fig. 6).This calculation was done using an "average year" calculated from the temperatures and precipitation data of the 9 hydrological years reported in Table 1.A 3 • C range in annual average temperature (+1.5 • C and −1.5 • C the present value) was explored to determine the corresponding precipitation conditions.This was done to evaluate (i) the uncertainty arising from the calibrated paramaters and (ii) the precipitation-temperature sensitivity of the models.The 4 obtained precipitation-temperature curves are characterized by different slopes, which reflects contrasted sensitivities of the models.Under warmer conditions (+1.5 • C), models including different ablation factors for snow and ice (Model 2 and 3) require lower precipitation (×3.2 for instead of ×3.7) to maintain the ELA at 5400 m a.s.l.than those having similar ablation factors for snow and ice.A comparison of models 1 and 2 with models 3 and 4 indicates that including solar radiation does not significantly change the precipitation temperature sensitivities.The error arising from the propagation of the parameters uncertainty is not significant.Indeed, even for models characterized by a high factor uncertainty (Model 3), the error envelop indicates that the temperature uncertainty remains below 1 • C, for a given precipitation.
The shape and slope of these obtained precipitation temperature curves can be compared with similar curves established by previous works.The curves yielded by the 4 models tested here match quite well with the curve established by Shi (2002) by plotting the summer temperature vs the annual precipitation measured at the ELA of several mid-latitudes glaciers from the Alps and Western China.Our models yield a dP/dT slope that is slightly lower than the average curve built by Ohmura et al. (1992) from a large worldwide dataset.The P vs. T curves predicted by our models remains however consistent with the precipitation-temperature sensitivity range reported by this author, suggesting that our approach does not yield unrealistic results.

Conclusions
The ablation parameters of several glacier-climate models were calibrated by using monthly time-scale meteorological and mass balance data from the Zongo glacier (Bolivia, 16 • S).The majority of the so-obtained degree-day ablation factors are within the range of other values obtained from other worldwide glaciers (Hock, 2003).The "short wave ablation factors" of models including solar radiation are quite different from the parameters calibrated from mid latitude glaciers (Hock, 1999;Pellicciotti et al., 2005).
The models presented here are based on a semi-empirical approach, which hampers precise and accurate physical interpretations.In particular, it is important to keep in mind that classical temperature-index models are not suited for high elevation (> 6000 m) tropical glaciers, where wind erosion and sublimation are the dominant ablation mechanisms (Wagnon et al., 2003).However, (i) our study shows that the inclusion of different melting parameters for snow and ice only induce a slight improvement of the model ability in fitting the observed mass-balance.Calibrated ablation parameters of snow are systematically smaller than those of ice, suggesting that the effect of albedo is, at least partially, accounted by the model structure.(ii) Models with short wave radiations offer the advantage to include a quantitative estimate of spatial and temporal changes in solar radiations.This property can be of major interest to quantify the effects of shading (e.g. in deeply incised valley) or the influence of orbital insolation changes.
Reconstructed precipitation-temperature curves able to fit a given ELA show that these models are able to yield temperature reconstructions with ∼ 1 • C uncertainty.
However, the temperature-precipitation sensitivity can be different between each model: models with different ablation parameters for snow and ice require lower precipitation at the same temperature.The performance of these empirical models could certainly be improved by further studies, notably by considering new and larger dataset acquired under different atmospheric and climatic conditions.Another way to improve these types of glacier-models would be to increase their sophistication.For example, it could be interesting to include a parameterization of cloudiness modulating the incoming solar radiation.Eventually, it could be appropriate to move toward energy balance models under certain favorable conditions.But it must be kept in mind that it is quite tough to increase the complexity of the model without adding extra parameters that cannot be constrained accurately for the past.Altitude (m) 1997-1998 1998-1999 1999-2000 2000-2001   Model 1      This calculation shows that the ablation arising from the a × I term becomes significant only for elevations higher than 5600 m.This calculation shows that the ablation arising from the a × I term becomes significant 8 only for elevations higher than 5600 m. 9 10 Used parameters are those calibrated from the present study (Table 4).Envelope shows the 1σ 6 error calculated by taking into account the parameters uncertainty given in Table 4. 7 Input temperature and precipitation are the mean of the 9 hydrological years.Used parameters are those calibrated from the present study (Table 4).Envelope shows the 1σ error calculated by taking into account the parameters uncertainty given in Table 4.

2121
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | , with a strong precipitation seasonality.The Zongo watershed has a total area of 3.3 km 2 .It is orientated toward the South East in the upper part (between 6000 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

2127
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | 1 ) is calculated for each month assuming precipitation falls as snow when temperature is below the rain-snow threshold T s .The model assumes 2129 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 3.3.2Model 1: A = MF × T p ) 2131 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | is the radiation-temperature melt factor and I i (W m −2 ) is the monthly average short wave solar radiation.In the present Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

2137
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

2139
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | considered with caution given the highest MF factors obtained in the present study are all determined from the 1997 to 2002 hydrological years.Indeed, these years of observations might be potentially affected by a lower reliability of the temperature data from the MEVIS stations.

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

2143
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

2149
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Table 2. Mass balance data (m yr −1 ) measured on the Zongo 1 glacier for the nine hydrological years between 1997 and 2006.Each hydrological year starts on 1 September and ends on 31 August.Data were collected by the ORE Glacioclim IRD program.
• C −1 d −1 ) Monthly precipitation a (mm W −1 m 2 d −1 ) Clear-sky shortwave radiations 2151 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Table 4. Calibrated parameters.R 2 (all years) are the correlation coefficient calculated using the whole set of modeled and observed mass balances of all hydrological years.R 2 pos and R 2 neg are calculated considering the positive and the negative mass balance, respectively.

Fig. 1 .Figure 2 .5
Fig. 1.Map of the Zongo glacier.Locations of the weather stations are shown, as well as those of the snow pits and the stakes used to measure the mass balance during the 2005-2006 year.2153

Fig. 2 .Figure 3 .Fig. 3 .
Fig. 2. Values of the efficiency criterion R 2 vs. ablation parameters (MF, a) for models 2 to 4 inferred for the 9 hydrological years between 1997 and 2006.Red star shows the best parameters, which yields the maximum R 2 value.2154

2155Fig. 4 .
Fig. 4. For each of the 4 models, plot of the observed vs. best-fit modeled annual mass balances for the 9 hydrological years between 1997 and 2006.
Figure 5. a) Mean annual positive temperature (Tp, in °C, calculated from the positive degree day approach) and mean annual clear sky incoming short-wave radiations on the Zongo glacier surface (I, in kW.m -2 , calculated taking into account the glacier slope and orientation) vs elevation (m).b) For Model 4, shown are the respective proportions of total ablation arising from the (a × I) and the (MF × Tp) terms vs elevation (m).The MF and a values used in the calculation are 11.8 mm.°C -1 .d - and 2.1 ×10 -4 mm.W -1 .m 2 .°C -1 .d -This calculation shows that the ablation arising from the a × I term becomes significant only for elevations higher than 5600 m.
annual positive temperature (Tp, in °C, calculated from the positive 2 degree day approach) and mean annual clear sky incoming short-wave radiations on the 3 Zongo glacier surface (I, in kW.m -2 , calculated taking into account the glacier slope and 4 orientation) vs elevation (m).b) For Model 4, shown are the respective proportions of 5 total ablation arising from the (a × I) and the (MF × Tp) terms vs elevation (m).The MF 6 and a values used in the calculation are 11.8 mm.°C -1 .d - and 2.1 ×10 -4 mm.W -1 .m 2 .°C- 7 1 .d -

Fig. 5 . 2 3
Fig. 5. (a) Mean annual positive temperature (T p , in • C, calculated from the positive degree day approach) and mean annual clear sky incoming short-wave radiations on the Zongo glacier surface (I, in kW m −2 , calculated taking into account the glacier slope and orientation) vs. elevation (m).(b) For Model 4, shown are the respective proportions of total ablation arising from the (a × I) and the (MF × T p ) terms vs. elevation (m).The MF and a values used in the calculation are 11.8 mm• C −1 d −1 and 2.1 ×10 −4 mm W −1 m 2 • C −1 d −1This calculation shows that the ablation arising from the a × I term becomes significant only for elevations higher than 5600 m.

Fig. 6 .
Fig. 6.Precipitation -Temperature curves able to maintain the ELA at 5400 m for models 1 to 4.

Table 1 )
were collected by two weather stations located in the vicinity of the Zongo glacier.Locations of the weather station are shown on Fig. 1.For the 6 hydrological years ranging between 1997 and 2003, monthly temperature and precipitation data were recorded by the MEVIS station, which is located at an elevation of 4750 m, at a distance of approximately 1 km from the glacier terminus (16 • 16.9 S, 68 • 7.3 W).The climatic data for the 3 hydrological years between 2003 and 2006 are from a weather station located at 5050 m (16• 17 S, 68 • 8.3 W), on the lateral moraine of the Zongo glacier (ORE GLACIOCLIM network), but far enough from the ice to avoid major perturbation due to change in the glacier surface conditions.Moreover, since the present model aims to infer regional temperatures from the inversion of past glacial extents, it is relevant to calibrate the model using local temperatures measured out of the Zongo glacier.Mean annual temperature is 0.4 ± 0.3 • C at 5050 m (ORE station -from 2003 to 2006) and 2.0 ± 0.7 • C at 4750 m (MEVIS station -from 1996 to 2003) (Table 1).The average annual amplitude is 3 ± 1 • C. Temperature standard deviations around monthly means

Table 1 .
Meteorological data near the Zongo glacier 1 for the 1997-2006 hydrological 2 yr.Hydrological years range from 1 September to 31 August.

Table 3 .
Models description, input data and 1 calibrated parameters.