In situ continuous visible and near-infrared spectroscopy of an alpine snowpack

Snow spectral albedo in the visible/near-infrared range has been continuously measured during a winter season at Col de Porte alpine site (French Alps, 45.30◦N, 5.77◦E, 1325 m a.s.l.). The evolution of such alpine snowpack is complex due to intensive precipitation, rapid melt events and Saharan dust deposition outbreaks. This study highlights that the resulting intricate variations of spectral albedo can be successfully explained by variations of the following snow surface variables : snow 5 specific surface area (SSA), effective light-absorbing impurities content, presence of liquid water and slope. The methodology developed in this study disentangles the effect of these variables on snow spectral albedo. The presence of liquid water at the snow surface results in a spectral shift of the albedo from which melt events can be identified with an occurrence of false detection rate lower than 3.5%. Snow SSA mostly impacts spectral albedo in the near-infrared range. Impurity deposition mostly impacts the albedo in the visible range but this impact is very dependent on snow SSA and surface slope. Our work 10 thus demonstrates that the SSA estimation from spectral albedo is affected by large uncertainties for a tilted snow surface and medium to high impurity contents and that the estimation of impurity content is also affected by large uncertainties, especially for low values below 50 ng g−1 black carbon equivalent. The proposed methodology opens routes for retrieval of SSA, impurity content, melt events and surface slope from spectral albedo. However, an exhaustive accuracy assessment of the snow properties retrieval would require more independent in situ measurements and is beyond the scope of the present study. Such time series 15 of snow spectral albedo nevertheless already provides a new insight into our understanding of the evolution of snow surface properties.


Introduction
Snow is among the most reflective materials on Earth (Dozier et al., 2009) and its albedo exhibits large spectral variations in the solar spectrum (Warren, 1982).Snow albedo is therefore a crucial variable of the Earth energy balance (Imbrie and Imbrie, 1980) through the surface energy budget of snowcovered surfaces.Snow spectral albedo varies as a function of many factors such as (i) the spectral and angular characteristics of the solar incident radiation and (ii) the physical and chemical properties of the snowpack (Warren, 1982;Gardner and Sharp, 2010).Since the absorption of solar energy affects in turn the physical and chemical properties of the snowpack, snow albedo is involved in several feedback loops (e.g.Flanner and Zender, 2006;Doherty et al., 2010) which generally enhance snow metamorphism and melt and are thus of crucial importance for the evolution of snow-covered area and more generally for the Earth climate.Measuring, understanding and modelling snow spectral albedo is therefore essential.
The variations of snow spectral albedo with snow microstructure can be well explained using the specific surface area (SSA) of snow, i.e. the ratio of the surface of the iceair interface to the mass of ice (e.g.Grenfell and Warren, 1999;Carmagnola et al., 2013).Snow spectral albedo de-Published by Copernicus Publications on behalf of the European Geosciences Union.
creases with SSA in the near-infrared so that an increased amount of solar energy is absorbed, which further accelerates snow metamorphism (Picard et al., 2012).Ice crystal habit influences the spectral albedo but this effect is not fully understood yet (e.g.Picard et al., 2009;Libois et al., 2013;Kokhanovsky and Zege, 2004;Malinka, 2014;Liou et al., 2014;He et al., 2014).
Light-absorbing impurities (LAI) such as mineral dust, soot or algae cause a decrease in snow albedo in the visible wavelengths (Warren, 1982).This impact is enhanced at low SSA, which results in complex interdependencies and gives raise to an additional snow albedo feedback loop (Doherty et al., 2010).LAI also tend to concentrate at the surface during melt, further lowering the albedo (Sterle et al., 2013).Modelling the effect of LAI on snow spectral albedo is challenging since large uncertainties are associated with their nature, refractive indices (e.g.Dang et al., 2015) and physical properties.The size distributions and exact location with respect to the ice matrix of LAI in snow also induce large uncertainties in the simulated spectral albedo (Flanner et al., 2012;Liou et al., 2014;He et al., 2014).
Because the contrast between ice and water refractive indices is small, the effect of the presence of liquid water on snow spectral albedo is subtle.Warren (1982) explained that it simply increases the effective optical radius.As detailed in Green et al. (2006), Dozier et al. (2009) and Gallet et al. (2014), the liquid water absorption features are shifted towards shorter wavelengths.Accurately predicting the amplitude of the shift for a given liquid water content (LWC) is somewhat difficult since it would require knowledge of the exact location of the liquid water with respect to the ice matrix (e.g.Gallet et al., 2014).Although the retrieval of LWC from spectral albedo seems challenging, the presence of liquid water in snow can be detected using hyperspectral measurement by searching for absorption features minimum around 1000 nm (Green et al., 2006).
Slope and roughness of the snow surface also largely affect snow spectral albedo.Increased roughness generally leads to a decrease in albedo (e.g.Zhuravleva and Kokhanovsky, 2011).Slope modifies the effective incident zenith angle of solar radiation.The effect on the albedo thus depends on the ratio of diffuse to total incoming radiation (e.g.Wang et al., 2016;Dumont et al., 2011).This implies that the surface slope effect on measured snow spectral albedo varies with the wavelength along with the spectral ratio of diffuse to total solar irradiance.
Snow spectral albedo is therefore extremely informative on the state of the snow "surface".Such spectral observations from either in situ or airborne sensors are, however, sparse in time and space.Dozier et al. (2009) provided an overview of the capability offered by hyperspectral remote sensing to study the evolution of optical radius, surface liquid water and LAI radiative forcing.Painter et al. (2013) and Seidel et al. (2016) developed algorithms to retrieve optical radius (which can be related to SSA via r opt 3 ρ ice SSA ; e.g.Gallet et al., 2009), LAI radiative forcing and broadband albedo from airborne hyperspectral images.These data offered a new insight into spatial and temporal evolution of snow surface parameters over large areas.Carmagnola et al. (2013) and Negi and Kokhanovsky (2011) used in situ measured spectral albedo and related its variations to SSA and LAI content.These studies all focused on spectral measurements sparse in time.More recently, Libois et al. (2015) and Picard et al. (2016a) presented a 3-year time series of spectral albedo acquired at Dome C, Antarctica, using an automatic spectral albedometer called Autosolexs.They developed a methodology to correct for the instrument artefacts, assess the measurements uncertainties and retrieve near-surface snow SSA.
Snow surface albedo variations at Dome C are essentially due to SSA variations because the LAI content is too low to significantly affect the albedo (Warren et al., 2006) and the snowpack remains dry throughout the year.In contrast, in alpine conditions the presence of LAI and liquid water can significantly affect snow albedo (e.g.Di Mauro et al., 2015).This motivates the deployment of a second Autosolexs at the Col de Porte site in the French Alps to monitor the evolution of snow spectral albedo during one snow season.
The main goal of the paper is to investigate how alpine snow spectral albedo variations can be attributed to variations of surface and near-surface snow properties, namely SSA, effective impurity content, presence of liquid water and surface slope.We first investigate how an analytical formulation of spectral albedo as a function of SSA, effective impurity content and slope can be used to simulate the measured albedo.We then investigate the theoretical uncertainties associated to the optimal snow characteristics.Finally, the consistency of the time series of optimal SSA, slope, effective impurity content and presence of liquid water with respect to independent snow and meteorological measurements is investigated.Section 2 provides an overview of the study site and instruments.Section 3 describes the method used to analyse albedo variations.Section 4 provides results and discussion.

Data
The measurements were taken at Col de Porte, which is located at 1325 m altitude in the Chartreuse mountain range, France.An exhaustive description of the study site and the long-term measurement instruments is provided in Morin et al. (2012).The long-term measurements used in this study are 2 m air temperature, surface temperature, direct and diffuse broadband incident short-wave radiations, snow depth and snow water equivalent (SWE), snowfall and rainfall rates.

Spectral albedo measurements
From January to May 2014, the Autosolexs instrument was installed at Col de Porte site (Fig. 1) to automatically mea-sure spectral albedo.This instrument is similar to the one described in Libois et al. (2015) and Picard et al. (2016a), except that it only features one albedo head.The upward (downward) optic is set up 2.4 (2.1) m above bare soil.The height of the albedo heads and, consequently, the field of view of the sensor are thus varying with snow depth.The device acquired an upward and downward spectrum every 12 min over the 350-1050 nm range with an effective spectral resolution, i.e. spectral bandwidth, of 3 nm.Two automatic cameras provided qualitative information on the weather, snow and device state during the entire season.
Contrary to Dome C site, the Col de Porte site features large snow precipitation events which cover the upwardlooking channel.Consequently, the device was cleaned up manually after each snowfall although both the upward-and downward-looking domes are heated and ventilated.The snow surface below the measurements head was not flat due to local topography.Figure 1 indeed shows that there is a slight north-facing slope below the measurement head.Note that the tilt of the sensor was recorded during the measurements and remained smaller than 0.5 • .
Following Picard et al. (2016a), every acquired spectrum was corrected for (i) dark current and stray light, (ii) integration time scaling, (iii) calibration and (iv) collector angular response.The correction of the collector angular response was applied in two steps: (i) cross calibration of the two entrance optics under the same illumination conditions and (ii) cosine response correction on the direct component of the incident radiation as detailed in Picard et al. (2016a) (Sect. 3.3.3 and 3.3.4).The corrected spectra are used to compute the bi-hemispherical reflectance (Schaepman-Strub et al., 2006), simply called spectral albedo in the following.

Additional in situ measurements
During this snow season, additional in situ measurements were carried out: -The surface SSA was measured on 22 January 2014 using the ASSSAP instrument (Arnaud et al., 2011) between 12:00 and 13:00, a few metres away from Autosolexs.The experimental protocol for surface SSA measurements is described in Libois et al. (2015) and Carmagnola et al. (2014).
-The vertical profile of impurity content was measured on 11 February 2014.Refractory black carbon (BC) was measured using an SP2 instrument, following the procedure described in Lim et al. (2014).Insoluble dust measurements were performed using a microparticle counter (Coulter counter © Multisizer III) for particles with a diameter ranging from 1 to 30 µm, divided in 300 equivalent size channels.The total mass of dust was calculated from the volume size distribution, assuming a density of 2.5 g cm −3 (Delmonte et al., 2002).
-The vertical LWC profile was measured for the whole snowpack on a weekly basis at Col de Porte site using a dielectric probe at 13 MHz (Brun et al., 1989).Fifteen measurements were available within the period of observation of the albedo (four in January, four in February, four in March and three in April).

Atmospheric model outputs
Knowledge of the spectral ratio of diffuse to total solar irradiance and of atmospheric aerosol deposition fluxes is crucial to understand the evolution of the measured snow spectral albedo.In this respect, we used outputs from the atmospheric model ALADIN-Climate at Col de Porte site to calculate diffuse and direct spectral solar irradiance and to investigate the temporal evolution of dry and wet aerosols deposition fluxes.ALADIN-Climate is a regional climate model based on a bi-spectral semi-implicit semi-Lagrangian scheme.The version 5.3 (Nabat et al., 2015) is used in the present study with a 50 km horizontal resolution, 31 vertical levels and the ERA-Interim reanalysis (Dee et al., 2011) as lateral boundary forcing.This model includes a prognostic aerosol scheme, adapted from the GEMS/MACC aerosol scheme (Morcrette et al., 2009;Benedetti et al., 2011;Michou et al., 2015).The main aerosol species are represented: dust, sea salt, sulfate, black carbon and organic particles.Natural aerosols (dust and sea salt) are emitted from the surface depending on surface wind and soil characteristics, while anthropogenic emissions come from external inventories (Lamarque et al., 2010).The spatial domain of our simulations has been designed to include all the sources generating aerosols that can be transported over the French Alps, such as the Saharan desert and a large part of the northern Atlantic Ocean.Aerosols coming from longer-range transport (e.g.fires in America) are not considered.
More details about the aerosol scheme, as well as an evaluation showing the performance of the scheme, can be found in Nabat et al. (2015).All these aerosols interact with the short-wave and long-wave radiation scheme.Even though ALADIN-Climate is a regional climate model, the model has the ability to reproduce the observed weather chronology thanks to the spectral nudging method (Radu et al., 2008), which enables us to keep large-scale atmospheric conditions from the boundary forcing.The accuracy of the chronology of meteorological episodes is indeed essential to correctly represent the chronology of aerosol deposition on the snowpack.In this simulation, surface pressure, wind vorticity and divergence and specific humidity were nudged towards ERA-Interim.
The outputs of the ALADIN-Climate model used in this study are hourly total aerosols optical thickness, total water vapour column and total ozone column and dust and black carbon wet and dry deposition fluxes.In this study, we used the outputs of a single model cell, the closest to Col de Porte site.The horizontal distance of the model cell centre to the site is 22.6 km and the grid cell is located 800 m below Col de Porte site.

Estimation of the direct-to-diffuse solar irradiance ratio
The ratio of diffuse to direct irradiance is required to perform accurate correction of the measured spectrum (e.g.Picard et al., 2016a).The calculation of such ratio requires the knowledge of the atmospheric profiles and of the local cloud optical thickness, τ , at 550 nm.We thus estimate this local cloud optical thickness using both the atmospheric profiles from ALADIN-Climate and local meteorological observations at Col de Porte (Morin et al., 2012) to overcome the coarse resolution of ALADIN-Climate.The first step consists in estimating a relationship between τ and the broadband direct (SW dir ) to total (SW tot ) ratio.For this purpose, the SBDART detailed radiative model (Ricchiazzi et al., 1998) was used to calculate SW dir over SW tot from varying cloud optical thicknesses and mean atmospheric conditions at Col de Porte (aerosols optical thickness, total ozone column and total water vapour column).The mean atmospheric conditions were derived from ALADIN-Climate outputs and the measured 2 m air temperature average over the measurements period.A regression equation (Eq. 1) was derived from those results.
where µ s = cos θ s is the cosine of the solar zenith angle.As a second step, Eq. ( 1) was used to estimate τ from SW dir over SW tot ratio measured at Col de Porte.
Finally, the outputs of the ALADIN-Climate (namely hourly total aerosols optical thickness, total water vapour column and total ozone column) together with measured air temperature and estimated τ were used as inputs for the SB-DART (Ricchiazzi et al., 1998) model in order to compute the hourly diffuse and direct spectral irradiance.The ratio between these variables is used in our analysis.Note that the broadband SW dir and SW tot estimated from Col de Porte measurements and simulated with SBDART agree within ±10 W m −2 .The accuracy of the simulated spectral direct to diffuse solar irradiance ratio has not been evaluated due to the absence of measurements.The difference in elevation between the ALADIN-Climate grid cell and Col de Porte site might lead to an overestimation of the diffuse fraction in the visible wavelengths.

Spectral albedo dependence on snow SSA and impurity content
In order to relate albedo to snow properties, we use the theoretical formalism of Kokhanovsky and Zege (2004).Several assumptions are made: (i) the snowpack is horizontally and vertically homogeneous, which means only one bulk SSA and impurity content value are used to explain albedo variations; (ii) the surface is flat; (iii) snow phase function and single scattering albedo are implicitly described by the asymmetry factor, g, the absorption enhancement parameter, B, and SSA; (iv) the surface and the sensor are perfectly horizontal.The effect of surface slope is investigated in Sect.3.3.Under these conditions, following Libois et al. (2013) and Picard et al. (2016a), snow albedo can be written as where r diff (λ, θ ) is the diffuse-to-total-irradiance ratio, ρ = 917 kg m −3 is the ice density at 0 • C and n i (λ) is the imaginary part of the ice refractive index taken either from Warren and Brandt (2008) (default) or from Picard et al. (2016b).B = 1.6 and g = 0.85 are constant and taken from Libois et al. (2014).m imp is the complex refractive index of BC taken from Flanner et al. (2012).ρ imp is BC density and is set to 1270 kg m −3 according to Flanner et al. (2012).c imp is the effective impurity mass per unit of snow mass (kg kg −1 ) where effective stands for BC optically equivalent content.Following Picard et al. (2016a), we introduce a scaling factor A to account for several artefacts and shortcomings in the measurement technique.This scaling factor relates the measured albedo, α meas , and the theoretical albedo obtained from Eq. ( 2), α th .

Effect of the slope on measured albedo
The slope of the snow surface introduces a change in the solar irradiance with respect to the solar radiation incoming on a perfectly horizontal surface.This change is of crucial importance for our application since it is wavelength dependent.
In addition to the change in the effective sun incident angles, the surface slope modifies (i) the solid angle under which the sky is viewed from the surface and thus the incoming amount of diffuse solar radiation and (ii) the total amount of radiation received by the snow surface since it receives some of the radiation reflected by the adjacent slopes.The upward radiation measured by the horizontal sensor is also modified with respect to what would happen if the sensor were parallel to the surface since part of the field of view sees the atmosphere and not the snow surface.In the following, we assume that (i) both diffuse solar radiation and reflected radiation are isotropic and (ii) the surface slope is small and local enough not to modify significantly the solid angles under which the incoming and reflected radiations are measured with respect to the solid angles that would apply in the case of an horizontal surface.With these assumptions, the slope of   6) for an horizontal surface and Eq. ( 8) for tilted surface (φ n = 300 • , θ n = 5 • ) using estimated diffuse-to-total-irradiance ratio and φ s and θ s for 8 March 2014 at noon.The black lines correspond to SSA = 40 m 2 kg −1 and the blue lines to SSA = 5 m 2 kg −1 .(b) Example of measured, predicted albedo and theoretical albedo on an horizontal surface for three dates in the season.For the red lines (and blue and black), optimal SSA is 36 m 2 kg −1 (6.6, 2.6) and optimal c imp is 0.5 ng g −1 (15.1, 328).A is set to 0.943.
the surface only affects the effective sun zenith and azimuth angles and thus the direct solar irradiance (see details in Appendices A and B) .Note that for fully cloudy days, under these assumptions, the slope has consequently no influence on the measured albedo.
Let θ s and φ s be the direction of sun with respect to a perfectly horizontal surface, θ n and φ n the slope and aspect of the surface, respectively" and θ s and φ s the effective direction of the sun with respect to the tilted surface (Fig. 2).Then cos θ s = K cos θ s (e.g.Dumont et al., 2011), where K is the relative change in the cosine of the sun effective incident angle to the local tilt of the snow surface: This leads to the following relationship between the measured and the theoretical albedo on a tilted surface (see details in Appendices A and B): The variations of measured spectral albedo with SSA, c imp and surface slope are illustrated in Fig. 4a.The effect of the anisotropy of diffuse solar radiation is discussed in Bogren et al. (2016) while the effect of the anisotropy of the reflected radiation is discussed in Dumont et al. (2010) and Carmagnola et al. (2013).The anisotropy of the sky diffuse component and of the reflected radiation is second order as long as the cosine response correction is small (Picard et al., 2016a;Carmagnola et al., 2013).
Figure 3 illustrates the raw measured albedo diurnal cycles for a cloudy day (red lines) and for a clear-sky day (blue lines).As expected from Eq. ( 8), the diurnal cycle is more pronounced for the clear-sky day, the albedo evolution being non-symmetric with respect to solar noon probably due to both slope and changes in snow properties effects.

Albedo variations analysis method
In order to relate variations of spectral albedo to variations of surface snow properties, we apply the following methodology to the measured albedo.The main idea of the methodology is to fit Eq. ( 8) using optimal parameters for each spectrum.Equation (8) indeed contains four unknowns, namely A, SSA, c imp and K.For moderate to high amounts of impurities in snow, several (A, c imp , SSA) triplets can lead to approximately the same modelled spectrum.For this reason, we chose to first set the scaling factor A to a constant value for the whole season (step 1).Optimal SSA, c imp and K are then estimated for each spectrum (step 2).The diurnal cycles of parameter K are then used to estimate a daily value for surface slope and aspect (step 3).Although this step is not required for the estimation of optimal SSA and c imp , it provides a further verification of the physical consistency of the methodology.Finally, the measured spectra are analysed with respect to the presence of liquid water (step 4).
Step 1: estimate the scaling factor A A seasonal value of A is estimated using Eq. ( 8) with r diff = 1 (selected fully cloudy days from visual inspection of the photographs and θ s smaller than 65 • ) and a non-linear least squares method (provided by Python scipy.optimize.leastsqfunction).To avoid avoid multiple solutions (A and c imp ) in case of moderate to high amounts of impurities in snow, only the spectra from the beginning of season were used for this estimation.
Quantiles 25 and 75 of A distribution are used to propagate uncertainties on the near-surface and surface properties predicted from the spectra.
Step 2: estimate optimal SSA and c imp Once the value of A is set, optimal SSA and log 10 (c imp ) values are estimated from the spectrum within 400-1050 nm along with K using Eq. ( 8) and the non-linear least square optimization method.For cloudy days, the optimization is performed with K = 1.Clear-sky days are selected as days with τ < 0.01, τ being estimated using Eq.(1).Illustration of step 2 is provided for three spectra in Fig. 4b.
After these steps, spectra are filtered based on the root mean square deviation (RMSD) between the measured spectrum and the optimal spectrum calculated with Eq. ( 8) with a threshold of 0.022.This filtering ensures that the measurement artefacts are reasonably accounted for in the optimal snow surface properties estimation.The threshold value is higher than the one used in Picard et al. (2016a) and was set to account for the discrepancies between the modelled and measured albedo spectra due, among others, to the impurity type assumption in the model.Indeed, the model used in this study only includes BC.Alpine snowpacks are frequently affected by deposition of Saharan dust (e.g.Di Mauro et al., 2015) that has a reddish spectral signature.The presence of red dust in snow can induce discrepancies between the modelled spectrum with BC only and the measured spectrum (Fig. 4 in Warren and Wiscombe, 1980).This "dusty" pattern is clearly visible on the black measured spectrum of Fig. 4b in the visible wavelengths (400-500 nm).
Note that by using a seasonal A value for every spectra, we assume that the measurement artefacts are the same under cloudy and clear-sky conditions and for varying solar zenith angles.
Step 3: estimate daily optimal surface slope and aspect Section 3.3 explains how the measured albedo varies with the surface slope and aspect.The slope and aspect below the sensor evolves in time with respect to precipitation, melt and snow transportation by the wind.Unfortunately no measurement of surface slope and aspect is available during this winter season.Consequently, using K diurnal cycles and Eq. ( 7), we estimate daily optimal values of surface slope angle and aspect for fully clear-sky days.
This step is not required for the estimation of optimal SSA and c imp but it indirectly validates that K optimization has not compensated for other artefacts than slope.In other words, estimating physically consistent values of daily slope and aspect from K diurnal cycles further ensures the consistency of the spectra correction.
Step 4: detect if the snow surface is wet or dry Taking benefit of the spectral shift of the absorption feature at 1030 nm in presence of liquid water (Green et al., 2006), we apply the following method to distinguish wet from dry snow.The spectra are first averaged using a 20 nm moving window average in order to reduce noise before minimum calculation.We then compute the wavelength of minimum albedo in the 1000-1050 nm range and apply the following criteria to detect liquid water presence: λ water threshold has been set after studying the distribution of argmin(α) on the whole spectra dataset.0.943.This value is close to one, which indicates that the actual instrumental response deviates only slightly from the ideal one.The spread of A values is small (quantile 25, q 25 , of 0.920 and quantile 75, q 75 , of 0.964).This spread could probably be attributed to residual measurement artefacts such as deposition of precipitation particles on the measurement head, small direct incoming radiation and reproducibility errors due to the optical switch (Picard et al., 2016a).In the following, A is set to 0.943.Quantiles 25 and 75 are used to propagate uncertainties of the near-surface and surface properties predicted from the spectra, i.e. the optimal parameters.

Determination of slope angle and aspect
Figure 6a illustrates that the optimal parameter K (green diamonds) follows a diurnal cycle with lower values in the morning and higher values in the afternoon as predicted by Eq. ( 7).It also illustrates the good agreement between the optimal K and simulated K (blue crosses, step 3) for 11 March (perfectly clear day) and the poorest agreement for 10 March where clouds were detected, especially in the afternoon.
Figure 6b shows the daily estimates of θ n and φ n obtained after step 3 over the season.The estimated slope varies from 3 to 10 • while the aspect mainly ranges between 300 and 360 • (discarding obvious outlier values).The uncertainty range of estimated θ n is generally lower than 2 • while larger uncertainty ranges, up to 50 • , are estimated for aspect.The average slope and aspect are in agreement with what can be visually estimated from Fig. 1.The seasonal evolution of slope and aspect is related to snow evolution (the red line in Fig. 6b shows the measured snow depth).Slope is smaller just after a precipitation event and increases during the melt season.
For comparison, the same method has been applied to the albedo database measured in Dome C and described in Picard et al. (2016a).The snow surface is horizontal at large scale but the surface can be locally rough due to wind-drift effects.The method applied to Dome C data during the 2012-2013 season leads to slope angles between ±2 • , thus providing an insight into the accuracy of the method.

Determination of λ water
Figure 7b shows the distribution of the wavelength at which the minimum reflectance value is reached for all the measured spectra with RMSD lower than 0.022 over the 400-1050 nm range.The distribution is bimodal as predicted by Green et al. (2006) and exhibits two peaks, the first one at 1029 nm and the second one at 1034 nm.The first one is likely to correspond to wet snow, while the second is for dry snow.λ water in Eq. ( 9) is consequently set to 1032 nm in the following.Figure 7a provides an illustration for two spectra (a wet and a dry one) measured on 9 March 2014.
The Cryosphere, 11, 1091-1110, 2017  9).Spectra are discarded when the RMSD between predicted and measured spectra is lower than 0.022 over the 400-1050 nm range.λ water is indicated by the black vertical line.

Theoretical uncertainty and representativeness of the estimated snow parameters
Section 3.3.1 describes the methodology applied in this study to estimate optimal SSA, c imp and presence of liquid water from the measured spectra.The uncertainty related to the optimal SSA is discussed in detailed in Picard et al. (2016a) together with the vertical representativeness of the estimated SSA.However, our study differs from Picard et al. (2016a) since (i) the snowpack contains larger amount of LAI, (ii) the snow surface is not perfectly horizontal and (iii) the snow can be wet.These points might indirectly affect the estimated uncertainty of the optimal SSA and of the other snow parameters, which is investigated in the section below.Note that the results presented below in Sect.4.2.1 and 4.2.2 have been obtained under clear-sky conditions.The methodology applied in this study for cloudy sky intrinsically leads to larger uncertainties for cloudy-sky conditions.

Effect of slope on the optimal parameters
Figure 8a shows an example of measured spectrum on 10 March 2014 (grey solid lines) compared to the optimal spectrum predicted with slope (i.e. from step 2, red dotted lines), assuming a perfectly horizontal surface (blue dotted lines) and slightly varying slope zenith (and azimuth) angles from +2 • (+10 • ) (green dotted line).It illustrates that the (i) the shape of the spectrum computed with no slope does not agree with the measured albedo, (ii) the best agreement is obtained for the red dotted line (i.e. with slope) and (iii) a small variation of the slope angles induces a variation of 10 % for the optimal SSA and of 13 ng g −1 for optimal c imp .
More generally, Fig. 6b shows that the spread of A distribution is translated in an uncertainty of typically less than 2 • for slope zenith angle and of less than 10 • for the aspect.Using the simulated spectrum from Fig. 4 and adding slope and aspect variations within these ranges, we obtain retrieved SSA variations up to 20 % for 40 m 2 kg −1 and up to 10 % for 5 m 2 kg −1 .The larger uncertainty associated with the higher SSA value is because tilt effect is proportional to the albedo value (higher for higher SSA).Optimal c imp variations are up to 25 (and 40) ng g −1 for SSA = 40 m 2 kg −1 , c imp = 0 (100) ng g −1 .For lower SSA (5 m 2 kg −1 ), c imp variations range from 4 to 40 ng g −1 for initial impurity content ranging from 0 to 1000 ng g −1 .
www.the-cryosphere.net/11/1091/2017/The Cryosphere, 11, 1091-1110, 2017 In addition, the surface slope also affects the value of argmin(α).Using the spectral albedo from Fig. 4 and slope angle varying between −20 and +20 • , we nevertheless found that argmin(α) does not vary for slope angle lower than 10 • and that the maximum variation for steeper slope is less than ±0.5 nm.

Coupled effect of SSA and c imp on spectral albedo
Figure 8b shows an example of measured spectrum on 3 April 2014 (grey solid lines) compared to the optimal spectrum predicted with slope (i.e. from step 2, red dotted lines).The blue dotted line corresponds to the optimal albedo obtained while increasing the original optimal value of c imp by 50 ng g −1 and the green dotted line to the optimal albedo obtained while adding +15 % to the original optimal value of SSA.It illustrates that changes in the SSA (and c imp ) induce changes in the optimal values of c imp (SSA).
More generally, using again the simulated spectrum presented in Fig. 4, we first investigate the change in estimated SSA while varying c imp by ±50 ng g −1 .The larger variations, up to 15 %, are obtained for very low c imp (1 ng g −1 ) and small SSA (5 m 2 kg −1 ).As soon as c imp is higher than 100 ng g −1 the SSA variations are smaller than 13 % for SSA = 5 m 2 kg −1 and than 5 % for higher SSA (40 m 2 kg −1 ).
Then, we investigate the change in optimal c imp while varying SSA by ±15 %.It shows that a ±15 % uncertainty on SSA leads to a ±20 % uncertainty on c imp .

Optimal SSA and liquid water
The presence of liquid water modifies the spectrum.Thus it necessarily affects the value of optimal SSA estimated from the spectrum.Gallet et al. (2014) provide an overview on how it affects the SSA estimated from 1310 nm reflectance.The conclusion is different in our case where the whole spectrum is used.Although a detailed modelling study is beyond the scope of the study, using the relative difference in optimal SSA, SSA = SSA wet −SSA dry SSA wet , between two consecutive spectra (absolute time difference of 12 min), we were able to investigate this effect.The distribution of SSA for the whole measurements period is peaked.Quantiles 25 (and 50 and 75) are −4.15% (0.025 and 6.99 %).Assuming the change in SSA between the two spectra is only due to refreezing and not to metamorphism, this indicates that the SSA of the wet spectrum is not generally higher or smaller than the SSA of the dry spectrum and that an upper bound of the uncertainty on the optimal SSA due to the presence of liquid water is 7 %.

Vertical representativeness of optimal impurity content
Real snowpacks are not vertically homogeneous as assumed in our method and often displays high vertical gradient of SSA or impurity content (e.g.Sterle et al., 2013).However only one bulk SSA and c imp values are estimated is order to stabilize the optimization.The vertical representativeness of SSA is discussed in Picard et al. (2016a).Using the same methodology, we analyse here the vertical representativeness of c imp .We use the twostream radiative model TARTES (Libois et al., 2013) to compute the albedo of semi-infinite medium with a fixed SSA, density and c imp .A layer with variable SWE h is added on top of it with a slightly different value of c imp (c imp +δc).The relative contribution of the uppermost layer to the albedo is then defined as (α(h) − α(0))/(α(∞) − α(0)), where α(h) is the albedo computed with TARTES averaged over 400-1050 nm for a two-layer snowpack with the uppermost layer of h.The value of δc does not change the results as long as it is small (e.g.+10 % in Fig. 9).Thus, the vertical representativeness calculated here is only valid for slightly inhomogeneous snowpack.
Figure 9 shows the relative contribution of the uppermost layer as a function of its SWE for two SSA values (5 m 2 kg −1 plain lines, 40 m 2 kg −1 crosses) and varying values of c imp (from 1 to 500 ng g −1 ).It shows that the higher the SSA and c imp , the higher the contribution of the uppermost cen-timetres of the snowpack to the albedo value.For instance for low SSA (5 m 2 kg −1 ) and high c imp (500 ng g −1 ), the top 10 kg m −2 of the snowpack contributes to more than 80 % of the signal.For a density of 400 kg m −3 , this corresponds to the uppermost 5 cm.For higher SSA (40 m 2 kg 1 ) and lower c imp (10 ng g −1 ), the top 30 kg m −2 of the snowpack contributes to more than 80 % of the signal.For a density of 200 kg m −3 , this corresponds to the uppermost 15 cm.This can be explained since (i) the higher the SSA the lower the light penetration depth in the snowpack and (ii) the presence of impurities shortens the light penetration depth in the snowpack (e.g.Libois et al., 2013).

Time evolution of the optimal SSA
Section 4.3 and 4.4 focus on of the optimal SSA and effective impurity content time series predicted from the measured spectra and their consistency with snow and meteorological conditions.

Seasonal evolution of "noon" SSA
Figure 10a, b present the seasonal evolution of all the optimal SSA and the optical radius, r opt , predicted from the measured spectra between 12:00 and 13:00 each day.Clearsky conditions are represented in blue and cloudy conditions in red.Vertical bars represent the uncertainties derived from q 25 and q 75 of A distribution.The SSA ranges from more than 70 m 2 kg −1 after snowfall down to 2 m 2 kg −1 .For clearsky days uncertainties are small (typically less than 10 %) whereas cloudy days feature larger uncertainties.The range of variations is consistent with measured values of surface SSA presented in Morin et al. (2013) using the DUFISSS instrument (Gallet et al., 2009) in 2010.The SSA decay is also consistent with current understanding and parameterizations of SSA evolution (e.g.Carmagnola et al., 2014;Flanner and Zender, 2006;Schleef et al., 2014).
On 22 January 2014, surface SSA was measured between 12:00 and 13:00 with ASSSAP at 44 ± 4 m 2 kg −1 .The sky was partially clear during the day and the sensor was covered by snow in the morning due to recent precipitation leading to valid spectra only after 13:00.The optimal SSA value at 13:00 is (42.6, 43.1, 43.4) m 2 kg −1 (q 25 , q 50 , q 75 ), which is only slightly lower than ASSSAP value and within ASSSAP uncertainty range.
Note that the small SSA increase (mostly visible in Fig. 10b) by 2 April 2014 is an artefact induced by the appearance of grass in the field of view of the sensor, leading to higher reflectance in near-infrared wavelengths and thus higher optimal SSA.

SSA diurnal cycles
Figure 11a, b show the diurnal evolution of SSA inferred from albedo variations.Figure 11a zooms on the details of the diurnal cycle from 6 to 9 March 2014.The diurnal cycle evolves from higher SSA in the morning to lower SSA in the afternoon in the absence of precipitation, which can be explained by snow metamorphism and by the snow intrinsic albedo feedback (Dumont et al., 2014).
Additionally, Fig. 11c presents the daily SSA decay rate in % per day inferred from the diurnal cycle.The grey area corresponds to the 15 % uncertainty estimated in Picard et al. (2016a) for the optimal SSA retrieval in Antarctica.SSA decay rates are larger after snowfall.Such variations are also described by Seidel et al. (2016) for the investigation of r opt variations from airborne hyperspectral imagery.the Alps (e.g.Di Mauro et al., 2015).Low impurity content (typically smaller than 50 ng g −1 ) is affected by large uncertainties.On 11February, the equivalent BC content measured in the first 2 cm of the snowpack was 17 ± 6 ng g −1 .The day was partially cloudy.Only two spectra before 12:00 are valid and measured during clear-sky conditions (from the camera images).For these spectra, the optimal SSA values are [22,23,24] m 2 kg −1 and c imp = [0.5, 2., 6.6] ng g −1 (q 25 , q 50 , q 75 ).The measured c imp is higher than the optimal one though lying in the range of impurity content values for which the retrieval is highly uncertain.

Time evolution of the optimal effective impurity content
Diurnal cycle of c imp is not investigated in this study since the uncertainties associated with the c imp values are too high with respect to the diurnal variations.
Figure 12b shows the seasonal evolution of c imp but distinguishes the prevailing colour of impurities (black indicated by black crosses or red indicated by red dots).This distinction is based on the RMSD calculated between the modelled and the measured albedo within the 400-500 nm wavelengths range.If this RMSD in the visible wavelengths (400-500 nm) is larger than the RMSD over the whole spectrum, c imp is represented in red.This indeed indicates that the predicted and modelled spectrum agrees well except in the 400-500 nm wavelengths as illustrated by the black spectrum in Fig. 2b and that the prevailing impurities at that date may have a reddish spectral signature.The blue vertical bar represents days for which melt was detected on the SWE measurements (Fig. 12a).Figure 12c shows the wet and dry deposition predicted by ALADIN-Climate for black carbon and mineral dust at Col de Porte site.
High values of c imp in January and February may be related to dust deposition events.Low values are found beginning of March after a precipitation period.From March to mid-April, the exponential increase in c imp can be related to melt during which part of the impurities concentrate on the surface (Sterle et al., 2013) and to dust deposition events beginning of April.Note that the coarse resolution of ALADIN-Climate limits the accuracy of the precipitation events.Consequently the wet deposition event predicted on 30 March 2014 has no impact on the snowpack because no precipitation was observed during that day at Col de Porte.
Finally, the uncertainty associated with the imaginary part of the refractive index of ice (e.g. Warren and Brandt, 2008;Carmagnola et al., 2014;Picard et al., 2016b) strongly af-fects the relationship between spectral albedo and c imp (see Eq. 5). Figure 12b shows the values of c imp estimated using the values of the ice refractive index proposed in Picard et al. (2016b) instead of Warren and Brandt (2008) (green dots).c imp is significantly lower using Picard et al. (2016b), especially for low impurity content.

Presence of liquid water
Figure 13a shows the comparison between the wavelength of minimum reflectance and the surface LWC measured in the field for 15 snow pits.The agreement between our method of detection and field observations is perfect for these 15 dates.
Figure 13b, c show the results of the liquid water detection method over the whole season.The detection is in good agreement with measured surface temperature.The accuracy of the measured surface temperature is estimated equal to 1 K. False detection cases, i.e. when the snow is detected as wet but surface temperature is negative, i.e. less than 272.15K to account for the measurement accuracy, are less than 3.5 %.False detection cases may originate from measurement artefacts since the signal-to-noise ratio for the considered wavelengths is high and differences in the field of view of the spectrometer and of the long-wave downlooking sensor used to calculate surface temperature.Figure 13c shows that melt-freeze diurnal cycles are well depicted with dry snow in the morning and late afternoon and wet snow around noon.

Concluding remarks
This study introduces a continuous dataset of spectral albedo acquired from January to May 2014 from 400 to 1050 nm at the Col de Porte alpine site as a follow-up of the study conducted using the same instrument at Dome C in Antarctica during 3 years (Picard et al., 2016a).The alpine snowpack radically differs from inner Antarctica snowpack and the spectral albedo variations cannot be only explained by the evolution of near-surface SSA and solar zenith angle.We investigated how the variations of measured spectral albedo can be explained by the evolution of other factors such as ef-fective impurity content, surface slope and aspect and presence of liquid water.For this, we used analytical formulation of spectral albedo as a function of these near-surface parameters.Results show that the measured spectra are accurately simulated using this formulation.This formulation disentangles the effect of slope, SSA and effective impurity content on spectral albedo and consequently provides optimal values of these parameters for each measured spectra.This study also demonstrates that wet and dry snow surfaces can be distinguished due to the spectral shift between the refractive index of ice and water around 1000 nm.
The uncertainty associated to the optimal SSA was discussed in detail in Picard et al. (2016a) and theoretically estimated to be better than 15 %.In our case, in the presence of impurities and surface slope, the accuracy is degraded.The detection of wet snow surface (only binary information) at the snow surface is less affected by measurement artefacts and surface tilt.Finally, the estimation of an effective impurity content is subject to large uncertainties, making c imp values lower than 50 ng g −1 challenging to estimate as pointed out by Warren (2013).The study emphasizes the need for  measuring snow albedo over snow surface as flat as possible in order to be able to infer accurately snow surface variables.
The influence of surface roughness, incident and reflected radiations and of the location of the impurities with respect to the ice matrix deserves future work.Further refinements on the calculation of the spectral ratio of direct to total irradiance or simultaneous measurements of this ratio and albedo are required to improve the accuracy of the method.It must be underlined that the detailed study of the accuracy of the estimation of near-surface snow parameters from spectral albedo would require more extensive validation measurements than presented in this paper, i.e. micro-tomography evaluation of SSA and systematic measurements of surface impurities content.However, optimal near-surface SSA predicted from the spectra exhibits strong seasonal and diurnal cycles that can be related to snow and meteorological conditions.Near-surface optimal effective impurity content also exhibits strong variations along the season that can be related to surface enrichment process due to melt and Saharan dust deposition events that frequently occur in the French Alps (Di Mauro et al., 2015).
Such time series of spectral albedo are required to understand and quantify the evolution of snow albedo in relationship with snow surface variables such as SSA and surface impurity content.They are also a unique opportunity to better understand the evolution of near-surface SSA and effective impurity content during the snow season.They provide a unique dataset to evaluate and refine detailed snowpack model such as Crocus (Vionnet et al., 2012) or SNOWPACK (Bartelt and Lehning, 2002).Furthermore, without the need of a retrieval methodology for snow properties, such measured spectral albedo time series can be assimilated in the snowpack model in order to improve simulation of the snowpack structure (Charrois et al., 2016).This study thus emphasizes the usefulness of hyperspectral optical observations of snow.
Data availability.Time series of snow spectral albedo and superficial snow-specific surface area and impurity content are available through the PANGAEA database (doi:10.1594/PANGAEA.874272).

Figure 1 .
Figure 1.Part of the col de Porte field measurements site on 17 February 2014, 14:35.North direction is the left of the picture.The two measurement heads of Autosolexs are circled in blue.

Figure 2 .
Figure 2. Schematic drawing of the tilted snow surface (blue plane) and associated angles.The grey plane corresponds to the horizontal plane.The black reference frame is the one attached to the tilted surface.The grey arrow represents the vertical with respect to the horizontal plane.H (φ) is the horizon line in the tilted reference frame.θ s and φ s are the zenith and azimuth sun angles in the tilted reference frame.θ n is the slope elevation and φ n the aspect.

Figure 4 .
Figure 4. (a) Spectral albedo simulated using Eq.(6) for an horizontal surface and Eq.(8) for tilted surface (φ n = 300 • , θ n = 5 • ) using estimated diffuse-to-total-irradiance ratio and φ s and θ s for 8 March 2014 at noon.The black lines correspond to SSA = 40 m 2 kg −1 and the blue lines to SSA = 5 m 2 kg −1 .(b) Example of measured, predicted albedo and theoretical albedo on an horizontal surface for three dates in the season.For the red lines (and blue and black), optimal SSA is 36 m 2 kg −1 (6.6, 2.6) and optimal c imp is 0.5 ng g −1 (15.1, 328).A is set to 0.943.

Figure 5 .
Figure5.Distribution of the scaling factor, where A is for overcast days.Spectra used for this analysis corresponds to θ s smaller than 65 • and RMSD between predicted and measured spectrum lower than 0.02 over the 400-1050 nm range.

4
Figure 5 describes the distribution of A for fully cloudy days before 5 March 2014.The distribution has a median value

Figure 7 .
Figure 7. Effect of liquid water presence on measured spectra.(a) Example of two measured spectra on 9 March at 09:00 in blue and 13:12 in green and the associated wavelength of minimum reflectance represented by the vertical dotted lines.The vertical black line corresponds to λ water .(b) Distribution of wavelength of minimum albedo in Eq. (9).Spectra are discarded when the RMSD between predicted and measured spectra is lower than 0.022 over the 400-1050 nm range.λ water is indicated by the black vertical line.

Figure 9 .
Figure 9. Relative contribution of the top layer of a semi-infinite snowpack to albedo averaged over 400-1050 nm as a function of top layer SWE (h) for low SSA (5 m 2 kg −1 , plain lines) and high SSA (40 m 2 kg −1 , dotted lines).The relative contribution is defined as (α(h) − α(0))/(α(∞) − α(0)), where h is the SWE of the top layer of same SSA and c imp modified by +10 % with respect to the bottom layer.α(h) is the albedo of the snowpack constituted of these two layers.

Figure 10 .
Figure 10.Seasonal evolution of near-surface SSA (a), r opt (b) and c imp (c) predicted from measured spectra.Only spectra with RMSD between measured and predicted spectrum lower than 0.022 and measured between 12:00 and 13:00 UTC are indicated.Vertical error bars are obtained using quantiles 25 and 75 of A distribution.Clear-sky days are indicated in blue and cloudy days in red.The measured snowfall rates are indicated in grey in panel (a).

Figure
Figure 10b describes the evolution of the effective impurity content, c imp , inferred from spectral albedo values.The values are affected by large uncertainties for cloudy days as already discussed for SSA.c imp values range between 0 and 1000 ng g −1 .Values up to 400 ng g −1 are found in January and February.The highest values are found at the end of the season after long periods without precipitation.c imp variations range is in agreement with other studies conducted in

Figure 11 .
Figure 11.(a) Diurnal evolution of estimated near-surface SSA between 6 March and 9 March 2014.Blue circles are the values kept for decay range analysis (c).The red circles (and diamonds) are the estimated morning (afternoon) SSA values used to estimate the daily decay rate (c).(b) Seasonal evolution of estimate SSA for every spectra with RMSD lower than 0.022.(c) Estimated daily change rate in % of daily mean SSA values.The grey area corresponds to the 15 % accuracy given in Picard et al. (2016a).

Figure 12 .
Figure 12.(a) Measured snow water equivalent (black cross) in kg m −2 and snow and rainfall rates (grey line) in kg m −2 d −1 .(b) Nearsurface c imp estimated between 11:00 and 13:00.Red circles correspond to c imp values for which the RMSD within 400-500 nm is larger than RMSD within 400-1050 nm.Only values with RMSD lower than 0.022 over 400-1050 are reported.The vertical blue bars indicated melting days estimated from measured SWE. Green dots are the values of near-surface c imp estimated using the value of ice refractive index proposed in Picard et al. (2016b) instead of Warren and Brandt (2008).(c) Wet and dry deposition fluxes simulated by ALADIN-Climate.Dust is in red and black carbon in black.

Figure 13
Figure 13.(a) Wavelength of minimum reflectance as a function of measured volumetric liquid water content (LWC).(b) Measured surface temperature (black crosses).Spectra detected as wet are represented by the vertical blue bars and spectra detected as dry by the vertical yellow bar.Values presented in this figure are only for spectra with RMSD lower than 0.022.(c) Zoom of (b) between 6 March and 9 March 2014.

Figure 14 .
Figure 14.Comparison of Eq. (8) (dotted lines) and Eq.(B6) (solid lines) for the calculation of measured albedo for varying slope angles.The azimuth is set to 0 • and the ratio of diffuse to total solar irradiance and solar zenith angle are taken from 7 March 2014 noon data.