Implementation and evaluation of prognostic representations of the optical diameter of snow in the SURFEX / ISBA-Crocus detailed snowpack model

In the SURFEX/ISBA-Crocus multi-layer snowpack model, the snow microstructure has up to now been characterised by the grain size and by semi-empirical shape variables which cannot be measured easily in the field or linked to other relevant snow properties. In this work we introduce a new formulation of snow metamorphism directly based on equations describing the rate of change of the optical diameter ( dopt). This variable is considered here to be equal to the equivalent sphere optical diameter, which is inversely proportional to the specific surface area (SSA). dopt thus represents quantitatively some of the geometric characteristics of a porous medium. Different prognostic rate equations ofdopt, including a re-formulation of the original Crocus scheme and the parameterisations from Taillandier et al. (2007) and Flanner and Zender (2006), were evaluated by comparing their predictions to field measurements carried out at Summit Camp (Greenland) in May and June 2011 and at Col de Porte (French Alps) during the 2009/10 and 2011/12 winter seasons. We focused especially on results in terms of SSA. In addition, we tested the impact of the different formulations on the simulated density profile, the total snow height, the snow water equivalent (SWE) and the surface albedo. Results indicate that all formulations perform well, with median values of the RMSD between measured and simulated SSA lower than 10 m 2 kg−1. Incorporating the optical diameter as a fully fledged prognostic variable is an important step forward in the quantitative description of the snow microstructure within snowpack models, because it opens the way to data assimilation of various electromagnetic observations.


Introduction
Snow is a dynamic medium that undergoes continuous thermodynamical and mechanical processes leading to changes in its microstructure (Colbeck, 1983;Dominé and Shepson, 2002;Flin et al., 2004;Schneebeli and Sokratov, 2004).A good representation of this so-called "snow metamorphism" within snowpack models is crucial, since the snow microstructure has a significant impact on many macroscopic properties of the snowpack itself.For instance, optical properties such as the transmission of light trough the snowpack or the surface albedo and chemical exchanges between air and snow are strongly affected by the snow microstructure (Warren, 1982;Meirold-Mautner and Lehning, 2004;Domine et al., 2008).Therefore, designing snow physical models in such a way that they can accurately simulate the metamorphic processes can improve the calculation of the energy and mass budget of snow-covered surfaces (Flanner and Zender, 2006;Taillandier et al., 2007), avalanche forecasting (Brun et al., 1989(Brun et al., , 1992;;Durand et al., 1999) and the modelling of climate and of air-snow exchanges of reactive chemical species (Barret et al., 2011).
Few snowpack models incorporate an explicit representation of snow metamorphism.In the simplest physically based models, snow layers are described by their thickness, Published by Copernicus Publications on behalf of the European Geosciences Union.
temperature, density and liquid water content (Essery et al., 2013).More detailed models, in order to simulate snow aging and albedo evolution, include a representation of grain size growth.In the Community Land Model (CLM, Oleson et al., 2010), for instance, dry snow aging is represented as an evolution of the ice effective grain size, which evolves following the microphysical model described by Flanner and Zender (2006).More complex models also incorporate a notion of grain shape (Brun et al., 1989(Brun et al., , 1992;;Lehning et al., 2002).This notion can be important if the aim is to predict an estimate of the avalanche hazard (Durand et al., 1999).In the SURFEX/ISBA-Crocus (Crocus hereinafter) multi-layer snow model (Brun et al., 1989(Brun et al., , 1992;;Vionnet et al., 2012), metamorphism is currently represented by equations describing the evolution of the snow grain size, the dendricity and the sphericity.These variables, however, are semi-empirical and can neither be measured easily in the field nor directly linked to other relevant snow properties.In addition, this representation makes data assimilation efforts operating on physical properties of snow particularly cumbersome (Toure et al., 2011;Dumont et al., 2012a).
Among scalar variables describing the microstructure of snow and which can be derived from the 3-D geometry of this porous medium (Flin et al., 2003;Löwe et al., 2011), the specific surface area (SSA) has gained increased attention over the past decade.SSA is defined as the total area at the ice/air interface in a given snow sample per unit mass.This variable is impacted by snow aging in a potentially predictable way (Flin et al., 2004;Legagneux and Dominé, 2005;Flanner and Zender, 2006) and generally decreases over time, with values ranging from 224 m 2 kg −1 for diamond dust crystals (Domine et al., 2012) to less than 2 m 2 kg −1 for melt-freeze crusts (Domine et al., 2007).Moreover, SSA is a practical metric for relating the microphysical state of the snowpack to snow electromagnetic characteristics, such as snow albedo and penetration depth (Wiscombe and Warren, 1980;Flanner and Zender, 2006) and microwave behaviour (Brucker et al., 2011;Roy et al., 2013a).A rich data set of SSA measurements is now available, since SSA can be retrieved from remote sensing (Kokhanovsky and Schreier, 2009;Dumont et al., 2012b;Mary et al., 2013), computed by microtomography under both isothermal (Flin et al., 2004;Löwe et al., 2011;Schleef and Loewe, 2013) and temperature gradient conditions (Schneebeli and Sokratov, 2004;Calonne et al., 2014;Riche et al., 2013) and measured in the field using optical methods (Matzl and Schneebeli, 2006;Gallet et al., 2009;Arnaud et al., 2011).
SSA is inversely proportional to the equivalent spheres' diameter, which is the diameter of a monodisperse collection of disconnected spheres featuring the same surface area/mass ratio.The equivalent spheres' diameter is often used interchangeably with the snow optical diameter (d opt ).This allows for writing where ρ ice is the density of ice (917 kg m −3 ).Unlike grain size, a quantity with an ambiguous definition and for which it is difficult to obtain accurate values from visual inspections (Fily et al., 1997;Grenfell and Warren, 1999;Fierz et al., 2009), d opt is a well-defined variable representing some geometric characteristics of a porous medium (Giddings and LaChapelle, 1961;Warren, 1982;Grenfell et al., 1994).If the purpose is to model air-snow exchanges of chemical species, it is more convenient to use SSA, whereas d opt is more adequate for describing snow optical properties.In the following we will present results in terms of both d opt and, alternatively, SSA.
Crocus already included a description of the optical diameter, because of its impact on near-infrared albedo.However, this variable was only diagnosed from the dendricity (δ), the sphericity (s) and the snow grain size (g s ).This generates errors in the parameterisation of d opt which add up to the errors in the model internal variables.The main consequence, as stated by Morin et al. (2013), is that "this formalism hampers direct improvements of the model performance, because improving the optical diameter prediction would require improving either the relationships between δ, s and g s or the metamorphism laws acting on them".In this work we introduce an alternative approach, in which snow metamorphism within the Crocus model is now described by equations formulated in terms of the rate of change of two state variables, sphericity and optical diameter.In other words, we replaced two of the primary Crocus variables (dendricity and size) with optical diameter, turning the latter into a prognostic variable.We decided to represent the evolution of the optical diameter instead of SSA because small errors in low SSA values produce large differences in d opt (see Eq. 1), with a large impact on the optical properties.This new formalism allows for simplifying of the model by reducing from 3 to 2 the number of variables which evolve over time.Moreover, it makes it easier to implement in Crocus different parameterisations of the rate of increase of d opt , which is one of the purposes of this paper.In particular, we tested four evolution laws of dry metamorphism: the original Crocus formulation (Brun et al., 1992;Vionnet et al., 2012, B92 hereinafter) using δ, s and g s ; our new formulation (C13) using s and d opt , in which the rate of change of the optical diameter is deduced, making some simplifications, from the same equations as B92; the parameterisation proposed by Taillandier et al. (2007) (T07); and the model developed by Flanner and Zender (2006) (F06).All the above-mentioned formulations were evaluated by comparing them to field measurements.Two instruments which retrieve snow specific surface areas from infrared reflectance measurements at 1310 nm, DUFISSS (Gallet et al., 2009) and ASSSAP (a light version of POSSSUM, Arnaud et al., 2011), were used at Summit Camp (Greenland) in May and June 2011 and during the 2009/2010 and 2011/2012 winter seasons at Col de Porte (French Alps).During these field campaigns, SSA data were acquired with high vertical resolution (about 1 cm), allowing for testing of the accuracy of the different representations of dry metamorphism.
2 Metamorphism in SURFEX/ISBA-Crocus 2.1 Model overview SURFEX (Surface Externalisée, Masson et al., 2013) is the surface modelling platform developed by Météo-France.It has been designed to be coupled with atmospheric and hydrological models and contains independent physical schemes like ISBA for natural land surface, TEB (Town Energy Balance model) for urbanised areas and FLake for lakes.ISBA (Interaction Sol-Biosphère-Atmosphère) includes, in turn, several sub-modules simulating the exchanges of energy and water between the soil-vegetation-snow continuum and the atmosphere above.For snow, different scheme options are available within ISBA, the most detailed of them being Crocus (Brun et al., 1989(Brun et al., , 1992;;Vionnet et al., 2012).
Crocus is a unidimensional model able to compute the energy and mass balance of the snowpack.To this end, the vertical profile of the physical properties of snow is represented by a large number of numerical layers (Vionnet et al., 2012).The model includes a detailed description of the time evolution of the snow microstructure.To implement the metamorphism laws, a semi-quantitative formalism describing snow as a function of continuous parameters has been introduced into Crocus (Brun et al., 1992).These parameters are the dendricity δ (dimensionless, varying between 0 and 1), the sphericity s (dimensionless, also varying between 0 and 1) and the grain size g s (corresponding to the diameter of the grain, in m).Two main classes of snow types are considered by Crocus.Initial δ and s values are prescribed to every freshly fallen snow layer, depending on wind speed and temperature.In the case of low wind and cold temperature conditions, δ and s are set to 1 and 0.5, respectively.A layer is considered dendritic as long as precipitated snow crystals are still recognisable.When δ, which always decreases with snow aging, reaches 0, snow enters the non-dendritic state and the variables describing its microstructure are switched to s and g s (Vionnet et al., 2012;Morin et al., 2013).The time evolution of dendricity, sphericity and grain size follows empirical laws whose parameters were adjusted through experimental investigations.In the case of dry snow, these laws depend mostly on temperature T and temperature gradient G (Brun et al., 1989(Brun et al., , 1992;;Vionnet et al., 2012).In particular, three regimes are distinguished: weak temperature gradient (G ≤ 5 K m −1 ), middle temperature gradient (5 K m −1 < G ≤15 K m −1 ) and strong temperature gradient (G > 15 K m −1 ).In the latter case, g s increases over time fol-lowing the parameterisation described by Marbouty (1980).For wet snow metamorphism, evolution laws only depend on liquid water content (Brun, 1989).

Impact of the snow microstructure within the model
The snow microstructural variables (δ, s and g s ) are computed for every time step and for each layer within the metamorphism routine of Crocus.Other processes, however, can modify these variables.In addition, the microstructural properties are used in turn to calculate different physical quantities.In this section we describe briefly where and how, apart from the metamorphism routine, the microstructural variables are used or modified within the model.

Snow compaction
Snow layer settling due to the combined effect of the metamorphism and the weight of the overlaying layers is computed using a Newtonian viscosity law, which leaves the mass of each layer unchanged and reduces the layer thickness in proportion to density increase.Snow viscosity depends on snow density, temperature and liquid water content, but also on microstructural properties (depth hoar, for instance, has a lower compaction rate).

Snow surface albedo and solar radiation transmission through the snowpack
In order to calculate the surface albedo and the transmission of short-wave radiation through the snowpack, the optical diameter of snow was empirically derived from the snow microstructural properties based on the experimental work by Sergent et al. (unpublished).Laboratory measurements of the optical diameter using an optical method and of δ, s and g s using 2-D image analysis (Lesaffre et al., 1998) allowed the formulation of two different equations, described in Vionnet et al. (2012): in the dendritic case the optical diameter is computed as a function of δ and s, whereas in the non-dendritic case it is computed as a function of s and g s (see Fig. 1).
The albedo is then calculated from the snow properties of the two upper numerical layers by splitting the solar radiation into three separate spectral bands ([0.3-0.8],[0.8-1.5] and [1.5-2.8]µm).In the UV and visible range ([0.3-0.8]µm), albedo depends on the optical diameter and on the amount of light absorbing impurities, the latter being parameterised from the age of snow.In the infrared bands ([0.8-1.5] and [1.5-2.8]µm), albedo depends only on the optical diameter of snow (Vionnet et al., 2012).For instance, increasing the optical diameter from 1.09 • 10 −4 m (corresponding to an SSA of 60 m 2 kg −1 ) to 3.27 • 10 −4 m (corresponding to an SSA of 20 m 2 kg −1 ) leads to an albedo decrease from 0.79 to 0.67 in the range 0.8-1.5 µm.

Snowfall and grid resizing
Crocus can modify the discretisation of the vertical grid, in order to keep the number of layers below a predetermined value (typically 50).When a new snowfall layer is added to an existing snowpack, the model first prescribes specific values accounting for its microstructure.Then, if the freshly fallen snow layer and the existing top layer have similar characteristics, they are merged.The similarity between both layers is determined from the value of the sum of their differences in terms of δ, s and g s , each weighted with an appropriate coefficient ranging from 0 to 200: 0 corresponds to the case in which the same snow type is present in both layers and 200 to very different snow types.In other words, merging is only possible for layers which are similar enough in terms of snow types.If a new numerical snow layer is built from two older layers, its characteristics are calculated in order to conserve the averaged weighted optical diameter of the former layers.This ensures a strong consistency in the evolution of surface albedo.When, instead, merging is not possible, a new numerical layer is added to the existing snowpack.
A complete description of the grid resizing in Crocus can be found in Vionnet et al. (2012).

Wind drifting
Snow compaction and metamorphism due to wind drift are taken into account by Crocus (Brun et al., 1997;Guyomarc'h and Merindol, 1998;Vionnet et al., 2013).The impact of wind on surface layer properties is evaluated in three steps.First, a mobility index is calculated for each layer from its snow type and density.This index, which describes the potential for snow transport, is highest for fresh snow and tends to decrease with sintering and compaction.Mobility is then combined with wind speed to compute a driftability index.Finally, density and snow types are modified in the case of transport: in each layer, δ and g s are reduced and s is increased according to the driftability index.

Time evolution of the optical diameter within the model
In Sect.2.1 we have presented the original Crocus formulation of snow metamorphism laws, as developed by Brun et al. (1992).That formulation describing snow grains by their dendricity, sphericity and size and treating the optical diameter as a diagnostic variable of the model is called B92 hereinafter.Here we introduce a new formulation (called C13) in which the original evolution laws are re-formulated in terms of d opt and s.Indeed, both dendricity and grain size, plotted as a function of time, are monotonic functions, the former always decreasing when snow is in a dendritic state and the latter always increasing in a non-dendritic state.Thus, it is possible to replace δ and g s by a single variable, the optical diameter, which is always increasing with time and whose rate of change can be easily deduced from the same equations as B92.d opt is then turned into a prognostic variable of the model.However, since this variable alone is not sufficient to describe uniquely the snow microstructure, in C13 sphericity still remains, in order to take into account the shape of the crystals (rounding and facets).
In addition to the original formulation B92 and the new formulation C13, two other representations of dry snow metamorphism were implemented in Crocus: the parameterisation of the d opt rate of change from Taillandier et al. (2007) (T07) and that from Flanner and Zender (2006) (F06).For C13, T07 and F06, in the case of wet snow we used the original B92 equations, reformulated in terms of d opt (Table 1).Formulations C13, T07 and F06 are described in detail in the following sections.

Formulation C13
In formulation C13, all original equations of B92 including δ, s and g s were re-written in terms of d opt and s.In practice, δ and g s were replaced by expressions that are functions of only d opt and s.Based on the original relationships linking δ, s and g s to d opt described in Sect.2.2.2, we replaced the dendricity by δ = and the grain size by where d opt and g s are expressed in m and α is equal to 10 −4 m.Since δ and s vary between 0 and 1, Eq. ( 2) establishes that d opt values for dendritic snow range from 0.1 mm, corresponding to an SSA value of 65 m 2 kg −1 , to 0.4 mm, corresponding to an SSA of 16 m 2 kg −1 (see Fig. 1).The transition from a dendritic to a non-dendritic regime, which occurs when δ reaches 0 in B92, was also re-formulated in terms of d opt and s: Thus, snow enters its non-dendritic state if d opt grows beyond a certain threshold.The higher the sphericity value, the easier it will be for d opt to exceed this threshold.The snow optical diameter was introduced in all Crocus routines described in Sect.2.2.In most cases, this leads to a significant simplification of the equations.For example, using the optical diameter as a prognostic variable allows for computing of the snow albedo directly, which depends explicitly on d opt , and the microstructural characteristics of merged layers, which are calculated in order to conserve the averaged weighted d opt of the former layers.Elsewhere, the change of variable is not as trivial.This is the case, for instance, for wind drifting.When a new fresh snow layer is created, its optical diameter is set by default to 10 −4 m.This value corresponds to the combination of d = 1 and s = 0.5 in B92, as can easily be seen using Eq.(2).In the case of wind drifting, however, the d opt value is modified depending on the wind speed.Therefore, the evolution rates of the snow microstructural properties caused by the snow drifting (reported in Table 3 of Vionnet et al. (2012) in their B92 formulation) here have to be re-written for the optical diameter as follows: where t is time expressed in hours, τ is computed from the driftability index (which depends on the wind speed and the microstructural properties of snow, see Sect.2.2.4) and represents the time characteristic for snow grain changes under wind transport, and δ can be derived through Eq. ( 2).In the case of no drifting, τ tends to infinity and Eq. ( 5) leads to unchanged values for d opt (Brun et al., 1992(Brun et al., , 1997;;Vionnet et al., 2012).
The most significant differences between C13 and B92 concern the re-formulation of the dry metamorphism laws.Combining the original Crocus rate equations for the snow microphysical variables (reported in Table 1 of Vionnet et al., 2012) and Eqs.(2-4), we obtained new parametric laws describing the rate of change of d opt .These equations are reported in Table 2 along with the rate equations for s, which are identical to those of B92.The rate of change of d opt (always increasing) and s (either increasing or decreasing) are a function of the vertical temperature gradient (G, in K m −1 ), the temperature (T , in K) and the time (t, in days).Six cases are distinguished, depending on the G value (weak, middle and strong temperature gradient) and the regime (dendritic and non-dendritic snow).When G > 15 K m −1 and s = 0, d opt of non-dendritic snow increases over time following the parameterisation of Marbouty (1980), which allows for predicting of depth hoar growth rate.In the original B92 formulation, non-dendritic d opt became a function of g s only when the latter exceeded an empirical threshold set to 8 • 10 −4 m.In C13 we do not have any information about g s and therefore we have removed this threshold.This means that in the case of a strong temperature gradient metamorphism, C13 will lead to d opt values higher than those of B92.Aside from that difference, formulations C13 and B92 are supposed to give the same results (see Sect. 2.3.4), as C13 consists in a reformulation, in terms of optical diameter, of the original metamorphism laws of B92 expressed in terms of dendricity and snow grain size.

Formulation T07
Taillandier et al. ( 2007) performed several experiments during which the SSA was measured under isothermal and temperature gradient conditions.Based on this data set, they proposed a parameterisation of the rate of decay of SSA (formulation T07), whose implementation in Crocus is described in this section.
and ∂s ∂t = 0 where A, B and C are fitted coefficients that depend on the actual temperature of the layer (T ) and its initial snow SSA (SSA 0 ).No density effect is predicted by those equations, since experiments performed by Taillandier et al. (2007) failed to detect a significant influence of that variable.For SSA 0 , in order to be consistent with the maximum SSA value allowed by formulation B92, we used: Equation (6a and b) can be merged, using a hyperbolic tangent form, into a single equation where SSA evolution depends continuously on temperature gradient.For this aim we defined two new coefficients as follows: This allowed a unified equation to be written which formulates SSA as a function of the temperature gradient and the temperature: Equation ( 9) computes SSA using Eq.(6a) when ∇T < 5 K m −1 and Eq.(6b) when ∇T > 15 K m −1 .When the temperature gradient is between these two values, a weighted average of SSA ET and SSA TG is used.From Eq. ( 9), with a second-order Taylor series development we can finally obtain the SSA decay: Two main issues had to be overcome in order to implement formulation T07 in Crocus.The first one concerned the temperature appearing in Eq. (6a and b), which should be, according to Taillandier et al. (2007), the mean temperature of evolution of the snow layer during the period of interest.Since we do not know at each time step the full temperature history of the numerical layers, we decided to use their actual temperature.This approach makes sense, since the actual temperature is used to compute changes in SSA at each model time step, whereas the full temperature history of the layers is already contained in the current SSA value (SSA(t) in Eq. 10).The second issue was the fact that in T07 the SSA rate of change depends, amongst other variables, on the age of the layer.This does not constitute a problem for one-layer snow models such as CLASS (Canadian   (Taillandier et al., 2007).Layer 1 and layer 2 are created at different times to take into account new snowfalls.At t = 30 days, these two layers aggregate and give rise to layer 3. The age of the newly formed layer has to be chosen arbitrarily and this has an impact on the SSA rate of change.
LAnd Surface Scheme), in which Roy et al. (2013b) have already implemented parameterisation T07 successfully.In multi-layer snow models, however, when two or more layers merge the age of the new-formed layer has to be chosen arbitrarily.For instance, we can decide to keep the age of either the younger or the older of the former layers, or a mean of them.All choices may be equally valid and lead to different SSA rates of change (Fig. 2).In our simulations, we computed the age of the new layer by making a weighted average of the ages of the former layers.The weight chosen was the same as that used to compute the new optical diameter (see Sect. 2.2.3).Formulation T07 uses a logarithmic function to fit the decreasing trend of SSA.Legagneux et al. (2004) showed that this approach leads to a good approximation of a more general evolution law of SSA that can be easily derived from grain coarsening theories.This approximation, however, only applies for times less than 150 days and predicts that the SSA will become negative after about a year (Taillandier et al., 2007).For this reason, in order to avoid unrealistic values, in our simulations the minimum SSA was forced to 8 m 2 kg −1 .This allows the parametric laws of T07 to be adequate for most applications to seasonal snowpacks, in which SSA do not go, in general, below that value (Taillandier et al., 2007).In contrast, formulation T07 was not designed to simulate the decreasing trend of SSA over time periods longer than a few months.In light of these considerations, we decided to evaluate this formulation only on the seasonal snowpack at Col de Porte and not on the permanent snowpack at Summit.

Formulation F06
Flanner and Zender (2006) developed a model (referred to as F06 here) able to simulate diffusive vapour flux amongst collections of disconnected ice spheres of different radii.In this model, limited to dry snow, the optical diameter can be computed for any combination of snow density, temperature, temperature gradient and initial size distribution.A parametric equation, which fits model output and is computationally inexpensive, was included in the Community Land Model by Oleson et al. (CLM, 2010).Here we followed the same approach by incorporating this equation into Crocus.
The optical diameter was computed as follows: where: r is the optical radius (in µm), r 0 is the initial optical radius for fresh snow (in µm) and t is time (in hours).The initial rate of change of optical radius dr dt 0 and the coefficients τ and κ depend on snow density, temperature and temperature gradient and can be retrieved from look-up tables www.the-cryosphere.net/8/417/2014/The Cryosphere, 8, 417-437, 2014 (Flanner and Zender, 2006).The domain covered by these tables includes densities ranging from 50 to 400 kg m −3 , temperatures ranging from 223.15 to 273.15 K and temperature gradients ranging from 0 to 300 K m −1 .The initial SSA can be set to 60, 80 or 100 m 2 kg −1 .In order to be consistent with the maximum SSA value allowed by formulation B92, we chose an initial SSA of 60 m 2 kg −1 , which corresponds to r 0 = 50 µm in Eq. ( 12).

Comparison of the different formulations
The SSA rates of decay predicted by the above-described formulations are shown in Fig. 3.All curves were plotted for a constant temperature of −20 • C and a constant density of 200 kg m −3 .Different temperature gradient values ranging from 0 to 50 K m −1 were tested; in the lower right panel, ∇T was set to 0 K m −1 until t = 22.5 days and to 50 K m −1 thereafter.In order to have a similar initial SSA for all schemes, about 60 m 2 kg −1 , B92 was initialised with d = 0.964 and s = 0.5 and C13 with d opt = 109 µm and s = 0.5.Regardless of the formulation chosen, the most rapid snow aging is always produced by the combination of low density, high temperature and large temperature gradient.However, on closer inspection the comparison of the four metamorphism formulations shows some differences.B92 and C13, for instance, display a discontinuous derivative when snow enters the non-dendritic state.This discontinuity, which has no physical meaning and comes from the empirical parameterisation of the rate equations, is not present in T07 and F06.Indeed, in terms of SSA rate of change there is no distinction for T07 and F06 between the dendritic and nondendritic regimes; this distinction remains true only for the time evolution of sphericity, which is identical to that of B92 for all formulations.In addition, it is easy to see that B92 and C13 differ only when ∇T is high (> 15 K m −1 ).In this case, in the C13 formulation SSA starts decreasing following the parametrization of Marbouty (1980) as soon as the nondendritic state is reached; B92, instead, makes SSA decrease only when g s >8 • 10 −4 m.In Fig. 4, the formulations of the SSA rate of change were compared to cold room measurements.Flin et al. (2004)  ).Results reveal that, albeit with some differences, all formulations perform similarly.Formulation F06 reproduces the experimental data slightly better than the other formulations, whereas B92 and C13 generally perform not quite as well, mostly because of their change in regime between dendritic and non-dendritic snow, which does not match the observed behaviour of SSA decrease.Analogous results are expected to be found when these different representations of the SSA rate of change are implemented directly in Crocus.

Materials and methods
In this section we describe the sites where measurements were carried out, the meteorological forcing data used to run the simulations, the snow data set acquired to evaluate the accuracy of the metamorphism formulations, the model runs and the statistical metrics used to compare Crocus outputs to observations.

Summit Camp
A two-month field campaign took place in May and June 2011 at Summit Camp (72 • 36 N, 38 • 25 W).This research station is located at the peak of the Greenlandic ice cap, at 3210 m a.s.l.(http://www.summitcamp.org/).In this area, snowfall can occur in all seasons (Albert and Hawley, 2000) and the accumulation rate is not seasonally uniform.Wind conditions can change dramatically during the year (Albert and Shultz, 2002) and have a strong impact on snow surface layers.During our campaign, the air temperature was always negative and no liquid water was ever found within the snowpack.Several wind drifting events occurred, with wind speeds up to 15 m s −1 .

Col de Porte
Snow stratigraphies and SSA measurements were performed at Col de Porte on a weekly basis during the winter seasons of 2009/2010 and 2011/2012.This experimental site, located in the French Alps at an elevation of 1325 m a.s.l., has been used for snow research for over 50 years, with a continuous record of snow and meteorological variables at an hourly time resolution since 1993 (Morin et al., 2012).Due to the relatively low altitude of the site, snow melting and rain can occur any time in the winter.In addition, even if strong wind events are generally sporadic, spatial heterogeneities in terms of snow height, snow water equivalent (SWE) and grain properties are commonly observed.The great time variability of its seasonal snowpack makes the mid-altitude Col de Porte site particularly well suited to test the different metamorphism laws implemented in Crocus.

Meteorological data
In order to run SURFEX/ISBA-Crocus simulations, a complete meteorological forcing has to be provided to the model.In its basic configuration, Crocus requires the following driving variables: air temperature and specific humidity at 2 m above ground, wind speed at 10 m above ground, snowfall and rainfall rates and incoming short-wave and long-wave radiations.
For simulations at Summit, input variables were provided by the ERA-Interim reanalysis (Dee et al., 2011) with a 3 h time step.This data set extends back to 1979 and has a spatial resolution of about 80 km, with 60 vertical levels up to 0.1 hPa (about 64 km).Since meteorological conditions over the upper Greenlandic ice cap are quite homogeneous, we used the 80 km grid cell nearest to the Summit Camp location, without downscaling the data to a smaller numerical grid.For the period of our field campaign, ERA-Interim data were replaced by in situ measurements acquired on an hourly basis by the Summit Automatic Weather Station (AWS, Steffen et al., 1996).Air temperature and relative humidity were measured using a Type-E thermocouple (estimated accuracy of 0.1 • C) and a Campbell Sci.CS-500 (accuracy of 10 %), respectively.Wind speed was recorded with a RM Young propeller-type vane (estimated accuracy of 0.1 m s −1 ) and short-wave incoming radiation was measured using a Li Cor photodiode (nominal accuracy of 15 %).No measurement of precipitation rates and incoming long-wave radiation were performed by AWS. Figure S1 shows a comparison between AWS measurements and ERA-Interim reanalysis for four different variables during the period April-July 2011.
For simulations at Col de Porte, we used in situ meteorological observations with a time step of 1 h, spanning the period from 1 August 1993 to 31 July 2012.A detailed description of the procedure followed to generate such an interrupted driving data set is presented in Morin et al. (2012).Briefly, in situ observations were used in the period of the year concerned with snow on the ground (i.e. from mid-September to mid-June, approximately).A manual procedure for data inspection and correction was followed, including filling of episodic data gaps using redundant sensors.The summer period may be used for site maintenance and not all sensors operate continuously during this time.Hence, during such periods, to which the snowpack conditions in the following winter season are marginally sensitive, meteorological data were provided by the SAFRAN meteorological analysis (Durand et al., 1993).

Snow data
Measurements of several snow properties were performed at Summit Camp and Col de Porte in order to evaluate the different representations of the optical diameter of snow implemented in SURFEX/ISBA-Crocus (see Table 3).
At Summit, density was measured every day, down to about 80 cm and with a vertical resolution better than 4 cm, by sampling the snow with a 250 cm 3 rectangular steel cutter (Fierz et al., 2009;Conger and McClung, 2009).At the same time, snow specific surface area measurements were performed using DUFISSS (DUal Frequency Integrating Sphere for Snow SSA measurement, Gallet et al., 2009), an integrating sphere which allows for retrieving of SSA from infrared reflectance measurements.In the top 10 cm we kept a vertical resolution of 1 cm, whereas below that depth and down to about 80 cm a 2 to 4 cm resolution was used.The snow height was monitored by the AWS, using a Campbell SR-50 with 1 mm precision.A more in-depth description of the snow data set acquired at Summit can be found in Carmagnola et al. (2013).
At Col de Porte, 14 and 16 vertical profiles of density and SSA were acquired during the winter of 2009/2010 and 2011/2012, respectively.Density was recorded for each snow layer with cylindrical steel cutters and SSA was measured in 2009/2010 using DUFISSS, whereas in 2011/2012 highresolution measurements of SSA were obtained using ASS-SAP (Alpine Snowpack Specific Surface Area Profiler).This instrument is the handy and lightweight version of POSS-SUM (Arnaud et al., 2011).Its working principle is similar to that of DUFISSS (retrieval of SSA from infrared reflectance measurements), the main difference being that it allows the continuous acquisition of SSA profiles down to about 90 cm, with a vertical resolution of 1 cm.In addition to density and SSA measurements, manual weekly observations of snow height and SWE were carried out.The daily integrated surface albedo was computed from the ratio between incoming and reflected short-wave fluxes, measured using radiation sensors (Morin et al., 2012).Rather than computing an average of the hourly albedo estimates, which would be impacted by noise induced by low radiation conditions early and late in the day, the ratio of the sums of total reflected and incoming radiations, respectively, was taken for each day, featuring enough incoming radiation.This approach mitigates the impact of low radiation conditions and puts more weight on the time of the day when maximum radiation occurs.These data can be safely compared to albedo simulations at noon since diurnal variations of snow albedo, induced by variations of the solar zenith angle, are currently not represented in Crocus.

Model runs
At Summit, the simulations were run from 1 February 1979 to 1 July 2012.The output time step was set to 6 h and the total initial snow height was set to 10 m.The snowpack was initialised with eight layers, prescribing thinner and less dense layers near the surface.The choice of the profile used to initialise the simulations has a very small impact on the final results of the model more than 30 yr later.The impact of replacing ERA-Interim forcing data by in situ measurements during our field campaign period is also negligible.
At Col de Porte, the simulations were run from 1 August 1993 to 31 September 2012, with an output time step of 6 h.We used, for the whole time period, the same initialisation procedure, consisting of assigning to all ground layers the mean annual temperature at Col de Porte.This procedure ensures that the temperature data of the simulated uppermost ground layers attain reasonable equilibrium with environmental conditions after 16 and 18 yr of simulation (for the 2009-2010 and 2011-2012 snow seasons, respectively).

Definitions of metrics
In order to evaluate the agreement between observations and simulations, we define the quantity SSA , which represents the root mean square deviation (RMSD) between measured and simulated SSA for a given date: where SSA obs,h and SSA sim,h are the measured and simulated SSA, respectively, and N h is the total number of considered SSA values for a given vertical profile.In practice, to compute SSA both measurements and simulations were interpolated on a 1 mm vertical grid and the difference SSA obs,h −SSA sim,h was calculated for every SSA obs,h value.
The same quantity can be expressed in terms of d opt .Every single measured and simulated SSA value was then converted into respectively d opt obs,h and d opt sim,h , in order to compute the RMSD for the optical diameter: 4 Results

Field measurements
At Summit, the upper metre of the snowpack was dominated by faceting rounded grains with density values ranging from 280 to 350 kg m −3 , interspersed with hard wind slabs with densities up to 400 kg m −3 .At the surface, in the case of fresh snow precipitation or surface hoar formation, density could be significantly lower (120-140 kg m −3 ).SSA reached 65-75 m 2 kg −1 (corresponding to optical diameters of 109 to 93 µm) at the surface when rime, fresh snow or surface hoar were present.For layers deeper than 60 cm, SSA was about 20 m 2 kg −1 (optical diameter of 327 µm).Even if the cold temperatures of the Arctic ice cap make the metamorphic processes slow, during the campaign we observed a decrease in SSA over time, along with a slight density increase.Indeed, if we consider the depth between 2 and 15 cm, our first eight vertical profiles (from 5 May to 19 May) have a mean density value of 315 ± 13 kg m −3 and a mean SSA value of 37 ± 3 m 2 kg −1 .These values become, respectively, 335 ± 20 kg m −3 and 23 ± 3 m 2 kg −1 for the last eight profiles (11-25 June).In terms of snow height, results show that during our field campaign, despite frequent small amounts of precipitation, no significant accumulation occurred.At Col de Porte, density spread over a wide range of values, from 65 kg m −3 for fresh snow to 450 kg m −3 for melt clusters.In the same way, measured SSA values ranged from 5 m 2 kg −1 (corresponding to an optical diameter of 1.3 mm) for melt-refrozen forms up to 80 m 2 kg −1 (optical diameter of 73 µm) for precipitation particles fallen in weak-wind conditions and at low temperatures.The maximum snow height measured during the 2009/2010 winter season was just above 1 m and the maximum measured SWE was about 300 kg m −2 .In 2011/2012, these values were respectively 1.24 m and 390 kg m −2 .The daily integrated surface albedo values ranged from 0.95 in the presence of freshly fallen snow to less than 0.5 for dirty and old snow during melting periods.

Summit camp
Figure 5 shows the time series of snow height and SWE over the whole simulated period.After about 33 yr of simulations, the different metamorphism formulations (see Sect. 2.3) give similar results, with differences less then 4 % in snow height and less than 3 % in SWE.A comparison with the measured data for the period January-July 2011 (Fig. 5c) shows that the simulated snow height is not always able to capture the observed variability: the general trend is reproduced by the model, but the variations at the daily scale are not.This can be mostly explained by the spatial variability of the snowpack due to the effect of wind, and in particular to the eventdriven deposition of snow, which is frequent in polar regions (Groot Zwaaftink et al., 2013).
The vertical profiles of density and SSA simulated at Summit from 1 April to 31 July 2011 using the different metamorphism formulations are shown in Fig. 6.In all cases it is possible to follow the density increase and the SSA decrease over time within the different layers, which tend to aggregate when their properties become similar.Changing the evolution over time of the optical diameter can in principle have an impact on the simulated density (see Sect. 2.2.1).This impact, however, is generally low and all density profiles look similar.In terms of SSA the discrepancies are more important.In Fig. 7 we show the difference between the SSA profiles obtained with formulations C13 and F06 and the SSA profile obtained with formulation B92.Since every formulation leads to a different snow height and snowpack layering, in order to calculate this difference all SSA profiles were previously interpolated over a 1 mm vertical grid; then, starting from the snow surface, SSA values from B92 were subtracted from those from C13 and F06 at every 1 mm depth.
Results from C13 differ from those from B92 by less than 5 m 2 kg −1 over the entire profile considered, which is not surprising since C13 is designed after B92.Also, SSA val- ues computed with F06 look consistent with those obtained with B92, except for the upper layers.Contrary to the other formulations, in F06 the rate of change of the optical diameter also depends on the density and makes light layers evolve more rapidly.
A comparison between simulated and measured density and SSA down to 0.8 m is shown in Fig. 8. Profiles refer to 10 May 2011, when surface hoar was present in the upper 1 cm and the rest of the snowpack was made up of a layered system of hard wind slabs interspersed with faceting rounded grains.Since their rules of aggregation depend on the snow properties (see Sect. 2.2.3), the numerical layers simulated using B92, C13 and F06 have different thicknesses.Simulated density values range from 180 kg m −3 at the surface to more than 300 kg m −3 for deep layers.Differences between simulated densities are generally lower (within 10 %) for the upper layers, down to about 40 cm.The overall features of the observed density profile are captured, whereas the observed SSA profile is not completely reproduced: measured and simulated SSA decrease with depth, but at different rates.In addition, it can be noticed that the differences between parison in terms of SSA between the metamorphism formulations implemented into order to compare results obtained using different formulations, all SSA profiles were over a 1 mm vertical grid.The plots show the difference between the new formulations d F06 (bottom) and the original formulation B92, for simulations run at Summit from 1 uly 2011.different simulations are significantly smaller than those between simulations and observations.
The time evolution of the surface SSA at Summit is presented in Fig. 9.In order to obtain the SSA of the top 1 cm, all simulations were interpolated over a 1 mm vertical grid; then, an exponentially weighted average of the SSA values within the top 1 cm of snow was made (Mary et al., 2013).The general behaviour of the observed SSA is well captured by the simulations.Sometimes, however, the measured SSA increased over time, whereas the simulated SSA tended to decrease.This always occurred in conjunction with strong wind events (grey bands represent periods during which the AWS measured wind speed > 9 m s −1 ).These events were only strong enough to reduce the rate of decrease of the simulated SSA but not to make it increase.Observed surface SSA values, on the contrary, show an increase, since measurements were probably performed on wind-drifted snow.Moreover, it is clear that the simulated SSA values underestimate the observations.In the original B92 formulation, indeed, the maximum allowed SSA value was set to 65 m 2 kg −1 (see Eq. 2).This constitutes a limit for the simulations at Summit, where the surface SSA can easily reach 75 m 2 kg −1 .The new metamorphism formulations using the optical diameter as a prognostic variable allows for reduction of the discrepancies between simulations and observations.In Fig. 9b, for instance, formulations C13 and F06 were run setting the maximum SSA value to 80 m 2 kg −1 and show an improved agreement with the observations.Even in this case, however, the observations are not perfectly reproduced, meaning that our understanding of the processes involved in the SSA decrease over time is not complete.

2009/2010 winter season
Figure S2 shows the snow height, SWE and surface broadband albedo simulated by Crocus over the 2009/2010 winter season.The observed values are reported as well.Model results are generally similar, except during the melting period, when T07 gives higher height, SWE and albedo values.RMSD between simulations and observations are larger than those between simulations themselves, with values reaching 0.10 m in snow height, 40 kg m −2 in SWE and 0.12 in albedo.These results in terms of RMSD are reasonable for this site and consistent with previous studies (Vionnet et al., 2012;Essery et al., 2013;Morin et al., 2013).
The vertical profiles of density and SSA simulated at Col de Porte in 2009/2010 using the four metamorphism formulations are presented in Fig. S3.As for Summit, differences   in density are negligible, whereas those in SSA are more significant.The difference in terms of SSA between results obtained with formulations C13, T07 and F06 and results obtained with formulation B92 was computed by interpolating all SSA profiles over a 1 cm vertical grid (Fig. S4).Formulation C13 gives, as expected, results similar to those obtained with B92, with differences lower than 8 m 2 kg −1 .Right after the snowfall, SSA from T07 are generally lower than the corresponding values from B92; then, within a few days, they become higher (up to 20 m 2 kg −1 ) during several weeks and finally match the SSA from B92 after the sharp transition at the beginning of the melt period.It can also be noticed that for the dry and warm conditions that occurred in the middle of March 2010 the surface SSA values obtained using T07 remain larger than the corresponding values simulated using other formulations.This results in higher albedo, which leads in turn to higher snow height and SWE (see Fig. S2).As at Summit, SSA computed with F06 is consistent with that obtained with B92, except for the recent upper layers, whose SSA evolves more rapidly using the F06 parameterisation.Regardless of the formulations considered, the smaller differences are found during the melt period, not only because in this case the SSA values are generally lower, but also because all representations include the same formulation for wet metamorphism (see Table 1).
A comparison between simulated and observed density and SSA profiles at Col de Porte on 11 February 2010 is presented in Fig. S5.The total snow height measured that day was 1.04 m, the record for the 2009/2010 winter season.A 1 cm-thick layer of precipitation particles was present at the surface and was sitting on about 35 cm of decomposing and fragmented snow, the rest of the snowpack being made up of melt forms.Even if all formulations underestimate the snow height, the simulated density and SSA values follow the same pattern as the observations.For the density, the simulated values range from 80 kg m −3 at the surface to about 350 kg m −3 at the bottom and a quantitative match with the observed values is attained.For the SSA, the range of variations and the vertical layering are not well reproduced.
Figure S6 shows that all formulations perform similarly in terms of snow type and liquid water content.The only significant differences appear when the temperature gradient is between 5 and 15 K m −1 .In this case, B92 and C13 make the SSA decrease faster than T07 and F06 in a time period ranging from about 10 and 60 days since snowfall (see Fig. 3b).This is the main reason why B92 and C13 simulate the presence of faceted crystals between January and February, whereas T07 and F06 indicate, for the same period, the presence of decomposing and fragmented snow.Since visual observations between January and February revealed the coexistence, within deep layers, of both faceted and decomposing crystals, it is not easy to determine which formulation matches the observed snow profiles better in terms of snow type.In fact, since the very notion of grain types represents a discontinuous evolution of a continuous process, it is inevitable that thresholds between types vary slightly between formulations, and even different observers would place the limit between types differently.Other comparisons performed at Summit and at Col de Porte during 2011/2012, however, showed that F06 seems to reproduce more accurately the observed snow types.The differences during the melt period, as well as those in liquid water content, are negligible, since the laws of metamorphism are the same for all formulations in the case of wet snow.melting period.If we consider snow height and SWE, RMSD between simulations and observations are slightly lower than those for 2009/2010 (0.07 m and 30 kg m −2 , respectively); for the albedo, instead, RMSD are higher (between 0.14 and 0.15).Figures 11-12 present, respectively, the simulated vertical profiles of density and SSA and the difference in terms of SSA between results obtained with different metamorphism formulations.All the considerations we made for the 2009/2010 winter season still apply here.In particular, in B92 (Brun et al., 1992), C13 (this work), T07 (Taillandier et al., 2007) andF06 (Flanner andZender, 2006).39 B92 (Brun et al., 1992), C13 (this work), T07 (Taillandier et al., 2007) and F06 (Flanner and Zender, 2006).
terms of density the discrepancies are negligible and in terms of SSA the results obtained using formulations C13 and B92 are similar.Moreover, SSA values computed with F06, even if generally consistent with those obtained with B92, tend to be lower in the upper layers.Lastly, the fact that during the melt period the snow height, SWE and albedo simulated using T07 remain higher than the other model results (see Fig. 10) can be explained by the large SSA values given by T07 during the dry and warm period between March and April 2012.A comparison between simulated and observed profiles at Col de Porte on 6 February 2012 is presented in Fig. 13.That day, the measured snow height was 1.09 m and the snowpack was mainly made up of decomposed forms in the top 5 cm, faceted crystals between 5 and 20 cm and melt forms further   down.All formulations estimate the snow height well and the simulated density values are in agreement with the observations.The simulated SSA profiles are almost flat within the melt form layers and are able to capture the separation, at about 0.9 m above the ground, between older and more recent snow layers.However, some differences between model results appear at the surface: in the decomposed form layer SSA values from F06 are lower, whereas in the faceted crystal layer SSA values from T07 are higher.

Quantitative comparison between simulations and observations
The visual comparison between observed and simulated SSA profiles was complemented and consolidated by computing the quantity SSA (see Sect. 3.5).SSA values for Col de Porte were computed with and without stretching vertically the simulated profiles, in order to match the measured snow height (in practice, according to Morin et al. (2013) the thickness of each numerical layer was linearly scaled with the ratio between observed and simulated total snow height).In the same way, correlation statistics were also performed in terms of optical diameter by computing the quantity d opt .
Due to the inverse relationship between SSA and optical diameter (Eq.1), values computed with one of these two related variables are sensitive to the model performance in a different way.Good predictions for high SSA values, for instance, have more impact on SSA than on d opt ; in contrast, statistics expressed in terms of optical diameter will favour the model performances for the lower range of SSA values (Mary et al., 2013;Morin et al., 2013).
The results for SSA and d opt are presented in Fig. 14.The performance in terms of SSA of the simulations using B92, C13 and F06 is similar for Summit, with median values ranging from 8.7 to 10.1 m 2 kg −1 , and for Col de Porte in 2009/2010, with median values ranging from 7.5 to 10.1 m 2 kg −1 .For the 2011/2012 Col de Porte winter season, median values of SSA are lower, between 4.8 and 6.5 m 2 kg −1 .In terms of optical diameter, B92, C13 and F06 give median values of d opt between 0.09 and 0.12 mm for Summit, between 0.36 and 0.42 mm for Col de Porte in 2009/2010 and between 0.55 and 0.71 mm for Col de Porte in 2011/2012.Summarising, we can state that, regardless of the formulation chosen, the high SSA are estimated better for the 2011/2012 season at Col de Porte than for the other data sets, whereas the low SSA are better simulated at Summit.SSA simulated using formulation T07 are consistent with those of the other formulations and d opt are even lower, with median values of 0.35 and 0.34 mm for the 2009/2010 and 2011/2012 Col de Porte winter seasons, respectively.Lastly, it could be noticed that the vertical stretching of the simulated profiles does not improve significantly the agreement between the measured and observed values of SSA and the optical diameter.

Discussion
In the SURFEX/ISBA-Crocus multi-layer snowpack model, the snow metamorphism can now be described in terms of the rate of change of d opt , the optical diameter of snow.This makes Crocus different from any other existing detailed snowpack model.In SNOWPACK (Lehning et al., 2002), for instance, the snow microstructure is characterised by the grain geometric radius, whose growth is driven by the water vapour gradient between ice matrix and pore space, and by the dendricity and sphericity, which follow the same evolution laws used by Crocus; in addition, SNOWPACK also incorporates a representation of bonds between grains.The SMAP model (Niwano et al., 2012) employs the same variables and the same equations as those used in SNOW-PACK.The implementation of the optical diameter as a completely prognostic variable constitutes a different approach compared to that followed in previous works, in which SSA was estimated from other prognostic variables (Jacobi et al., 2010;Morin et al., 2013;Roy et al., 2013b).This enables Crocus initialisations with an observed SSA profile and makes it possible to run the model using various parameterisations of SSA decay.
The original Crocus metamorphism formulation (B92), based on empirical evolution laws of grain dendricity, sphericity and size, already led to a satisfactory agreement between simulated and observed SSA values, considering that the accuracy of the SSA measurements is estimated to be about 10 % (Gallet et al., 2009;Arnaud et al., 2011).Using B92, we found median values of SSA lower than 10 m 2 kg −1 both at Summit and at Col de Porte.The new formulation using sphericity and optical diameter (C13), in which the rate of change of d opt is deduced from the same equations as B92, also leads to maximum SSA median values of less than 10 m 2 kg −1 .The fact that B92 and C13 give quite similar results in terms of SSA means that the optical diameter has been integrated successfully into Crocus as a prognostic variable.In terms of SSA , the model devel- oped by Flanner and Zender (2006) (F06) performs equally to the other above-mentioned formulations, with results particularly close to those of C13.The only significant difference stands out for the low-density layers at the very near surface.Indeed, F06 makes the SSA decrease faster in the case of low-density layers with non-zero temperature gradients (Fig. 3d in Flanner and Zender, 2006).This dependence on snow density does not appear in the other formulations.Under isothermal conditions, the rate of SSA change is independent of snow density in all formulations.This www.the-cryosphere.net/8/417/2014/The Cryosphere, 8, 417-437, 2014 is consistent with the recent study of Schleef and Loewe (2013), who found that the rate of decrease of the SSA under isothermal conditions, computed on micro-tomographic images, was independent of the density.At Col de Porte, the parameterisation proposed by Taillandier et al. (2007) (T07) gives median values of SSA that never exceed 8 m 2 kg −1 and are consistent with those of the other formulations.This means that, regardless of the formulation chosen, high SSA values are simulated well.During melt periods, SSA values from T07 remain relatively higher than those obtained using B92, C13 or F06.This leads to higher snow height, SWE and albedo values and also results in a better agreement between measured and simulated optical diameters.
Apart from some minor differences, formulations B92, C13, T07 and F06 lead to similar results in terms of simulated SSA.The statistical agreement between measured and simulated SSA profiles is rather satisfactory and comparable to previous studies.Morin et al. (2013), for instance, carried out simulations using Crocus and two parameterisations of snow SSA, one determined from density and snow type and the other one simply derived from the internal computation of the optical radius in Crocus.In both cases, comparing simulations with SSA measurements in an Alpine environment, they found RMSD values of 10 m 2 kg −1 .Jacobi et al. ( 2010) implemented a parameterisation in Crocus in which the SSA was estimated from other snow variables using the equations of Taillandier et al. (2007).Comparing their results against measurements in a dry subarctic snowpack at Fairbanks (Alaska), they observed generally good agreement between simulated and measured SSA, with RMSD ranging from 8 to 11 m 2 kg −1 .These values are close to our findings.
Density and optical diameter are not sufficient to characterise the snow microstructure uniquely.A notion of shape has to be introduced to account for different degrees of roundness from the entirely angular crystals, such as depth hoar, to the mostly rounded grains, such as melt forms.This need for a third variable to describe the snow microstructure completely can also be made clearer by looking at the two-point correlation function of the microscopic density.Indeed, the first orders of the expansion of this function at the origin are linked not only to the volume fraction and the surface area per unit volume, but also to the curvature (Torquato, 2002;Löwe et al., 2011).Moreover, a good characterisation of the crystal shapes in snowpack models has proven to be important for several applications, such as the evaluation of the stability of the snow layers for avalanche forecasting (Durand et al., 1999) and the study of the transmission of light through snow (Meirold-Mautner and Lehning, 2004;Libois et al., 2013).For this reason, the notion of sphericity present in the original Crocus formulation still remains in the current version of the model.However, in the same way in which we replaced two semi-empirical quantities (dendricity and size) with the optical diameter, in the future it is desirable to replace the sphericity with some other fully fledged physical variable easily measurable in the field.A possible candidate to discriminate between different grain shapes seems to be the ratio of the vertical to the horizontal component of the thermal conductivity.This anisotropy of the thermal conductivity ranges from 0.7 for rounded grains up to 1.5 for faceted crystals (Calonne et al., 2011).Unfortunately, for now this quantity can only be computed by micro-tomography and only a very limited data set is available.Recently, Libois et al. (2013) studied the impact of the grain shape on the macroscopic optical properties of snow.In their approach, the snow particle shape is fully defined by two quantities, the absorption enhancement parameter, quantifying the enhancement of absorption due to lengthening of the photon paths within the grain, and the geometric diffusion term, accounting for the angular anisotropy in diffused light.These quantities are both independent of the size and the former can be retrieved from reflectance and extinction measurements.They might constitute a possible alternative to the sphericity to describe the grain shape.When the sphericity will be replaced by one or more quantities which can be measured readily in the field and linked objectively to other relevant snow properties, then the snow microstructure in Crocus will be characterised in an entirely physical manner.

Conclusions
A new approach to describe snow metamorphism has been implemented in the SURFEX/ISBA-Crocus model.The optical diameter of snow, so far estimated indirectly from other state variables, has been turned into a prognostic variable and different parameterisations of its rate of change have been tested, by comparing results of simulations to field measurements.Characteristics and limits of such parameterisations have been discussed.Results indicate that all metamorphism formulations perform well in terms of simulated SSA, with median values of the RMSD between observed and simulated SSA lower than 10 m 2 kg −1 .
Compared to the previous description of snow metamorphism based on semi-empirical variables, this new representation in terms of optical diameter does not reduce significantly the discrepancies between simulated and measured SSA profiles.Its interest rather consists in using a physical and easily measurable prognostic variable to characterise the snow microstructure.This approach opens the way to several future improvements.Firstly, it will make it easier to revise parametric laws (such as those of the layer viscosity, the mobility index for wind transport and the albedo) which are directly related to the snow microstructure and potentially depend on the optical diameter.Moreover, this new metamorphism scheme will simplify the data assimilation of various electromagnetic observations, which may lead to a reduced dependency of Crocus simulations on the quality of the meteorological forcing.Lastly, a more complete representation of the properties of the surface layer, including for instance a

Fig. 1 :Fig. 1 .
Fig. 1: Relationship between SSA and grain variables in the original Crocus version.Empirical equations (Vionnet et al., 2012) were used to estimate the optical diameter (a) from dendricity and sphericity in the dendritic case and (b) from sphericity and size in the non-dendritic case.Then, the optical diameter was converted into SSA using Eq. 1.

Fig. 2 .
Fig. 2: Time evolution of layer age (a) and SSA (b) in a multi-layer model, simulated using the parametrization T07(Taillandier et al., 2007).Layer 1 and layer 2 are created at different times to take into account new snowfalls.At t = 30 days, these two layers aggregate and give rise to layer 3.The age of the new-formed layer has to be chosen arbitrarily and this has an impact on the SSA rate of change.
computed snow SSA from micro-tomography during an isothermal experiment at T = −2 • C; Schleef and Loewe (2013) analysed micro-tomographic images acquired during an isothermal experiment at T = −20 • C; Taillandier et al. (2007) determined SSA by measuring the adsorption of methane at 77 K on snow that had evolved under temperature gradient conditions (T mean = −10 • C, ∇T = 33 K m −1 ); Calonne et al. (2014) used micro-tomographic images to compute SSA under temperature gradient conditions (T mean = −4 • C, ∇T = 43 K m −1

Fig. 3 :Fig. 3 .
Fig. 3: Time evolution of snow SSA, according to different metamorphism formulations.B92 stands for Brun et al. (1992), T07 for Taillandier et al. (2007) and F06 for Flanner and Zender (2006), whereas C13 refers to this work.Temperature was set to -20 • C and density to 200 kg m −3 .The four panels correspond to different temperature gradient values.

Fig. 4 :
Fig. 4: Comparison between the SSA rate of decay predicted by different metamorphism formulations and measured during cold room experiments.The isothermal experiments of Flin et al. (2004) and Schleef and Loewe (2013) were performed at T = -2 • C and T = -20 • C, respectively.In the temperature gradient experiment of Taillandier et al. (2007), ∇T was 33 K m −1 and the temperature was -10 • C.These values were respectively 43 K m −1 and -4 • C for the temperature gradient experiment ofCalonne et al. (in prep.).In (b), we plotted both the data (black dots) and the fit (black dotted line) according to Eq. 7 ofSchleef and Loewe (2013). 32

Fig. 4 .
Fig. 4. Comparison between the SSA rate of decay predicted by different metamorphism formulations and measured during cold room experiments.The isothermal experiments of Flin et al. (2004) and Schleef and Loewe (2013) were performed at T = −2 • C and T = −20 • C, respectively.In the temperature gradient experiment of Taillandier et al. (2007), ∇T was 33 K m −1 and the temperature was −10 • C.These values were respectively 43 K m −1 and −4 • C for the temperature gradient experiment of Calonne et al. (2014).In (b), we plotted both the data (black dots) and the fit (black dotted line) according to Eq. (7) of Schleef and Loewe (2013).

Fig. 5 .
Fig. 5. Time evolution of SWE and snow height at Summit.The metamorphism formulations described in Sect.2.3 are represented by different colours: B92 (Brun et al., 1992) in red, C13 (this work) in blue and F06 (Flanner and Zender, 2006) in green.Observations from AWS are in solid black.(a) Simulated SWE since 1 February 1979.(b) Simulated snow height since 1 February 1979.(c) Close up on simulated and measured snow height during our field campaign (all curves are reset to zero on 1 January 2011, in order to better compare them).

Fig. 7 .
Fig. 7. Comparison in terms of SSA between the metamorphism formulations implemented in Crocus.In order to compare results obtained using different formulations, all SSA profiles were interpolated over a 1 mm vertical grid.The plots show the difference between the new formulations C13 (top) and F06 (bottom) and the original formulation B92, for simulations run at Summit from 1 April to 31 July 2011.

Fig. 9 :
Fig. 9: Time evolution of the SSA of the top 1 cm at Summit.The simulated SSA were obtained by making a weighted average of the SSA values within the top 1 cm of snow.Grey vertical bands represent periods of strong wind (wind speed > 9 m s −1 ).Formulations C13 and F06 were run setting the maximum SSA value to 65 m 2 kg −1 in (a) and to 80 m 2 kg −1 in (b). 37

Fig. 9 .
Fig. 9. Time evolution of the SSA of the top 1 cm at Summit.The simulated SSA were obtained by making a weighted average of the SSA values within the top 1 cm of snow.Grey vertical bands represent periods of strong wind (wind speed > 9 m s −1 ).Formulations C13 and F06 were run setting the maximum SSA value to 65 m 2 kg −1 in (a) and to 80 m 2 kg −1 in (b).

Fig. 11 :
Fig. 11: Simulated density (left column) and SSA (right column) at Col de Porte during winter 2011/2012.Simulations were run using different metamorphism formulations.From top to bottom:

C
. M. Carmagnola et al.: Prognostic representations of the optical diameter of snow Comparison in terms of SSA between the metamorphism formulations implemented into In order to compare results obtained using different formulations, all SSA profiles were ted over a 1 cm vertical grid.The three plots show the difference between the new formu-13 (top), T07 (middle) and F06 (bottom) and the original formulation B92, for simulations ol de Porte during winter 2011/2012.40

Fig. 12 .
Fig. 12.Comparison in terms of SSA between the metamorphism formulations implemented in Crocus.In order to compare results obtained using different formulations, all SSA profiles were interpolated over a 1 cm vertical grid.The three plots show the difference between the new formulations C13 (top), T07 (middle) and F06 (bottom) and the original formulation B92, for simulations run at Col de Porte during the winter of 2011/2012.

Fig. 14 .
Fig. 14.Comparison between simulations and observations in terms of SSA (left column) and d opt (right column).(a, b) values computed for Summit.(c, d) values computed for Col de Porte during the winter of 2009/2010.(e, f) values computed for Col de Porte during the winter of 2009/2010, in which simulations were stretched vertically in order to match the measured snow height.(g, h) Same as (c, d) for the winter of 2011/2012.(i, j) Same as (e, f) for the winter of 2011/2012.

Table 2 .
use optical diameter (d opt ) and s.See text for more details.Metamorphism laws for formulation C13 in the case of dry snow.G is the vertical temperature gradient (∇T = |∂T /∂z|, in K m −1 ), T the temperature (in K) and t the time expressed in days.f , g, h and are empirical functions to predict depth hoar growth rate in the case of high G values

Table 3 .
Overview of the measurements acquired at Summit Camp in 2011 and at Col de Porte during the winters of 2009/2010 and 2011/2012.These data were used to evaluate the different representations of the snow optical diameter growth rate implemented in SURFEX/ISBA-Crocus.