Snow specific surface area simulation using the one-layer snow model in the Canadian LAnd Surface Scheme ( CLASS )

Snow grain size is a key parameter for modeling microwave snow emission properties and the surface energy balance because of its influence on the snow albedo, thermal conductivity and diffusivity. A model of the specific surface area (SSA) of snow was implemented in the one-layer snow model in the Canadian LAnd Surface Scheme (CLASS) version 3.4. This offline multilayer model (CLASS-SSA) simulates the decrease of SSA based on snow age, snow temperature and the temperature gradient under dry snow conditions, while it considers the liquid water content of the snowpack for wet snow metamorphism. We compare the model with ground-based measurements from several sites (alpine, arctic and subarctic) with different types of snow. The model provides simulated SSA in good agreement with measurements with an overall point-to-point comparison RMSE of 8.0 m2kg−1, and a root mean square error (RMSE) of 5.1 m2kg−1 for the snowpack average SSA. The model, however, is limited under wet conditions due to the single-layer nature of the CLASS model, leading to a single liquid water content value for the whole snowpack. The SSA simulations are of great interest for satellite passive microwave brightness temperature assimilations, snow mass balance retrievals and surface energy balance calculations with associated climate feedbacks.


Introduction
Snow grain size is of particular interest for microwave snow emission models, the surface energy balance (albedo and turbulent fluxes) and atmospheric-snow chemistry interactions (Domine et al., 2008).At high microwave frequencies (generally measured at 19 GHz and 37 GHz), snow grain size is an important variable affecting snowpack extinction and scattering properties (Kontu and Pulliainen, 2010;Grody, 2008;Durand et al., 2008;Roy et al., 2004).Thus, snow grain size must be considered in microwave snow emission models (MSEM) for the retrieval of snow properties from satellite passive microwave observations (Langlois et al., 2012;Huang et al., 2012;Pardé et al., 2007).Hence, in passive microwave applications, prior information such as snow grain size from a snowpack physical model is required for snow water equivalent (SWE) estimates (Durand and Liu, 2012).The surface albedo is also sensitive to the snow grain size and its vertical profile (Wiscombe and Warren 1980;Jin et al., 2008;Lyapustin et al., 2009;Aoki et al., 2011).Gardner and Sharp (2010) found that the broadband albedos of snowpacks show a logarithmic relationship with specific surface area (SSA).The thermal properties of snow, such as snow conductivity and diffusivity, are also related to snow microstructure (Domine et al., 2008;Adams and Sato, 1993).Surface albedo and snow conductivity are thus key parameters for modeling the surface energy balance in order to understand the impact of snow cover on global and regional climate dynamics (Armstrong and Brun, 2008).They also have a major impact on the prediction of the snow water equivalent as well as the timing of melt onset (Franz et al., 2010).
However, many snow evolution models do not take into account snow grain size.The Canadian LAnd Surface Scheme (CLASS: Verseghy, 1991;Verseghy et al., 1993) is used in the Canadian global circulation models (Scinocca et al., 2008) and the Canadian Regional Climate Model (CRCM: A. Roy et al.: Snow specific surface area simulation Music and Caya, 2007;Caya and Laprise, 1999); it includes a one-layer snow model that does not simulate snow grain metamorphism.This is a major limitation for the assimilation of passive microwave brightness temperature (T B ) data for the improvement of snow simulations.In the context of data assimilation, where physical and emission models of snow are coupled, estimates of snow grain size are needed (Durand et al., 2009;Toure et al., 2011;Langlois et al., 2012).The implementation of snow grain metamorphism within CLASS is thus of particular interest for assimilation purposes.This implementation is not, however, straightforward in a one-layer snow model because snow metamorphism depends on many variables, such as snow age and the temperature gradient, which lead to a stratification of snow layers with different grain sizes.Thus, a major difficulty is that the vertical stratification is not considered in single-layer physical snow models.This study aims to address this issue, as simply as possible, using the CLASS one-layer snow evolution model.Even if one layer snow models are less physically correct then multi-layered models (Brun et al., 1992;Bartlelt et al., 2002;Bougamont et al., 2005;Ettema et al., 2010;Niwano et al., 2012), the SnowMIP experiments have shown that CLASS performs relatively well (Brown et al., 2006;Rutter et al., 2009).Furthermore, in climate and meteorological models, the errors in snow simulations are often related to the precipitation inputs.Hence, in these contexts, a more complex multi-layer model would not necessarily produce better results.
Grain size is a parameter that is difficult to characterize accurately and measure in the field.The specific surface area (SSA), which represents the ratio of the surface area per unit of mass, is a well-defined parameter representing the geometric characteristics of a porous medium, such as snow (Dominé et al., 2001).Methods based on snow reflectance in the shortwave infrared (SWIR) can now provide rapid and reproducible field measurements of SSA (Gallet et al., 2009;Arnaud et al., 2011;Montpetit et al., 2012), which can be related theoretically to grain size.SSA can be related to the radius of a monodisperse collection of ice spheres, each having the same surface area to volume ratio, called the optical radius (R opt ): (1) Recent studies have shown that SSA offers a reliable representation of snow grain size in the context of microwave emission snow modeling (MESM) (Roy et al., 2013;Montpetit et al., 2013;Brucker et al., 2011).These studies showed that a scaling factor on R opt derived from SSA is required to simulate brightness temperatures in order oversimplification of snow grain representation in models.From a good representation of snow grain size in the snowpack for microwave emission simulations, it is possible to determinate which part of the signal is attributable to snow grain and which is attributable to other snow characteristics of interest like SWE.
Considering the importance of snow grain size and the advances made in snow microstructure characterization with the SSA metric, many studies have developed approaches for modeling the evolution of SSA throughout the winter season.Cabanes et al. (2003) first proposed an empirical exponential decay function of time and temperature for snow SSA.Legagneux et al. (2003) showed, using laboratory experiments under isothermal conditions, that the decreasing trend of SSA is best fitted using a logarithmic function.That trend has also been confirmed with X-ray microtomography measurements (Flin et al., 2004;Kaempfer and Schneebeli, 2007;Chen and Baker, 2010).Taillandier et al. (2007), using methane adsorption SSA measurements (Domine et al., 2001) in a taiga environment, proposed empirical relationships for the decrease of SSA as a function of time based on the snow age, snow temperature and the temperature gradient within the snowpack.A similar approach relating SSA to snow type (fresh snow, recognizable particles, aged and rounded crystals, aged and faceted crystals, and depth hoar) and snow density was developed by Domine et al. (2007).Jacobi et al. (2010) implemented these last two approaches in the Crocus multi-layer snow model (Brun et al., 1992).With the model based on snow type and density (Domine et al., 2007), SSA was overestimated in surface snow, but this was mainly because Crocus underestimated density, as this model does not take into account the upward water vapor flux induced by the large temperature gradient in the subarctic snowpack (Taillandier et al., 2006); however, a generally good agreement between SSA simulations and measurements (methane adsorption) was observed when the SSA was calculated based on prognostic equations using snow age (Taillandier et al., 2007).Flanner and Zender (2006) developed a physically-based model to predict the evolution of dry snow SSA.The model considers the snow temperature, temperature gradient and snow density and uses two adjustable parameters for the distribution of crystal sizes and for the irregularity in particle spacing.A weakness of most of these previous approaches is that the wet snow metamorphism is not taken into account, whereas water within the snowpack leads to a drastic decrease of SSA (wet snow metamorphism) due to rapid rounding and an increase in the size of snow grains (Brun, 1989).However, Flanner et al. (2007) implemented wet snow metamorphism following Brun (1989) in the model of Flanner and Zender (2006).Note that wind can also have complex effect on snow grains by enhancing the rate of SSA decrease (Cabanes et al., 2003) or, on the contrary, leading to an increase in SSA (Domine et al., 2009).Morin et al. (2013) compared SSA deduced from the R opt values simulated by Crocus with SSA measured from SWIR reflectance (Gallet et al., 2009) in an alpine environment.They showed qualitative agreement between measured and simulated SSA and that its simulation is difficult under wet snow conditions mostly because of the difficulty in simulating adequately the vertical profile of liquid water content in the snowpack.
The main objective of this study is to evaluate an offline SSA model implemented in the one-layer CLASS snow model for different northern climate environments.More specifically, the SSA model is a multi-layer snow model driven by CLASS outputs to simulate the evolution of SSA in the different snow layers.The evolution of SSA is computed for dry snow using the model of Taillandier et al. (2007) based on snow aging, and the equation of Brun (1989) for wet snow metamorphism.The simulated SSA values are compared with measured SSA derived from SWIR reflectance (Montpetit et al., 2012;Gallet et al., 2009) for five different sites (two northern mid-latitude, arctic tundra, taiga and alpine) throughout the winter season.The model is developed in a perspective of passive microwave applications for SWE retrievals at a large scale, but could be used for other applications such as snow albedo estimates.The study also provides an additional validation of the Taillandier et al. (2007) equations using new sets of accurate in situ SSA measurements for different environments.

CLASS-SSA model
The CLASS-SSA model operates in an "offline" mode, meaning that it uses the CLASS simulated state variables to simulate the SSA evolution, but without feedback on the snowpack evolution.The CLASS snow model is a one-layer model (a detailed description of the snow model in CLASS is given in Bartlett et al., 2006;Brown et al., 2006).Version 3.4 of the standalone driver for CLASS (Verseghy, 2009), which allows running the model using meteorological data, was used in this study.CLASS has been designed to run at a time step of 30 min or less, to ensure numerical stability of the modeled prognostic variables (Verseghy, 2009).In this study, CLASS is run at a time step of 30 min.In our case, the meteorological data used to drive the CLASS model (precipitation rate, air temperature, wind speed, air humidity, and incoming shortwave and longwave radiation) were derived from in situ measurements or from the North American Regional Reanalysis (NARR) data (Mesinger et al., 2006) (more details on driving data are provided in Sect.2.2).The use of NARR data is motivated by the necessity to run the model at a large scale in the perspective of passive microwave space-borne applications.The thermal conductivity of snow was calculated from snow density using the empirical relationship described in Sturm et al. (1997).
The offline SSA model is a multilayer model, where layer evolution is constrained by snow density, snow depth and SWE from CLASS simulations.Figure 1 shows the flowchart of CLASS-SSA.The SSA evolution of dry snow is based on the logarithmic relationship for snow aging developed by Taillandier et al. (2007).The CLASS-SSA model adds snow layers when snowfall occurs.Consecutive precipitation (precipitation occurring during two or more consecutive time steps) is however considered as one precipitation event and contributes to the same layer.The initial SSA (SSA initial ) was set to 73.0 m 2 kg −1 (the mean SSA value for fresh snow measured by Domine et al., 2007; the SSA initial value is discussed in Sect.4), and the density of each new snow layer was set to the fresh snow density calculated by CLASS (Hedstrom and Pomeroy, 1998).However, because we want the CLASS-SSA model to be coherent with the CLASS snow model, the snow parameters (SWE, snow depth and snow density) of the CLASS-SSA model are corrected during the snowpack evolution in order to match the corresponding values simulated by CLASS.Thus, prior to each time step, a correction factor is applied to the SWE value of every snow layer to fit the multilayer model SWE with the CLASS simulation.A densification routine is then implemented, mostly to estimate the position and thickness of each layer within the snowpack.The same densification model as the one used in CLASS (Bartlett et al., 2006) is applied to every layer.After compaction is applied to each layer, if the summed multilayer snow depth is lower than the snow depth simulated by CLASS, a correction factor is applied to the thickness of the www.the-cryosphere.net/7/961/2013/The Cryosphere, 7, 961-975, 2013 top snow layer so that the summed multilayer snow depth corresponds to the snow depth simulated by CLASS (Fig. 1).However, in this context, the top layer cannot be less than 100 kg m −3 .If it reaches 100 kg m −3 , densification is applied to the second layer and so on.In this case, the correction is applied only to the top layers to avoid an unrealistic thickening of the dense bottom layers.On the other hand, if the sum of the snow depths for all the layers is higher than the snow depth simulated by CLASS, a correction is applied to all the layers, but the density of any layer cannot exceed the maximum snow density estimated by CLASS.Thus, when the snowpack melts, the density of every layer increases due to the decreasing thickness, leading to wet densification.
The SSA evolution for each snow layer is then calculated considering the model of Taillandier et al. (2007) (Fig. 1).The model parameterizations for SSA evolution are based on snow age and snow temperature (T snow ).Two algorithms are available, depending on the temperature gradient regime: one for equi-temperature (ET) metamorphism, (2) and the other for strong temperature gradient (TG) metamorphism, where t is the age of the snow layer in hours.Note that in Eqs. ( 2) and (3), SSA is in cm 2 g −1 .T snow is the snow layer temperature ( • C) calculated by linearly interpolating the air temperature and the CLASS simulated snow-soil interface temperature.Air temperature appears more accurate and representative than the surface (skin) temperature for estimating the snowpack temperature gradient.Figure 2 shows a rapid decrease in the SSA over the first few days, which is related to destructive metamorphism when snow crystals lose most of their complicated shape and break up into smaller grains with less total surface area (Sommerfeld and Lachapelle, 1970).This metamorphic process is faster in warmer snow (higher T snow ) (Colbeck, 1983).After a few days, the decrease in SSA slows down earlier in the ET regime when compared with the TG regime.The process of constructive metamorphism is dominant when the temperature gradient induces water vapor transport from warm to cold temperatures, causing rapid grain growth from vapor deposition at the bottom of the snow grains (Colbeck, 1983).Hence, in the absence of that mechanism in ET conditions, the decrease in SSA slows and SSA rapidly reaches its minimum value.
According to Eqs. ( 2) and ( 3), the rate of SSA decrease for a given time step ( SSA) depends on snow age, snow temperature, temperature gradient and SSA initial .Based on  2) and ( 3), as well as a function of liquid water content (LWC) from Eq. ( 5).Jacobi et al. (2010), we calculate the SSA from Eqs. ( 2) and (3) according to the following: where t corresponds to the time step (0.5 h).The SSA is then subtracted from the model's previous SSA value.Jacobi et al. (2010) used a temperature gradient threshold (TG threshold ) of 15 K m −1 to distinguish between ET and TG conditions, which will be evaluated in this study.Taillandier et al. (2007) also suggested a minimum value for SSA, because the logarithmic equation for SSA can lead to unrealistic low values.The minimal SSA value is set to 8.0 m 2 kg −1 , as proposed by Taillandier et al. (2007).
Nevertheless, the parameterization reported by Taillandier et al. (2007) does not take into account metamorphism during wet snow conditions.The equation of Brun (1989), derived from experimental measurements, provides a way to simulate the evolution of snow grain volume under wet snow conditions with respect to the liquid water content of the snowpack.The equation of Brun (1989) can be expressed with optical radius growth ( R opt ): where C 1 and C 2 are empirical coefficients (C 1 = 1.1 × 10 −3 mm day −1 , C 2 = 3.7 × 10 −5 mm day −1 ) and LWC is the liquid water content in mass percentage.Note that in the experiment of Brun (1989), the empirical relationship was based on the volume equivalent sphere deduced from the measured mean convex snow grain radius, which is a definition closely related to the SSA. Figure 2 shows that the SSA decrease is more pronounced when LWC increases in the snowpack.In this study, when the CLASS liquid water content is greater than zero, the model-derived SSA value is converted into its equivalent R opt using Eq. ( 1) in order to apply Eq. ( 5), and then reconverted to SSA.Furthermore, because the LWC distribution is not homogeneous within the snowpack, even though CLASS uses a single LWC value for the whole snowpack, here the LWC is first distributed in the top 10 cm.If the LWC in the top 10 cm is greater than 10 % in mass, the excess water is distributed to the rest of the snowpack.The 10 % limit can thus be considered as the water retention capacity.This value was chosen because 10 % in mass is the value where, in the experiment of Brun (1989), LWC reaches the irreducible water content and percolation occurs, leading to a saturation of grain growth rate increase for high LWC.However, the liquid water retention capacity of CLASS for the whole snowpack was kept at 4 %.

Sites and data
Snowpit measurements were conducted at five sites.Measurements were taken during the winter of 2010-2011 at the first two sites, which were located in an open mid-latitude northern environment.The sites were at the Site interdisciplinaire de recherche en environnement extérieur (SIRENE) experimental station at the Université de Sherbrooke (45.37At the first four sites, SSA profiles were taken at a vertical resolution of 5 cm.The SSA was measured with the shortwave InfraRed Integrating Sphere (IRIS) system (Montpetit et al., 2012), based on the principle described by Gallet et al. (2009), which exploits the relationship between the SWIR snow reflectance and the SSA (Kokhanovsky and Zege, 2004).The density was measured with a 185 cm 3 density cutter, and the samples were weighed with a 100 g Pesola light series scale with an accuracy of 1 g.The temperature was measured with a Traceable 2000 digital temperature probe.Liquid water content of snow was also measured with the Snow Fork (Toikka Engineering Ltd., Espoo, Finland) at the Churchill arctic fen site during wet conditions on 13 and 16 April 2010.At the Col de Porte site (the fifth site), 16 SSA profiles were taken using the Dual Frequency Integrating Sphere for Snow SSA instrument (DUFISSS: Gallet et al., 2009), also based on the relationship between the SWIR reflectance and the SSA.Note that from late February onwards, warm conditions led to several snowmelt events, which caused a significant decrease in the snow SSA values.Hence, a distinction was made between dry snow conditions at Col de Porte (7 sets of data from 6 January to 16 February) and wet snow conditions (9 sets of data from 25 February to 20 April).The total snow depth and snow density profiles were also measured (Morin et al., 2013) and ultrasonic snow depth observations were acquired at SIRENE and Col de Porte.
NARR data (Mesinger et al., 2006) (2 m air temperature and air humidity, precipitation, 10 m wind speed, surface shortwave and longwave radiation) were used to force the CLASS model at the first four sites.Langlois et al. (2009) show that the NARR product delivers reliable input data for snowpack modeling.Forcing data from the NARR nearest neighbor pixel of each site was employed.As NARR provide data on a three-hour time step, the variables were interpolated to a 30 min time step, except for precipitation, which maintained a three-hour interval.To initialize the starting conditions, the CLASS model was run starting the year prior to the winter in this study:

CLASS snow parameter evaluation
First, an analysis of the one-layer CLASS snow model simulations was conducted.Simulated snow density, total snow depth and SWE were compared with all measurements taken in snow pits where SSA profiles were measured.Figure 3 shows that the simulation accuracy varies from one site to another.Snow density is generally accurate except at Col de Porte, where the density was overestimated.The overestimation is probably due to the high densification of snow under wet conditions with CLASS.For snow depth, there is an underestimation for the Churchill sites.Since there were underestimates at both forest and fen sites, NARR precipitation is probably the main cause.In fact, the cumulated NARR precipitation from the beginning of the snow season to the first snowpit measurement in February at the dry fen is lower (97.1 mm) than the snowpit measured SWE (157.3 mm).However, other phenomena such as blowing snow and interception by vegetation could lead to differences between the simulated and measured SWE (consequently snow depth as well).Comparisons between continuous ultrasonic snow depth observations at SIRENE and Col de Porte also show that errors in diagnosing precipitation phase at the beginning of the snow season lead to an offset of snow depth (overestimation at SIRENE) (Fig. 4a).This sensitivity to precipitation phase in CLASS is also demonstrated in Langlois et al. (2013).Figure 4b shows that underestimation of melt events at the beginning of the season also lead to positive offset in the snow depth.However, the snow depth RMSE is comparable to what was found with CLASS 3.1 in the Snow Model Intercomparison Project (Brown et al., 2006).The overall SWE RMSE is 64.7 mm, which is close to what was found in Langlois et al. (2013) between modeled (CRCM) and observed SWE values for northern Québec.There is a good correlation between the measured and simulated SWE for the SIRENE and Col de Porte sites, where there is, however, a consistent overestimation.However, the SWE is underestimated at St-Romain and both Churchill sites.

CLASS-SSA model evaluation and validation
In this study, SSA is considered for the evaluation and validation because measurements of shortwave infrared reflectance of snow are related to SSA (see Sect. 2.2).To evaluate the CLASS-SSA simulations, an analysis of the TG threshold was first conducted.However, differences between the simulated and measured snow depths (Figs. 3 and 4), caused problems when relating measured SSA to its corresponding simulated SSA for a given snow depth in the snowpack.A correction was applied to the measured snow depth in order to match the simulated snow depth; this caused the measured profiles to be stretched or compressed.The root mean square error values between the simulated and measured SSA (RMSE SSA ) were calculated at the first four sites for different TG threshold values.The Col de Porte site where wet snow metamorphism dominated was excluded because wet metamorphism has a strong influence on SSA evolution that is not related to the TG threshold .Figure 5 shows that, for TG threshold values between 10 K m −1 and 20 K m −1 , the RMSE SSA for the 5 sites is relatively constant.The minimum RMSE SSA value (7.3 m 2 kg −1 ) for the 5 sites is TG threshold = 16 K m −1 , which is close to the 15 K m −1 value used by Jacobi et al. (2010).This value is also consistent with Taillandier et al. (2007), who proposed that the TG threshold should be between 9 K m −1 and 20 K m −1 .The SIRENE site reached a minimum RMSE at TG threshold = 14 K m −1 , whereas the RMSE at both Churchill sites slightly increased with TG threshold .The RMSE at the St-Romain site and the dry snow condition data for the Col de Porte site slightly decreased from 10 K m −1 to 25 K m −1 before reaching a constant value (Fig. 5).
The minimum RMSE was at TG threshold = 16 K m −1 and was thus used to simulate SSA with CLASS-SSA at the 5 sites (Fig. 6).Previous studies have generally defined the TG threshold for depth hoar formation between 10 K m −1 and 20 K m −1 (Taillandier et al., 2007;Colbeck, 1983;Marbouty, 1980).However, it should be noted that the formation of faceted snow crystals has been observed at low growth rates under low gradient thresholds (Domine et al., 2003;Flin and Brzoska, 2008).Pinzer and Schneebeli (2009) proposed that alternating temperature gradients also leads to formation of rounded grains, similar to those observed in ET metamorphism.However, overall comparisons show good agreement between simulated and measured SSA (Fig. 6).More specifically, the SIRENE and St-Romain results show similar patterns with a gradient from low SSA at the bottom to higher SSA coming from fresh precipitation at the top.Nevertheless, there is a low SSA layer that appeared in mid-December 2010 caused by a melt event.This layer was observed as a melt ice-crust layer of 3 cm with low SSA (measured with a SWIR camera: Montpetit et al., 2012) during the snowpit measurements, but SSA was not measured with IRIS because it is difficult to extract this kind of snow (i.e., crusts) with the IRIS instrument.For the Churchill sites, both measured and simulated SSA are low near the bottom (≈ 25 cm), which is related to the formation of depth hoar in the presence of a high temperature gradient.However, the simulated SSA values in the top layers are generally higher than the measurements.This may be due to the underestimation of the snow depth at the beginning of the season causing an underestimation of the relative thickness of the bottom layers with low SSA within the snowpack, which leads to an overestimation of the top layer thickness (Fig. 6c and d).Underestimation of the April 2010 measurements at both sites should, however, be attributed to an underestimation of snow LWC by CLASS during the spring melt, which limited the decrease of simulated SSA by wet metamorphism.In fact, LWC as measured with the Snow Fork instrument on 13 and 16 April at the Churchill arctic fen site suggests a strong underestimation by CLASS (0.2 % vs. 3.8 % on 13 April and 1.0 % vs. 15.5 % on 16 April, for CLASS and the Snow Fork, respectively).The issue with LWC is discussed in Sect. 4. The Col de Porte site illustrates the difference between the first seven dry sets of data showing good agreement, and the second period, starting on 25 February, when wet snow becomes predominant in the snowpack, giving a systematic overestimation of the SSA.
The comparison of the simulated SSA values to their corresponding measurements gives a RMSE of 8.0 m 2 kg −1 , which represents an error of 42.3 % (Fig. 7).Part of the error could be attributed to the fact that we did not necessarily compare the same snow layers due to different positions between the simulated and measured points.The correction applied to the simulated snow depth profile might be a factor, but the high variability within a SSA profile might also be a source of error.The simulated SSA variations are also strong within the snowpack, mainly for high SSA, where the evolution is faster (see Fig. 2).Considering the mean depthaveraged SSA weighted by the snow layer thickness, the RMSE decreases significantly to 5.1 m 2 kg −1 , representing an error of 25.7 %.Furthermore, the coefficient of determination (R 2 ) increases from 0.60 to 0.84.As mentioned previously, another major source of error corresponds to the influence of wet conditions, as observed at the Col de Porte site www.the-cryosphere.net/7/961/2013/The Cryosphere, 7, 961-975, 2013 after the mid-season.In fact, by removing data for this wet period, the depth-averaged RMSE decreases to 4.1 m 2 kg −1 (17.5 %).As mentioned in Morin et al. (2013), even with a multi-layer model, limitations on the precision of LWC simulations exacerbate the difficulty of modeling snow grain evolution under wet conditions.The weakness of the model under wet snow conditions will be analyzed below in the Discussion (Sect.4).The SSA at St-Romain and at Col de Porte (dry snow period) are underestimated by the model, while at the Churchill sites SSA is slightly overestimated due to the high simulated SSA in the top layers (Fig. 6c and d).

Discussion
The simulation of a stratified phenomenon such as SSA using a one-layer snow model such as CLASS requires certain assumptions and simplifications of the physics within the snowpack.These assumptions may induce errors in estimates of the SSA evolution.Here, we discuss the different elements that may impact the precision of the model and how they may influence the estimates.It is thus possible to identify the conditions under which CLASS-SSA is more limited and propose possible improvements.We then consider the proposed model's application, mainly in the context of passive microwave simulations.

Sources of errors
An assumption made in the CLASS-SSA approach is that the temperature profile of the snowpack is linear.
In general, the temperature variations will be larger in the top layers that are responding to the variations in air temperature, while the bottom layers are less affected as the air temperature fluctuations do not reach these layers because of the low snow thermal conductivity ( Armstrong and Brun, 2008;Vionnet et al., 2012).Hence, in cold weather like in Churchill, the linearity of the temperature profiles is likely to induce underestimated snow layer temperatures.This phenomenon could also partly explain why the SSA of top layers at both Churchill sites is overestimated, considering that the SSA decrease is more pronounced with higher snow temperatures.Furthermore, the linearity of the temperature gradient would generally underestimate the local temperature gradient in the top layers and overestimate the local temperature gradient in the bottom layers.However, in the north during winter, this diurnal temperature cycle is generally in most cases less pronounced than over temperate or mountainous regions (Leathers et al., 1998).Thus, using a linear gradient throughout a dry and relatively shallow snowpack (below about 1 to 1.5 m depth) appears as a satisfactorily hypothesis in most cases over northern areas.Kondo and Yamazaki (1990) demonstrated that a linear temperature profile can be successfully employed in a snowmelt model.As shown through the wet metamorphism simulation, CLASS-SSA is limited by the modeling of snow parameters in CLASS.Hence, the use of a one-layer model giving a LWC for the entire snowpack becomes a limitation.Furthermore, there might be an underestimation of LWC by CLASS.Measurements of LWC with the Snow Fork instrument at the Churchill arctic fen site suggests a strong underestimation by CLASS.Moreover, raising the limit of the simulated snowpack water retention capacity from 4 % to 10 % did not improve the simulated LWC and the SSA calculation under wet conditions in Col de Porte because the CLASS LWC rarely reaches the retention capacity.Part of the problem was, however, resolved by distributing the LWC mostly in the top layers, but the SSA evolution under wet conditions remains a weakness.Table 1 shows that the bias is significantly reduced when wet metamorphism is modeled with snowpack liquid water distributed in the top 10 cm at the Col de Porte wet sites compared to a wet metamorphism considering a uniform LWC or compared to no wet metamorphism.We also tested simulations by drastically increasing the total LWC artificially in CLASS (if LWC > 0 then LWC = 10 %); such conditions significantly reduced the simulated SSA, as expected, with a bias of −2.6 m 2 kg −1 for Col de Porte wet snow condition data (Table 1).This last case confirms that the problem comes from an underestimation of LWC in CLASS under warm conditions.Snow depth errors from CLASS might also impact CLASS-SSA simulations.In fact, as shown for both Churchill sites (Fig. 6c and d), a bias in snow precipitation can impact the representation of the thickness of a given snow layer.Thus, in this study, part of the SSA error could be related to uncertainties in the NARR precipitation data (Langlois et al., 2009).
Other phenomena not parameterized in the CLASS snow model, such as blowing snow, could influence the simulated snow depth (Liston and Hiemstra, 2011).In open areas (four of our five sites), strong wind shear stress could have exceeded the snow particle resistance to dislocation (Li and Pomeroy, 1997) arctic region like Churchill (Baggaley and Hanesiak, 2005).Furthermore, the snow thermal conductivity strongly varies between the tundra, where the snowpack has a high conductivity due to hard wind slabs, compared to taiga and forest snowpacks, which have three to four times lower thermal conductivity due to lower wind compaction and depth hoar development (see Gouttevin et al., 2012).These differences impact the snow temperature and temperature gradient but are not represented in CLASS-SSA.
In the CLASS-SSA model, the SSA initial value is fixed at 73.0 m 2 kg −1 .This value was chosen based on freshly fallen snow SSA measurements (sampled maximum 24 h after snowfall) from methane adsorption by Domine et al. (2007).However, this study shows a range of 33.1 to 155.8 m 2 kg −1 with a standard deviation of ±26.2 m 2 kg −1 based on 63 samples.Freshly fallen snow SSA is rarely modeled as it depends on the type of solid precipitation, which depends on the meteorological conditions (air temperature, wind, type of clouds, atmospheric stratification) when the snowflake is formed.Domine et al. (2007), however, proposed freshly fallen snow SSA values based on four types of fresh snow that can be related to density.As CLASS calculates the fresh snow density from the air temperature using the equation from Hedstrom and Pomeroy (1998), we implemented SSA initial values based on the Domine et al. (2007) relationship.However, this implementation did not change the results significantly: a slight increase in RMSE from 8.0 to 8.3 m 2 kg −1 was found.Figure 8 shows that SSA initial has a relatively low impact on simulations.The sensitivity to SSA initial values appears to be more important for the snowpack where measurements were taken mostly at the beginning of the season (SIRENE and St-Romain).A precise dynamic parameterization of freshly fallen SSA could probably improve the results, mostly for snow with high SSA at the beginning of the snow season.

Comparison with other models
Despite the above simplifications, the CLASS-SSA model simulates SSA with a reasonable accuracy for a wide range of snow types.Our RMSE of 8.1 m 2 kg −1 (Fig. 7) is comparable to the result obtained at Col de Porte by Morin et al. (2013) from internal computation of the optical radius in Crocus (6.37 m 2 kg −1 ) and the method of Domine et al. (2007) based on the density and snow type (8.08 m 2 kg −1 ).Snow data from the 2010 winter season at Col de Porte provide a unique and very accurate time series of SSA measurements (Morin et al., 2013).Figure 9 shows a comparison of temporal snowpack averaged SSA values at Col de Porte for CLASS-SSA, the Crocus model (Morin et al., 2013), and the measurements.When the snowpack is dry, both models underestimate the SSA.On 25 February and after, when wet conditions occurred, CLASS-SSA overestimates the SSA due an underestimation of the snowpack LWC, while Crocus still underestimates the SSA.For this dataset, CLASS-SSA simulations seem comparable to or better than Crocus in dry conditions.However, in wet conditions, Crocus better simulates the decrease in SSA as LWC increases (Fig. 9).Hence, Crocus seems to better capture the dynamics of the SSA evolution.Note that Jacobi et al. (2010) obtained, from 162 snow SSA measurements at a taiga site, an RMSE of 8.6 m 2 kg −1 with the implementation of the Taillandier et al. (2007) approach within Crocus, whereas the implementation of Domine et al. (2007) resulted in a RMSE of 16.2 m 2 kg −1 (the results were highly affected by the underestimation of snow density by Crocus).

R opt analysis for MESM
In the context of using the simulated SSA to assimilate microwave brightness temperatures (T B ), we now examine the impact of errors generated by the proposed approach in terms of T B error.As mentioned in Sect. 1, the DMRT MESM calculates T B from R opt derived from SSA. Figure 10 shows the comparison between R opt derived from measured SSA and R opt derived from simulated SSA.As the relationship between SSA and R opt is not linear (Eq.1), we see that for individual points (Fig. 10, left panel), the differences between simulated and measured R opt are more important for larger grains.This is caused by the fact that for low SSA, a given variation of SSA leads to a larger change in the simulated R opt (a change of SSA from 10 to 8 m 2 kg −1 leads to a change in R opt from 0.327 to 0.409 mm).R opt is then much more sensitive to error in SSA for larger grains.This also partly explains why the error for wet snow condition data in Col de Porte is large.For mean R opt values over the snowpack, there is a RMSE of 0.043 mm and R 2 of 0.84, comparable to SSA results, if the wet snowpacks are excluded.There is, however, a small positive bias of 0.034 mm, mainly caused by the strong influence of large grains in the averaging.
Moreover, the density stratification is another parameter that should be considered when modeling the radiative transfer within the snowpack (e.g., Picard et al., 2012).In our case, the densification routine in CLASS-SSA is only used to calculate the depth and the position of every snow layer, and not necessarily to calculate a precise density.The simulated densities might differ from real densities.An example is the decrease of densification at the bottom snow layer when depth hoar formation occurs that is not taken into account in the SSA offline model.The correction of densities in the top layers (see Sect. 2.1) might also lead to low densities in the top layers.Here, to attenuate the effects of these simplifications in the T B simulations, we consider a bulk snowpack characterized by an effective snow grain effective radius R eff calculated from the averaged R opt following the form suggested by Kontu and Pulliainen (2010): where a = 1.3 and b = 4.7.These empirical parameters were fitted from simulations using averaged R opt derived from measured SSA compared to ground-based radiometric measurements, except for the bulk snowpack (see Roy et al., 2013).Simulations with the DMRT-ML (Picard et al., 2013) were conducted to analyze the effect of the R opt error on T B simulations.The T B simulations were conducted employing single-layer averaged R eff values (Eq.6) ± the derived RMSE of R opt (±0.043 mm).Considering a snowpack of 0.5 m with a bulk density of 250 kg m −3 , an error of 0.043 mm in R opt leads to maximum variations of T B of the order of ±23.3 K at 36.5 GHz and ±2.7 K at 18.7 GHz.We thus see the high sensitivity of 36.5 GHz to grain size, given an error that could be significant in some cases (high depth hoar layer), while the proposed simple approach can be applied for T B simulations at 18.7 GHz with acceptable accuracy.Such a sensitivity analysis would benefit from further development as many combinations of snowpack parameters and conditions could occur.

Model applications
The simulation of snowpack parameters, such as SWE, at individual sites using an operational land-surface scheme designed for use in large-scale climate models could include large errors, as illustrated in Fig. 3.These errors could result from uncertainties in the meteorological forcing data, model parameters, as well as the nonlinearity and scaling effects of the processes modeled (e.g., Andreadis et al., 2008).The proposed model opens opportunities to couple CLASS with MESM for improving SWE estimates.Data assimilation offers the potential to merge information on snow variables from satellite observations and land-surface model simulations.CLASS-SSA was developed mainly for passive microwave T B assimilation in CLASS to improve estimates of snow parameters.The model employed in this study provides a good estimate or "first guess" of the snow grain size and a description of the snow type at a given time during the snow season.Inversion approaches, where parameters (snow depth, snow density) are retrieved by minimizing the differences between simulated and measured brightness temperatures (Langlois et al., 2012;Vachon et al., 2010;Pardé et al, 2007) will benefit from SSA simulations by taking into account the important effect of snow metamorphism on the microwave signal.This "first guess" methodology could also be used as a state initial condition in more complex data assimilation system approaches (Toure et al., 2011;Durand et al., 2009;Reichle, 2008) because the grain-size parameterization is no longer the dominant source of uncertainty.Grain size could still be considered as one of many sources of uncertainty, but with known likely error or variation.Hence, the CLASS-SSA model could be applied to improve SWE estimates at large scales from satellite-borne passive microwave information.From this perspective, as mentioned previously (Sect.4.3), attention needs to be paid to the effect of the conversion of SSA to R opt on the uncertainty related to the grain size simulation that depends on the type of grains.In fact, considering Eq. ( 1), errors in SSA with large grains (low SSA, such as depth hoar) will lead to higher variation of R opt than for smaller grains (high SSA) (see Morin et al., 2013).
The proposed methodology could also be implemented for a hydrology land-surface scheme (HLSS), such as the one developed within the framework of Environment Canada's community environmental modeling system, MESH.MESH evolved from the WATCLASS model that links a hydrological routing model (WATFLOOD) (Pietroniro et al., 2006) to the Canadian LAnd Surface Scheme (CLASS) discussed in this paper.It is used as a basis for coupling horizontal surface hydrology (river routing) with both weather and climate atmospheric models (see discussion by Teutschbein and Seibert, 2010).
Furthermore, snow surface albedo (mostly in the infrared) is driven by snow grain size.Hence, the use of SSA estimated with CLASS-SSA could lead to improved estimates of snowpack albedo, which are derived from a physically-based model.Pure snow albedo (no impurities) could be related to the SSA using a simple optical equation model suggested by Kokhanovsky and Zege (2004).Based on the simple radiative transfer model of Gardner and Sharp (2010), an error of 8.0 m 2 kg −1 in SSA leads to an uncertainty in the broadband albedo calculation of around ±3 % for small grains (30 m 2 kg −1 ) to around 6 % for large grains (10 m 2 kg −1 ).It should be noted that the grains at the top of the snowpack that drive the broadband surface albedo are generally smaller (mostly in dry conditions) and thus less affected by grain uncertainty.

Conclusions
This study analyzes the coupling of a SSA evolution model with a one-layer snow model from the Canadian LAnd Surface Scheme (CLASS).The simulated SSA values were compared with a unique SSA database for five different sites, representing four different climatic environments, including a wide range of snow types.Based on the SSA decrease due to snow aging in snow layers (Taillandier et al., 2007), the CLASS-SSA model is an offline multi-layer parameterization driven by CLASS single-layer snow model outputs.The CLASS-SSA model also considers wet metamorphism, using the equation of Brun (1989) based on the liquid water content of snow.
Despite the limits of a simple one-layer snow model, it provides a SSA estimate with an overall RMSE of 8.0 m 2 kg −1 for individual layers, and a depth-averaged RMSE of 5.1 m 2 kg −1 for the snowpack SSA.The model, however, shows weaknesses in the wet snow metamorphism regime, which is mostly due to a low bias in the snow model simulations of LWC within the snowpack.
The proposed implementation of the SSA model in an offline mode and driven by a one-layer snow model offers a simple, computationally efficient and versatile approach.It would not be difficult to implement for other models as it only needs six basic inputs that are normally available from snow models (snow depth, SWE, snow density, LWC, soilsnow interface temperature, and air temperature).This approach is thus applicable for other one-layer snow models (Turcotte et al., 2007;Bélair et al., 2003), and also for multilayer models where SSA is not explicitly modeled, such as the Snow-Atmosphere-Soil Transfer (SAST) energy balance snow physics model (Sun et al., 1999).
Future work will evaluate the use of these SSA simulations for satellite passive microwave brightness temperature assimilations and surface snow albedo calculations.
CLASS-SSA model flow chart.

Fig. 2 .
Fig. 2. SSA evolution as a function of the temperature gradient regime (ET or TG with LWC = 0) and snow temperature (T snow ) from Eqs. (2) and (3), as well as a function of liquid water content (LWC) from Eq. (5).
from 1 October 2009 to 1 June 2011 at SIRENE and St-Romain; from 1 October 2008 to 1 June 2010 at the two Churchill sites.At the Col de Porte site, meteorological variables (air temperature, humidity, windspeed, precipitation and incoming shortwave and longwave radiation) recorded with an hourly time resolution throughout the snow season of 2009-2010 (from 20 September 2009 to 10 May 2010) were interpolated to a 30 min time-step and used to drive the CLASS model (see Morin et al., 2012 for more details on the Col de Porte meteorological data).

Fig. 3 .
Fig. 3. Comparison of CLASS simulated snow properties with field measurements at different sites.Dry and wet sites at Col de Porte are separately indicated (see legend).

1112 Fig. 4 .
Fig. 4. CLASS snow depth simulations compared with Ultrasonic measurements at (a) SIRENE and (b) Col de Porte.Black dots are the snow depths from snowpit measurements.

Fig. 5 .
Fig. 5. RMSE SSA between the measured SSA and the simulated SSA using CLASS-SSA as a function of the temperature gradient threshold (TG treshold ).

Fig. 6 .
Fig. 6.Seasonal profile of simulated SSA at the five sites, (a) SIRENE, (b) St-Romain, (c) Churchill Arctic fen, (d) Churchill Forest, and (e) Col de Porte), compared to SSA measurements (squares: where the measured SSA profiles were adjusted to the simulated snow depth, see text).The measured profiles are stretched or compressed to fit with simulated snow depth.Black dots correspond to the measured snow depths.

Fig. 7 .
Fig. 7. Measured vs. simulated SSA comparison (left panel shows all points; right panel shows the one-layer average).The RMSE in averaged SSA (right panel) without the wet points is 4.1 m 2 kg −1 .

Fig. 8 .
Fig. 8. Bias between the measured SSA and the simulated SSA using CLASS-SSA as a function SSA initial (the vertical dotted blue line represents the SSA initial set in CLASS-SSA at 73.0 m 2 kg −1 ).

Fig. 9 .
Fig. 9. Snowpack-averaged SSA evolution with time at Col de Porte for CLASS-SSA, Crocus (Morin et al., 2013) and the measurements.The last Crocus value was excluded because the simulations give no snow on ground.Error bars on measurements are the measurements accuracy (12 %: Gallet et al., 2009).

Fig. 10 .
Fig. 10.Measured vs. simulated R opt derived from SSA (Eq. 1) comparison (left panel shows all points; right panel shows the one-layer average).The RMSE in averaged R opt (right panel) without the wet points is 0.043 mm.

Table 1 .
SSA bias and RMSE for Col de Porte wet snow condition data for different configurations of CLASS-SSA wet metamorphism.