The Energy and Mass Balance of Peruvian Glaciers

Peruvian glaciers are important contributors to dry season runoff for agriculture and hydropower, but they are at risk of disappearing due to climate change. We applied a physically based, energy balance melt model at five on‐glacier sites within the Peruvian Cordilleras Blanca and Vilcanota. Net shortwave radiation dominates the energy balance, and despite this flux being higher in the dry season, melt rates are lower due to losses from net longwave radiation and the latent heat flux. The sensible heat flux is a relatively small contributor to melt energy. At three of the sites the wet season snowpack was discontinuous, forming and melting within a daily to weekly timescale, and resulting in highly variable melt rates closely related to precipitation dynamics. Cold air temperatures due to a strong La Niña year at Shallap Glacier (Cordillera Blanca) resulted in a continuous wet season snowpack, significantly reducing wet season ablation. Sublimation was most important at the highest site in the accumulation zone of the Quelccaya Ice Cap (Cordillera Vilcanota), accounting for 81% of ablation, compared to 2%–4% for the other sites. Air temperature and precipitation inputs were perturbed to investigate the climate sensitivity of the five glaciers. At the lower sites warmer air temperatures resulted in a switch from snowfall to rain, so that ablation was increased via the decrease in albedo and increase in net shortwave radiation. At the top of Quelccaya Ice Cap warming caused melting to replace sublimation so that ablation increased nonlinearly with air temperature.

To elucidate the mechanisms of mass change of Peruvian glaciers and their sensitivity to climatic forcing, it is necessary to understand their energy balance and ablation characteristics. In the Cordillera Blanca distributed mass balance modeling over Shallap Glacier revealed that net shortwave radiation dominates the energy balance, with melt energy higher in the dry season due to greater dry season net radiation (Gurgiser et al., 2013). The phase and quantity of precipitation below the equilibrium-line altitude (ELA) and its influence on the snowpack was the main control on glacier wide mass balance, with sublimation more importantly in the dry season (Gurgiser et al., 2013). The importance of dry season sublimation was also found on Artesonraju Glacier, with measured sublimation increased under conditions of low specific humidity, moderate wind speeds, and rough surface conditions (Winkler et al., 2009). On Cuchillacocha Glacier, Aubry-Wake et al. (2018) identified the importance of the longwave flux from nearby exposed rock in increasing daily melt near the glacier margins (by 106 W m −2 on average).
In the Cordillera Vilcanota, analysis of the energy balance conditions on Quisoquipina Glacier revealed that the net shortwave radiation also dominated the energy balance, primarily controlled by albedo and hence by snowfall. The energy available for melt was highest between September and December (Suarez at al., 2015). Perry et al. (2017) investigated the character of the precipitation on Quelccaya Ice Cap, finding that the diurnal precipitation distribution had a main afternoon and secondary nighttime peak, and was nearly exclusively solid phase. Graupel (heavily rimed snow crystals associated with convective precipitation) occurred frequently, comprising more than 30% of precipitation hours in all months except May (Perry et al., 2017).
Despite the advances of these individual studies, they all focused on a single site, occasionally comparing with data from Bolivia and Ecuador (Perry et al., 2017;Suarez et al., 2015). A better understanding of glacier contribution to regional water supply necessitates the assessment of the energy balance and ablation characteristics across Peruvian glaciers covering a variety of glacier and climate settings.
To address this, we focused on the Cordillera Blanca and Cordillera Vilcanota (the largest and second largest tropical glaciated mountain ranges in the world, respectively [Drenkhan et al., 2015]), which have never previously been the subject of a comparative study. We collated data from five glaciers, three in the Cordillera Blanca and two in the Cordillera Vilcanota ( Figure 1) and applied the energy and mass balance module of the land surface model Tethys-Chloris at the point scale to all five records. To contextualize the results from these five sites we also investigate the climatology and glaciology of their catchments: the upper Rio Santa catchment in the Cordillera Blanca (4,953 km 2 , draining to the Callejon de Huaylas), and a catchment containing the Cordilleras Vilcabamba, Urubamba and Vilcanota, and the Quelccaya Ice Cap (11,048 km 2 , draining to the Rio Urubamba at Santa Maria, henceforth the Vilcanota catchment) (Figure 2). This study aims to understand the key drivers of the energy and mass balance of these five glaciers and to determine the differences between the Cordillera Blanca and Cordillera Vilcanota. To do this we will address the following objectives: 1. Determine the relative importance of each of the energy fluxes on melt rates of the five glaciers. 2. Determine the temporal variability and spatial differences between sites in the magnitude of ice melt, snow melt, and sublimation.  Table S3 in Supporting Information S1. Panel (d) gives an overview of Peru and the location of the Cordilleras Blanca and Vilcanota, with panels (a) to (c) showing the study sites in the Cordillera Blanca and panel (e) showing the study sites in the Cordillera Vilcanota. The elevation data and underlying hillshade are from ASTER GDEM v.3, which is a product of METI and NASA. The glacier and lake outlines are from the National Peruvian Glacier and Lake Inventories, respectively. Garreaud et al., 2009;Marengo et al., 2012;Vera et al., 2006). By contrast, the annual temperature range is small, unlike the large diurnal temperature range (Juen et al., 2007;Kaser et al., 2003).
The climate of the Peruvian Andes is considerably influenced on interannual timescales by the El Niño-Southern Oscillation (ENSO) (e.g., Vuille et al., 2003).  Thompson et al., 2017). However, the effect of ENSO fluctuations on precipitation is spatially variable across Peru, with the sign of precipitation change varying with latitude, El Niño regime, and in southern Peru, with distance to the Amazon (Jonaitis et al., 2021;Lavado-Casimiro & Espinoza, 2014;Sanabria et al., 2018) There is strong evidence that air temperatures in the tropical Andes will increase in the future, with a warming of 4.5-5°C predicted by the end of the 21st century under the SRES A2 scenario, which will be enhanced over the Cordillera Blanca (Urrutia & Vuille, 2009;Vuille, Francou, et al., 2008). Yarleque et al. (2018) projected a warming over the Quelccaya Ice Cap of 2.4°C under RCP 5.4 and 5.4°C under RCP 8.5 by the end of the 21st century. However, there is less certainty in the sign of precipitation changes. Enhanced December-January-February (DJF) precipitation is expected over southern tropical South America, with a precipitation increase more likely north of ∼12°S, but a decrease further south (Urrutia & Vuille, 2009). Strong drying is suggested over the central Andes with Neukom et al. (2015) expecting a DJF precipitation decrease under RCP 8.5 of 19%-33% and Minvielle and Garreaud (2011) a 10%-30% reduction under moderate to strong emissions. This precipitation reduction will likely lead to both reduced accumulation and increased energy for ablation over the glaciers of the Cordillera Vilcanota (Kronenberg et al., 2016).
To investigate the climate of the two study regions we analyze long-term, bias-corrected outputs from the Weather Research and Forecasting (WRF) model version 3.8.1, which was run from 1980 to 2018 over each of the study catchments at a 4 km horizontal resolution (Text S1.1 in Supporting Information S1). Glaciers in the Rio Santa catchment are only found on the mountains on the eastern side (Cordillera Blanca), whereas the glaciers in the Vilcanota catchment are found in the north (Cordilleras Vilcabamba and Urubamba) and south-east (Cordillera Vilcanota) (Figures 2a and 2b). The spatial variation in air temperature across both regions is primarily controlled by elevation (Figures 2a and 2b), with slightly cooler dry season air temperatures and a greater seasonal difference in Vilcanota than in the Rio Santa ( Figure 3a). Due to the higher elevations of the Vilcanota compared to the Rio Santa glaciers (Figures 3e and 3f), dry season air temperatures over the glacier areas are lower in the Vilcanota region, but in the wet season are similar ( Figure 3a).
Spatial patterns in precipitation in the Rio Santa are primarily determined by elevation, with the highest precipitation over the Cordillera Blanca ( Figure 2c). In the Vilcanota catchment precipitation totals are much higher over the Cordilleras Vilcabamba and Urubamba than over the Cordillera Vilcanota toward the south ( Figure 2d). Here, precipitation comes most commonly from the northwest, with moisture originating from deep moist convection over the Amazon lowlands (Perry et al., 2014). The north-south spatial precipitation gradient therefore results in a lack of an elevation dependent precipitation gradient ( Figure 3d).
The lower total and wet season precipitation within the Vilcanota catchment ( Figure 3b) results in a comparably less favorable climate for the existence of glaciers compared to the Rio Santa. This is also reflected in the higher off-glacier wet season snowlines in the Vilcanota ( Figure 3e). Consequently, the minimum and median glacier elevations are higher in the Vilcanota than in the Rio Santa (Figures 3e and 3f).

Study Sites and Data
In this study, we model the energy and mass balance at point locations on five glaciers across the two glacierized regions described above: the Artesonraju (AG), Shallap (SG), and Cuchillacocha (CG) Glaciers in the Rio Santa region, and the Quelccaya Ice Cap (QIC) and Quisoquipina (QQG) Glacier in the Vilcanota region ( Figure 1 and Table S3 in Supporting Information S1). All the Rio Santa stations and the Quisoquipina station are near the glacier snout, whereas the Quelccaya Ice Cap station is at the top of the accumulation area ( Figure 3f). Figure 3f shows that the Vilcanota study glaciers do not reach down to the elevation of the Rio Santa stations. Hourly data from on-glacier meteorological stations were available at all five glaciers, but not all on-glacier records were complete or included all necessary variables to run the model. All data were carefully quality checked, and then gaps in the on-glacier record were filled using data from near-by off-glacier stations whose data were transformed to be relevant to the on-glacier site. Conversely, precipitation data were measured off-glacier at all sites except Quelccaya. Where data gaps remained, they were either filled by a day of data on the same date of a randomly selected year (if the records were long enough) or data were derived from the WRF model outputs. Full details of the station instrumentation and data preprocessing are given in Text S2 in Supporting Information S1. 3, which is a product of METI and NASA. See Text S1.2 in Supporting Information S1 for details of the methodology to derive the snowlines and Table S2 in Supporting Information S1 for the data.
M sno E = 450. These values are at the high end of the literature ranges, but resulted in a good comparison with the SR50 data and were chosen based on the understanding that initial snow densities in the ablation area of tropical glaciers are likely to be relatively high (Sicart et al., 2002). At Cuchillacocha, Quelccaya and for 8% of the Artesonraju record, measured albedo and incoming longwave radiation were not available, and we used parameterizations that introduced additional empirical parameters. Albedo was parameterized following Brock et al. (2000), introducing the parameters describing the albedo decay and t E Pr , the threshold precipitation required to reset the snow albedo to that of new snow. Incoming longwave radiation was parameterized with the Dilley and O'Brian (1998) and Unsworth and Monteith (1975) equations, introducing further empirical coefficients (Text S3.1 in Supporting Information S1). For Cuchillacocha and Artesonraju we used the standard albedo and longwave parameter values as given in the literature. At Cuchillacocha we checked the suitability of these parameterizations by comparing modeled net radiation with the 2 weeks of measured net radiation available (Text S3.9 and Figure S4 in Supporting Information S1). Calibration was necessary at Quelccaya since the deep snowpack meant results were sensitive to the snow density and initial runs suggested modeled mass loss was too high. We therefore calibrated t E Pr (best value = 5 mm) and both the snow density parameters against SR50 snow depth data (best value for both = 510 kg m −3 ). Although physically  1 the best results at Quelccaya were when they were equal, which may be due to the rarity of melting conditions or the deep snowpack.
To understand the model uncertainty due to variations in model parameters we performed a multiparameter Monte-Carlo analysis at all five sites (the parameters and their ranges are given in Table S8 in Supporting Information S1). We performed 1,000 model runs with randomly selected parameters from within assumed ranges. Note, it was not possible to include variations in the Unsworth and Monteith (1975) parameters since this caused instabilities in the calculation of surface temperature. The results of the Monte-Carlo analysis are given in Text S4.2 and Figure S6 in Supporting Information S1.

Quantifying the Impact of ENSO and Sensitivity to Climate
To determine the impact of ENSO on the energy and mass balance at the study sites, the NOAA ONI (Oceanic El Niño Index) was correlated against the Artesonraju and Quisoquipina model outputs, since these sites had long records (almost 7 years for Artesonraju, ONI range −1.7°C-1.6°C, and almost 5 years for Quisoquipina, ONI range −1.1°C-2.6°C). The variables were averaged over rolling 3-month periods, to match the temporal resolution of the NOAA ONI index, in which positive values are associated with El Niño conditions and vice versa. The analysis was conducted for all years in each rolling 3-month period. Linear regression was conducted if the correlation was significant at p < 0.05.
To understand the climate sensitivity of each of the sites, we reran the model with changes to the air temperature (0-6°C) and precipitation (−30% to +30%) inputs (both separately and combined). Unlike in the standard runs, at all sites we used modeled albedo, so that changes in albedo due to changes to the snow cover duration could be accounted for. The vapor pressure and dew point temperature were recalculated after changes were made to the air temperature, to make sure that there was consistency between them. This results in an increase in the vapor pressure with an increase in air temperature and vice versa. This will have an impact on sublimation rates and the precipitation partition, since the Ding et al. (2014) precipitation partition implemented in T&C depends on the vapor pressure (Text S3.7 in Supporting Information S1). At Quelccaya, the initial snow pack depth of 1,000 mm w.e. was used for all the climate sensitivity runs. During warm conditions the model snowpack was completely melted, initiating the melting of the ice pack underneath. There is 20-30 m of firn beneath the Quelccaya station (Salzmann et al., 2013), so to replicate this the ice albedo was set to 0.53, a reasonable value for firn (Oerlemans & Knap, 1998).

Model Validation
To assess model performance the outputs were compared against surface height change measured using an SR50 sensor (for Shallap, Artesonraju and Quelccaya), ablation stakes (Shallap, Artesonraju, Cuchillacocha) and the surface cover (ice or snow) as derived from albedo measurements (Shallap, Artesonraju, Cuchillacocha, Quisoquipina), and satellite imagery (Cuchillacocha). Full details of the validation are given in Text S4.1 and Tables S9-S11 in Supporting Information S1.
The modeled surface height matches the SR50 records well ( Figure S5 in Supporting Information S1), with RMSE values in the order of tens of centimeters in most periods. The average RMSE across all periods is 0.360 m at Shallap, 0.220 m at Artesonraju, and 0.141 m at Quelccaya. There is a tendency for the model to underestimate surface height (overestimate mass loss) at Shallap (mean ME across periods is −0.251 m) and to overestimate the surface height at Artesonraju (mean ME is 0.169 m), whereas the overall bias at Quelccaya is very small (mean ME is 0.004 m).
Stake data were compared against cumulative ice mass loss (due to melt and sublimation) since stakes were measured when the ice was snow-free. The stake data also confirm the good model performance ( Figure S5 and Table S10 in Supporting Information S1). At Shallap, modeled mass loss was between the values of the two closest stakes in 2011, although melt was underestimated compared to the closest stake to the station (S7), whereas in 2012 errors were small, with the mean absolute difference equating to 0.077 m w.e. At Artesonraju, modeled mass loss was underestimated compared to the A14 stake in 2011, whereas modeled mass loss fell between the range of the 2012-2013 stake measurements. At Cuchillacocha, modeled mass loss is high compared to the mean of the 12 stakes measured over a 2-week period (the average difference is 0.055 m w.e.), although it still fits within the range. The reason for the large range in stake values is likely observed ice weathering, resulting in cryoconite holes, which had not collapsed at all stakes. T&C cannot account for ice weathering processes and small-scale spatial variability in ablation.
Comparisons of measured and modeled surface type showed that the surface was mainly modeled correctly at Shallap (85.5% agreement), Artesonraju (65.7% agreement), and Cuchillacocha (76.6% agreement using 2 weeks of measured albedo and 68.4% agreement when the surface was compared to Planet Labs imagery). At Quisoquipina, the agreement between simulated and observed surface was lower (54.3%), with a bias toward modeling ice when the surface should be snow. This could be explained by the evening snowfall melting through the day to leave bare ice in the late afternoon, which would affect the modeled surface type but have little influence on measured albedo.

Micro-Meteorology
Here, we use the model simulations and input meteorological data to characterize the micro-meteorology at the sites. Wet season air temperatures are warmer than in the dry season, with the difference between wet and dry season air temperatures greater at the Vilcanota sites than the Rio Santa sites ( Figure S7 in Supporting Information S1). The air temperature at all Rio Santa sites (Shallap, Artesonraju, and Cuchillacocha) is similar, although shading delays the morning rise in air temperature at Cuchillacocha. Air temperatures are colder at Quisoquipina, and colder still at Quelccaya due to their higher elevations. Peak wind speeds are much higher at Cuchillacocha (7.7 ms −1 ) than at the other sites (3.0-4.5 ms −1 ), with wind speeds in general higher in the dry compared to wet season.
The wet season is characterized by much higher precipitation amounts compared to the dry season, with mean wet season precipitation ∼3 to 6 times dry season precipitation. There is a clear afternoon peak in wet season precipitation, although the peak is earlier (∼15:00-16:00) and smaller in magnitude in the Vilcanota compared to the Rio Santa (at 18:00). The wet season is also characterized by much higher relative humidity (wet season diurnal minimum ranges from 66%-77% compared to 54%-61% in the dry season) with a larger diurnal amplitude compared to the dry season. The higher atmospheric water content and cloud cover in the wet season also results in higher incoming longwave radiation, which has a lower diurnal amplitude than in the dry season at the sites where it has been measured (Shallap, Artesonraju, and Quisoquipina). Incoming longwave radiation is lowest at Quelccaya due to its higher elevation.
The diurnal pattern of incoming shortwave radiation is similar across all the sites, except at Cuchillacocha where shading delays the morning rise in incoming shortwave radiation. Incoming shortwave radiation is either similar (Vilcanota) or higher (Rio Santa) in the dry season compared to the wet season, suggesting that the greater cloudiness in the wet season counteracts the higher potential radiation receipt.

Contribution of Energy Balance Fluxes
Net shortwave radiation is the most important flux for providing energy to the surface, in all months and seasons (Figures 4 and 5). The net shortwave radiation is higher in the dry season at all sites (Figures 4 and 5), suggesting that increased wet season cloudiness and a higher surface albedo due to snow cover counteracts the higher wet season potential radiation receipt. This is especially true at Shallap, which has the largest difference in wet and dry season net shortwave radiation. The month with the highest net shortwave radiation tends to be in late winter/ spring (August at Shallap and Cuchillacocha, September at Quelccaya, and October at Quisoquipina), when low albedo is combined with an increasing radiation receipt (Figure 4). The exception is Artesonraju where net shortwave radiation is highest in January.
The sensible heat flux provides energy to the surface at the lower elevation Rio Santa sites, although it is much smaller in magnitude than net shortwave radiation. Values are greatest at Cuchillacocha due to its high wind speeds and here the wet season sensible heat flux is greater than in the dry season, with the highest monthly mean in November. At Shallap and Artesonraju, the sensible heat flux is slightly higher in the dry season, with the highest monthly mean in May. At Quisoquipina and Quelccaya, the sensible heat flux is rather small (0.7 W m −2 and -1.0 W m −2 , respectively, averaged over the full record).
The energy flux leading to the greatest loss of energy from the surface is net longwave radiation. At all sites it is more negative in the dry season (Figures 4 and 5), reflecting less cloud cover in the dry season reducing incoming longwave radiation. However, it is also a positive flux at Quisoquipina in December, January, and February (Figure 4), because of increased incoming longwave radiation counteracting the outgoing longwave radiation, which does not vary much during the year. The net longwave radiation is more negative at Quelccaya compared to Quisoquipina due to a combination of colder surface temperatures and lower incoming longwave radiation than Quisoquipina in the wet season. Diurnal values tend to be most negative in the early morning and less negative or even slightly positive in the middle of the day ( Figure S8 in Supporting Information S1).
The latent heat flux results in the next greatest loss of energy from the surface. It is a negative flux for all sites and in both seasons, except at Cuchillacocha in the wet season. At all sites, it is more negative in the dry season (Figures 4 and 5). It tends to be most negative in the morning hours and decreases in magnitude during the day ( Figure S8 in Supporting Information S1).
The remaining fluxes are all small contributions to the surface energy balance once averaged at a monthly resolution. The heat flux due to precipitation ( v E Q ) is less than 1 W m −2 , with usually higher values in the wet than in the dry season, and peaks in the afternoon. The monthly average heat flux due to melting or freezing of water in the snowpack ( fm E Q ) is always less than 0.1 W m −2 . The sign of this flux varies during the day, being negative in the morning due to melting and positive in the evening due to refreezing ( Figure S8 in Supporting Information S1). The ground heat flux ( E G ) is slightly larger but the monthly average is always less than 2 W m −2 in magnitude.
The energy available for melt is usually greater in the wet than the dry season ( Figure 5), with the highest mean monthly melt in January for Artesonraju and Quelccaya, February for Cuchillacocha, and November for Quisoquipina ( Figure 4). The exception is Shallap where the melt is greatest in August. This demonstrates that the greater dry season net shortwave radiation is balanced by the increase in the net longwave radiation (due to cloudless skies) and the latent heat flux (due to low relative humidity). Conversely, during the wet season, the net longwave radiation and latent heat flux are small (at least at Shallap, Artesonraju, and Quisoquipina) so the net shortwave radiation alone drives variations in the energy for melt, and variations in albedo (driven by the occurrence and magnitude of snowfall) determine variability in wet season melt. At Shallap, melt is lower in the wet season than in the dry season, because its continuous snowpack resulted in a high albedo, reducing net shortwave radiation, and its annual energy for melt.

Patterns of Melt and Sublimation
Ice and snow melt occurs in all months of the year at Artesonraju, Cuchillacocha, and Quisoquipina, with the fraction of ablation that is snow melt in the wet season usually less than 50% ( Figure 6). This is a result of the snowpack at these sites not being continuous throughout the wet season, but instead being continually melted and reformed over periods of days to weeks. The snowpack tends to be quite shallow, with the maximum mean monthly snow depths being only 0.13 m at Artesonraju (in March), 0.07 m at Cuchillacocha (in November), and 0.02 m at Quisoquipina (in February).
Despite the shallow snow depths, the snow cover is still very important for the overall energy and mass balance at these sites. Unlike the sites above, Shallap had a continuous snowpack over the 2010/2011 wet season (mean wet season percentage snow melt of ablation was 91%), with a maximum mean monthly snow depth of 0.61 m in January. The higher surface albedo resulted in lower wet season melt rates than the dry season, in contrast to the four other sites (Figures 5 and 6). This is likely an exceptional situation for Shallap caused by below average wet Mass loss is only from snow melt and sublimation at our highest site, Quelccaya. Snow melt is very small, and occurs in only 7 months of the year, with the maximum mean monthly snow melt only 0.66 mm w.e. day −1 in January. Sublimation, on the other hand accounts for 42%-100% of mass loss in each month, with the maximum mean monthly sublimation of 0.96 mm w.e. day −1 occurring in July.
At the lower elevation sites sublimation is a much smaller contributor to mass loss than ice or snow melt, but it does occur in all months at Shallap, Artesonraju, and Quisoquipina, and at Cuchillacocha in January and May to November. Sublimation is highest in the dry season, with the proportion of ablation which is sublimation Figure 6. Ablation from T&C split into mean monthly ice melt, snow melt, and sublimation for each month. The averages are applied over the whole record. The modeled record lengths for each weather station are given in Table S3 in Supporting Information S1. Also shown is the fraction of snowmelt of total ablation, the fraction of sublimation of total ablation, the albedo (monthly mean of daily minimum values), and on the right y-axis the mean monthly snow depth. Note that the actual snow depth at Quelccaya Ice Cap will be much deeper than the values shown, since these are a function of the snow depth used to initiate the model run, and over the full model run the snow depth increases. However, the relative seasonal variation of the snowpack depth is correct. reaching a maximum in July, with values of 6% for Shallap, 8% for Artesonraju, 16% for Cuchillacocha, and 27% for Quisoquipina. Although the quantity of mass loss directly attributed to sublimation is small, it reduces overall mass loss because the strongly negative dry season latent heat flux reduces the energy available for melt (Section 5.3).
Mean ablation decreases with elevation, which is especially clear in the dry season when the snowpack at Shallap does not influence results (Figure 7a). The gradient of change in ablation with elevation is steeper in the dry season at the Rio Santa sites, compared with the Vilcanota sites, suggesting that ablation rates in Vilcanota are higher for a given elevation than in the Rio Santa. However, it would be beneficial to have more sites with which to compare given the small elevation range of the Rio Santa sites and the possibility of site-specific factors influencing results. The decrease in ablation with elevation across all sites in the full time series is driven by a decrease in melt rates. Mean sublimation rates are similar at all the sites.
Interestingly, the proportion of ablation comprised of ice and snow melt is rather similar between the Rio Santa sites and Quisoquipina (Figure 7b). The mean ice (snow) melt percentage for the Rio Santa sites is 72% (26%), compared to 68% (29%) for Quisoquipina. This suggests that despite Quisoquipina being ∼400 m higher in elevation than the Rio Santa sites, and having lower overall melt rates, the snowpack conditions are similar. This is likely caused by the lower precipitation rates at Quisoquipina compared to the Rio Santa sites.
The 479 m increase in elevation between Quisoquipina and Quelccaya is associated with a decrease in mean snow melt rates (from 2.57 mm w.e. day −1 to 0.15 mm w.e. day −1 ), with the percentage snow melt of total ablation falling from 29% to 19%. However, the key difference is in the shift in the importance of sublimation to total ablation, which increases from 4% at Quisoquipina to 81% at Quelccaya (Figure 7).

Impact of El Niño/La Niña
At Artesonraju, positive ONI values (associated with El Niño conditions) are associated with increased wet season total melt (significant relationships in DJF and JFM) ( Figure S9 in Supporting Information S1). This is caused by an increased wet season air temperature (significant positive relationship with ONI from September to May), which drives an increase in the sensible heat flux (significant December to April). There is also an increase in rainfall (significant November to March), which reduces the albedo (significant in DJF and JFM) and increases net shortwave radiation (significant in JFM). The correlation of precipitation with ONI in the wet season is usually positive but not significant.
ONI values at Quisoquipina are also positively correlated with increased total melt (significant in MAM), but the positive correlation of air temperature with ONI is not significant and the sign of the relationship of precipitation with ONI is less clear with no months with significant correlations ( Figure S9 in Supporting Information S1). ONI values are negatively correlated with albedo, which increases net shortwave radiation (both significant in DJF) and net longwave radiation is negatively correlated with ONI (significant in SON). There is also a significantly negative relationship between wind speed and ONI in DJF, AMJ and NDJ, which decreases sublimation from snow (significant in NDJ).

Sensitivity to Air Temperature and Precipitation
Results of the model sensitivity experiments to changes in the air temperature and precipitation show that, as expected, an increase in air temperature increases ablation (Figure 8). However, the drivers of this differ with elevation. At Shallap, Artesonraju, Cuchillacocha, and Quisoquipina ice melt is increased at the expense of snowmelt because of a change in the precipitation phase from snow to rain. A 2°C increase in air temperature resulted in a decrease of snowfall of total precipitation by 44% at Shallap, 53% at Artesonraju, 41% at Cuchillacocha, 17% at Quisoquipina, and 0.8% at Quelccaya, resulting in lower albedos by 0.09, 0.11, 0.11, and 0.03 for the four lower sites, while albedo remains similar at Quelccaya. Increased net shortwave radiation therefore drives the increase in ablation with air temperature at the lower sites ( Figure S10 in Supporting Information S1).
At the highest elevation site, Quelccaya, the rate of increase in snowmelt with air temperature steepens dramatically between 1 and 2°C. This corresponds with when sublimation switches from increasing to decreasing and marks the air temperature threshold at which sublimation is replaced by melting (Figure 8a). This results in ablation increasing nonlinearly with air temperature at Quelccaya (a 1°C air temperature rise enhances ablation by 0.03 mm w.e. h −1 , but an additional 1°C rise increased ablation by an additional 0.20 mm w.e. h −1 ). This corresponds to the percentage sublimation of total ablation decreasing from 81% under the standard run to 9.7% with a 2°C air temperature rise. Since the energy requirements for melt are much less than for sublimation, this results in the rapid increase in mass loss.
Higher precipitation reduces ablation and vice versa ( Figure 8). However, the impact of increasing precipitation to decrease ablation is reduced at higher air temperatures because the percentage snowfall of total precipitation is reduced. However, our sensitivity experiments do not account for the likely change in cloudiness and humidity associated with increased precipitation, which would impact net radiation and the latent heat flux. The difference in the specific mass balance under each of the experiments (Figure 8c) shows that under a 30% increase in precipitation enhanced specific mass loss compared to the standard run occurs when the air temperature increase is at least 1°C at most of the sites and at least 2°C at Quelccaya.

What are the Key Drivers of the Energy and Mass Balance of Peruvian Glaciers?
The key driver of the energy and mass balance of Peruvian glaciers is the precipitation dynamics. This is because net shortwave radiation is the most important contributor to energy available for melt ( Figure 5) and its magnitude is determined by albedo and cloudiness: despite higher top-of-the-atmosphere incoming shortwave radiation in the wet season, net shortwave radiation was lower than in the dry season due to increased cloudiness and higher wet season albedo ( Figure 5). Obviously, there is a close link between precipitation and cloudiness, and increased snowfalls maintain a high snow albedo and reduce the possibility that the snowpack is melted completely. This is encapsulated by the situation at Shallap where a continuous snowpack over the 2010/2011 wet Figure 8. T&C sensitivity experiments showing the influence of the increase in air temperature (row a) and precipitation (row b) on snow melt, ice melt, total melt, sublimation, and ablation. Note that the scenarios here involve change in only air temperature or precipitation, the other variable remains unchanged. For Quelccaya, ice melt corresponds to the melt of firn, which was not represented as a separate compartment in the model. Row (c) shows difference in the specific mass balance (in m w.e. per year) between the standard runs and each of the air temperature and precipitation scenarios. The mass balance in this case is the difference in the water equivalent mass between the first and last time step of the model run, which may not be equivalent to over a hydrological year due to differences in the model run time periods.
season dramatically reduced net shortwave radiation and energy available for melt. Meanwhile, at Artesonraju, Cuchillacocha, and Quisoquipina the difference between the wet and dry season net shortwave radiation (and energy available for melt) was much smaller because of their ephemeral wet season snowpack, which allowed ice to be exposed in the wet season (Figures 5 and 6). The importance of precipitation gleaned from an energy balance perspective supports 's findings that glacier mass balance fluctuations in the Cordillera Blanca are driven by interannual precipitation variability.
The dynamics of the wet season precipitation are therefore likely very important in determining annual variations in total melt at a given location. The wet season precipitation in the Peruvian Andes is driven by the rapid southward shift of the ITCZ in October , associated with easterly airflow bringing moisture from the Amazon across to the west of South America (Mourre et al., 2016). Within the wet season, the weekly variations in the snowpack depth are likely driven by the episodic nature of precipitation, which in the central Andes alternates weekly between wet and dry conditions (Garreaud, 1999;Guy et al., 2019).
The afternoon precipitation in the Cordillera Vilcanota is due to convective precipitation resulting from diurnal heating of the planetary boundary layer. But nocturnal precipitation is stratiform and instigated by cooling of the troposphere as the atmosphere stabilizes after convective storms so that the moisture bearing layer reaches saturation, forming thickening clouds (Perry et al., 2014(Perry et al., , 2017. Later in the evening thick cirriform clouds spreading from Amazonian convection precipitate into the moisture bearing layer, resulting in snow over the mountains and rain at lower levels. Most precipitation events in Cusco (95%) originate in the Amazon basin (trajectories from the north or northwest) -indicating the importance of Amazonian moisture to the Cordillera Vilcanota and central Andes (Perry et al., 2014).
Fluctuations in ENSO between warm (El Niño) and cold (La Niña) phases are the main drivers of interannual precipitation and air temperature variability across South America ). However, the influence of ENSO on precipitation in Peru varies between the north and south of the country, with distance from the Amazon and according to the El Niño regime (Jonaitis et al., 2021;Lavado-Casimiro & Espinoza, 2014;Sanabria et al., 2018). At Shallap, Maussion et al. (2015) found a strong anticorrelation between glacier mass balance and ENSO, with El Niño years associated with reduced snowfall and increased net shortwave radiation, which generally agrees with the results for Artesonraju (Section 5.5). The strong La Niña conditions during the 2010/2011 wet season in our study period at Shallap are also likely the cause of the continuous snowpack and low wet season melt rates which contrast with the mean conditions at the other sites (Section 5.4). The more complex relationship of ENSO with precipitation in southern Peru (Jonaitis et al., 2021) may explain the lack of a clear relationship between precipitation and ONI at Quisoquipina in the Cordillera Vilcanota (Section 5.5).
Climate change will determine the long-term evolution of Peruvian glaciers, and here the positive trend in air temperatures is rather certain, compared to the spatial variability in precipitation trends (Urrutia & Vuille, 2009). The direct impact of air temperature via the sensible heat flux is a less important determinant of melt rates due to the relatively cold air temperatures found on Peruvian glaciers ( Figure 5). However, our sensitivity experiments suggest that warmer air temperatures will change the precipitation phase from snow to rain, especially at lower elevations, acting to decrease albedo and increase net shortwave radiation, reducing the snow cover that protects the glacier ice from melting ( Figure S10 in Supporting Information S1). Importantly, at locations currently in the accumulation zone (e.g., the Quelccaya site) warmer air temperatures will allow melting to replace sublimation, resulting in a nonlinear increase in ablation with air temperature (Section 5.6). This reduction in sublimation with increased air temperature is seen at all the stations, but at a 2°C air temperature increase it only becomes more important annually than the increase in net shortwave radiation (due to the change in precipitation phase) at Quelccaya ( Figure S10 in Supporting Information S1).
Both longwave radiation and the latent heat flux (resulting in sublimation) are important in reducing dry season ablation on Peruvian glaciers ( Figure 5). Cloudless conditions in the dry season limit incoming longwave radiation, resulting in negative net longwave radiation, with values most negative overnight and at higher elevations ( Figure S8 in Supporting Information S1). High dry season sublimation rates are driven by low relative humidity, but between-site variations are also influenced by local wind fields and surface roughness variations. Wind speeds are higher over Cuchillacocha and Quelccaya compared to their neighbors, increasing dry season sublimation rates ( Figure S7 in Supporting Information S1 and Figure 4). At Quelccaya this is likely due to its exposed position at the summit of the Quelccaya Ice Cap, whereas at Cuchillacocha funneling of the prevailing synoptic winds into the north-south oriented valley are a possible cause. The surface roughness also controls sublimation rates, and it is known to be higher in the dry season over penitentes and lower over smooth fresh snow surfaces (Favier, Wagnon, Chazarin, et al., 2004;Wagnon et al., 1999;Winkler et al., 2009). Results of the Monte-Carlo analysis (Text S4.2 and Figure S6 in Supporting Information S1) suggest that varying the surface roughness between 0.001 and 0.05 m (covering the range of values mentioned in Winkler et al. (2009) (without penitentes), Favier, Wagnon, Chazarin, et al., 2004, Sicart et al., 2011, and Gurgiser et al., 2013 varied modeled total melt by 17% on average. Incorporating the spatial and temporal variability of surface roughness into future modeling would be a useful advance.

What Are the Differences Between the Two Regions?
The key glaciological difference between the two regions, as shown in Section 2, is that glaciers in the Rio Santa catchment extend to lower elevations than those in the Vilcanota catchment. This is possible primarily because of the greater wet season precipitation at the elevation of the glaciers in the Rio Santa, compared to the Vilcanota. This is a consequence of the northwesterly precipitation trajectory in the Vilcanota region (Perry et al., 2014), which means that the glaciers further south receive less precipitation, and explains the lack of a clear dependence of precipitation on elevation (Figures 2 and 3). The higher Vilcanota glacier elevations result in mean wet season air temperatures over the glaciers being similar to the Rio Santa (for the same elevation they would be warmer), but they are colder in the dry season due to the greater amplitude of the seasonal air temperature cycle (Figure 3a).
Contrasting the conditions at Quisoquipina with those at the Rio Santa sites allows the condition in the ablation zone of both regions to be compared, since all these stations are near the snouts of their respective glaciers (Figure 3ef). Many of the meteorological inputs are similar at Quisoquipina to the Rio Santa sites ( Figure S7 in Supporting Information S1). However, the cooler dry season air temperature and lower dry season net shortwave radiation reduce dry season ablation compared to in the Rio Santa (Figure 7). Both ice and snowmelt occur within the wet season at Quisoquipina, similar to Artesonraju and Cuchillacocha but at lower rates than those stations ( Figure 6). This results in a similar wet season percentage of ice and snowmelt ( Figure 6) and is a function of the lower precipitation rates in the Vilcanota region.
Based on our analysis there is more energy available for melt at a given elevation for the Vilcanota glaciers compared to the Rio Santa glaciers. This is evident from Figure 7 where ablation is very similar in the wet season between the Rio Santa sites and Quisoquipina, and in the dry season, where the gradient of the change in ablation with elevation is much steeper in the Rio Santa sites than the Vilcanota sites. There is of course uncertainty in this due to the small number of stations with which to compare, but the overall lower precipitation rates and warmer wet season air temperatures in the Vilcanota catchment support this pattern.

Peruvian Glaciers in a South American Context
The energy balance characteristics of the Peruvian study glaciers can be compared to those in Ecuador, Bolivia, and Chile in Table S12 in Supporting Information S1 and Figure 9. Some caution in direct comparison is required given the differences in the record lengths and periods, energy flux calculation approaches (see Ayala et al., 2017;Favier, Wagnon, Chazarin, et al., 2004;Favier, Wagnon, & Ribstein, 2004;Schaefer et al., 2020;Wagnon et al., 1999), and the relative station elevation compared to the glacier ELA. Nevertheless, the Peruvian glaciers' average energy balance characteristics are similar to Antizana Glacier 15 in Ecuador, Zongo Glacier in Bolivia, and the Chilean glaciers from Guanaco to Bello. At these sites the net shortwave radiation dominates the energy balance (average across the glaciers Antizana to Bello 169 W m −2 ), the net longwave radiation is strongly negative (average −63 W m −2 ), the sensible heat flux is small (average 25 W m −2 ), and the latent heat flux is negative (average −41 W m −2 ). Mean air temperatures are also close to or below 0°C (average −0.8°C), except for Juncal Norte. Due to the importance of the net shortwave radiation, albedo variations drive temporal variations in melt on the Antizana 15 and Zongo Glaciers (Favier, Wagnon, Chazarin, et al., 2004;Favier, Wagnon, & Ribstein, 2004;Sicart et al., 2011), mirroring our findings on the Peruvian glaciers (Section 5.3).
These similarities reflect the climate of these glaciers, which exist either in "cold and dry" or "intermediate" conditions, where precipitation rates in both classes are less than ∼1,100 mm year −1 (Sagredo & Lowell, 2012), with this also following the classification of the Andes into dry regions north of 35°S and wet regions further south (Lliboutry, 1998;Schaefer et al., 2020). These glaciers are therefore restricted to high elevations where temperatures are low. Further south, the air temperatures on the Mocho, Exploradores, and Tyndall glaciers in Chile are much warmer, with the sensible heat flux more important, net longwave radiation less negative, and the latent heat flux small or positive (Figure 9 and Schaefer et al., 2020). These differences result in higher melt rates, although of course this is countered by higher precipitation rates at these southern sites (Sagredo & Lowell, 2012).
Despite the similarities of the mean conditions on the Peruvian glaciers to those in Bolivia and Ecuador, there are important seasonal differences. For instance, on Zongo Glacier in Bolivia, the main accumulation period (January to March) is shorter than in Peru (Sicart et al., 2011) and on Antizana Glacier 15 in Ecuador the seasonal variation in the meteorological conditions is very small, with precipitation throughout the year (Favier, Wagnon, Chazarin, et al., 2004).

Conclusion
This work provides the first comprehensive understanding of the mass and energy fluxes of glaciers across Peru.
Our key findings are as follows: Figure 9. Comparing South American glaciers, with sites ordered by latitude (north-south) and data from on-glacier meteorological stations. Panel (a) station elevation (*) in context of the glacier maximum and minimum elevation (error bars), and the mean on-glacier air temperature (red points). Panel (b) radiation and melt characteristics. S* is net shortwave radiation, L* is net longwave radiation, H is the sensible heat flux, λE is the latent heat flux, and 'Melt' is ice and snow melt. Data for sites marked * were averaged over a few months or less, the data for the other sites were averaged over at least a year. Data from Favier, Wagnon, and Ribstein (2004) Barcaza et al. (2017). See Table S12 in Supporting Information S1 for station details. Panel (c) location of the on-glacier meteorological stations, shown on top of all South American glaciers classified following Sagredo and Lowell (2012). The data on the South American glaciers were provided by Esteban Sagredo.
1. The Rio Santa glaciers reach lower minimum elevations and have a greater elevation range than the glaciers in the Vilcanota catchment. This is primarily due to higher precipitation in the wet season in the Rio Santa catchment. There is also a clear precipitation dependence on elevation in the Rio Santa catchment, which is lacking in the Vilcanota catchment; here the northwest prevailing precipitation trajectories lead to lower precipitation over the glaciers to the south. 2. The net shortwave radiation is the greatest contributor to the energy available for melt, with the sensible heat flux being much smaller in comparison due to low air temperatures. The net shortwave radiation is higher in the dry season due to a lower albedo, but cloudless skies and low humidity mean it is counter-balanced by a strongly negative net longwave radiation and latent heat flux. The energy available for melt is therefore highest in the wet season, except at Shallap Glacier where strong La Niña conditions resulting in below average air temperatures meant the snowpack remained continuous throughout the wet season. 3. At three of the sites (Artesonraju Glacier, Cuchillacocha Glacier, and Quisoquipina Glacier) both ice and snow melt occurred throughout the wet season. The wet season snowpack remained thin and was continually formed and melted over periods of days to weeks. Sublimation was more important in the dry season at all sites and became a greater proportion of total ablation with increasing elevation, reaching 81% of total ablation at the summit of Quelccaya Ice Cap. 4. The key drivers of the mass and energy balance of Peruvian glaciers are the dynamics and phase of wet season precipitation since they control the surface albedo and hence net shortwave radiation. Air temperature is still important in determining the phase of precipitation and the elevation at which sublimation becomes more important than melting for ablation. 5. Sensitivity experiments revealed that increased air temperature led to enhanced ablation, but the mechanisms differed by elevation. At the lower elevation Rio Santa sites increased ablation was driven by a change in precipitation phase from snow to rain, reducing albedo and increasing net shortwave radiation, alongside an increase in the sensible heat flux. However, at higher elevations (specifically at the top of Quelccaya Ice Cap), the increase in ablation with temperature was nonlinear and driven by a switch from sublimation to melt dominated ablation.

Data Availability Statement
The Artesonraju and Shallap meteorological data (both on and off glacier) were collected by the University of Innsbruck, with raw data available here: https://acinn-data.uibk.ac.at/. Meteorological data used to bias-correct the WRF outputs from SENHAMI are available here: https://www.senamhi.gob.pe/servicios/?p=estaciones. All of the data associated with this paper are available from the NERC EDS Environmental Data Centre, with the averaged WRF data at: https://doi.org/10.5285/7dbb2d72-7032-4cfa-bc9b-aa02bebe8df5, the meteorological input data and T&C code available at: https://doi.org/10.5285/b69b8849-6897-47eb-a820-f488f8bca437 and the analyzed outputs and associated code available at: https://doi.org/10.5285/5f6661e4-1d34-4b01-8f3a-9fc86c546f73. used to analyze the WRF data. We are grateful to Evan Miles for supporting the derivation of the snowline elevations for both catchments. Esteban Sagredo kindly provided the data on South American glaciers used within Figure 9 and we thank Álvaro Ayala for providing meteorological data for the Guanaco, Yeso, and Juncal Norte Glaciers given in Table S12 in Supporting Information S1 and Figure 9.