Intercomparison of retrieval algorithms for the specific surface area of snow from near-infrared satellite data in mountainous terrain, and comparison with the output of a semi-distributed snowpack model

This study compares different methods to retrieve the specific surface area (SSA) of snow from satellite radiance measurements in mountainous terrain. It aims at addressing the effect on the retrieval of topographic corrections of reflectance, namely slope and aspect of terrain, multiple reflections on neighbouring slopes and accounting (or not) for the anisotropy of snow reflectance. Using MODerate resolution Imaging Spectrometer (MODIS) data for six different clear sky scenes spanning a wide range of snow conditions during the winter season 2008–2009 over a domain of 46× 50 km in the French Alps, we compared SSA retrievals with and without topographic correction, with a spherical or non-spherical snow reflectance model and, in spherical case, with or without anisotropy corrections. The retrieved SSA values were compared to field measurements and to the results of the detailed snowpack model Crocus, fed by driving data from the SAFRAN meteorological analysis. It was found that the difference in terms of surface SSA between retrieved values and SAFRAN-Crocus output was minimal when the topographic correction was taken into account, when using a retrieval method assuming disconnected spherical snow grains. In this case, the root mean square deviation was 9.4 m2 kg−1 and the mean difference was 0.1 m 2 kg−1, based on 3170 pairs of observation and simulated values. The added-value of the anisotropy correction was not significant in our case, which may be explained by the presence of mixed pixels and surface roughness. MODIS retrieved data show SSA variations with elevation and aspect which are physically consistent and in good agreement with SAFRAN-Crocus outputs. The variability of the MODIS retrieved SSA within the topographic classes of the model was found to be relatively small (3.9 m 2 kg−1). This indicates that semi-distributed snowpack simulations in mountainous terrain with a sufficiently large number of classes provides a representation of the snowpack variability consistent with the scale of MODIS 500 m pixels.


Introduction
Snow on the ground is both a resource and a hazard in mountain regions.As a temporary reservoir of water that is Published by Copernicus Publications on behalf of the European Geosciences Union.

A. Mary et al.: Mountain snow SSA from MODIS
released during the melt season with implications for natural ecosystems, agriculture or hydropower, or as an essential component of winter sports and tourism, it has a significant impact on socio-economic activities.Various types of snow avalanches, and debris flows and floods occurring when snow melts, represent the darker side of the impact of snow on the ground on mountain communities.Not only the presence of snow on the ground, but also its amount and physical properties govern its positive or negative impact on human activities.For example, the vertical profile of the mechanical properties of snow plays a major role in determining its susceptibility to the occurrence of avalanches.During the melt season, knowing the amount of snow available for melting and the melt rate are key factors to assess the melt water flux into rivers and aquifers.These variables not only depend on the amount of snow accumulated, but also on the time evolution of the surface energy balance of the snowpack and its internal physical properties.All of these variables vary greatly in space and time in mountain areas due to complex interactions between meteorological conditions (in particular, wind and precipitation), surface snow conditions and topography.
A wide number of numerical models have been developed to simulate and predict the physical state of snow on the ground, ranging from conceptual representations of snow mass and melt rate (so-called "positive degree-day" approach -PDD) to physically based models explicitly resolving the surface energy balance (Armstrong and Brun, 2008).Coupling such models with appropriate methods to estimate the meteorological conditions in mountain areas leads to integrated modelling systems which are helpful to quantify the amount of snow on the ground and predict its behaviour in the near future.Point-scale ground based data can be used to evaluate model results but due to the large heterogeneity of the snow cover their relevance is often questionable.Alternatively, satellite observations of the snow cover extent can be used as an independent evaluation of the model results and provide indications on the level of confidence that the model can be given.Assimilation of satellite observations into the snow component of land surface models is a promising avenue where the added value of observations and model results should be optimally combined.This requires that the satellite observations are relevant to the state variables used in the snow model.For example, a PDD-based model in which snow mass is the only state variable for snow would not be able to assimilate snow surface temperature data.
Satellite monitoring capabilities with ever increasing spatial resolution on the ground, revisit frequency and spectral resolution, provide information that is increasingly relevant to multi-layered physically-based snow models.Using such models in mountainous areas has been initiated several decades ago with numerous successful applications such as avalanche hazard forecasting (e.g.Durand et al., 1999;Bellaire et al., 2011), mountain hydrology (e.g.Braun et al., 1994;Magnusson et al., 2011) and glacier mass bal-ance (Obleitner and Lehning, 2004;Gerbaux et al., 2005).However, there have been few examples of direct assimilation of satellite data into models of such a degree of complexity (Toure et al., 2011).Dumont et al. (2012) have recently demonstrated how remotely sensed spectral albedo of snow can be assimilated into the detailed snowpack model Crocus, with positive impact on the simulated mass balance of an alpine glacier.Regardless of the frequency range (microwave, visible/near infrared, thermal infrared, . ..) or the acquisition type (active vs. passive), satellite data obtained over mountainous terrain are strongly affected by the topography.This can either occur because a too coarse resolution of the acquired images results in a blending of multiple types of surface conditions into the same pixel (this effect is also present in flatter terrain, but the presence of different types of slopes dramatically increases the mixing), or because of issues inherently associated to the sloping conditions, i.e. multiple reflections on neighbouring slopes, etc.As an example of how such effects impact the retrieval of snow properties from satellite images in mountainous terrain, here we describe and compare different methods used to post-process near-infrared (NIR) reflectance data from the MODerate resolution Imaging Spectrometer (MODIS) on board TERRA and AQUA.Such data feature a high temporal and spectral resolution (daily coverage and about 20-30 nm wide bands) along with a moderate spatial resolution (i.e., 500 m for the NIR bands).Such data are highly relevant for the surface energy balance of the snowpack.
The spectral albedo of snow is determined by the spectral and angular characteristics of the solar irradiance but also by the physical and chemical properties of the snowpack.In the visible part of the spectrum (VIS), snow albedo is mainly influenced by the impurities content.In the NIR, the effect of the snow grain size is predominent (Warren, 1982).Among numerous and sometimes ambiguous definitions of the grain size of snow (Wiscombe and Warren, 1980;Aoki et al., 2000;Fierz et al., 2009), for optical radiative transfer applications the notion of an "optical grain radius" is commonly employed.It corresponds to the radius of an optically semi-infinite collection of monodisperse disconnected spheres featuring the same NIR hemispherical reflectance.In practice, as far as hemispherical albedo is concerned, the optical radius is not significantly different from the equivalent spheres radius, r eq , i.e. the radius of a collection of spheres having the same surface to volume ratio also called specific surface area of snow (SSA).
r opt ∼ = r eq , r eq = 3 ρ i SSA , where ρ i is the density of ice (917 kg m −3 at 0 • C).In the following, the terms optical radius and equivalent spheres radius are used interchangeably.Different methods have been tested and implemented for grain size retrieval (Tedesco and Kokhanovsky, 2007).Early work by Dozier et al. (1981b) and Dozier and Marks (1987) demonstrated the potential for grain size retrieval from remote sensing.First successful demonstrations of grain size retrieval were presented in Nolin andDozier (1993, 2000) and Bourdelles and Fily (1993).Recently, Painter et al. (2009) introduced the MODSCAG algorithm to retrieve snow covered area, grain size and albedo on a sub-pixel basis from MODIS data.This method is based on the spectral unmixing of MODIS ground reflectance product (MxD09) whereby a range of snow spectra with various grain sizes are considered as potential end-members.However, although MxD09 data product addresses atmospheric effects (Vermote and Vermeulen, 1999;King et al., 2004), these methods ignore topographic effects, such as the illumination angle or the reflected terrain irradiance, despite their importance in mountainous terrain (Proy et al., 1989;Fily et al., 2000;Sirguey, 2009).In addition, the spectral endmembers representing snow used in MODSCAG are based on theoretical spectra in which snow grains are assumed to be spherical, the effect of soot on reflectance is ignored, and the effect of the anisotropy of snow reflection is addressed using an assumption of a constant viewing geometry.Recent updates of MODSCAG make it possible to evaluate the surface radiative forcing of light absorbing impurities in the snowpack (Painter et al., 2012).Alternatively, the semianalytical snow retrieval algorithm (ART) was developed by Kokhanovsky and Zege (2004) and applied to MODIS data (Tedesco and Kokhanovsky, 2007;Lyapustin et al., 2009;Negi and Kokhanovsky, 2011a,b;Zege et al., 2011).The various applications differ in the way the ART model was used, although all assume the distribution of grain shapes to be a mix of plates and columns instead of solely spherical grains (Zege et al., 2011).Each study departs regarding the number of MODIS spectral bands that were used to retrieve grain size (see Negi and Kokhanovsky, 2011a) as well as in terms of the way atmospheric corrections were addressed.For example, Tedesco and Kokhanovsky (2007) relied on the standard correction associated with the MxD09 data product, while Lyapustin et al. (2009) or Zege et al. (2011) used custom corrections. Only Zege et al. (2011) accounted for the effect of the presence of soot in the snowpack.Yet, none of these studies addressed the impact of multiple reflections in mountainous areas, and few addressed the impact of the local illumination change due to the slope and the aspect of the studied surface.
The goal of this study is twofold: first of all, it aims at evaluating the effect of (1) the local topography, (2) the anisotropy of snow and ice reflection, (3) the assumption on the shape of snow grains, to retrieve snow SSA from MODIS data in mountainous areas.The second objective is to compare the MODIS retrieved SSA with the output of the detailed snowpack model Crocus (Brun et al., 1992;Vionnet et al., 2012).Following a description of the data and models used, this study introduces in detail the various methods used to retrieve SSA from MODIS data.In the next section, MODIS retrieved SSA from the different methods are com-pared one to each other, to the field measurements and to the outputs of the detailed snow model Crocus.In addition to a detailed discussion on the comparison of the retrieval methods, this study allows for the discussion of the variability of the SSA of surface snow in mountainous terrain in both MODIS and model output data 2 Material and methods

Application site
The study was carried out in a geographical domain located in the centre of the French Alps East of Grenoble covering an area of 46 × 50 km (see Fig. 1a).The terrain is highly rugged with an elevation ranging from 224 to 3983 m a.s.l.The mean elevation is 1860 m a.s.l.

Digital elevation model (DEM)
The DEM used for computing topographic parameters originates from the Shuttle Radar Topography Mission (SRTM, Farr et al., 2007).This 3 arc-sec DEM was projected and downgraded using a cubic interpolation from 90 to 125 m (for shadows computation) and 500 m spatial resolution to match the spatial resolution of MODIS data.This downgrading implies a smoothing of the orography, and ignore subgrid local variations in slope and aspect.The relatively loose resolution of 500 m (compared to the very rough nature of the studied area) should be considered as one of the main source of uncertainty in the retrievals.Figure 1b illustrates the DEM on the studied area.SRTM data have an absolute geolocation error of 8.8 m and a relative elevation error of 8.7 m over Europe (Rodriguez et al., 2005).

Satellite data
The main features of MODIS data making it an appealing source of information for remote sensing of snow from space are outlined below.First of all, it has a suitable spectral resolution (see Table 1), including a detection channel (band 5) centred at 1.24 µm, which is highly sensitive to SSA and barely affected by light absorbing impurities (Warren, 1982;Kokhanovsky et al., 2011).It offers a good compromise between spatial and spectral resolution (respectively 500 m and 0.02 µm for band 5).Last, it provides at least daily global coverage since 2000, which makes it a valuable dataset for multi-year studies (e.g.Box et al., 2012) and, despite a primary goal of being a research-oriented data source, many www.the-cryosphere.net/7/741/2013/The Cryosphere, 7, 741-761, 2013 operational applications makes increasing use of these data (e.g.Pettinato et al., 2009).
The region of interest is represented by 92 × 100 pixels at 500 m resolution.Note that 500 m is the nominal resolution, which is effective at nadir of the sensor, but data with large sensor zenith angle (on the side of the swath) are affected by a loss of effective pixels resolution.For the studied dates, the absolute sensor zenith angle ranged from 2 to 27 • , with less than a 4 • variation within each image.
The Goddard Space Flight Center (MODIS project, http: //modis.gsfc.nasa.gov/)provides several products processed from raw sensor acquisitions.For this study, we used two of these products in order to evaluate the impact of the topographic correction on the retrieved SSA: the MOD02 swath data product (level 1B, obtained with MOD03 geolocation file for reprojection), which contains calibrated top of atmosphere (TOA) radiances; and the MOD09 data product (level 2), which directly provides atmospherically corrected ground reflectances computed by the LSRSCF (Land Surface Reflectance Science Computing Facility: http://modis-sr.ltdri.org/) (Vermote and Vermeulen, 1999).

Numerical simulation of snow reflectance
Throughout the study, numerical simulations of snow reflectance were carried out using the DIScrete Ordinate Radiative Transfer (DISORT) model (Stamnes et al., 1988).We assumed grains to be spherical and disconnected and computed their single scattering parameters using the Mie theory.The impact of the grain shape on the reflectance can be significant (Picard et al., 2009).However, Carmagnola et al. (2012) showed that using field measurements of SSA and density together with DISORT using the same assumption leads to NIR spectral albedo simulation in very good agreement with co-located field measurements.
Besides input from near-surface vertical profiles of the physical properties of snow, DISORT requires the knowledge of a spectrally resolved optical index of ice.Recent work (Carmagnola et al., 2012) has shown that a recent compilation of the ice refractive index (Warren and Brandt, 2008) may not be the more relevant for certain near-infrared wavelengths.However, in the spectral range concerned with MODIS bands (Table 1), the agreement between observations and numerical simulations was found extremely satisfactory, so that values from Warren and Brandt (2008) are used here.Note that the relevant database is available at http://www.atmos.washington.edu/iceoptical constants/.

SAFRAN-Crocus raw data
Numerical simulations of the physical properties of snow on the ground were carried out using meteorological data from the SAFRAN downscaling system (Durand et al., 1993) and the snowpack scheme Crocus coupled to the land surface scheme ISBA (Boone et al., 2000;Decharme et al., 2011) within the SURFEX interface (Masson et al., 2012).
The meteorological forcing provided by SAFRAN includes air temperature and humidity, wind speed, cloud cover, precipitation, and downward longwave and shortwave direct and diffuse incident radiation.Assessing such meteorological conditions on mountain slopes is challenging due to the scarcity of ground-based measurements and the coarse resolution of numerical weather prediction models (NWP).The SAFRAN meteorological downscaling system was designed to overcome this issue and provide relevant meteorological forcing, combining remotely sensed cloud cover information and ground-based and radiosondes observations with an a priori estimate of meteorological conditions obtained from the ARPEGE NWP model (Courtier et al., 1991).The downscaling is performed within geographical areas (ca.400 km 2 ) assumed to be meteorologically homogeneous and referred to as "massifs".In each massif, the surface analysis is carried out on a vertical profile and provides meteorological data typically on a vertical grid of 300 m.In the case of our simulations and for each altitude within a given massif, this meteorological forcing was provided to the land surface model including modifications of the direct component of the solar radiation accounting for five slope angle values (0, 10, 20, 30 and 40 • ) and for eight aspects (N, NE, E, SE, S, SW, W, NW).Outputs of the land surface model runs thus depend on time and on four variables defining a topographic class: geographical zone ("massif"), altitude, aspect, and slope.Table 2 details the discretisation of these variables in our study area.The classification chosen is arbitrary, and we note that more elaborated clustering of topographical parameters could also be employed such as recently proposed by (Fiddes and Gruber, 2012).The three massifs included in our study area are called Belledonne, Grandes-Rousses, and Oisans and are indicated in Fig. 1b.
Crocus is a 1-D detailed snowpack model simulating the energy and mass balance of the snowpack including a detailed description of internal processes such as snow settling, liquid water percolation and snow metamorphism (Brun et al., 1989(Brun et al., , 1992;;Vionnet et al., 2012).For each of the layers of its modelled stratigraphy, it simulates the evolution of grain characteristics, thickness, density, liquid water content and temperature.Numerical simulations were initialized with no snow on the ground on 2000-08-01 using an initial ground temperature corresponding to the mean air temperature (SAFRAN) for the period 2000-08-01-2010-07-01 for each simulation point.This procedure ensures that ground thermal conditions are initialized realistically, and the 9 yr spin-up time leads to consistent numerical simulations of ground-snow interactions.Crocus variables describing snow microstructure, namely sphericity, dendricity, and size, are semi-quantitative and depend on the snow metamorphism history.The sphericity represents the ratio of rounded grains to faceted grains.The dendricity is equal to 1.0 for fresh snow and then decreases to zero as the snow ages.From both grain shape variables, the Crocus model diagnoses an optical diameter (Brun et al., 1992;Vionnet et al., 2012;Morin et al., 2012) that we converted to an SSA value according to Eq. (1).Morin et al. (2012) demonstrated that they were in very good agreement with field measurements carried out at the research station Col de Porte over the course of one snow year (2009)(2010), with a root mean square deviation (RMSD) found to be on the order of 6 m 2 kg −1 .The snowpack model Crocus has been successfully evaluated at several sites and under different climatic conditions (e.g.Durand et al., 1999;Bouilloud and Martin, 2006;Wagnon et al., 2009;Brun et al., 2012).

Post-processing of SAFRAN-Crocus data and comparison with MODIS-retrieved data
To compare the SSA values simulated by SAFRAN-Crocus with those retrieved from MODIS data, we computed a ver- tical average of the SSA of the surface layers using Crocus outputs relevant to the the first centimeters of the snowpack.This average was performed to a single and appropriate surface SSA value to be compared to MODIS retrieved data.This was done considering an exponential decay to accomodate the attenuation of the solar radiation through the snowpack (Warren, 1982) as where s is a truncation of the penetration depth, and d denotes the e-folding depth.The latter varies with SSA, density, and wavelength (Warren, 1982).In order to simplify the processing steps of SAFRAN-Crocus outputs, we used a single value of d = 2 cm and s = 4 cm for the whole study.Indeed, the sensitivity study conducted by Kokhanovsky et al. (2011) revealed that NIR wavelengths are only sensitive to the optical radius on the very first centimeters of the snowpack.3 gives an illustration of the variation of the e-folding depth for two extreme cases (fresh light snow and aged dense snow) and for each of the MODIS bands used in this study.These values have been calculated with DISORT and assuming that grains are disconnected spheres.Note that, as noticed in Sergent et al. (1987), the e-folding depth computed using this assumption may be higher than for real snow.Following typical extremes values for these two cases, SAFRAN-Crocus SSA data were processed for the six studied dates using an e-folding depth of 1 mm and 40 mm.The root mean square deviation between all SSA values obtained in the two cases range from 0.7 to 8.0 m2 kg −1 depending on the date (mean 4.24 m 2 kg −1 ).These variations are sufficiently small to justify the use of a single value for d in this study.Further developments may consider using wavelength and perhaps even SSA and density specific d values for an even more detailed computation of a representative surface SSA value corresponding to the remotely-sensed variable.
According to the minimum and maximum values presented in Table 3, the surface SSA was processed for the six studied dates from the outputs of SAFRAN-Crocus using an e-folding depth of 1 mm and 40 mm respectively.The standard deviation of the SSA processed in these two cases is small (3.53 m 2 kg −1 ).This shows that a single value of the penetration depth can be used throughout the study to compute a relevant value for surface snow SSA from Crocus outputs to be compared to MODIS-derived data regardless the spectral bands upon which the retrieved values are computed.
Outputs of the SAFRAN-Crocus model corresponding to the six acquisition dates of MODIS data were processed using the procedure indicated above.The simulated SSA values for each topographic class were distributed throughout the area of interest using elevation, slope and aspect computed from the DEM.The value of the nearest topographic class was simply assigned to each pixel, without any interpolation of the simulated data.

Field measurements
To complement the comparison between MODIS-retrieved data and SAFRAN-Crocus outputs and introduce a slight dose of ground truthing in our analysis, field measurements of SSA were performed in winter 2012.Two field sites featuring a maximum homogeneity in terms of slope, aspect, absence of shading over a characteristic horizontal scale on the order of MODIS pixel sizes were selected for this purpose within the Belledonne massif.One of them is a flat terrain at 2000 m .a.s.l. and the second one is a south-westerly facing slope tilted around 15 • , at 2100 m a.s.l.Near-surface vertical profiles of snow SSA and density were measured using the optical instrument DUFISSS (Gallet et al., 2009).The estimated accuracy of each SSA measurement is ±10 %.Repeated measurements within the same areas were performed with the intent to quantify the in situ variability of surface snow SSA.

Detailed description of the satellite retrieval methods
Table 4 presents an overview of the acronyms used to refer to the different retrieval methods which are described in detail below.

Atmospheric and topographic correction: computation of ground reflectance from MOD02 radiance data
The radiance detected by a satellite sensor (L TOA ) at a given wavelength 1 is the sum of several contributions including (Fig. 2): the ground radiance L g ↑ (i.e. the radiance of the ground target) modulated by the ground-to-sensor atmospheric transmittance T v , the path radiance L p (i.e. the radiance emitted by the column of atmosphere between the sensor and the target) and the background radiance L k (i.e. the radiance emitted by pixels surrounding the target and diffused towards the sensor) (Sirguey et al., 2009).Consequently, the ground radiance in the direction 2 defined by the zenith angle θv and the azimuth φv can be written as (Dumont et al., 2011) where L g is the incident radiance received by the target, ρ is the bidirectionnal reflectance distribution function (BRDF), and (θ i , φ i ) the zenith and azimuth angles of the illumination, respectively.The radiance L g can be written as the sum of several components (Fig. 2) related to the direct solar irradiance E s , the diffuse solar irradiance E d , the diffuse environmental irradiance E m , and the reflected terrain irradiance E t .
For simplicity, we assume that the diffuse irradiances E d and E m are isotropic (Schaepman-Strub et al., 2006), as well as E t , although it has been shown that anisotropic treatment of the diffuse irradiance may be more appropriate (Wang et al., 2012).E m and E t differ in the size of the neighbourhood considered for their computation, 500 m of radius for E t and 1 km for E m .As such, all three components can be accounted for as a single diffuse irradiance component, where θs is the solar illumination angle on the tilted target defined as cos( θs where θ s is the solar zenith, φ s is the solar azimuth, and (θ n , φ n ) are the slope and aspect of the ground tilted pixel, respectively.The isotropic assumption of E diff allows Eq. ( 4) to be written as where δ 2 θs ,0 Combining Eqs. ( 3) and ( 5) yields where α is the hemispherical-directional reflectance in the direction θv .The Helmoltz reciprocity principle ρ(θ, θv , φ − φv ) = ρ( θv , θ, φv − φ) enables us to inverse Eq. ( 6) so that where R, the anisotropy factor, is defined as the ratio of the BRDF to the hemispherical-directional reflectance (see Dumont et al., 2010).
While allowing the topographic effects to be addressed, this formulation also enables the SSA to be retrieved under two scenarios: (1) snow is a Lambertian surface (R = 1); ( 2) snow has a marked anisotropic reflectance function measured by the anisotropy factor R.
Under the first scenario, Eq. ( 7) corresponds to Eq. ( 5) in Sirguey et al. (2009): where b varies between 0 if the pixel is considered shaded to 1 if it is entirely illuminated.T v and T s are, respectively, the ground-to-sensor and sun-to-ground atmospheric transmittances.d is the Earth-Sun distance, and E 0 the extraterrestrial irradiance, both provided within MOD03 products.The hemispherical-directional reflectance, α, is computed using Eq. ( 8) with the MODImLab algorithm (Sirguey et al., 2009;Dumont et al., 2011).It calculates E t , E m and α using an iterative method.The quantities E d , E s , T v , T s , E m , L p , and L k are obtained using the SPCTRAL2 radiative transfer model (Bird and Riordan, 1986).The shadowing factor b is computed within MODImLab using the 125 m DEM based on an implementation of the horizon line algorithm of Dozier et al. (1981a).
Under the second scenario, we retrieved SSA using a correction of the anisotropic reflectance of snow.The anisotropy factor is defined as R(θ i , φ, θ v ) = π ρ(θ i , φ, θ v )/α(θ i ).Mean anisotropy factors inferred from measurements in a cold room on three types of snow and under three different illumination angles by Dumont et al. (2010) were used in Eq. (7).

Computation of ground reflectance in MOD09 products
Ground reflectance provided in the MOD09GA V005 data product are computed using the atmospheric radiative www.the-cryosphere.net/7/741/2013/The Cryosphere, 7, 741-761, 2013 transfer model 6S (Vermote et al., 1997;Vermote and Vermeulen, 1999).Although anisotropic surface reflectance behaviour are taken into account in the MOD09 product, the effects of topography are ignored in this product, i.e. b = 1, E m = E t = L k = 0, and θs is replaced by θ s .

From reflectance to SSA
Maps of ground reflectance were obtained with MODImLab from MOD02 data (corrected for topographic effects, and corrected or not for anisotropic reflectance effects), or readily available from MOD09 data products.Two general classes of methods were then employed to convert spectral reflectances into surface SSA values.In one case, the spectral signature of snow was converted based on a look-up table built from a wide range of DISORT runs.Alternatively, the ART theory was employed.The two methods are described in detail below.An overview of the methods and options is provided in Table 4.

Retrieval methods based on look-up tables built using DISORT
The retrieval of SSA values was only carried out for pixels which are both non-shaded and identified as fully covered with snow (snow cover fraction = 1).We ignored shaded pixels because of the fact that limited signal to noise ratio introduced large uncertainties.Considering only fully snowcovered pixels further permitted to ignore the compounded effects arising from the mixture of spectral responses and bi-directional reflectance of different targets.The sub-pixel snow cover fraction was determined based on its linear relationship with the Normalized Difference Snow Index (NDSI) after the atmospheric and topographic corrections (Salomonson and Appel, 2006;Hall and Riggs, 2007).Thus, all nonshaded pixels whose NDSI was larger than 0.7 were considered to be fully covered.Although snow fraction was found to be variable for a given NDSI by Salomonson and Appel (2006), the use of a relatively high threshold and the fact that NDSI was computed from ground reflectance corrected for atmospheric and topographic effects were considered a conservative means to retain fully covered snow pixels.We used the DISORT model to establish a relationship between SSA and reflectance (Stamnes et al., 1988).We computed reflectance values for a set of incident angles ranging from 2 to 88 • by steps of 2 • , SSA values ranging from 2 to 160 m 2 kg −1 by steps of 1 m 2 kg −1 , and wavelengths corresponding to MODIS bands 1 to 7. DISORT calculations for high zenith angle are highly uncertain but these high zenith angles concern very few pixels of the MODIS data.Average viewing zenith angles vary between 23 and 32 • depending on the date.We ignored the presence of impurities in the snow because they are not expected to affect substantially the reflectance in the spectral bands that are most sensitive to SSA (Warren, 1982).A look-up table (LUT) was built to link SSA values to the spectral reflectance for each MODIS band as a function of the viewing angle.
In order to retrieve SSA from spectral reflectance values, we used several processing options: 1. using either MODIS band 5 only or MODIS bands 2, 5, 6, and 7 together although band 2 is sensitive to impurities; 2. relying on absolute reflectance values or on the relative shape of the snow spectrum (i.e. the ratio between the bands considered and band 4).This ratioing enables to account for the relative shape of the snow spectrum instead of absolute reflectance which could be affected by, for example, atmospheric perturbations or local shade.
Similarly to Nolin and Dozier (1993), for each pixel to be processed, the algorithm searches the LUT and selects the SSA whose computed reflectance spectrum resembles most that of the image according to the spectral distance D defined as where α M is the observed reflectance vector from MODIS and α D is the reflectance vector from DISORT.The parameters η(i) are weighting coefficients quantifying the sensitivity of each band to SSA.They were established by evaluating the quantity ρ D (SSA = 5, b) − ρ D (SSA = 160, b).This yielded the following values η(2) = 0.2, η(5) = 0.7, η(6) = 0.05, and η(7) = 0.05 when using four bands, or simply η(5) = 1 when only band 5 was used.The SSA value for a pixel was not assigned when D min was greater than a predefined threshold, denoted D 0 .The value of D 0 is discussed in Sect.4.3.1.

Retrieval method using non-spherical grains (ART)
One alternative method for obtaining SSA was considered that takes into account the grain shape in a different method than the previously described method.Indeed, in the ART theory (Kokhanovsky and Zege, 2004), snow grains are considered to be fractal rather than spherical as described in Negi and Kokhanovsky (2011a,b).The fractal particle model was introduced by Macke et al. (1996): the ice grain is modelled as a tetrahedron with small tetrahedrons attached to each plane of the initial tetrahedron.This model was found to be capable of describing the snow BRDF (Kokhanovksy et al., 2005).The method relies on one visible (band 4) and one NIR (band 5) channels.It accounts for the effect of slope on the incidence and viewing angles, as well as for the anisotropy of snow reflection.The method was applied only to MOD09 products for pixel whose NDSI exceeded 0.6, reflectance on band 4 exceeded 0.6 and incidence angle was lower than 70 • .This method is primarily to be compared with the DA method (see Table 4) which uses the same products, no topographic correction but uses the DISORT LUT-based retrieval method.

Results and discussion
According to Eq. ( 1), the SSA is inversely proportional to the optical radius.Studies reporting on the retrieval of snow properties from satellite data often refer to the optical radius (Fily et al., 1999;Painter et al., 2009;Negi and Kokhanovsky, 2011b), while some others tend to refer to SSA (Dumont et al., 2012).It is worth noting that the use of SSA is more appropriate when studying small grains with high albedo.Indeed, in terms of snow albedo, a given variation of SSA in the higher range of SSA values has a stronger impact in terms of snow albedo than a similar variation in the lower range of snow SSA values.Therefore, as long as radiative transfer and surface albedo are concerned, statistics computed based on SSA are a more appropriate metric than optical grain size (Morin et al., 2012).Nevertheless, here the results of this study are reported in terms of both SSA and optical radius.

Comparison to field measurements
Table 5 shows an overview of the comparison between simultaneous estimates the SSA of surface snow using in-situ measurements, MODIS-derived estimates and SAFRAN-Crocus simulated values.It shows that MODIS DTA SSA retrievals and SAFRAN-Crocus data are the closest to field measurements.SSA MODIS data retrieved without topographic corrections (DA and DR) exhibits significant biases compared to field measurements.The measured intrapixel variability is up to 5.1 m 2 kg −1 .This underlines the need for making numerous field measurements inside the pixel area to be as representative as possible of the snow properties within each MODIS pixel.For the two dates considered, the standard deviation of MODIS-derived data within a topographic class is low (2.4 and 1.7 m 2 kg −1 ).This indicates that on these dates, the semi-distributed approach followed by SAFRAN-Crocus seems sufficient to represent the variability of the SSA of surface snow.
Although it may seem satisfactory that the ground measurements are in good agreement with SAFRAN-Crocus output and one of the MODIS retrieval methods (DTA), we emphasize that this attempt to use ground measurements as a validation of the satellite retrieval methods and numerical simulation outputs is considered preliminary and should be complemented by many more ground-based in situ measurements in the future.Figure 4 shows the colour composite of MODIS images on the six dates of the study.Figure 5 illustrates the SSA obtained from DTA retrievals for all dates.

Comparison of the retrieval methods and options
In this section, we compare the values obtained from the various retrieval methods tested here using different options to SAFRAN-Crocus outputs.Figure 6 shows the SSA maps obtained with SAFRAN-Crocus and the DTA method for 2009-03-21.Tables 6 and 7 summarize the SSA and optical radius results obtained with each method on the six dates.We used four statistics, applied to both SSA and optical radius, to report on and assess the performance of each method: 1. the mean and standard deviation of the variable over the whole area, 2. the root mean square deviation (RMSD) and mean value of the pixel-to-pixel difference between each of the tested methods and SAFRAN-Crocus estimates (see for example  All statistics were computed based on pixels where all retrieving methods successfully provided an estimate of snow SSA.
Although SAFRAN-Crocus outputs can by no means be considered ground truth for the evaluation of the different retrieval methods employed, it was used here as a common benchmark to compare the different methods.It was considered that the overall realism of the meteorological driving data along with the detailed numerical simulations of the physical properties of the snowpack, evaluated independently in several previous studies, provides a physically consistent view of the general features of surface snow SSA in the domain and the time span considered.Whether the output of a given retrieval method is in agreement with the general features of the SAFRAN-Crocus output is hypothesized here to mean that the large scale features of the physical properties of the uppermost snow layers are well represented in the MODIS-derived surface snow SSA estimate.Figure 7 illustrates how the mean SSA and optical radius values compare with SAFRAN-Crocus estimates for each method.Figure 8 summarizes the corresponding RMSD of all methods for both SSA and optical radius.

Influence of selected spectral bands
Here we describe the impact of using either only one (band 5) or four (bands 2, 5, 6 and 7, see Table 1) MODIS spectral bands in the retrieval algorithm using DISORT-based LUTs.
In the case of the DTA method, we observed that when using only band 5 (1240 nm) the RMSD slightly increases with D 0 .In contrast, the performance of the method using four bands 2, 5, 6, 7 actually depends on the value of D 0 chosen.For a D 0 smaller than 0.2, the number of pixels that can be used falls and the mean difference strongly increases.On the other hand, the RMSD varies between 1 and 40 m 2 kg −1 for a threshold lower than 0.15 and then increases regularly with D 0 .An optimal value of D 0 minimizing both RMSD and mean difference in average would be about 0.14.With this threshold, the method using four bands provides data for about 87% of the pixels; in this case, the RMSD between the MODIS-derived and SAFRAN-Crocus values were approximately the same whether four or one band were used, with 14.0 and 14.8 m 2 kg −1 , respectively.The overall SSA mean difference on the subset was reduced from 1.7 m 2 kg −1 when using band 5 alone to −0.1 m 2 kg −1 with four bands.
Using four bands generally yields smaller SSA values: this result is consistent with Li et al. (2001) who reported the fact that the SSA usually decreases with depth in the first centimeters, consistent with the generally observed decrease of SSA with time (Domine et al., 2008), and the fact that band 2 (860 nm) penetrates deeper in the snowpack, as shown in Table 3.In the following, due to the fact that the obtained results are not significantly different, only computations using only Table 5. Measured, MODIS retrieved and SAFRAN-Crocus simulations of SSA (m 2 kg −1 ) on two different locations of the Belledonne massif during winter 2012.For field measurements, σ and n are the standard deviation and number of measurements carried out in each pixel.DTA, DA and DR refer to the MODIS-derived SSA at the location of measurements.MODIS DTA Topographic Avg is the mean SSA within the retrieved pixels of the Belledonne massif belonging to the same topographic class (massif, elevation, slope, aspect).In this case, σ is the standard deviation within these pixels and n their number.band 5 with a D 0 value of 0.02 are reported and discussed here.

Influence of topographic correction
The retrieval methods ignoring topographic correction (i.e., DA and DR) yielded higher SSA (smaller optical radius) estimates compared to DTA and SAFRAN-Crocus.This was largely marked for dates with high SSA/small optical radius (Table 7).The biases and RMSD in terms of SSA and optical radius were also significantly higher for DA and DR than for any other method (Table 7).
Figure 7a shows that DTA and SAFRAN-Crocus were in good agreement for relatively small SSA values, while higher SSA exhibited more dispersion.The overall mean difference associated with the DTA method over the six dates was 0.1 m 2 kg −1 , with the largest mean difference being −12.0 m 2 kg −1 ≡ 30 % relative mean difference (Table 7) on 2009-01-09.The RMSD computed for each date ranged from 2.9 to 25.9 m 2 kg −1 with lower RMSD associated with lower SSA.The average RMSD for the DTA method over the six dates was the smallest of all methods with 9.4 m 2 kg −1 .However, over the six dates 54 % of pixels had an absolute difference lower than 5 m 2 kg −1 , so that the RMSD is largely due to a relatively small number of pixels showing a large deviation to SAFRAN-Crocus values.In terms of optical radius estimates, the DTA method was the second in best consistency with SAFRAN-Crocus after the ART method with an overall mean difference of −7.8 µm and RMSD of 129.6 µm (see Table 7).
Therefore, this comparison suggested that the correction of topographic effects in the reflectance calculation results in smaller SSA (higher optical radius) with better agreement to those simulated by SAFRAN-Crocus.
This observation can be explored analytically by studying the sensitivity of Eq. ( 8) to topographic effects (TE).TE is a variable quantifying the amount of local topographic effects taken into account.It increases with the number of topographic corrections used, e.g.slope, multiple reflections, etc.The reflectance variation with regards to TE can be writ-ten with partial derivatives as follow: The shadowing factor b, although included in TE, is omitted in this equation since only non-shaded pixels (i.e.b = 1) were considered in our retrievals (cf.Sect.3.3.1).The signs of ∂α/∂L k and ∂α/∂E d are known from Eq. ( 8).L k / TE, E t / TE and E m / TE are positive since E m , E t and L k are equal to 0 without topographic correction and positive otherwise.E d / TE has a negative sign since the topographic correction restricts E d to the part of the atmosphere above the horizon line (Sirguey et al., 2009).E t / TE represents the irradiance directly reflected by the surrounding pixels in the portion of the atmosphere below the horizon line.In a first approximation and given the relatively high reflectance of snow, we can assume that E d is decreasing less than E t is increasing when taking into account topographic effects.Under this assumption [ E t / TE + E t / TE] has a positive sign.∂α/∂ θs has a positive sign since the reflectance is increasing with the solar zenith angle.Finally, θs / TE represents the variation from θ s to θs .Consequently, it only depends on the slope and aspect of the ground pixel, according to the definition of θs in section 3.1.Thus, its sign is pixel-dependent and is determined by the interaction between the sun, the topography and the surface reflectance.Nevertheless, Löwe and Helbig (2012) showed in the simplified case of a Lambertian BRDF, atmospheric effects being neglected and with single terrain reflections, that the higher the subgrid orography, the lower the effective albedo.This yields that α/ TE has probably a negative sign in a first approximation, and this confirms that including topographic effects in reflectance calculations results in a lower apparent reflectance, which is converted to a smaller SSA by the retrieval algorithms.In addition, E m , E t , and L k increase with reflectance and SSA, and hence so do | α/ TE|.The impact of correcting topographic effects is therefore more pronounced on dates with relatively higher SSA.

Influence of normalization
The use of band ratios (as defined in Sect.3.3) instead of absolute reflectance values proved to have an impact on the retrieved optical radius and SSA.When considering topographically corrected data only, the use of band ratios yielded generally higher SSA.SSA mean values were twice higher  75) 439 ( 106) 452 ( 108) 424 ( 185) 263 ( 173) 503 ( 130) (and the other way around for the optical radius) and largely departed from SAFRAN-Crocus estimates.
Band ratioing involves a division by the reflectance in a visible band (typically high for snow).In fact, reflectances at visible wavelengths are influenced by other factors, such as impurities, which can result in reflectances lower than simulated ones.This underestimation of visible reflectances induces overestimated ratios.This leads to an overestimation of SSA values when using band ratioing.
When considering data retrieved without taking into account topographic effects, the impact of band ratios differ whether looking at SSA or optical radius.SSA values still show higher SSA for band ratioing (DR) than for absolute reflectance method (DA) on dates with small mean SSA, but smaller SSA on dates with higher SSA values (Table 6).Optical radius values are always lower when using the band ratio method, except on 2010-01-12 where both methods give very similar values (Table 6).
The use of band ratio is a rather simple approach which can be used to overcome the error due to ignoring topographic effects on absolute reflectance.Indeed, since all spectral bands are similarly affected by local illumination angle, using a band ratio can be viewed as a form of topographic normalization.However, the amplitude of topographic effects remains dependent on wavelength due to the importance of other topographic contributions such as L k , E m , and E t .We demonstrated in Sect.4.3.2 that α/ TE < 0. Calling α the MOD09 reflectance (i.e., without topographic correction) and α * the topographically corrected reflectance, it comes that α = α * + δα, with δα > 0. This overestimated reflectance would lead to overestimates of SSA as confirmed by higher mean difference values for DTA than for DA in Table 7.Alternatively, comparing the DR reflectance ratio to the theoretical reflectance ratio can be formulated as where α nir and α v , respectively, represent NIR and visible reflectances.Since visible wavelengths are barely sensitive to SSA, α * v and δα v can be assumed independent of SSA.In contrast, α nir varies with more than a factor three when SSA varies from 5 to 60 m 2 kg −1 .Over the same range, δα nir varies as a first order of α nir with SSA.Thus, the second term α * nir δα v of the numerator has larger variance than the first term, making > 0 when SSA and α v are small, and < 0 when SSA and α v are large.This explains why DR overestimates SSA when the SSA is small, and underestimates it when the SSA is high.This effect, together with overestimations due to impurities can explain the sign of the difference between DA and SAFRAN-Crocus, and DR and SAFRAN-Crocus, respectively.

Influence of grain shape and anisotropy
Figure 7a shows that DTAA method retrievals leads to higher mean SSA on dates with high average SSA values (about +50 %), whereas it has very little impact on dates with lower mean SSA values.The mean difference and RMSD to SAFRAN-Crocus are higher for DTAA than for DTA.This is consistent with the fact that this correction has very little impact in terms of optical radius (Figs.7b and 8).The anisotropy correction implemented here is quite simple, as it uses the anisotropy factor inferred from measurements by Dumont et al. (2010) for three types of snow, featuring relatively low SSA.This correction seems to be too strong for snow with a high SSA, whose anisotropy is generally less pronounced than small SSA snow (Dumont et al., 2010).This www.the-cryosphere.net/7/741/2013/The Cryosphere, 7, 741-761, 2013 is also the reason why the correction has a significant effect in SSA but not in optical radius.
As for DTAA, the ART method also gives SSA mean values significantly higher than SAFRAN-Crocus on dates with high SSA values, and values close to the model results on dates with low mean SSA values (Fig. 7a).The ART retrieved optical radius are the closest to SAFRAN-Crocus estimates.Compared to the outputs of the DA method, the mean difference and RMSD are reduced (Table 7).This indicates that the impact of taking into account the grain shape is as significant as the impact of addressing multiple reflections on the retrieved SSA values at least for low SSA values.

Mixed pixels and surface roughness
In mountainous terrain, 500 m × 500 m pixels fully covered by snow are sparse.The presence of rock or vegetation modifies the directional reflectance distribution of the pixel which in this case largely departs from the snow anisotropy factor or Lambertian assumption applied here.To address this problem, spectral unmixing including snow, soil, vegetation and rock spectra could be applied as described in Painter et al. (2009) and the anisotropy corrections would be applied only on the part of the pixel reflectance attributed to snow.
In addition, it is likely that within a pixel area the surface of the snowpack largely differs from a flat surface due for example to wind drifted snow.Surface roughness also induces changes in BRDF (e.g.Leroux and Fily, 1998;Hudson and Warren, 2007;Zhuravleva and Kokhanovsky, 2011;Kuchiki et al., 2011).This point is not taken into account in our retrieval and consequently, on highly rugged pix-els the accuracy of retrieved SSA is likely to decrease.Consequently, application of the algorithm in polar regions where sastrugi are numerous should be treated with caution.Using data from other sensors (e.g.multi-angle sensors (Nolin and Payne, 2007) or laser altimeter) could help in the future to better characterizing surface roughness and improve the algorithm.
Mixed pixels and surface roughness can thus explain the fact that the simple anisotropy correction applied in this study leads to potentially overestimated SSA values when using the DTAA method.

Detailed comparison of the DTA retrieval method with SAFRAN-Crocus model outputs
In this section we compare SSA values obtained from DTA retrievals (MOD02 band 5 topographically corrected radiances, using absolute reflectance) for comparison with estimates from SAFRAN-Crocus.As shown in Sect.4.3 this method is the most consistent with of SAFRAN-Crocus surface SSA values.The following analysis takes into account pixels for which estimates existed both for DTA and SAFRAN-Crocus.Many more pixels were available for this method (more than 2000 on average) than in the subset for which all methods overlapped, thus making the statistics reported below slightly different than those presented earlier.
In the Crocus model, the representation of snow microstructure leads to SSA values bounded between 1 and 65 m 2 kg −1 .The SSA retrieved from MODIS data can range between 2 and 160 m 2 kg −1 .This discrepancy between the two SSA ranges may introduce an asymmetry in their statistics, despite the limited number of pixels for which the SSA was estimated to exceed 65 m 2 kg −1 based on MODIS retrieved data.

Intraclass SSA variability
Figure 9 shows the standard deviations of SSA values from SAFRAN-Crocus and MODIS data on the six dates.The standard deviation is systematically higher for satelliteretrieved data.Maps of SSA illustrated in Fig. 6 also illustrate the greater dispersion of DTA values.Both data sources have very different signal entropy, i.e. even on the same number of pixels, the two datasets feature different degrees of freedom.Indeed, while SAFRAN-Crocus provides estimates based on a discrete subset of topographic parameters (classes defined by massif, elevation, aspect and slope) representing 891 different classes altogether (cf.Table 2), each MODIS pixel is considered independant from each other.
In fact, the mean standard deviation for DTA within classes represented by at least 10 pixels is 10.4 m 2 kg −1 (whereas this standard deviation is, by definition, zero for SAFRAN-Crocus).We computed "smoothed" SSA maps, where the SSA on each pixel was replaced by the mean SSA of all pixels corresponding to the same SAFRAN-Crocus topographic class.The standard deviation of the resulting map was significantly reduced as shown in Fig. 9.It demonstrates that most of the difference of variability between SAFRAN-Crocus and DTA originates from the inherent variability of DTA estimates within the same SAFRAN-Crocus class, which could not be captured by the model.Based on the "smoothed" SSA map, the RMSD between DTA and SAFRAN-Crocus decreased from 9.4 to 7.0 m 2 kg −1 .This suggests that about 25 % of the RMSD between MODIS and SAFRAN-Crocus SSA values originates from the variability of the SSA retrieved from satellite measurements for a given constant altitude, slope and aspect.This intra-classes variability in the MODIS retrieved values may stem from local topographic effects such as sun direct light shadowing and wind drift, from mixed pixels, and also from surface roughness.
Figure 10 represents the standard deviation of the difference bewteen DTA and SAFRAN-Crocus SSA values within each topographic class of the SAFRAN-Crocus model represented by more than 10 pixels.For this figure we used the configuration of the model described in Sect.2.5 but also a configuration with a reduced number of classes (only 2 slopes and flat and 6 aspects).This corresponds to the operational configuration of the SAFRAN-Crocus-MEPRA model chain used for avalanche warning activities at Météo-France (Durand et al., 1999).In this case, 287 topographic classes are represented by more than 10 pixels, whereas 423 are present in the configuration used throughout this article.The standard deviations are higher for the reduced number of topographic classes (median 4.9 and mean 8.0 m 2 kg −1 ) than for the usual configuration of the model (median 3.9 and mean 6.0 m 2 kg −1 ).This indicates that the variability of   retrieved SSA is relatively low for the studied dates when using a larger number of topographic classes.Thus the use of a semi-distributed approach seems reasonable to represent the variability of the snowpack consistent with the scale of the MODIS pixels.

Interclass SSA variability
We analyse here the variability of SAFRAN-Crocus and DTA retrievals by means of a multiple linear regression of log(SSA) over topographic parameters and date.We used log(SSA) rather than SSA to overcome limitations of the linear regression technique due to the asymmetry of SSA distributions between the two datasets compared.The predictors used are date, elevation, slope, aspect, and massif.We also took into account predictor interactions, because topog- raphy effects tend to depend strongly on the date.In addition, effects related to the position of the sun are related to a combination of slope and aspect.For DTA, standard regression techniques (not shown here) show the significance of all predictors, as well as of the interactions between date and massif, and date, slope and aspect.This means that the influence of elevation, massif, slope and aspect varies significantly from date to date.All listed factors, variables and interactions are found significant at the 0.05 confidence level.For SAFRAN-Crocus, all predictors are found significant, except slope whose influence appears only through interactions.The percentage of variance explained by these statistical models reaches 68 % for DTA with 11 957 degrees of freedom.It raises to 88 % for "smoothed" DTA, which shows that the signal retrieved from the sensor is deeply linked to topographic parameters once the inter-pixel variability has been removed within a given class.For SAFRAN-Crocus,   the percentage of variance explained reaches 93 %, consistent with the fact that SAFRAN-Crocus simulations are explicitly linked to these predictors in a physically deterministic manner.
The influence of date is predominant on two items: the mean SSA, as Table 6 shows, and the magnitude of topographic gradients of SSA.Indeed, topographic parameters have a variable correlation with SSA depending on the date.This result is consistent with the meteorological conditions: just after a snow fall, or late in the season, the snowpack becomes more homogeneous so that gradients tend to vanish.
DTA and SAFRAN-Crocus SSA data show a significant increase of SSA with elevation (Fig. 11).The magnitude of the gradient varies from date to date, with smaller altitude gradients shortly after a snowfall or in the season.This result is consistent with the general understanding of the physical properties of snow in mountainous terrain (Armstrong and Brun, 2008).We note that the altitude gradient is often more pronounced for DTA retrievals than for SAFRAN-Crocus estimates.For the lowest elevations and lowest SSA values, we can link this observation to the fact that SAFRAN-Crocus seems to underestimate the rate and magnitude of SSA decrease in aged snow (Morin et al., 2012).For the highest elevations and the highest SSA values, the difference of SSA definition range between Crocus and DTA, and especially the upper limit at 65 m 2 kg −1 for  the SSA in Crocus, can explain the lower SSA observed for SAFRAN-Crocus.
Figure 12 shows SSA from DTA and SAFRAN-Crocus with regards to aspect for the six studied dates.It illustrates both for SAFRAN-Crocus and DTA the variations for each date upon the general influence of the aspect on SSA.The highest SSA are located on northern slopes, then western, eastern, and the smallest SSA on southern faces.This result is consistent with the amount of solar radiation received by these slopes.Again, it is noticeable on Fig. 12 that the variability due to the aspect is more pronounced for DTA than for SAFRAN-Crocus.This difference may come from the idealized relief of SAFRAN-Crocus, which ignores the possible presence of neighbouring slopes, bringing either shade during part of the day or additional illumination by reflection.
The influence of slope is less visible than that of aspect.SSA tends to decrease when the slope increases, but the influence of slope actually depends on the aspect and date.The influence of slope is complementary to the one of aspect, as the local solar incidence angle is a combination of both parameters (see Sect. 3.1).For sun facing slopes and when the sun elevation is low in winter, the solar illumination angle decreases as the slope increases.Hence, the SSA decreases faster with higher slopes, which explains the negative impact of slope.Though, this parameter has less influence than aspect on the SSA, as the slope range is smaller than the aspect range.Here again, this effect is more pronounced for DTA than for SAFRAN-Crocus.In addition, the forward scattering peak is more often sampled from MODIS for sun facing slopes.Consequently, the MODIS retrieved SSA with DTA www.the-cryosphere.net/7/741/2013/The Cryosphere, 7, 741-761, 2013 may be biased high.This suggests that the SSA decrease is even more pronounced in reality than the DTA method shows.Additional field measurements are needed to lend further support to this observation.The influence of the massif is briefly presented and discussed below.The Belledonne massif features the lowest SSA of the DTA record, but the highest SSA for SAFRAN-Crocus outputs.These trends are not monotonic and actually dependent on the date, although SAFRAN-Crocus and DTA also often differ in their temporal trends.This is the only predictor considered where the retrievals and the model disagree on the trends of variability.This behaviour needs a larger dataset to be analyzed in detail.
These results strengthen our confidence in MODIS data retrievals, as the technique show physically consistent behaviours of retrieved SSA.The relative differences between retrievals and SAFRAN-Crocus -mainly, gradient intensity or massifs heterogeneities -calls for further in-depth comparisons.A larger set of data over a whole snowy season, and validation against field measurements of SSA, are required to further analyze these results, buiding on the methodology outlined in the present study.

Conclusions
In this work, we assessed the effect of accounting for (1) the local topography and multiple reflections, (2) the anisotropy of snow and ice reflection, (3) the shape of snow grains, in snow grain size retrievals from MODIS data in mountainous areas.We compared different methods (with combinations of these effects) on six clear sky dates representative of the different states of the winter snowpack, in the French Alps.In terms of SSA, the method using a topographic correction on absolute reflectance (though assuming Lambertian surface and spherical disconnected particles) is the most consistent with the SSA provided by the estimates of SAFRAN-Crocus and with field measurements.In this case, the overall RMSD is 9.4 m 2 kg −1 and the mean difference is 0.1 m 2 kg −1 .The differences with the other methods are larger on date with high mean SSA values (high reflectance), where the other methods (including band ratioing) tend to give higher SSA than SAFRAN-Crocus.
Furthermore, the method shows significant SSA relationships to topography.Mainly, the retrieved SSA increases with elevation, and decreases with the amount of incident solar radiation, which is linked to the solar zenith angle.The comparison of these SSA variations to SAFRAN-Crocus reveals that they are consistent, although the MODIS retrieved values show more variability linked with topography than SAFRAN-Crocus estimates.
The results expressed in terms of optical radius are less pronounced.They still show a positive influence of the topographic correction, but with as much benefit as when only assuming non-spherical particles and correcting for the anisotropy.The respective improvements of using either topographic correction or non-spherical particles encourages the combination of the two methods in further work.Furthermore, the presence of mixed pixels could not be excluded and surface roughness is not taken into account.These two facts can partly explain that the anisotropy corrections, not adequate to mixed pixels and roughness, does not benefit to the retrieval.In future work, these problems should be addressed in order to reduce the uncertainties of retrieval implied by these features.
The mean variability of SSA within the topographic classes of SAFRAN-Crocus estimated using MODIS data varies from 8 to 6 m 2 kg −1 when the number of classes increases from roughly 300 to 400.Along with the results of the field measurements, this seems to indicate that semidistributed snowpack modelling with a sufficient number of classes gives a satisfying representation of the snowpack variability for these dates in the French Alps at the scale of the MODIS 500 m pixels.This is consistent with the results of Fiddes and Gruber (2012), who demonstrated that lumped land surface simulations could match the performance of fully distributed simulations provided that a minimum number of adequately chosen classes are considered.This study opens the way to the assimilation of MODIS-derived information over wide geographical areas in the detailed snowpack model Crocus.Our findings indicate that such an assimilation may be appropriate in the fixed altitude/slope/aspect/massif framework of the SAFRAN-Crocus-MEPRA chain used operationally for avalanche warning activities in France (Durand et al., 1999), and even potentially surpass direct assimilation in a fully distributed simulation where inter-pixel variability ("noise") from MODIS data could be detrimental to the model overall performance.

Fig. 1 .
Fig. 1.(a) Geographical location of the studied domain (small red rectangle).(b) Elevation from SRTM DEM on the studied domain.Stars show locations of field measurements.

Figure 3
Figure 3 shows the daily mean air temperature provided by SAFRAN at 2100 m, and snow height simulated by SAFRAN-Crocus at 1800 m and 2100 m on flat terrain for the season 2008-2009, on the Grandes-Rousses massif.The six dates of MODIS acquisitions in 2009 are marked as vertical lines.The scenes are characterized by recent snowfalls for 2009-01-09, 2009-01-25, settling periods for 2009-02-26 and 2009-03-21, a recent snowfall only above 2000 m for 2009-04-22, and a melting period for 2009-05-06.Figure4shows the colour composite of MODIS images on the six dates of the study.Figure5illustrates the SSA obtained from DTA retrievals for all dates.

FigFig. 3 .
Fig. 3. Snow height simulated by SAFRAN-Crocus for the season 2008-2009 on the Grandes-Rousses massif, at 1800 m (blue dashed line), and 2100 m (blue plain line).The green curve plots daily mean air temperature, issued from SAFRAN at 2100 m.Red vertical lines indicate MODIS acquisitions dates.

Fig. 3 .
Fig. 3. Snow height simulated by SAFRAN-Crocus for the season 2008-2009 on the Grandes-Rousses massif, at 1800 m (blue dashed line), and 2100 m (blue plain line).The green curve plots daily mean air temperature, issued from SAFRAN at 2100 m.Red vertical lines indicate MODIS acquisitions dates.

Fig. 5 .
Fig. 5. SSA maps obtained from DTA. Elevation lines are 1000 m spaced.Grey pixels are shaded areas.White pixels are not detected as snow.Delimitation of the SAFRAN massifs Belledonne, Grandes-Rousses and Oisans are also drawn.

Fig. 7 .
Fig. 7. Mean SSA (a) and optical radius (b) obtained from satellite retrievals versus SAFRAN-Crocus mean value, over the six studied dates.DTA in blue, DTAA in red, DA in green, DR in purple and ART in orange.
Fig. 10.Standard deviation of the difference between DTA and SAFRAN-Crocus SSA values within each topographic class of the model represented by more than 10 pixels.The values are computed including data for the six dates studied for two configurations of the model : the one presented in Section 2.5 and a reduced one using only 3 slopes (20 and 40 • ) and flat and 6 aspects (N, E, SE, S, SW, W).The thick black line is the median value.The box represents the first and the third quartiles.Minimum and maximum values are also indicated for each configuration. 39

Fig. 10 .
Fig. 10.Standard deviation of the difference between DTA and SAFRAN-Crocus SSA values within each topographic class of the model represented by more than 10 pixels.The values are computed including data for the six dates studied for two configurations of the model: the one presented in Section 2.5 and a reduced one using only 3 slopes (20 and 40 • and flat) and 6 aspects (N, E, SE, S, SW, W).The thick black line is the median value.The box represents the first and the third quartiles.Minimum and maximum values are also indicated for each configuration.

Table 1 .
Spectral definition of the seven MODIS channels at 500 m resolution.

Table 2 .
Discretisation of the SAFRAN-Crocus topographic parameters on our study area.

Table 3 .
E-folding depth for the MODIS bands used in this study for two extreme cases of snow: F -fresh light snow, r opt = 50 µm, density = 150 kg m −3 ; A -aged dense snow, r opt = 1000 µm, density = 350 kg m −3 .The values were calculated with DISORT and assuming spherical and disconnected snow grains.

Table 4 .
Nomenclature of the retrieval methods.

Table 6 .
Retrieved SSA, on first line (m 2 kg −1 ), and optical radius, on second line (µm), for the different methods over the six studied dates.The standard deviations are provided in parenthesis.

Table 7 .
Root mean square deviation (and mean difference in parenthesis) of the methods over the six studied dates.SSA on first line (m 2 kg −1 ), optical radius on second line (µm).The number of pixels on which the computations were carried out is similar to Table6.