Influence of grain shape on light penetration in snow

The energy budget and the photochemistry of a snowpack depend greatly on the penetration of solar radiation in snow. Below the snow surface, spectral irradiance decreases exponentially with depth with a decay constant called the asymptotic flux extinction coefficient. As with the albedo of the snowpack, the asymptotic flux extinction coefficient depends on snow grain shape. While representing snow by a collection of spherical particles has been successful in the numerical computation of albedo, such a description poorly explains the decrease of irradiance in snow with depth. Here we explore the limits of the spherical representation. Under the assumption of geometric optics and weak absorption by snow, the grain shape can be simply described by two parameters: the absorption enhancement parameterB and the geometric asymmetry factor gG. Theoretical calculations show that the albedo depends on the ratio B/(1−gG) and the asymptotic flux extinction coefficient depends on the product B(1−gG). To understand the influence of grain shape, the values of B andgG are calculated for a variety of simple geometric shapes using ray tracing simulations. The results show that B and(1− gG) generally covary so that the asymptotic flux extinction coefficient exhibits larger sensitivity to the grain shape than albedo. In particular it is found that spherical grains propagate light deeper than any other investigated shape. In a second step, we developed a method to estimate B from optical measurements in snow. A multi-layer, two-stream, radiative transfer model, with explicit grain shape dependence, is used to retrieve values of the B parameter of snow by comparing the model to joint measurements of reflectance and irradiance profiles. Such measurements were performed in Antarctica and in the Alps yielding estimates of B between 0.8 and 2.0. In addition, values ofB were estimated from various measurements found in the literature, leading to a wider range of values (1.0–9.9) which may be partially explained by the limited accuracy of he data. This work highlights the large variety of snow microstructure and experimentally demonstrates that spherical grains, withB = 1.25, are inappropriate to model irradiance profiles in snow, an important result that should be considered in further studies dedicated to subsurface absorption of short-wave radiation and snow photochemistry.


Introduction
Snow is involved in strong feedback loops in the climate system (Cess et al., 1991;Hall, 2004) because its optical properties depend on some snow physical characteristics (grain size, density, impurities) that themselves evolve with atmospheric conditions (Colbeck, 1982;Picard et al., 2012).Snow is highly reflective in the visible and near-infrared (near-IR) range (Warren, 1982) so that the amount of solar energy absorbed by the snowpack is small with respect to vegetation, bare soil or bare ocean for instance, but varies significantly with slight changes of grain size and impurity content (Wiscombe and Warren, 1980;Warren and Wiscombe, 1980).Snow is also translucent so that the solar energy penetrates within the snowpack and is partly absorbed in depth below the surface (Colbeck, 1989).Light penetration depends on density, grain size (Bohren and Barkstrom, 1974) and is sensitive to very low amounts of impurities (Warren et al., 2006).Although Brandt and Warren (1993) highlight that most of the absorption occurs in the IR and takes place in the very top centimeters of the snowpack, the penetration of light has Q.Libois et.al: Influence of grain shape on light penetration in snow a crucial impact on the thermal regime (Liston and Winther, 2005;Flanner, 2005;Kuipers Munneke et al., 2009;Picard et al., 2012) and on the availability of photons for photochemical reactions in the ultraviolet and visible (King and Simpson, 2001;Domine et al., 2008).The amount and vertical distribution of the absorbed energy are determined by the albedo and the asymptotic flux extinction coefficient, respectively.Both quantities are highly wavelength dependent (Warren, 1982).Modelling these quantities from given snow physical properties is crucial for general circulation models (Hall and Qu, 2006;Waliser et al., 2011) and photochemistry models (e.g.Lee-Taylor and Madronich, 2002).
The computation of snow optical properties from microstructural properties (often reduced to grain size and density) has been investigated in many studies using the radiative transfer theory (e.g.Kokhanovsky, 2004).According to this theory, snow is represented by a continuous medium with appropriate absorption and scattering properties.In order to calculate these quantities from the snow physical properties that can be measured in the field or predicted by snow evolution models (Brun et al., 1989(Brun et al., , 1992;;Vionnet et al., 2012), it is necessary to assume that snow is composed of individual grains and that the grains have a particular shape.Equivalent spheres have often been used (Bohren and Barkstrom, 1974;Wiscombe and Warren, 1980;Brandt and Warren, 1993), firstly because Lorenz-Mie theory (e.g.van de Hulst , 1981) provides an analytical and rigorous formulation and, secondly, because the albedo calculated with other grain shapes can be close to the albedo calculated with spheres if the grain size is appropriately chosen (Warren, 1982).In particular, Grenfell and Warren (1999) emphasized that any grain type can be replaced by a collection of spheres of equal specific surface area as for irregular grains albedo calculations.However subsequent theoretical studies show that replacing snow grains by spheres of equal specific surface area is not entirely satisfying since single-scattering properties, including single-scattering albedo and asymmetry factor, are different from one shape to another (Neshyba et al., 2003;Kokhanovsky and Zege, 2004;Grenfell et al., 2005;Picard et al., 2009).
Also, the spherical assumption seems inadequate for the prediction of the bi-directional reflectance function of snow (Sergent et al., 1998;Aoki et al., 2000;Dumont et al., 2010), that is critical for the interpretation of remote sensing data.Negi and Kokhanovsky (2011) use fractal particles for satellite snow grain size retrieval.Zege et al. (2008) point out that for snow composed of dendrite or faceted crystals the snow grain size retrieval algorithm they use leads to errors of 50 % if grains are assumed spherical.On the contrary, the experimental study on metamorphosed snow by Gallet et al. (2009) shows that spheres are adequate to relate bi-hemispherical albedo and specific surface area within ±12 %.Although snow grains are not spherical, the spherical assumption remains the simplest way of representation and is largely used for meteorological and climate purposes (van den Broeke et al., 2008;Kuipers Munneke et al., 2012), as well as for photochemistry applications (Zatko et al., 2013).
Quantitative studies of irradiance profiles in snow are scarce compared to studies of snow reflectance (Bohren and Barkstrom, 1974;Perovich, 2007).Performing accurate irradiance measurements is difficult since the introduction of sensors perturbs the medium (Dunkle and Bevans, 1956;Giddings and LaChapelle, 1961) and repeat experiments demonstrate uncertainties of ∼ 20 % as a result of measurement error and natural variability of snow (France et al., 2011b;France and King, 2012).The few attempts to reproduce irradiance measurements, using radiative transfer models with the spherical assumption, failed (Sergent et al., 1987;Meirold-Mautner and Lehning, 2004).In particular, Sergent et al. (1987) measured asymptotic flux extinction coefficients twice greater than those predicted by a model with spherical grains (their Fig. 6).Haussener et al. (2012) have compared the transmittance of snow slabs obtained by Monte Carlo ray tracing performed on computed tomography images, to the transmittance calculated with the radiative transfer code DISORT (Stamnes et al., 1988) assuming equivalent spheres.They found that transmittance computed with DIS-ORT is always much larger than that obtained with the ray tracing method (their Fig. 9b and d).They attribute part of the difference to the simplification of snow morphology in DISORT.
The objective of the present paper is to investigate the performance of the spherical assumption to compute albedo and irradiance profiles.To this end, this work is based on an analytical formulation of the radiative transfer theory, where snow grain shape is fully defined by two parameters, the absorption enhancement parameter B and the geometric asymmetry factor g G (Kokhanovsky and Zege, 2004).In Sect.2, we derive analytical expressions of the two measurable optical properties (the albedo and the asymptotic flux extinction coefficient) in terms of B, g G , specific surface area, density, and optical indices.These expressions show that the albedo and the asymptotic flux extinction coefficient depend on grain shape in different ways.In order to quantify the impact of the grain shape on snow optical properties, ranges of variations of B and g G are needed.To the authors' knowledge, these shape parameters have never been measured for snow.Hence the following sections are aimed at determining the values of B and g G for snow by using two approaches.The first approach is numerical and uses Monte Carlo ray tracing methods (Sect.3) to determine the values of B and g G for a set of simple geometric shapes (Kokhanovsky and Macke, 1997;Picard et al., 2009).The second approach is experimental and uses optical measurements in snow (Sect.4).It allows the determination of B only.Section 5 discusses the values obtained from these two approaches.

A two-stream radiative transfer model with explicit shape dependence to compute irradiance profiles in snow
Diffusion of light in snow is modelled using the radiative transfer theory.Scattering and absorbing properties of snow depend on the single-scattering properties of individual snow grains.The latter depend on grain size and shape (Kokhanovsky, 2004), as detailed in Sect.2.1.In the field, only macroscopic optical properties of the whole snowpack or of individual snow layers (albedo and asymptotic flux extinction coefficient) are measurable.The asymptotic analytical radiative transfer theory (Kokhanovsky and Zege, 2004) provides analytical expressions for the albedo α and the asymptotic flux extinction coefficient k e of a homogeneous snow layer, as a function of snow grains single-scattering properties (Sect.2.2).As natural snowpack is usually stratified, these analytical expressions are used in a multi-layer model based on the Eddington approximation (Jiménez-Aquino and Varela, 2005), to calculate the albedo and irradiance profiles of a stratified snowpack (Sect.2.3).Snow is never pure, so that light-absorbing impurities have to be taken into account in the model in order to interpret experimental data in the visible range (Sect.2.4).

Single-scattering properties of pure snow
Diffusion of light in snow is often modelled using the radiative transfer theory, which assumes that snow is a continuous medium (Bohren and Barkstrom, 1974;Wiscombe and Warren, 1980).Each layer of the snowpack is characterized by the phase function, extinction coefficient σ e and singlescattering albedo ω which are determined from snow grain size and shape, and snow density (Bohren and Barkstrom, 1974).σ e is related to the scattering and absorption coefficients by σ e = σ s + σ a and ω = σ s /σ e .In the Two-stream Analytical Radiative TransfEr in Snow (TARTES) model that we developed, the snow grains in a particular layer have the same size and shape and their size is much larger than the wavelength.In such a case, calculations can be performed with the geometric optics approximation.In addition, we limit our study to the spectral range 300-1350 nm where ice is globally weakly absorbing (Warren and Brandt, 2008), that is 1−ω 1, even though ice is considerably more absorptive at 1350 nm than at 300 nm.Snow grains are treated as incoherent scatterers (Bohren and Barkstrom, 1974;Wiscombe and Warren, 1980), so that the extinction and absorption coefficients σ e and σ a (m −1 ) are proportional to the extinction and absorption cross sections C ext and C abs (m 2 ) of individual snow grains: where n is the grain concentration (m −3 ) related to the volume V of an individual grain and to snow density ρ by C ext is twice the average projection area of the grain (Zege et al., 2008): Scattering occurs each time a photon reaches the grain surface while absorption is due to the photon path within the grain.C abs is thus where γ = 4π χ λ is the ice absorption coefficient and the wavelength dependent imaginary part of the ice refractive index χ can be found in Warren and Brandt (2008).The coefficient B depends on the shape of the particle and on the difference between the real parts of the refractive indices of ice and air.B is called the "absorption enhancement parameter" (Kokhanovsky and Zege, 2004) since it quantifies the enhancement of absorption (compared to a straight trajectory through the grain) due to lengthening of the photon paths within the grain owing to multiple internal reflections.B is independent of the size of the particle.By definition of ω the single scattering co-albedo is given by The average cosine of the scattering angle, determined from the phase function, is called the asymmetry factor and is denoted g.In the model TARTES, the phase function is entirely characterized by g, which is the average of the geometric term g G and the diffraction term g D .For particles large compared to the wavelength, diffraction is essentially forward so that g D 1 ( Kokhanovsky and Zege, 2004) and The geometric asymmetry factor g G depends on the shape and real part of the ice refractive index.It is independent of the grain size, like B. In the spectral range 300-1350 nm, ice is weakly absorbing and the real part of the ice refractive index is nearly constant (according to Warren and Brandt (2008) it varies between 1.30 and 1.33), so that variations of B and g G are less than 3 % (Kokhanovsky, 2004).Hence in TARTES, B and g G are assumed independent of the wavelength.

The asymptotic analytical radiative transfer theory
The asymptotic analytical radiative transfer theory presented in Kokhanovsky (2004) assumes that the bi-hemispherical Q. Libois et.al: Influence of grain shape on light penetration in snow albedo α (albedo of a semi-infinite scattering medium illuminated by a diffuse source, hereafter referred as albedo) and the asymptotic flux extinction coefficient k e are expressed in terms of g and ω for weakly absorbing wavelengths by = k −1 e is the e-folding depth: the depth at which diffuse irradiance has decreased by a factor e. Similar expressions have been derived or used in former studies (Bohren and Barkstrom, 1974;Wiscombe and Warren, 1980).With the formalism presented in Sect.2.1, Eqs. ( 8) and ( 9) become In Eqs. ( 8)-( 11) the quantities k e and α are wavelength dependent since γ is.From the combination of Eqs. ( 10) and ( 11), the value B of a semi-infinite homegeneous slab of snowpack can be estimated from the albedo at wavelength λ α and the asymptotic flux extinction coefficient at wavelength λ k e : where λ k e and λ α are wavelengths in the range of validity of the theory.They can be different or equal.Eqs. ( 10) and ( 11) are valid for concave and convex particles.Although snow grains are partly convex and partly concave, in the specific case of convex particles Eqs. ( 10) and ( 11) can be simplified.The specific surface area SSA (m 2 kg −1 ) is the total surface area per unit mass.For convex particles, is proportional to the surface area S of the grain: S = 4 (Vouk, 1948;Zege et al., 2008), so that Hence for convex particles, α and k e are eventually expressed in terms of ρ, SSA, B and g G by: In Eqs. ( 14) and ( 15), (1−g G ) is multiplied by SSA, which shows that (1 − g G ) cannot be optically measured, unless an independent measure of SSA is available (e.g.Domine et al., 2001;Flin et al., 2004).

Extension to a multi-layer snowpack using the Eddington approximation
The two-stream formulation of the radiative transfer theory (e.g.Choudhury, 1981) has been widely used in snow optics to describe diffusion of light in a layered plane-parallel snowpack (Dunkle and Bevans, 1956;Schlatter, 1972;Bohren and Barkstrom, 1974;Wiscombe and Warren, 1980;King and Simpson, 2001).The formulation is accurate when photons experience a large number of scattering events before escaping the snowpack or being absorbed (i.e. at weakly absorbing wavelengths).The model TARTES uses the Eddington approximation (e.g.Jiménez-Aquino and Varela, 2005) to solve the two-stream radiative transfer equation.In that case, the direct irradiance I dir and the downward (I ↓ n ) and upward (I ↑ n ) diffuse irradiance in a horizontally homogeneous layer take the form (e.g.Jiménez-Aquino and Varela, 2005): where τ (z) is the optical depth at depth z, I the direct incident irradiance and µ 0 = cos θ 0 , with θ 0 the local zenith angle.α n and k en are the albedo and asymptotic flux extinction coefficient of layer n if the layer were infinite.The terms proportional to e −τ (z)/µ 0 in Eqs. ( 17) and ( 18) correspond to the contribution of direct light that has been scattered only once, and quickly vanish with depth (detailed expressions of G ↓ n and G ↑ n are given in Jiménez-Aquino and Varela ( 2005)).Equations ( 17) and ( 18) are applied to each layer of the Nlayer snowpack.The continuity conditions state that the diffuse irradiance I ↓ n and I ↑ n be continuous at each interface at depth L n : It gives 2(N −1) equations.The upper boundary condition is given by the diffuse incident irradiance I 0 and the lower one by the albedo α b of the underlying surface, which must be prescribed: This provides 2 more equations.These 2N equations form a linear system with 2N unknowns A n and B n .The coefficients A n and B n are determined by solving the system, which allows explicit calculation of the diffuse irradiance at any depth.The two-stream method has been extensively described (Shettle and Weinman, 1970;Wiscombe, 1977;Toon et al., 1989).The total irradiance is obtained by summing the diffuse and direct irradiance.The albedo is deduced as the ratio of the total upwelling to the total downwelling irradiance at the surface.

Impurities
Equations ( 10) and ( 11) have been derived assuming pure snow.However, light-absorbing impurities, such as black carbon (BC) (e.g. Warren and Wiscombe, 1980;Flanner et al., 2012), dust (Warren and Wiscombe, 1980;Painter et al., 2007) and humic-like substances (Hoffer et al., 2006;France et al., 2012), impact snow optical properties, especially in the visible range (Warren and Wiscombe, 1980;Reay et al., 2012).These impurities are taken into account in TARTES.Whether impurities are inside or external to ice grains is critical to model their impact.However little is known about it (Flanner et al., 2012).We assume that they are external to ice grains.This hypothesis has been questioned in several papers (e.g.Chýlek et al., 1983;Bohren, 1986).We also assume that these particles are small compared to the wavelength and that they do not affect the total phase function.In contrast, they have a strong influence on the absorption coefficient.The absorption cross section C i abs of a small spherical impurity of type i is proportional to the particle volume V i and depends on the complex refractive index m i (Kokhanovsky, 2004): The total absorption coefficient is where n i is the number density of impurity i particles (m −3 ) and can be expressed as , where ρ 0 i is the bulk density and ρ i the impurity concentration (kg m −3 ).The latter is related to the impurity content by c i = ρ i ρ snow (kg kg −1 ).For low impurity contents, we assume that scattering is essentially due to ice grains, and that σ e is unaltered.In that case, Eq. ( 6) is replaced by Similar expressions have been used by Zege et al. (2008) and Kokhanovsky (2013) to take into account analytically the influence of impurities on snow albedo.
To summarize, TARTES takes as inputs for each layer the density, depth, B, (1−g G ) /V (or SSA(1−g G ) if the particles are convex), the type and content of impurities.In addition, the direct and diffuse incident irradiance and the albedo of the underlying surface are required.The model computes the internal irradiance at any depth, from which the absorbed energy can be determined for each layer.The e-folding depth (λ) between two depths z 1 and z 2 (with z 2 > z 1 ) is calculated at any wavelength as To evaluate TARTES, we compared it to DISORT simulations run with a monodispersion of spherical particles (Wiscombe, 1980).To this end we chose B = 1.25 and g G = 0.79 in our model.For a large set of snowpack with black carbon content varying from 1 to 1000 ng g −1 and specific surface area varying from 5 to 100 m 2 kg −1 , we imposed various incident light conditions, mixing diffuse light and direct light with zenith angle ranging from 0 • to 85 • .The spectral relative error between albedos and irradiance calculated with TARTES and those predicted by DISORT was lower than 3 %.

Theoretical calculations of B and g G for various grain shapes
In the theory described in Sect.2, snow grain shape is fully described by two parameters: B and g G .Values of these parameters for geometric shapes are scarce in the literature.We gather values from the literature and from our calculations (for an extended set of geometric shapes), and we analyze the range of variations.Monte Carlo ray tracing methods have been used by Takano and Liou (1989); Mishchenko et al. (1996); Macke et al. (1996); and Xie et al. (2006) to determine singlescattering properties of ice particles in clouds.These studies showed that the phase function, single-scattering albedo and asymmetry factor largely depend on grain shape.The parameters B and g G were calculated for spheroids, hexagonal plates and fractal particles by Kokhanovsky and Macke (1997) and Kokhanovsky and Zege (2004) using ray tracing.In addition, we have performed calculations for cubes, cuboids, cylinders and depth hoar using the ray tracing model SnowRat (Picard et al., 2009) as follows.Photons are launched from different locations on one grain and their trajectories within the grain are calculated assuming geometric optics laws with optical index of ice at 900 nm.The process is repeated for 10 6 photons and 10 4 different grain orientations.g G is then computed by averaging the cosine of the deviation angle of photons that have escaped the grain.ω is the ratio of photons that have escaped the grain with respect to the number of launched photons.and V are directly computed in SnowRat and B is deduced using Eq. ( 6).We studied shapes with aspect ratio between 0.25 and 4, as defined in Picard et al. (2009).The values of B and g G for these shapes are summarized in Table 1 and reported in Fig. 1a.All numerical calculations were performed on crystals smaller than 1 mm so that absorption within a single grain is very low and the calculated values of B and g G do not depend on grain size.In Fig. 1a, (1 − g G ) is plotted as a function of B for all shapes (the quantity (1−g G ) rather than g G is used, to reflect the dependence in Eqs. ( 10) and ( 11)).For the tested shapes, B varies from 1.25 to 2.09 and is minimum for spheres.g G varies from 0.50 to 0.80.The relative variations of (1 − g G ) are larger than those of B. For spheres, g G = 0.79 and B = 1.25, highlighting great forward-scattering and weak absorption enhancement.For most of the hexagonal plates they studied, Grenfell et al. (2005) noticed that the asymmetry factor was lower than that of spheres.Jin et al. (2008) also point out that the spherical assumption overestimates the forward reflected radiances compared to measurements performed in Antarctica by Hudson et al. (2006).They attribute this to the forward peak of the scattering phase function of spheres.
We examine then the influence of grain shape on snow macroscopic optical properties.Since the albedo depends on the quotient B/(1 − g G ) and the asymptotic flux extinction coefficient on the product B(1 − g G ), Fig. 1b shows an alternative representation of Fig. 1a.In this representation, two snowpack composed of grains with shapes that have the same abscissa (respectively ordinate) have the same asymptotic flux extinction coefficient (respectively albedo).The vertical (respectively horizontal) line shows the shapes that have the same asymptotic flux extinction coefficient (respectively albedo) as spheres.These lines are "iso-e-folding depth" and "iso-albedo" and are also plotted in Fig. 1a. Figure 1b shows that snow grains with equal specific surface area but different shape can yield very different asymptotic flux extinction coefficients (given Eq. ( 14), the relative difference is 85 % between fractals and spheres of same specific surface area).As for the albedo, at λ = 1310 nm and SSA = 20 m 2 kg −1 for instance, the maximum relative difference is 37 % (see Eq. 15).Therefore the spherical assumption may be adequate for The Cryosphere, 7, 1803-1818, 2013  Kokhanovsky and Macke (1997) and were obtained with a complex refractive index m = 1.33.The value for fractal is from Kokhanovsky and Zege (2004) and was obtained at m = 1.31 − 10 −7 i.The values for the other shapes were obtained with SnowRat (Picard et al., 2009) at m = 1.30− 4.2 × 10 −7 i (which is the ice refractive index at 900 nm).For hexagonal plates, cuboids and cylinders, the aspect ratio is the ratio of the height to, the length of the hexagonal side, the length of the square side and the radius of the circle, respectively.For spheroids, it is the ratio of the largest semi-axis to the shortest one.

Shape
Aspect computing the albedo of a snowpack, but for the same snowpack may be inadequate to model irradiance profiles.For instance, spheroids with aspect ratio 0.5 (labeled in Fig. 1a and  b) are very similar to spheres in terms of albedo, but not in terms of e-folding depth.Figure 1b also highlights the position of spheres among all other shapes.The median value for B/(1 − g G ) is 4.7, in agreement with the scaling constant determined by Kokhanovsky (2006) for snow (that corresponds to a ratio of 4.6).For spheres, B/(1−g G ) = 5.7 and there are shapes above and below the horizontal line, suggesting that spheres may be a fairly good approximation to a mixture of various geometric shapes.This might explain the success of the spherical assumption for albedo calculations (e.g.Grenfell et al., 1994).In those calculations, the quotient B/(1 − g G ) in snow is assumed equal to that of spheres.In contrast, for irradiance profiles calculations, the spheres appear to be an extreme case, since the product B(1 − g G ) is minimum for spheres, almost two times smaller than the median value.Although settled snow that has experienced metamorphism is unlikely to look like geometric crystals (Bohren, 1983), this suggests that the spherical assumption is a shortcoming to model irradiance profiles in snow.

Experimental determination of B from optical measurements
In this section we propose two methods to retrieve B from optical measurements of snow.The first method assumes a single homogeneous layer.It is relevant for data from the literature (Sect.4.3.2) for which a detailed stratigraphy of the snowpack is absent.The second method treats the case of a multi-layer snowpack and uses field measurements we specifically performed for the B determination, under experimental conditions chosen to maximize the quality of the retrieval method.

Method
Method 1: Knowledge of the asymptotic flux extinction coefficient k e , and the albedo α, of a homogeneous snow layer allows calculation of the value of B by using Eq. ( 12) for pure snow, or Eqs. ( 8), ( 9) and ( 25  are measured.Using TARTES, an iterative procedure is completed to find the values of B and BC content that produce the best match between the modelled and measured irradiance profiles at two different wavelengths λ 1 k e and λ 2 k e (where λ 2 k e > λ 1 k e ), or similarly between the modelled and measured spectral e-folding depths within a particular layer.Values of λ 1 k e and λ 2 k e are chosen so that snow optical properties at λ 1 k e are much more sensitive to impurity content than at λ 2 k e .The greater the difference between ice absorption at λ 1 k e and λ 2 k e , the faster the convergence of the iterative process.In Sect.4.2 both wavelengths are in the range 500-780 nm.Each numerical layer of TARTES is 1 cm thick.For each layer, the density, (1 − g G ) /V , B, the type and content of impurities must be specified.Here we consider that density and reflectance at λ α are given (see Sect. 4.2 for measurements details, where λ α =1310 nm), while B and impurity contents are unknown parameters.Although B is likely to be different from one snow layer to another since snow type varies with depth, here B is assumed uniform within the snowpack.Regarding impurities, only black carbon is considered since it is the major contributor to light absorption by impurities above 500 nm, where optical measurements are taken (Sergent et al., 1993;France et al., 2011a).The optical properties of black carbon depend on the nature of its elementary constituents (Bond and Bergstrom, 2006).Here the refractive index of black carbon is taken from Chang and Charalampopoulos (1990) and its density is assigned to 1000 kg m −3 (Warren and Wiscombe, 1980).The black carbon content used in TARTES is denoted BC and is also assumed uniform within the snowpack.TARTES also requires the albedo of the underlying surface α b and the direct and diffuse incident irradiance I direct This assumption has no incidence on the irradiance calculations for the top 40 cm (maximum depth of our measurements) since we simulated at least 1 m-deep snowpacks.Irradiance measurements are taken at depths greater than 3 cm.At such depths, irradiance is entirely diffuse in the visiblenear-IR range (Barkstrom, 1972), so we can take I direct 0 = 0, I diffuse 0 = I 0 (λ) without loss of generality.I 0 (λ) is an unknown parameter.The iterative method to retrieve the optimal B opt , BC opt and I 0 (λ) follows a three-step procedure: 1.The measured irradiance profile at λ 2 k e is used to find a first guess for B, assuming BC 0 = 100 ng g −1 .
For each pair (B, I 0 (λ 2 k e )) in the range of optimization, (1 − g G ) /V is deduced for each layer from the reflectance profile using Eq. ( 11) and an irradiance profile is computed (B is varied in the range 0.5-2.5 with an increment of 0.1 and I 0 (λ 2 k e ) is varied by a factor 10 around the irradiance obtained from exponential extrapolation of the measurements at the surface).To measure the quality of the simulation, we consider two error metrics depending on the available measurements.
Case 1: The irradiance profile is measured.The difference 1 between the measured and modelled profiles is calculated as where z n is the depth where measurement n was taken and I (z n , λ) the modelled irradiance at that depth.
The pair (B, F 0 (λ 2 k e )) that minimizes 1 (λ 2 k e ) (or 2 (λ 2 k e )) provides a first estimate of B in the snowpack.2. The value of B found at step 1 is used to model the irradiance profile at λ 1 k e .BC is chosen by minimizing 1 (λ 1 k e ) (or 2 (λ 1 k e )). 3. The latter value of BC is then used at step 1 to refine B and so on.Iteration using λ 2 k e and λ 1 k e is performed until the optimal pair (B, BC) between two iterations does not change.In practice, less than 10 iterations are needed for convergence.
In both methods 1 and 2, the retrieved value of B is sensitive to the precision of the density, reflectance and irradiance measurements.In fact the derivation of Eq. ( 12) reads: It highlights the high sensitivity of B retrieval to the reflectance measurement if reflectance is close to 1.The larger λ α , the lower α(λ α ) and the greater the precision of B retrieval.Using λ α = 1310 nm for our measurements is thus more precise than using reflectance measurements in the visible range.

Materials
A first set of measurements were conducted at Dome C (75.10 • S, 123.33 • E; 3233 m a.s.l.), Antarctica, during the summer campaign 2009/2010.Spectral e-folding depths were measured in two different snowpits following the procedure detailed in France et al. (2011a).In each snowpit, e-folding depth was measured in a layer identified visually.Density at 3 cm resolution was measured with a rectangular cutting device of volume 300 cm 3 and a 0.3 g precision scale.Profiles of nadir-hemispherical reflectance at 1310 nm were obtained with the snow reflectance profiler POSSSUM (Arnaud et al., 2011).Following Eq. ( 1) of Picard et al. (2009), the conversion from nadir-hemispherical reflectance α nh (z) to diffuse reflectance α(z) is given by ln(α(z)) = 7 9 ln(α nh (z)). (30) Method 2 (case 2) described in Sect.4.1 is applied to each snowpit with λ 1 k e = 500 nm and λ 2 k e = 700 nm.The second set of measurements was collected during the winter 2012 at three different sites in the French Alps: Col de Porte (45.17 • N, 5.46 • E; 1326 m a.s.l.), Lacs Robert (45.08 • N, 5.55 • E; 2000 m a.s.l.), and Lac Poursollet (45.05 • N, 5.90 • E; 1658 m a.s.l.).For each site, a 1 m-deep pit was dug in a flat, horizontal and unaltered snow area.The spectral irradiance I (z, λ) is measured with fiber optics as follows.First, 60 cm long horizontal holes with the same diameter as the fibers are prepared at approximately 5 cm intervals from the surface to a maximum depth of 30 cm (a guiding support is used to ensure their horizontality).At greater depth, light entering the snowpack from the pit face perturbs the internal radiation field (Bohren and Barkstrom, 1974;Warren et al., 2006;France et al., 2011b).The holes are made in the sun direction in order to avoid shading by the operator.Then, one bare fiber optic with a 25 • field-of-view is successively inserted into the holes from the top to the bottom of the snowpack.The fiber is placed 2 cm away from the hole extremity to limit the perturbation due to the presence of the fiber.The spectral radiance is recorded at each depth with an Ocean Optics MAYA 2000 ® spectrophotometer with spectral resolution 0.5 nm.The measurement in the uppermost hole is duplicated at the beginning and at the end of the measurement series.Series with variations of the incident irradiance larger than 3 % were discarded.Since only the relative variations of irradiance with depth are of interest, and since the theoretical framework introduced in Sect. 2 ensures that the measurements taken with the fiber optic are directly proportional to irradiance, the measurements are hereafter referred as irradiance measurements.Profile of density with 5 cm vertical resolution is measured using a cylindrical cutting device of volume 100 cm 3 and a 0.1 g precision scale.Profile of reflectance at 1310 nm is measured at 1 cm resolution with the instrument ASSSAP (the Alpine version of POSSSUM), except for the 9 March and Lacs Robert measurements where diffuse reflectance was measured with the instrument DU-FISSS (Gallet et al., 2009) on samples extracted from the pit face.Method 2 (case 1) is used with λ 1 k e = 600 nm and λ 2 k e = 780 nm.

Application to measurements
The method 2 presented in Sect.4.1 is applied to the data collected at Dome C and in the Alps.Table 2 summarizes the values of B and the black carbon contents retrieved for each set of measurements, along with details on snow type and instruments used.At Dome C, the method yields B = 1.9 for the hard windpack layer and B = 2.0 for the depth-hoar layer.Figure 2a shows the irradiance profile at 600 nm measured at Col de Porte on 29 February 2012.Modelled irradiance profiles using the density and reflectance profiles shown in Fig. 2b are also plotted for various values of B and constant incident irradiance.Figure 2a illustrates the sensitivity of light e-folding depth to B. At 20 cm depth, irradiance at 600 nm varies by more than 50 % for B ranging from 1.25 to 2.0.With B = 1.25 (spheres), irradiance at depth is overestimated.The vertical variations of irradiance are well reproduced by the model, in particular the changes of slope at 10 and 20 cm.For this experiment, the best match is obtained for B = 1.9. Figure 3 shows the measured and modelled irradiance profiles at 780 nm for the other Alpine measurements.
The optimal modelled profile is symbolized by dark inverted triangles.B = 1.8 produces the best match for the 24 February 2012 measurements while B = 1.6 was the best for the measurements at Lacs Robert and Lac Poursollet.B = 0.8 is found for the fresh snow layer of the 9 March 2012 measurement.The seven values of B retrieved experimentally, ranging from 0.8 to 2.0, are shown in Table 2.This large range indicates that the asymptotic flux extinction coefficient can vary of a factor 1.6 solely due to the impact of grain shape on B.
The retrieved black carbon contents vary from 12 to 85 ng g −1 for the Alpine sites, which is in agreement with values between 34 and 247 ng g −1 reported by Sergent et al. (1998).Black carbon amounts retrieved for Dome C ( 5 ng g −1 ) are consistent with the 3 ng g −1 reported by Warren et al. (2006) for the upper 30 cm of the snowpack and with the values reported by France et al. (2011a).In fact, our method yields an upper limit of black carbon content since all impurities absorption is attributed to black carbon, which is probably not realistic (Zatko et al., 2013).Anyway, we do not aim here at retrieving a real black carbon content and this assumption does not affect the B estimation since such contents at Dome C have a negligible impact on snow optical properties at 700 nm (see Fig. 4).
Even though the retrieval method is based only on two wavelengths λ 1 k e and λ 2 k e for the irradiance profiles, we checked that modelled irradiance profiles at other wavelengths were in good agreement with the measured ones.
The measured and modelled e-folding depths are plotted for Dome C depth-hoar layer in Fig. 5a, and for the upper layer of the snowpack studied at Lacs Robert in Fig. 5b.Modelled e-folding depths for various values of B are also shown.For the best fit, the relative error is less than 5 % for the whole visible spectrum.In Fig. 5a, the discrepancies between the model and the measurement from 400 nm to 500 nm are probably due to the presence of other types of impurities (France et al., 2011a).

Application to data in the literature
Method 1 was applied to retrieve B using published data.There have been few simultaneous measurements of asymptotic flux extinction coefficient and reflectance of snow reported in the literature.Earlier studies dealt with broadband irradiance profiles (O'Neill and Gray, 1972;Fukami et al., 1985) so that general features were retrieved (such as the increase of broadband e-folding depth with depth) but no quantitative information on grain shape can be inferred.Studies with spectrally resolved measurements include the works of Grenfell and Maykut (1977) in the Arctic, Beaglehole et al. (1998) in Antarctica, Perovich (2007) on mid-latitude seasonal snow, France et al. (2011a) in Antarctica, France et al. (2011b) in Svalbard and France et al. (2012) in Alaska.Spectral diffuse reflectance and asymptotic flux extinction coefficient were measured or deduced from measurements.B values are calculated using Eq. ( 12) for all these cases, that is impurities are not taken into account.Experimental details and B values for each study are summarized in Table 3. B values are also reported in Fig. 6 The curve is deduced from Eq. ( 25).The dashed lines highlight the values at 600 nm and 780 nm (wavelengths used for the Alpine measurements) and the continuous lines highlight the values at 500 nm and 700 nm (wavelengths used for Dome C measurements).
retrieved value of B for each experiment, the largest available wavelengths are used in Eq. ( 12).Those wavelengths are specified in Table 3 since they depend on available data.
A very wide range of B is obtained, from 1.0 to 9.9.The largest theoretical value of B found in the literature was obtained by Kokhanovsky and Macke (1997) for spheroids with shape parameter ξ = 0.3 and equals 3.16.As a consequence, the values obtained from the literature that exceed 3.16 are questionable.Intuitively, high B values, i.e. high lengthening of the photon path within the grain, are only possible with specific shapes that trap the light, like fiber optics.Uncertainty in the measurements is a potential cause of these high values.In fact, Eq. ( 12) is applied with 700 nm < λ α < 900 nm, so that α(λ α ) is generally close to 1 and the uncertainty of the retrieved value of B is large (see Eq. 29).In Table 3, the five B values larger than 3.16 are in bold.They are obtained for the lowest values of α.This may indicate that Eq. ( 12) cannot be applied to wet or old coarse snow with low albedo, or that impurities amounts in these studies were too large for Eq. ( 12) to be be applied.This may also indicate that accurate reflectance measurements are more difficult than e-folding depth measurements.The spectral resolution of the instruments and the quality of the diffuse incident irradiance may also be accounted for.Eventually, the snow physical properties for those experiments (density, grain size) may be outside the range of validity of the theoretical framework presented in Sect. 2 and partly explain these high B values.

Discussion
According to Koh (1989), "until the phase function of snow is known, the use of a spherical equivalent phase function may be the most appropriate".Although the spherical assumption has proved successful for albedo modelling, it is much less appropriate to model irradiance profiles in the snowpack.The following discussion will explore the failure of the spherical grain assumption for irradiance data.
To model irradiance in snow, we developed a radiative transfer model, called TARTES, based on the asymptotic analytical radiative transfer theory.In this model snow grains are characterized by their specific surface area, and two shape parameters B and g G .The absorption enhancement parameter B quantifies the lengthening of photon paths inside a snow grain due to internal multiple reflections, and the geometric asymmetry factor g G quantifies the forward-scattering of light by the snow grain.Theoretical calculations performed in Sect. 3 show that spheres are highly forward-scattering (large g G ) and weakly absorbing (small B).Spheres are very singular among different other geometric grain shapes, explaining the difficulty to match irradiance measurements with models using spherical grains (e.g.Sergent et al., 1987;Meirold-Mautner and Lehning, 2004).Qualitatively, when B is large in a medium, photons are internally reflected several times in a grain before escaping, and are then likely to exit the grain in no particular direction.Conversely, when B is small, photons travel a minimal distance in the grain (i.e. they tend to exit at the opposite point of their entry, and tend to be scattered preferentially in the forward direction, which leads to large g G ). Thus B and (1 − g G ) are correlated.An important consequence is that the quotient B/(1 − g G ) -that determines the albedo -should be relatively constant with respect to grain morphology while the product B(1 − g G ) -that is important for the asymptotic flux extinction coefficient -varies widely.The numerical calculations presented in Fig. 1 confirm this analysis.This can be further explained by the fact that in snow with a large value of B, photons tend (i) to travel a long way within the ice grains (which increases likelihood of absorption) and (ii) to be scattered in all directions, including the backward direction.The penetration depth of the photons is limited because absorption is enhanced and because forward-scattering is reduced.Regarding the albedo, as the probability that a photon escapes the snowpack before being absorbed, photons in snow with large B are scattered in all directions and are likely to escape the medium after fewer scattering events than in snow with small B. However, this is counterbalanced by the higher probability of absorption between scattering events in snow with large B, because the path in the grains is longer.In short, the increase of photons paths within the grains for larger B is counterbalanced by the shortening of total path lengths in the snowpack due to lower g G .This explains the limited sensitivity of albedo to the shape of the grains.The low sensitivity was noticed experimentally by Domine et al. (2006)   studied snow samples with various crystal shapes and confirmed later by Gallet et al. (2009).They demonstrate indeed that the reflectance of snow samples of known SSA (Domine et al., 2001) is well reproduced using DISORT-Mie code (Wiscombe, 1980;Stamnes et al., 1988) (i.e. for a medium of independent spheres).
In this paper, we exploit the sensitivity of light e-folding depth to grain shape to estimate the parameter B using measurements of reflectance, density and irradiance in the snowpack, and TARTES model.Even though the optical measurements are taken in the visible-near-IR range, the retrieved values of B hold for the whole range 300-1350 nm (see Sect. 2).These B values -obtained by minimizing the difference between measured and modelled irradiance profilesare generally greater than the value for spheres (B = 1.25).It means that assuming spherical grains for these experiments would lead to an overestimation of the e-folding depth.B values calculated for various geometric shapes and deduced from data in the literature are also essentially larger than 1.25, which suggests that the spherical assumption is inadequate to model irradiance profiles in snow.By sampling different snow types, we found that old snow with large grains has a B close to that of spheroids with aspect ratio 0.7 (B = 1.6, see Table 1).Snow in Antarctica, that is faceted, has values around 1.9, which is also the case of hexagonal plates.Fine grains have intermediate values (1.6 ≤ B ≤ 1.9), prob-ably because they are composed of a mixture of rounded and regular faces.The result for fresh snow is more surprising.Indeed, the value B = 0.8 is incompatible with the definition of B (Kokhanovsky, 2004) that constrains its value to be larger than 1.Until more measurements are performed, this value should not be considered as typical of fresh snow.Inaccurate measurements, or the small size of fresh snow crystals, close to the wavelength for which geometric optics do not apply, may be the cause of this invalid estimation.Shadowing effects and inter-particle interference (Wiscombe and Warren, 1980;Warren, 1982;Kokhanovsky, 2004), which are not taken into account in the model, may also explain this low value of B. Except for fresh snow, the values of B obtained from our measurements are coherent with the values predicted by theoretical calculations (Fig. 6).In contrast, the values deduced from the literature vary to a large extent.In fact, the sensitivity of the retrieval method to the measurements, illustrated by Eq. ( 29), points out that measurements should be taken carefully, and at appropriate wavelengths.We recommend not to consider these latter B values in the future.
Although there is no mention of experimental determination of B in the literature to compare with our results, some studies explicitly mention the asymmetry factor g. Figure 4 of Bohren and Barkstrom (1974) compares irradiance measurements of Liljequist (1956) to a two-stream model The Cryosphere, 7, 1803-1818, 2013 and shows that the e-folding depth modelled with spheres is greater than that measured, and that g = 0.84 (or g G = 0.68) should be used for better agreement.Meirold-Mautner and Lehning (2004) show that taking g = 0.86 instead of 0.89 (that is g G = 0.72 instead of 0.78) in their model allows a better fit with measurements.Kokhanovsky and Zege (2004) suggest g G = 0.5, based on measurements on ice particles in clouds (Garrett et al., 2001;Barkey et al., 2002), and calculations on fractal particles described in Macke et al. (1996).Under the assumption that B/(1 − g G ) is nearly independent of grain shape for snow, and equal to that for spheres, a rough estimate of g G can be obtained from B values.Based on our B measurements (excluding fresh snow which gives erroneous B), the range of g G would be 0.66-0.73.These values are in agreement with previous studies.However, in most of these studies, geometric grain size was estimated from visual observation of grains which can differ from the optical grain size that determines snow optical properties (Langlois et al., 2010).Our approach is more accurate in that it does not use snow grain size but reflectance at 1310 nm.Reflectance can be measured more accurately and less subjectively than grain size, especially for non-spherical particles.

Conclusions
The paper shows that the spherical assumption is inappropriate to model irradiance profiles in snow.Indeed, theoretical calculations as well as measurements point out that the absorption enhancement parameter B is greater for snow than for spheres.As a consequence, the geometric asymmetry factor g G in snow is probably lower than that of spheres, and modelled e-folding depth is overestimated when grains are assumed spherical, which is in agreement with observations from the literature.Grain shape is highly variable with time and location (Colbeck, 1982), hence snow cannot be systematically represented by a collection of spheres for optical calculations.To account for the grain shape in optical radiative transfer models, such as those used in weather forecast and climate models, we need to further explore the relationship between B and snow type or traditional snow grain shape (Fierz et al., 2009).To this end, future work should focus on the systematic determination of B in laboratory and in the field.In addition, a method to determine g G should be developed which will require independent measurements of specific surface area.

Fig. 1 .
Fig. 1.(a) The shapes presented in Table 1 are displayed in a B, (1 − g G ) diagram.(b) The same shapes are presented in a B(1 − g G ), B/(1 − g G ) diagram, to illustrate their characteristics in terms of e-folding depth and albedo (see Eqs. 10 and 11).

Fig. 2 .
Fig. 2. Measurements conducted at Col de Porte (29 February 2012).(a) The measured irradiance profile at λ 1 k e = 600 nm is shown with the optimal modelled profile symbolized by dark inverted triangles.The other profiles have been calculated for various values of B and incident flux equal to the optimal I 0 (600 nm).All profiles are normalized at the surface.(b) Vertical profiles of density and reflectance at 1310 nm used to compute the irradiance profiles in (a).

0 and I diffuse 0 .
It computes the downwelling and upwelling irradiance I ↓ model (z, λ) and I ↑ model (z, λ) within a plane-parallel multi-layer snowpack.α b is set to 1 here.

Fig. 4 .
Fig. 4. Black carbon content beyond which black carbon absorption represents more than 5 % of the snow absorption coefficient.The curve is deduced from Eq. (25).The dashed lines highlight the values at 600 nm and 780 nm (wavelengths used for the Alpine measurements) and the continuous lines highlight the values at 500 nm and 700 nm (wavelengths used for Dome C measurements).

Fig. 5 .Fig. 6 .
Fig. 5. (a) Spectral e-folding depth measured between 24 cm and 42 cm in the depth hoar layer at Dome C (29 Dec 2009).Modelled spectral e-folding depths are also shown for different values of B. (b) Spectral e-folding depth measured between 7.5 cm and 23 cm at Lac Poursollet (10 May 2012) in a snowpack composed of coarse grains.

Table 1 .
Values of B and g G for various geometric shapes, computed using ray tracing methods.Values for spheroids, hexagonal plates and fractals are taken from

Table 2 .
Summary of the measurements collected at Dome C and in the Alps.λ 1 wavelengths used in the algorithm described in Sect.4.1.The retrieved values of B and black carbon contents are given.

Libois et. al: Influence of grain shape on light penetration in snow 1811
Case 2: If only the spectral e-folding depth between two levels is available, the difference is 2 The Cryosphere, 7, 1803-1818, 2013www.the-cryosphere.net/7/1803/2013/Q.

Table 3 .
Determination of B from data in the literature, using Eq.(12).λ k e and λ α are in nm.R, T and I stand for reflectance, transmission and irradiance, respectively, and refer to the measurements performed in the corresponding study.The five values of B larger than 3.16 and the corresponding α values are bold.These values are questionable.