Measurements and modeling of snow albedo at Alerce Glacier, Argentina: effects of volcanic ash, snow grain size and cloudiness

The impact of volcanic ash on seasonal snow and glacier mass balance has been much less studied than that of carbonaceous particles and mineral dust. We present here the first field measurements on Argentinian Andes, combined with snow albedo and glacier mass balance modeling. Measured impurities content (1.1 mgkg−1 to 30000 mgkg−1) varied abruptly in snow pits and snow/firn cores, due to high surface enrichment during the ablation season and possibly local/regional wind driven resuspension and redeposition of dust and volcanic ash. In addition, we observed a high spatial heterogeneity, due to 5 glacier topography and prevailing wind direction. Microscopic characterization showed that the major component was ash from recent Calbuco (2015) and Cordón Caulle (2011) volcanic eruptions, with a minor presence of mineral dust and black carbon. We also found a wide range of measured snow albedo (0.26 to 0.81), which reflected mainly the impurities content and the snow/firn grain size (due to aging). We updated the SNICAR snow albedo model to account for the effect of cloudiness on incident radiation spectra, improving the match of modeled and measured values. We also ran sensitivity studies considering 10 the uncertainty of the main measured parameters (impurities content and composition, snow grain size, layer thickness, etc) to identify the field measurements that should be improved to facilitate the validation of the snow albedo model. Finally, we studied the impact of these albedo reductions in Alerce glacier using a spatially distributed surface mass-balance model. We found a large impact of albedo changes in glacier mass balance, and we estimated that the effect of observed ash concentrations can be as high as a 1.25 meter water equivalent decrease in the annual surface mass balance (due to a 34 % of increase in the 15 melt during the ablation season).


Introduction
Since glaciers are highly sensitive to climate fluctuations, their unprecedented retreating rates ::: rates ::: of ::::: retreat : observed during the last decades represent one of the most unambiguous signals of climate change (Zemp et al., 2015;IPCC, 2019). Along the Wet Andes (below 35°S latitude), both precipitation decrease and air surface temperature increase have been pointed out as the drivers of the shrinkage of glaciers in the last decades (Dussaillant et al., 2019). Although some processes, like sublimation at the high and cold Dry Andes (37°S to 20°S) or the calving at the outlet glaciers of the Patagonian Ice fields (south of 45°S), could contribute, or be even more critical than melt for the shrinkage of glaciers in some particular cases, ablation is mainly ruled by melt. Along the Southern Andes, melt is driven by shortwave radiation and sensible turbulent flux (Schaefer et al., 2020). The shortwave :::::::: Shortwave : radiation absorption increases significantly during summer, due to the exposure of low 25 albedo areas in their ablation zones, which causes strong, positive feedback that enhances surface melt significantly and shapes the spatial ablation pattern (Brock et al., 2000). Furthermore, deposition of light-absorbing impurities ::::::: particles (LAP: mineral dust, volcanic ash, and black carbon) have a fundamental impact on the melting of glacier and snow-covered areas Bond et al., 2013;Molina et al., 2015). LAP decrease snow albedo, increasing solar radiation absorption and thus producing a direct effect on snow melting. But, in addition, the snowpack temperature increase due to the direct effect 30 accelerates the growth of snow grains, which produces a further albedo decrease (and thus an additional, indirect impact on snow melting) (Bond et al., 2013;Flanner et al., 2007). While LAP control the snow albedo mainly in the visible wavelengths (since ice is relatively transparent in the visible band), the snow grain size affects the albedo in the near-infrared (e.g., Hadley and Kirchstetter, 2012;Pirazzini et al., 2015;He and Flanner, 2020). Recently it has been highlighted that the growth of glacier algae could also decrease the albedo (Williamson et al., 2019). 35 Atmospheric particulate matter (PM) is diverse in size, chemical composition and optical properties; while most PM reflect :::::: reflects a large fraction of the incoming radiation and thus have a cooling effect on the atmosphere, other particles absorb a relevant :::::::: significant : fraction of the visible radiation (depending on the ratio of their absorption and scattering coefficients) and have a heating effect (Bond et al., 2013). In snow, the term LAP is used to refer to black carbon (BC), mineral dust, volcanic ash and all other particles that totally or partially absorb incident light and hence increase the snow energy absorption. Different 40 snow albedo models have been developed to include the direct effect of BC and other LAP as well as several positive feedbacks (Flanner et al., 2007;Koch et al., 2009;Krinner et al., 2006), such as the increase in surface concentration of impurities due to enhanced snow melting, or the albedo reduction due to the growth of snow grains by accelerated snow aging (Bond et al., 2013). More recently, models have included the effects of non-spherical snow grains He et al., 2017), and external/internal mixing of impurities with snow grains (He et al., 2018). Although some snow albedo models have been 45 successfully validated for laboratory conditions (Brandt et al., 2011;Hadley and Kirchstetter, 2012), the prediction of snow spectral albedo in environmental conditions is still challenging. When the snow has been undergoing heavy metamorphosis processes, a single snow grain size distribution is not enough to reproduce the snow spectral albedo due to the fact that the largest particles and the thinnest protrusions of the irregular crystals have contributions to the snow reflectance that depend on the wavelength Pirazzini et al., 2015). Notably, it has been shown that taking into account the 50 amount of LAP in the snow reduces the difference between simulated and measured albedo, especially in the visible range (Zhang et al., 2018).
Different studies have considered the effect of LAP in snow and ice albedo and its impact on glaciers mass balance or seasonal snow cover, and estimated its radiative forcing (Qian et al., 2015;Skiles et al., 2018). Some studies used point measurements of LAP content (ice cores) together with a snow albedo model to estimate potential melting, using a radiative transfer model to calculate the additional absorbed energy by BC and mineral dust Zhang et al., 2018) or perturbing a glacier mass balance model to include BC forcing (Painter et al., 2013). "Online" coupling of snow albedo models in global or regional atmospheric chemistry models (where both models are run simultaneously allowing two-way feedback) has been applied to study snow and glaciers interaction with the climate around the globe (Hansen et al., 2005;Flanner, 2013;Ménégoz et al., 2014). Although these global or regional atmospheric studies are beneficial to identify LAP 60 sources and dispersion patterns and to compare snow-atmosphere feedback in different regions, the spatial resolution can be inadequate to obtain accurate results in mountain regions (Ménégoz et al., 2014;Qian et al., 2015).
In recent years there has been an increase of measurement and modeling of albedo along the Southern Andes (Rowe et al., 2019). A three-year study  showed that glaciers closer to population centers in the Cordillera Blanca, Peru, have higher surface content of equivalent black carbon (EBC, BC plus other LAP, especially dust in this case), up to 70 ng g −1 EBC, as compared with remote glaciers (with surface content as low as 2.0 ng g −1 EBC). A one-week study successfully 70 connected the decreases in snow broadband albedo with heavy traffic days in the nearby road that connects Argentina and Chile (Cereceda-Balic et al., 2018). A more recent study along the Southern Andes of Chile found a mean albedo reduction due to light-absorbing impurities :::::: particles : in the snow, with its corresponding mean radiative forcing increase (Rowe et al., 2019).
They conclude that in the north (dusty, vegetation-sparse Atacama Desert), BC plays a smaller role than non-BC, whereas near Santiago and in the south (vegetation-rich), the BC contribution is higher. For example, the albedo reduction due to BC 75 alone in the north was estimated to be only about 43 % of that for all light-absorbing impurities ::::::: particles : (assuming spherical 100 µm radii snow grains). By comparison, these albedo reductions are 53 % and 82 % near Santiago and in southern Chile, where a greater share of light absorption is due to BC. In the Southern Andes of Argentina, the only available information on snow albedo is due to remote sensing (Malmros et al., 2018), and up to now, the impact of volcanic ash and other LAP on Argentinian glaciers mass balance has not been evaluated either.

80
Here we present the results from two field campaigns developed in the Alerce glacier during April 2016 and April 2017 to assess the bounds of PM deposition impact in the Alerce glacier mass balance. We show in situ albedo measurements and PM concentration values measured on surface and sub-surface snow and firn :: firn : samples in accumulation and ablation zones of the glacier. Albedo in situ measurements are compared with results from SNICAR snow albedo model (Flanner et al., 2007;He et al., 2018), using measured snow properties and LAP content as input data. We present here an improvement of SNICAR's 85 incident radiation spectra (presented as SNICARv2.1), to take into account changes in direct and diffuse solar radiation for partly cloudy skies. We study the effect of volcanic events that occurred in recent years (Puyehue-Cordón Caulle in 2011 andCalbuco in 2015). Finally, the influence of LAP on snow/ice albedo on the annual surface mass balance of Alerce glacier is assessed using an enhanced temperature index melt model (Oerlemans, 2001). This study is not only the first field study of the impact of LAP in Argentinian glaciers, but also one of the few studies of the long-term impact of volcanic ash on snow albedo.

2 Site Description and Experimental Methods
Alerce is a small (2.2 km 2 ), debris-free, mountain glacier located at Monte Tronador (41.15°S, 71.88°W), in the Northern Patagonian Andes. The climate on this region is primarily modulated by the weather disturbance embedded in the mid-latitude westerlies (Garreaud et al., 2009). Weather disturbances and prevailing winds coming from the Pacific Ocean are more frequent and stronger in winter. However, associated frontal precipitation system move over the Patagonian Andes all year round. In this 95 region, the hydrological year begins on April 1 st with the accumulation season. The accumulation season lasts until October 31 st , which marks the beginning of the ablation season.
Alerce glacier has an elevation range between 1650 m to 2400 m a.s.l. (above sea level), a gentle slope (mean of 10°), and is exposed to the southeast. Since 2013 it has been the focus of a glacier mass balance monitoring program by the IANIGLA (Instituto Argentino de Nivología, Glaciología y Ciencias Ambientales; Ruiz et al., 2015Ruiz et al., , 2017. Seasonal mass balance has 100 been studied every year using the traditional glaciological method of stakes and snow pits. An enhanced temperature index mass balance model has been developed (Huss et al., 2008;Huss, 2010) to study the surface mass balance of the glacier. This model is used here to analyze the influence of LAP, through glacier albedo changes, over the mass balance of Alerce glacier.
In recent years Monte Tronador glaciers have been reached by volcanic ash derived from two volcanic events: (i) Puyehue-Cordón Caulle volcanic complex, which had a long eruption between June 2011 and January 2012 and (ii) Volcán Calbuco, 105 which commenced on April 23 rd 2015.

Fieldwork
In April 2016 and April 2017, besides mass balance measuring, we took snow and firn ::: firn samples and we measured surface albedo at Alerce glacier. Figure 1 shows the sampling sites at Alerce glacier. We sampled accumulation and ablation zones and looked for similar sampling sites in both campaigns. Otto Meiling mountain hut served both as a base camp for field trips and 110 as a field laboratory for initial processing of the snow samples. April 2016 served as a exploratory campaign. Albedo measurements were improved for the 2017 campaign. We lowered instrumental uncertainty and used an improved mounting stand for the pyranometer, which allowed us to evaluate the variability/uncertainty of albedo measurements by repeatedly measuring in the same site. We also improved the measurement of snow grain size distribution. More details are given below. However, the second campaign duration and number of sampling sites were shortened due to poor weather conditions. Nevertheless, relevant 115 results of PM concentration and albedo measurements are presented for the first time for Monte Tronador glaciers.

Snow samples. Filters treatment
Before collecting snow/firn ::: firn/ice samples, we performed an in situ stratigraphy at each site to identify and date layers. Many of the sampling sites corresponded to the accumulation zone of Alerce glacier or accumulation pockets in the ablation zone. In those sites, we dated seasonal layers of snow/firn :: firn. The main elements to attribute layers were PM content and hardness 120 of the layers. Figure 2 shows the results of the stratigraphy and PM gravimetry, which are described in detail in Sect. 3.1. In sampling sites located on the ablation zone, we distinguish glacier ice from recent snow covering the glacier ice.
Most of the samples were taken from snow/firn :: firn : pits. In the 2016 campaign we also used a snow/firn ::: firn : hand auger to sample a 2.5 m snow/firn ::: firn : core (site Acc2-2016, Fig. 2). Samples were melted and filtered in the base camp, and filters were taken to the laboratory for gravimetric determination of PM content and further analysis. Further details are given in Sect. S1.1 125 in the Supplement.
PM in the filters was described and photographed using a Leica S8APO stereo microscope equipped with a DFC 295 camera.
Some samples were also studied by Scanning Electron Microscopy (FEI Quanta 200, equipped with an Edax accessory for energy dispersive X-ray analysis).

Albedo: measurements and corrections
We performed in situ albedo measurements in some of the snow sampling sites in both field campaigns. Raw albedo values were corrected to account for the diffuse or reflected light blocked by the operator or the mounting stand and, for upwelling radiation, the effect of shadows of the sensor and the stand in the snow surface (Wright et al., 2014;Carmagnola et al., 2013). Further details are given in Sect. S1.2 in the Supplement.

Pyranometer mounting stands and cloudiness effect 140
In the 2016 campaign, we used a fixed mounting stand with three stainless steel legs ( Fig. 3 (a)). It was designed to provide a stable irradiation measurement, with a precise tilt angle (parallel to the snow surface), and to minimize the blocking of incident light. When measuring clear-sky downwelling radiation, this stand does not block light at all (operators stand 4 m away from the sensor, blocking less than 0.1 % of incoming diffuse radiation). For clear-sky upwelling radiation, the percentage of blocked light is below 0.8 %, and shadows from the equipment represent another 0.4 %. Hence, total correction for upwelling radiation 145 sum up around 1.2 %, affecting around 1 % measured albedo. For cloudy or overcast conditions, due to the sharp changes in cloud cover, incoming radiation varies more quickly than the time needed for assembling/disassembling the pyranometer stand. To proceed faster under these conditions, the measurements were made differently: the sensor was held by two operators, each 0.45 m away from the sensor, without using the stand legs. Under these conditions 12 % of diffuse downwelling and 9 % of upwelling radiation is blocked by the operators, resulting in an albedo correction of 3.5 %, significantly higher than those 150 obtained for clear-sky conditions.
To overcome the difficulties due to cloudiness, for the 2017 campaign a new mounting stand was designed. The new lighter design has only one arm and one leg and is carried by one operator, located 1.25 m away from the sensor, and leveled manually with the help of a bubble level ( Fig. 3 (b)). This design allows fast and easy alternate downwelling/upwelling radiation measurements, making it possible to assess the variability of albedo under the same sky conditions. For downwelling radiation 155 the operator blocks around 1.1 % of diffuse light. For upwelling radiation, the operator blocks around 1.9 % of light, which, together with shadows of the equipment, brings corrections to a maximum of 2.4 %. Overall albedo corrections vary between 0.8 % and 2.0 %. Additional details on the mounting stands are given in Fig. S1 in the Supplement.

Diffuse and direct radiation fraction
For albedo calculation, the upwelling radiation measurement is used directly from measurements. But for downwelling radia-160 tion, direct and diffuse fraction must be distinguished (see Eq. (S1) in the Supplement).

6
The calculation of the diffuse fraction of downwelling radiation requires to add another measurement with the pyranometer (total downwelling, diffuse downwelling, and total upwelling radiation), and the operation of the accessory to block direct radiation. Fast changes in cloudiness during measurements made it very difficult to assure that all three measurements were performed under the same sky conditions. Therefore, we decided to prioritize that measurements required for albedo calcu-165 lation (total downwelling and total upwelling radiation) were performed under the same conditions, and thus we dropped the diffuse downwelling radiation measurement. Hence, the diffuse to global radiation ratio I dif f ↓ /I glob↓ (needed for albedo measurements corrections and comparison with modeled albedo) had to be estimated differently. We used in situ observations of cloudiness (or pictures of the sky taken before and after albedo measurements) together with the relations found by Kasten and Czeplak (1980)[eq. 4] to estimate the diffuse radiation ratios, which are presented in Table 1. In the 2016 campaign, snow was placed in a crystal grid (with three different scales: 2 mm, 1.2 mm, and 0.6 mm) and average size was determined with a magnifying lens. In the 2017 campaign, a similar in-house developed grid was used (with two scales: 1 mm and 0.5 mm) in combination with a macro lens and a mobile phone digital camera. High-resolution pictures ( Fig. S3 in the Supplement) were analyzed later with ImageJ software (Schneider et al., 2012). Snow grains were manually 175 fitted with ellipses; the metric choice was the average of the minor and major axes of the ellipse. The new equipment and methodology introduced in the 2017 campaign allows a more detailed description of the snow samples and a more precise average radius value.

Albedo: modeling
To analyze the different factors affecting measured albedo at each sampling site, we modeled albedo for the same conditions 180 using SNICAR (Flanner et al., 2007;He et al., 2017He et al., , 2018. Snow density and layer thicknesses were taken as parameters from in situ stratigraphies. Average snow grain size and shape were obtained from in situ measurements. LAP content was obtained from filters gravimetry. Based on in situ observations and the analysis of microscopy images (Sect. 3.2), we assigned all recollected PM mass to volcanic ash (in a similar way as previously done in sites where mineral dust represents most of LAP; Krinner et al., 2006;Painter et al., 2012;Wittmann et al., 2017). Albedo of the underlying layers was calculated explicitly 185 within the same model, using the properties of those layers.
SNICARv2 (He et al., 2017(He et al., , 2018 supported only four incident solar spectra: two clear-sky direct solar spectra (one for Summit, Greenland, and one for mid-latitude), and two overcast diffuse spectra (for the same locations). These spectra are used to calculate direct radiation albedo and diffuse radiation albedo, respectively. These are good approximations for clearsky albedo (where most for the incident radiation is direct, clear-sky solar radiation) or for overcast sky albedo (where most of 190 the radiation is diffuse). In this updated version of SNICAR (referred as SNICARv2.1 throughout the article) we provided an alternative for these spectra for cases when latitude, longitude or altitude differ significantly from those of the provided spectra, or where the sky is partly cloudy.
7 First, we calculated the clear-sky spectra for the site location and time using SMARTS model (Gueymard, 2001). Then, we calculated the direct and diffuse spectra for overcast or partly cloudy sky following Gueymard (1986Gueymard ( , 1987 and Ernst et al. 195 (2016): I dir , I dif f , and I glob are clear-sky direct, diffuse and global solar irradiance (where I glob = I dir + I dif f ), as calculated from SMARTS model. F dir,S (λ) and F dif f,S (λ) are the spectral distributions of clear sky direct and diffuse solar irradiance, also 200 from SMARTS model. F dir,norm (λ) and F dif f,norm (λ) are the normalized spectral distributions of direct and diffuse solar irradiance thus calculated for our sites. The cloud opacity factor N pt is calculated following Ernst et al. (2016): where ρ and ρ S are the diffuse to global irradiance ratio for the site and from SMARTS model, respectively.
The clear-sky direct radiation spectra available in SNICARv2 matches reasonably well SMARTS clear sky direct radiation 205 spectra. On the other hand, SMARTS clear sky diffuse radiation spectra is very different from diffuse radiation spectra available in SNICARv2 (Fig. 4). The spectral distribution obtained for 95 % cloud fraction for SNICARv2.1 closely matches the diffuse radiation spectra available in SNICARv2, which confirms that the latter was prepared to represent an overcast sky condition.
On the other hand, the spectral distribution obtained for a 50 % cloud fraction differs significantly from both spectra available in SNICARv2, showing a larger contribution from clear sky diffuse radiation (Fig. 4).

210
Hence, we expect to find a larger impact of our improved incident sun spectra for intermediate cloud cover fractions. For clear sky conditions, direct radiation spectra were already well represented. Even though diffuse radiation spectra were not accounted for, this fact has little impact on the calculated albedo, due to the low diffuse radiation fraction for clear sky conditions.
Conversely, for overcast conditions, diffuse radiation spectra were already well represented and neglecting direct radiation fraction has a low impact on albedo calculations.

215
Using different incident radiation spectral distributions, we obtained the pure direct and diffuse albedo with SNICARv2 and SNICARv2.1 (α dir and α dif f ). For SNICARv2.1 we also calculated the weighted average albedo, which should be compared to the net measured albedo:

Alerce glacier surface mass balance model 220
To analyze the role of albedo decrease over the surface mass balance of Alerce glacier, we used a spatially distributed surface mass-balance model (spatial resolution 20 m) driven by daily temperature, precipitation, and potential direct solar radiation (Huss et al., 2008). The model was calibrated by surface mass balance measurements performed on a seasonal to annual basis through the year 2016 over Alerce glacier.
Here we summarize the most relevant model components. Snow accumulation C (x,y,t) for all grid cells (x, y) and all time 225 steps (t) was calculated based on precipitation P (t) occurring below a threshold air temperature of 1.5 • C (Hock, 1999).
Accumulation distribution D s(x,y) was inferred based on a spatial distribution pattern derived from winter snow measurements and topographic parameters (slope, curvature) to account for small-scale snow redistribution (Huss et al., 2008;Sold et al., 2016).
230 P (t) was the daily precipitation at Tepual weather station (90 m altitude, ID = 857990; https://www.ncei.noaa.gov/access/search/datasearch/global-summary-of-the-day). The factor C pre allows adjusting precipitation measured at the weather station to the conditions on the glacier.
Snow and ice melt were calculated based on a simplified energy-balance formulation proposed by Oerlemans (2001), where the energy available for melt Ψ d(x,y,t) was defined as follows: where I (x,y,t) is the potential direct solar radiation in W m −2 , τ is the atmospheric transmission to solar irradiance, T (x,y,t) the air temperature and c 0 and c 1 represent parameters. T (t) was taken from the air surface temperature at Bariloche airport weather station (846 m altitude, ID = 877650; https://www.ncei.noaa.gov/access/search/data-search/global-summary-of-theday). Potential direct solar radiation for all grid cells and days was calculated following Hock (1999). The local surface albedo 240 α (x,y,t) was taken to be constant for bare-ice surfaces (α ice = 0.34), using most commonly applied literature value (Oerlemans and Knap, 1998;Cuffey and Paterson, 2010), for snow surfaces, α snow was calculated based on the snow aging function proposed by Oerlemans and Knap (1998) with a maximum snow albedo (α max ) of 0.8 and a minimum snow albedo (α min adjusted during the calibration procedure. The model was calibrated in two steps using surface mass balance measurements of year 2016 in Alerce glacier ( Fig. S4 in 245 the Supplement). First, the model was run over the winter period with an initial set of constants (c 0 and c 1 ) and a guess for the precipitation correction factor C pre . As melt is of minor importance in winter, this run was used to calibrate C pre , that scales D s for every snow fall event. After a good agreement of measured and calculated winter accumulation was obtained, the model was run over the entire year and the remaining constants were calibrated so that the root-mean-square error between modelled and observed point annual balances were minimized and the average misfit was close to zero ( Fig. S5 and S6 in the 250 Supplement). A random set of snow accumulation and ablation stakes measurements performed through the year and not used to calibrate the model were left apart to validate the results of the surface mass balance model.
We studied surface mass balance changes for different values of α min (Table 2), which are indicative of the sensitivity of glacier mass balance to a change in albedo that might occur in response to the darkening of the glacier surface. The fresh snow at the top of Abl2-2016 shows slightly higher content of PM than fresh snow sampled on the accumulation zone ((21.9 ± 0.6) mg kg −1 ). In the case of fresh snow on site Abl5-2017 (with a higher PM content of (1410 ± 30) mg kg −1 ) 295 we could not discard, due to its thin thickness, some contamination with PM from the glacier ice. Glacier ice was highly heterogeneous (relatively pure ice mixed with debris and cryoconite holes), in consequence a substantial variability of PM content over the ice surface was retrieved ((200 ± 20) mg kg −1 to ((4300 ± 900) mg kg −1 ).

Results and Discussion
Figure 5 combines data from both field campaigns and groups PM concentrations according to the attributed date of the layers, but excludes glacier ice samples, which cannot be assigned to an specific year/season. It must be noted that PM content 300 varies over several orders of magnitude (1.3 mg kg −1 to 21.9 mg kg −1 on fresh snow, to up to (30 000 ± 5000) mg kg −1 in end-of-summer layers of the ablation zones). As discussed in Sect. 3.3, this is one of the main causes of the albedo values variation.
The alternation of thin and high PM concentration with thick and low PM concentration is partially due to seasonality, as explained above. But in addition to seasonality, there is a large spatial heterogeneity, especially during spring/summer 305 (in winter, abundant fresh snow covers the glacier and gives a more homogeneous PM content and albedo distribution, as observed in other glaciers, Brock et al., 2000). The spatial variation is not only between the ablation and accumulation zones of the glacier. The interaction between glacier topography and prevailing winds produce accumulation pockets and windswept ridges, which have contrasted snow accumulation values. These areas of higher and lower accumulation lead to a wide range of spectral albedos. The detailed variations in PM concentrations, and therefore in the albedo, need to be accounted for in a 310 detailed mass balance of the glacier (see Sect. 3.4).
Field observations on Monte Tronador in 2013 and 2014 confirmed the presence of volcanic ash in the atmosphere, derived from resuspension of volcanic ash. The magnitude of resuspension events in Andean Patagonia, a region with strong, persistent westerlies and a dry season with low relative humidity, is well known. These aeolian remobilization events may produce huge ash clouds that may be even confused with true volcanic plumes, they can remobilize ash tenths of kilometers away (Toyos et al.,315 2017). In particular, deposits of volcanic ash that are covered by snow during the winter in the high mountain usually become exposed to remobilization during the summer, traveling through the atmosphere and redepositing over different surfaces due to decrease of wind competence or by adherence of particles on humid surfaces, even at considerably high altitudes.
The 2011 Puyehue-Cordón Caulle eruption produced several ashfall events during the second semester of 2011, by January 2012 explosive activity had declined. As a consequence, thick deposits of tephra with different grain size covered an extended 320 area in Argentina (see Fig. 2, Alloway et al., 2015). Calbuco eruption (April 2015) was active during a shorter period, but due to its location and predominant wind direction also affected Monte Tronador Reckziegel et al., 2016).
Direct ash deposition and resuspension events can affect the glacier surface in different ways. Continuous, thick layers of ash (few millimeters to few centimeters) have :::: been : shown to behave as an isolating layer when deposited over snow, in a similar way as in debris-covered ::::: similar ::: to ::: the ::::: effect :: of ::::: debris ::: in ::::::: covered glaciers (Brock et al., 2007), which reduces ablation. But 325 on the other hand, a thinner or disperse deposit may have the opposite effect, lowering the surface albedo of the glacier and increasing its melting. The effect of ash (or other PM, for instance from biomass burning events) deposition during autumn or winter can extend a few days until the next snow event, which covers the dark surface with the highly reflecting surface of fresh snow (see Fig. 7, Córdoba et al., 2015). But during spring and summer, warmer temperatures and fewer snow events result in an increase of ablation processes over accumulation. Snow melting can flush some of the smaller, hydrophilic PM, but larger 330 particles (or less water-soluble small particles) are concentrated on the glacier surface (Conway et al., 1996;Xu et al., 2012;Doherty et al., 2013;Li et al., 2017;Skiles and Painter, 2017), producing up to two orders of magnitude of surface enrichment of PM content (Doherty et al., 2016). Resuspension and surface enrichment explain the observed alternating thin, high PM concentration layers and thick, low PM concentration layers. They also impact the spatial variability of albedo on the glacier surface during summer ( Fig. 1) (Brock et al., 2000).

PM characterization
Three main types of particles were identified in samples collected in the field: mineral dust, volcanic ash and crystals derived from ashfall events, and carbonaceous particles.
Based on glass morphology, SEM images, and energy dispersive spectroscopy (EDS) microanalysis performed on selected fragments, we were able to identify the presence of volcanic glass derived from Cordón Caulle 2011 (CC) and Calbuco 2015

340
(Cal) eruptions. Isopach maps for both eruptions (Alloway et al., 2015;Villarosa et al., 2016;Reckziegel et al., 2016) show that Monte Tronador was reached by different plumes from ashfall events marginally, further confirming that most of the volcanic ash identified in the filters derive from these two recent eruptions. Though both eruptions deposited pumiceous ash east of the Andes in Patagonia, they can be distinguished by petrographic and morphological characteristics of the glass fragments ( Fig. 6). CC glass is very fine-grained colourless glass (rhyolitic in composition) while Cal pumice is light to pale brown, 345 clear glass (dacitic to andesitic in composition). SEM images show the presence of irregular glass fragments, with evidence of bubble coalescence, flat or slightly curved platy glass shards that are most probably pieces of broken thin vesicle walls and triangular (in cross section) to Y-shaped particles, which are vesicle walls from the junction of three adjacent vesicles (Fig. 7).
EDS analysis of individual fragments of glass from these samples were performed, and were compared with the composition of volcanic glass from samples collected in nearby locations during direct ashfall events. Results confirm the presence of glass 350 shards from 2015 Cal and 2011 CC eruptions (Fig. 8). One of the samples described under microscope, corresponds to :: In a sub-surface sample from site Abl3-2017, which was interpreted as winter snow from 2014, previous to 2015 Cal eruption, and ::: we ::::: found ::: that : approximately 75 % of the observed particles correspond to fine-grained colourless pumiceous ash. EDS of individual fragments confirmed that ash on that sample correspond to CC eruption, as expected.
Another evidence of the presence of volcanic material within the PM collected in the study area are crystals from pyroclastic 355 origin. They are clearly identified as they are partially surrounded by or associated with patches of glass and they are irregular in shape. Crystals that are not directly derived from CC 2011 or Cal 2015 are more or less rounded due to erosion and transport and they exhibit a dull lustre, and they are identified as mineral dust.
Another identified PM component is charcoal, present as black, elongated, brittle fragments. In addition, some of the samples showed evidence of the presence of BC particles, identified by their characteristic shape (carbon spherules of 100 nm to 200 nm 360 in aggregates of different morphology). Carbon content by EDS could not be used to confirm the identity of BC particles due to the usage of carbon tape to fix the particles for SEM imaging.
The predominance of volcanic glass in the collected PM indicates the need to take into account the effect of volcanic ash in the albedo of seasonal snow and glaciers of the region, which can be frequently affected by volcanic eruptions. It must be emphasized that ash from CC and Cal eruptions was observed in most of the samples, not only in layers dated immediately after 365 the eruptions, but also many years after direct deposition. Field stratigraphy together with these microscopy results suggest that we can study the effect of LAP on snow albedo considering that all PM content can be attributed to LAP (and more specifically, to volcanic ash). Further chemical studies will be performed on the PM samples to refine the representation of LAP in the snow albedo model, since optical properties can be very different for BC, mineral dust, volcanic ash, etc. (the ratio of light absorption to light scattering at different wavelengths depends on particle size, shape, and chemical composition). Table 1 shows measured and modeled albedo values for six sites (two from the first field campaign, 2016, and four from the latter, 2017), together with different measured properties of the snow topmost layer and site.

Albedo: measurements and models
Reported values of measured albedo include shadow corrections, although these corrections were quite small in all cases (below 3.5 % for worst conditions in the 2016 campaign and below 2 % for the 2017 campaign). In some cases (site Acc3-375 2016) the corrections in the measured incoming and reflected radiation are higher (10 to 14 %), but they largely balance out.
For the 2016 campaign, the reported measured albedo is a single measurement (registered after voltage reached a stable value) and is informed together with its instrumental uncertainty. It must be noted that for this campaign the reported uncertainty reached values as high as 15 % for worst conditions (low incident radiation and low albedo, as in Acc3-2016) or around 2 % for best conditions (clear sky, high albedo). For the 2017 campaign the instrumental uncertainty was lowered by improving the 380 accuracy of the digital multimeter used with the pyranometer, achieving uncertainties lower than 3.5 % (worst conditions) or lower than 1.2 % (best conditions).
Results from the 2017 campaign, obtained using the improved mounting stand, shed light on the reproducibility of albedo measurements. For this campaign, we found that repeated albedo measurements in the same site have a standard deviation corresponding around 5 to 10 % of the average values. This range could be partly due to the leveling of the stand, or to inherent 385 variability of the measurement at these sites (especially differences in the solar irradiance for situations with rapid changes in cloudiness).
Regarding snow grain size, it is relevant to notice the range of observed average radius. In fresh snow samples from the accumulation zone (sites Acc5-2017 and Acc6-2017) we found an average snow grain radius of (151 ± 41) µm, whereas in samples of older firn :: firn : in the ablation zone (or sub-surface snow/firn ::: firn in the accumulation zone) we measured values Table 1. usually around (1000 ± 200) µm. Pirazzini et al. (2015) also used 2D photos, but with a different metric. They suggest that SSK (shortest skeleton branch) is a proxy for "half the width of the shortest particle dimension", which they claim is a better approximation of the optically equivalent snow grain radius. Our metric (see Sect. 2.2.3) would probably give higher results than SSK, and hence we might have overestimated the optically equivalent snow grain radius. Nevertheless, as we show below in this section, our grain size measurements seem to be good enough to reproduce the measured albedo for fine and coarse 395 snow in SNICARv2.1 snow albedo model. It must be emphasized here that the developed method for characterizing snow grains for the 2017 campaign allows us to measure the size distribution and to assess the relevance of different grain shapes (when necessary). It has been shown that the shape of the snow grains can significantly affect snow albedo He et al., 2017). Except for fresh snow (snow less than one day old), where it is possible to still distinguish crystal fragments, in both campaigns the observed snow/firn :: firn : grains were rounded. This is related with the temperate climate at Monte Tronador, 400 where snow temperature is above −5 • C and the temperature gradient is low. Also, the presence of meltwater within the snow layers enhance the rate at which grains become rounded, because the grains melt first at their extremities. Finally, the average grain size increases because the smaller grains tend to melt before the larger ones (Flanner and Zender, 2006). Hence, we assumed spherical grains for all modeled albedo calculations. but the fraction of diffuse radiation is very low, and hence its contribution to net albedo is also low. For overcast conditions (Acc3-2016, Abl3-2017 and Abl4-2017), the pure diffuse albedo from both models is also similar, and weighted average albedo from SNICARv2.1 is coincident with the pure diffuse albedo. For both models, the diffuse radiation spectrum for overcast conditions is coincident with global solar radiation spectrum (see Fig. 4), which explains the similar results. It must be noticed that for site Abl4-2017, we observed rapid cloud movements, and we decided to register two sets of albedo measurements, The 415 average albedo of the second set is similar to the modeled weighted average albedo and to the measurement for site Abl3-2017.
We suggest that this coincidence means that the pictures of the sky above the site (taken after the two sets of measurements) and the estimate of cloud cover based on those pictures represent more accurately the sky conditions during the second set of measurements. Finally, partly cloudy skies (sites Acc5-2017 and Acc6-2017) are the main reason for the development of SNICARv2.1. For these cases, pure direct and pure diffuse albedo differ much more than the associated uncertainties, and 420 pure diffuse albedo from SNICARv2.1 also differs from that from SNICARv2. These differences are also evident from the comparison between the diffuse radiation spectra for partly cloudy skies developed for SNICARv2.1 and the diffuse spectra for overcast skies used in SNICARv2 (Fig. 4). For these sites, SNICARv2 cannot give a good approximation. For Acc5-2017 SNICARv2.1 weighted average albedo seems a good approximation of the measured albedo. For Acc6-2017, measured albedo is lower than pure direct and pure diffuse albedo, so both models give higher estimates for this site. As discussed below in this section, the effect of the diffuse radiation fraction does not seem to be the main source of this disagreement.
The updated model reproduces quite well the main features of the measured albedo (with a larger discrepancy for sampling site Acc3-2016). One of the most important parameters affecting albedo is PM content: the measurements with lower albedo values ( α meas < 0.4 ) correspond to sites with the highest PM content (Acc3-2016, Abl3-2017 and Abl4-2017), whereas the remaining sites have much lower PM content (fresh snow) and α meas > 0.6. It must be noted that for high PM con- On the other hand, comparison between sites with low PM content shows that snow grain size has a remarkable effect, as previously reported Hadley and Kirchstetter, 2012). Fresh snow with small grain size presents α meas ≈ 0.8 (sites Acc5-2017 and Acc6-2017), but snow with similar PM content that has aged a few days presents α meas ≈ 0.6 (site Acc2-2016). Spectral albedo measurements (not available in our field campaigns) would allow :: us : to study separately 445 the effect of grain size and LAP content (see for instance measurements of snow specific surface area, SSA, in Carmagnola et al., 2013), to confirm that our grain size measurements are a good estimate of the optically equivalent grain radius.
The last column in Table 1 reports the results of sensitivity studies to evaluate the impact on the calculated albedo of the uncertainty in key input parameters. We define the sensitivities as the modeled albedo changes increasing or decreasing one parameter in the same magnitude of its reported uncertainty (identified in Table 1 with a "+" or a "-" sign, respectively),

450
while keeping all other parameters unchanged. For each site, we studied PM content and grain size impact, together with other parameters that could be relevant at each site. We highlighted (with bold characters) the higher sensitivities for each site.
Concerning grain size uncertainty (the standard deviation of snow grain radii in each sample), it is clear that the impact on albedo is much larger when PM content is low (sites Acc2-2016, Acc5-2017 and Acc6-2017). For low PM content sites, the effect is comparable to experimental uncertainty, and is relevant both for sites with finer and coarser grain size snow. For sites 455 with high content of PM the uncertainty of grain size does not have an appreciable effect. Pirazzini et al. (2015) determined 11 % uncertainty in the grain size measurements from 2D photos (due to the subjectivity of the software operators). Although we did not determine such kind of uncertainty in our measurements, we suggest that the reported standard deviation (between 16 % and 26 % of the average value) is probably larger than the uncertainty of the method. The sensitivity studies showed that the effect on the modeled albedo is lower than 4.5 % for clean snow and lower than 0.8 % for dirty snow. We believe that this explains the fact that we can reproduce the measured albedo using the estimated grain size together with other snow properties (especially LAP content), even though our grain size estimates might not be as accurate as that ::::: those obtained by other methods.
The uncertainty of volcanic ash content does not have a relevant impact for any of the sites, although it is larger for site Abl4-2017. However, as previously mentioned, the presence of BC (not yet quantified in these samples) could have a more relevant 465 impact on albedo. For instance, it could explain the difference between measured and modeled albedo for site Acc6-2017, and the difference with site Acc5-2017.
Regarding the impact of the uncertainty of layer thickness, the results show that several factors determine the relevance of this parameter. The impact is maximum for very thin layers, especially when the underlying layer has a significantly different albedo (site Abl4-2017, 0.1 cm thick), and its minimum for the thicker layers (sites Acc5-2017 or Acc6-2017, 9 cm thick), or 470 for intermediate thicknesses with high PM content (i.e., low penetration of incident light, site Abl3-2017, 0.3 cm thick). The impact of uncertainty of snow density was not studied in detail, but the impact is inverse to that of the thickness of the layer.
Hence, we report only the moderate impact of snow density uncertainty for site Abl4-2017.
The impact of the uncertainty of the diffuse to global irradiance ratio is moderate but appreciable, which emphasizes the relevance of measuring the ratio on the field. Finally, the impact of the uncertainty of the incidence angle is low, and not 475 appreciable for this range of experimental albedo uncertainty.
Another possible reason for disagreement between modeled and measured albedo, especially for aged snow, is surface roughness. Millimeter scale surface roughness due to snow aging has :::: been : shown to reduce albedo, especially in the infrared region, due to multiple reflections in the cavities (Pirazzini et al., 2015). Computer simulations have studied the parameters that determine the magnitude of the effect of sastrugi (centimeter-scale roughness) on albedo (Zhuravleva and Kokhanovsky,480 2011). Quantification of the impact of surface roughness of snow in measured albedo is out of the scope of this work, but it must be remarked that in sites with higher PM content, where there has been longer snow metamorphosis processes (Acc3-2016, Abl3-2017 and Abl4-2017), we observed higher surface roughness.
Literature values of snow albedo mainly depend on the PM content. Two other studies that found snow albedo ranges similar to our measurements are connected with local/regional transport of dust (Painter et al., 2012;Wittmann et al., 2017). Young Recent studies in Chilean Andes measured or modeled small reductions on snow albedo, due to traffic related BC (Cereceda-Balic et al., 2018) or to a combination of urban BC and dust from desert regions (Rowe et al., 2019). Similarly, studies on Mera Glacier, Nepal , and at several sites on the Tibetan Plateau (Zhang et al., 2018) found small albedo reductions due to BC and dust, and almost negligible effects of impurities in Greenland Wright et al., 2014). Table 2 shows the modeled annual and winter surface mass balance, Equilibrium Line Altitude (ELA) and Accumulation Area Ratio (AAR) for different values of old snow albedo (α min ). Figure 9 shows the change in cumulative surface mass balance 495 and ablation and the annual mass balance elevation gradient for the different values of α min . The mass balance sensitivity to albedo change, defined as the change in surface mass balance per 0.1 of α min decrease is around of −0.6 mw e /yr and −0.07 mw e /yr, for annual and winter mass balance, respectively (Table 2). Aged snow albedo has a considerable effect on the surface mass balance of Alerce glacier (Fig. 9 A) increasing the amount of melt during the ablation period, from almost 2.4 m w.e. to more than 4.6 m w.e. when α min is decreased from 0.7 to 0.3 (Fig. 9 B). Although the accumulation of the glacier does 500 not change (the amount of precipitation for the different run test is the same) there is a decrease in the winter (accumulation) mass balance due to the albedo effect over ablation episodes at the begging of the accumulation season (Fig. 9, Table 2). The decrease in the aged snow albedo α min has an impact all over the glacier, decreasing the surface mass balance at all elevation range. Other glaciological parameters related to the surface mass balance of the glacier, like the ELA or AAR also seems to be profoundly impacted with the decrease of albedo, with a total increase of ELA of 250 m and a decrease of AAR of more than 505 50 % when the old snow albedo changes from 0.7 to 0.3. Nevertheless, since both ELA and AAR depend on the hypsometry of the glacier the changes do not increase constantly. To give physical meaning to the albedo values presented in Fig. 9 and  Our α min analysis allows us to estimate the impact of volcanic ash on the surface mass balance of Alerce glacier. In absence 520 of volcanic eruptions, if we assume that other local or regional PM sources (mineral dust, biomass burning, etc.) do not affect significantly fresh snow albedo, it is expected that the summer α min over glacier surface is similar to α min = 0.6 scenario.

Albedo and modeled impact on glacier mass balance
Although we could not sample summer firn ::: firn : layers previous to 2015 Cal eruption to test this hypothesis, this first order assumption would mean, that volcanic ash is responsible for a 1.25 mw e decrease in the annual surface mass balance (or a 36 % increase in summer ablation), if we compared the α min = 0.6 and α min = 0.4 scenarios.

525
Although more sampling of firn :: firn/snow layer and further chemical analysis on the samples are needed to confirm that the decrease of albedo is only due to the effect of volcanic ash, we have shown that PM content (and hence α min ) varies largely over the glacier surface. Taking into account these spatiotemporal changes in albedo for glacier mass balance models is a defying task. Defining a low number of representative regions over the glacier surface is not an easy task, due to the already mentioned high heterogeneity. In addition, it would be difficult to regularly measure PM content (and/or albedo) on those 530 regions, due to the distances and path conditions on the glacier. Regional atmospheric models could be of help in predicting deposition of volcanic ash, mineral dust, BC and other PM. But the spatial scale of those models (≥1 km) is too coarse to capture to reproduce the spatial variation of the albedo over the glacier.
These challenges have been acknowledged in literature, and several approaches have been followed to estimate snow/ice melting. The simplest approaches have used measured or modeled albedo changes together with measured or modeled solar 535 radiation to estimate melting, without taking into account spatial heterogeneity (in surface temperature, PM concentration, etc) Zhang et al., 2018). For Mera glacier, Ginot et al. (2014) calculate that BC and dust are responsible of approximately 26 % of total melting. Zhang et al. (2018) do not report the effect on melt rates but only the impact on seasonal snow cover duration, and hence the results are not easy to compare with ours. Painter et al. (2013) used a glacier mass balance model similar to ours, but introducing temperature anomalies (due to BC radiative forcing) to estimate mass balance changes.

540
They used several approximations to postulate BC concentrations over the glaciers based on limited ice cores. Their results are difficult to compare to ours due to the different approach, they analyze general mass balance trends over two centuries. Flanner et al. (2007) and Ménégoz et al. (2014) used emission inventories and general circulation models to study deposition of BC (and mineral dust, in the latter work) and its radiative forcing. The spatial resolution of their simulations make difficult the comparison with field PM concentration measurements, and hinder the accuracy of quantitative mass balance calculations 545 (Ménégoz et al., 2014;Qian et al., 2015). Young et al. (2014) used modeled ash deposition, SNICAR and a restricted degreeday radiation balance. They found melt rates between 140 % and 320 % higher than for pure snow, although the low spatial resolution of the simulations (≈ 18 km) may affect the precision of the results. Vionnet et al. (2012) used the detailed snow model CROCUS implemented on the soil model SURFEX to study the snowpack on the Grandes Rousses mountain range in the French Alps . They used a high resolution DEM (150 m) together with meteorological forcing from interpolation of 550 SAFRAN atmospheric reanalysis. They main weakness is that at that moment CROCUS did not explicitly treated PM in snow (it was only implicitly included in their parametrization of snow albedo changes with snow aging).
There are also some examples on literature that studied the coupling of meteorological models with glacier or snowpack models. Different authors studied climate feedback effects on Karakoram glaciers (Collier et al., 2013) and in the Svalbard glaciers (Aas et al., 2016), and the snowpack in Antarctica (Vionnet et al., 2012). The authors suggest that the next steps would 555 be to couple a regional atmospheric model with the ability of prognosis of PM deposition (such as Ménégoz et al., 2014) with a high resolution glacier mass balance model (such as ours or CROCUS implementation on SURFEX, Vionnet et al., 2012), and including explicit treatment of PM effect on snow albedo (such as SNICAR or recent CROCUS implementations Tuzet et al., 2017).

560
Our study combines field measurements and modeling to analyze the role of PM over the albedo of Alerce glacier in Monte Tronador. PM content of the samples varied in a wide range, from lowest to highest: fresh snow (1.1 mg kg −1 to 21.9 mg kg −1 ), old winter snow/firn ::: firn (4.9 mg kg −1 to 51 mg kg −1 , except for some samples from ablation zone), and thin, darker layers with contribution of local/regional resuspension of dust/ash (365 mg kg −1 to 410 mg kg −1 ) or with high PM enrichment due to spring and summer ablation (339 mg kg −1 to 9040 mg kg −1 , reaching even 12 250 mg kg −1 to 30 000 mg kg −1 in the ablation 565 zone). Microscopic characterization of PM showed that the major component on snow and firn :: firn : layers after 2014 and also glacier ice surface is volcanic ash, not only from the recent Calbuco eruption (2015), but also from the Cordón Caulle eruption (2011). Minor contributions of mineral dust and black carbon were also detected.
The fact that volcanic ash represents the largest fraction of the collected PM in all studied samples indicates that the effect of volcanic eruptions are expected not only immediately after direct deposition, but also many years later, due to surface 570 enrichment and wind resuspension and redeposition. The spatial and temporal distribution of PM is highly heterogeneous, due both to seasonality and to the combination of glacier topography and prevailing wind direction. These facts need to be accounted for when studying the effect of snow albedo on glacier mass balance. While the albedo parametrization used in the mass balance model partially accounts for the spatial heterogeneity of PM surface concentration (implicitly), we suggest that in the future it would be useful to couple our mass balance model with an atmospheric model which provides prognosis of PM 575 content and a snow albedo model that includes LAP effect explicitly.
The measured snow albedo also varied in a wide range (0.26 to 0.81), similar to that of other glaciers with dust or volcanic ash concentration in the same order of magnitude. We found that for our set-up (where the pyranometer must be inverted sequentially to measure upwelling and downwelling radiation) rapid changes in cloudiness hinder the repeatability of albedo measurements and may degrade the comparison with modeled albedo. Nevertheless, comparison of measured and modeled 580 snow albedo showed a good match, and illustrates the effect of PM content and composition (i.e., BC versus dust or volcanic ash), snow grain size, layer thickness, and cloudiness on snow albedo. To evaluate the latter, we updated the SNICAR snow albedo model to accurately represent the effect of cloudiness on direct and diffuse solar spectra (SNICARv2.1). This update improved considerably the match of measured and modeled albedo for partially cloudy sky conditions. The effect of uncertainties of field measurements was evaluated for different types of samples, suggesting strategies to reduce uncertainty in snow 585 albedo modeling or retrieval of snow properties from measured albedo. We found that snow grain size must be measured more carefully in samples with low volcanic ash content and that the accuracy of layer thickness can be relevant not only for very thin layers (0.1 cm) but also for thicker layers (6 cm) with low ash content. The accuracy of ash content was found to be good enough for reproducing our albedo measurements. However, it was remarked that the presence of small amounts of BC can affect the albedo significantly in samples with low ash content. 590 We showed that surface mass balance is highly sensitive to the parametrization of aged snow albedo. We find a glacierwide albedo change sensitivity of around −0.6 mw e /yr, mostly due to higher ablation during spring and summer. Finally, we suggest that the effect of volcanic ash in Alerce glacier can be as high as a 1.25 mw e decrease in the glacier annual mass balance or a 34 % of increase in the melt during the ablation season, considering a surface volcanic ash content compatible with that measured in sites Acc3-2016, Abl3-2017 and Abl4-2017. Nevertheless, a more accurate calculation of volcanic ash 595 impact would take into account the amount of other regional or local sources of PM present on the glacier in absence of such volcanic eruptions, which cannot be estimated with the results of the field campaigns reported in this article.
To the best of our knowledge, this work is the first study of PM content and snow albedo on Argentinian glaciers. Our results highlight the need of appropriately considering the effect of volcanic eruptions on snow albedo and glacier mass balance even years after the eruption events. We suggest possible future steps to improve prognosis ability and mass balance accuracy, using 600 a combination of measurements and modeling.
Code and data availability. The complete set of field measurements are available from the corresponding author on reasonable request.  The presence of the stand and the observer is taken into account to correct the albedo measurement through the angles θ and φ and Eq. (S1) and (S2)   Notice that for seasonal layers with only two measurements, the box represents those two values (coincident with the definition of standard deviation for N = 2). The plot includes data from both field campaigns, and excludes ablation ice samples, which cannot be assigned to a specific year/season. Fresh snow represent snow fallen a few days before field campaigns of 2016 or 2017.