Simulation of the specific surface area of snow using a one-dimensional physical snowpack model: implementation and evaluation for subarctic snow in Alaska

The specific surface area (SSA) of the snow con- stitutes a powerful parameter to quantify the exchange of matter and energy between the snow and the atmosphere. However, currently no snow physics model can simulate the SSA. Therefore, two different types of empirical parameter- izations of the specific surface area (SSA) of snow are im- plemented into the existing one-dimensional snow physics model CROCUS. The parameterizations are either based on diagnostic equations relating the SSA to parameters like snow type and density or on prognostic equations that de- scribe the change of SSA depending on snow age, snowpack temperature, and the temperature gradient within the snow- pack. Simulations with the upgraded CROCUS model were performed for a subarctic snowpack, for which an extensive data set including SSA measurements is available at Fair- banks, Alaska for the winter season 2003/2004. While a rea- sonable agreement between simulated and observed SSA val- ues is obtained using both parameterizations, the model tends to overestimate the SSA. This overestimation is more pro- nounced using the diagnostic equations compared to the re- sults of the prognostic equations. Parts of the SSA deviations using both parameterizations can be attributed to differences between simulated and observed snow heights, densities, and temperatures. Therefore, further sensitivity studies regarding the thermal budget of the snowpack were performed. They revealed that reducing the thermal conductivity of the snow or increasing the turbulent fluxes at the snow surfaces leads


Introduction
For large areas of the Earth, the snow cover constitutes either permanently or seasonally the interface between the Earth's surface and the atmosphere. Therefore, exchange processes involving the snow cover are of great importance for many scientific fields like meteorology, climatology, atmospheric chemistry, and biology. These exchanges govern the transfer of momentum, energy, and matter including water vapor and trace compounds between the snow and the atmosphere and vice versa, which depend on the structural and microphysical properties of the snow cover.
One important parameter that can be used to quantify the interaction between snow and the atmosphere is the specific surface area (SSA) (Domine et al., 2008). It represents the surface area of snow crystals that is accessible to gases per unit mass of snow (Hanot and Dominé, 1999). Therefore, it can be measured directly by the gas adsorption technique (Hanot and Dominé, 1999;Legagneux et al., 2002) or stereology (Arnaud et al., 1998). However, both techniques are time-consuming and the number of measurements in the field or in the laboratory obtained with these techniques Published by Copernicus Publications on behalf of the European Geosciences Union. will remain limited. Since the adsorption of gas phase constituents is related to the SSA, it provides a tool to quantify the uptake and release of volatile and semi-volatile impurities by the snow if the temperature-dependent adsorption coefficients of the concerned compounds are known (e.g. Roth et al., 2004;Domine at al., 2007a).
The snow albedo describing the fraction of reflected incoming solar radiation impacts the radiation budget in snowcovered regions. Since the albedo of pure snow and the snow grain size are related , the SSA can be utilized to simplify the calculation of the snow albedo (e.g. Wiscombe and Warren, 1980). This is based on the assumption that the snow can be represented by a collection of spheres (e.g. Flanner and Zender, 2006) with a so-called effective radius r eff related to the SSA according to: where ρ ice is the ice density. Therefore, a snow model that is capable of predicting the SSA of a snowpack is a powerful tool to investigate the snow albedo and its changes over time.
In real snow the further absorption by light-absorbing impurities at wavelengths below 900 nm must also be considered . In principle, the SSA and its evolution can be calculated precisely using a full three-dimensional snow structure model including information about the spatial distribution, the shape, and the size of each snow crystal. Techniques to measure such snow structure properties like x-ray (absorption) (micro-)tomography have recently become available (Flin et al., 2003;Kaempfer and Schneebeli, 2007;Kerbrat et al., 2008). However, snowpack models simulating three-dimensional snow structural properties and its development for an entire snowpack are still unavailable. The most comprehensive one-dimensional physical snowpack models (e.g. Brun et al., 1989Brun et al., , 1992Jordan, 1991;Yamazaki et al., 1993;Bartelt and Lehning, 2002;Lehning et al., 2002a) currently available have mainly been developed for operational avalanche forecasting or to deliver macroscopic information like density or snow type for different horizontal layers of the snowpack, but not for SSA. Fortunately, two different types of empirical relationships have been proposed for the SSA relying on parameters that can be delivered by such snow models. The first relationship calculates the SSA based on a diagnostic equation, in which the SSA is computed using actual snow properties like the snow type and density (Domine et al., 2007b). The second relationship is described by a prognostic equation that uses the snow age, the snow temperature, and the temperature gradient in the snowpack to predict the decrease of the SSA as a function of time (Taillandier et al., 2007).
The parameterizations based on the snow age use snow temperatures and temperature gradients within the snowpack and, thus, rely directly on the thermal budget of the snowpack. The development of specific snow types and the evolu-tion of snow density during the metamorphism of dry snow both depend on the temperature gradient within the snowpack (Colbeck, 1986;Sturm and Benson, 1997). Therefore, it may be expected that for both parameterizations a successful simulation of the SSA depends on a reasonable representation of the thermal budget of the snowpack in the model.
Here, both empirical relationships between SSA and snow density and age are implemented in the one-dimensional snow physics model CROCUS (Brun et al., 1989(Brun et al., , 1992. The model is then applied to the subarctic snowpack at Fairbanks, Alaska. The results for this snowpack are compared to observations of snowpack properties including SSA measurements performed during the 2003/2004 winter. Results of sensitivity studies are presented regarding the parameters influencing the thermal budget of the snowpack and the SSA simulations. Finally, further results of the simulations regarding the snow area index and the depth hoar formation are shown and compared to observational data.

Field measurements
Measurements made during the winter 2003/2004 at the Large Animal Research Station (LARS) of the University of Alaska Fairbanks (64.87 • N, 147.73 • W) were used as input data for the snowpack modeling as well as for validation purposes. Details of the intensive investigation of the snowpack properties were described previously (Taillandier et al., 2006).
Routine meteorological measurements of wind speed, air temperature, relative humidity, and down-welling shortwave radiation were recorded at the field site as 15-min averages. They have been used to calculate 1-h averages as input for the snowpack modeling. The modeling further requires information about the cloudiness, the precipitation rate and phase, and the long-wave radiation, which are, however, not available or only available to a limited extent from the field site. The cloudiness was extracted from hourly meteorological observations performed 7 km away at the Fairbanks International Airport (64.80 • N, 147.88 • W) available at the US National Climatic Data Center (www.ncdc.noaa.gov). The observations include six different categories of cloudiness between clear sky and overcast. Since the snowpack model requires a fractional cloud cover, each observational category was converted into a value between 0 and 1. The precipitation rate and its phase were extracted from daily measurements of the precipitation at Fairbanks International Airport. The daily measurements were used together with visual observations at the field site indicating snowfall events to extract hourly averages of the precipitation. In these cases the daily values have been equally distributed over the period of the snowfall events. The simulations start only with the beginning of the continuous measurements at LARS The Cryosphere, 4, 35-51, 2010 www.the-cryosphere.net/4/35/2010/ commencing at the end of October. However, snowfall was recorded at the airport during the entire month with a total precipitation of 6.6 mm snow water equivalent (SWE) for the period between 1 and 28 October. To correct for this missing precipitation, we added an artificial snow event at the beginning of the simulations (30 October) with the same total precipitation. After the modifications, the total precipitation for the entire winter used as input for the calculations corresponds to 91.2 mm SWE. The phase of the precipitation was set to a value of 1 corresponding to snow throughout the winter period. Mixed rain-snow events were observed only early in November (1, 2, 7, and 8) and for these days the precipitation phase was reduced to 0.9 (=10% rain). For the down-welling long-wave radiation we used 3h averages from the North American Regional Reanalysis (NARR) data (Mesinger et al., 2006) downloaded from the NARR homepage at the National Centers for Environmental Prediction (wwwt.emc.ncep.noaa.gov/mmb/rreanl/index. html). The horizontal resolution is 0.3 degrees (approximately 32 km) and data from the grid including the field site were selected. Previous investigations revealed that the NARR re-analysis delivers reliable input data for snowpack modeling (Langlois et al., 2009). The snowpack observations included at least weekly measurements of snowpack profiles regarding snow type, snow density, temperature, and SSA (Taillandier et al., 2006). The examination of the snowpack properties included continuous recordings of the temperature using a fixed vertical string of eleven white temperature sensors 75 mm apart. The snowpack height was observed using four different stakes distributed over the sampling site. The results of the observations of the air temperatures at 2 m, the snowpack heights, and the snowpack temperatures for the entire winter season are shown in Fig. 1. The observations demonstrate the insulating property of the snowpack: temperatures at the soilsnow interface always remain well above −10 • C even with air temperatures plummeting to values below −40 • C. Figure 1 indicates the occurrence of trends in the snowpack temperature lasting 5 to 10 days induced by low air temperatures and leading to decreases of snowpack temperatures. Such cooling periods occurred several times: the first one between 9 and 17 November. However, later in the winter season when temperatures dropped to −20 • C, the effect of the cooling of the snow was subsequently visible throughout the entire snowpack. This is most obvious for the period from 5 to 10 March, when a substantial cooling extended to the lowest layers of the snowpack (Fig. 1), although air temperatures were even cooler in periods before.
It must be noted that the snowpack temperatures in Fig. 1 are not corrected for the radiative heating of the temperature sensors during sunlit times that can have an impact on the measurements in the top 20 cm of the snow and can reach several Kelvin close to the snow surface (Brandt et al., 1991). Even covered sensors easily register temperatures higher than snow temperatures in the top few centimeters (Brandt and Warren, 1997). Therefore, temperatures measured close to the snow surface were overestimated. This concerns mainly the measurements in the second half of the winter season, when short-term variations in the snowpack temperature become visible, that were induced by the diurnal cycle of the solar radiation. For example, in February direct shortwave radiation regularly exceeded values of 200 W m −2 .

Snowpack modeling
The physical snowpack model CROCUS (Brun et al., 1989(Brun et al., , 1992 was applied to simulate the snowpack properties. It is the operational tool of the French National Weather Service (Météo France) used to forecast avalanche risks in the French Alps and other mountain regions in France (Durand et al., 1999). It is a one-dimensional multi-layer physical model of the snow cover that explicitly evaluates at hourly intervals the surface mass and energy budgets of the snowpack. Variables used to calculate these budgets include outgoing radiation, turbulent heat fluxes, and the exchange of moisture with the atmosphere. The simulated snow cover consists of multiple homogeneous layers according to snowfall events taken parallel to the slope surface. Mass and energy is exchanged between adjacent layers to account for physical processes such as heat diffusion, transfer of radiation, and latent heat due to sublimation and condensation, and liquid water percolation. Phase changes are taken into account and snow densification and metamorphism are parameterized, affecting mass and energy transfer and changing surface albedo. The model requires topographical input parameters like latitude, longitude, altitude, inclination, exposition and shading by other mountains. It runs with the following meteorological input parameters: air temperature, wind speed, relative humidity, precipitation quantity and phase, incoming solar direct and diffuse radiation, incoming long-wave radiation, and cloud cover. It generates information about the physics and properties of the modeled snowpack including the thickness, the snow temperature, the snow density, the liquid water content, the snow type, the age of the snow for each separate snow layer and the snow surface temperature. The model further calculates budgets of the snowpack like the total height of the snowpack, the run-off, the latent and sensible heat fluxes, and the fluxes of infrared and short-wave radiation.
In the CROCUS model the ground heat flux is imposed depending on the geographic location, the elevation, and the season. The heat fluxes gradually decrease from fall to winter and slowly increase again afterwards. The model also takes into account that additional heat is transferred to the snow due to the release of the latent heat if the soil moisture freezes. Therefore, the simulated total heat fluxes at the soilsnow interface are increased if the temperature at the interface decreases below 0 • C. Besides the heat flux at the soilsnow interface, the thermal budget of the snowpack depends on two further parameters: the thermal conductivity of the www.the-cryosphere.net/4/35/2010/ The Cryosphere, 4, 35-51, 2010 snowpack and the turbulent fluxes of heat and water vapor between the snowpack and the atmosphere. In the model, the thermal conductivity of the snowpack is related to the density of the snow according to the equation reported by Yen (1981): where k eff is the effective thermal conductivity of the snow in W m −1 K −1 , k ice is the thermal conductivity of ice (=2.22 W m −1 K −1 ), ρ snow is the density of the snow in g cm −3 , and ρ water is the density of water (=1 g cm −3 ). However, this relationship is used only for snow densities larger than 0.28 g cm −3 . For snow densities smaller than 0.1 g cm −3 the value of k eff is fixed to 0.1254 W m −1 K −1 . Between 0.1 and 0.28 g cm −3 the thermal conductivity increases linearly with the snow density. Turbulent fluxes of latent and sensible heat between the snow surface and the atmosphere are the third important parameter for the thermal budget of the snowpack. In the CRO-CUS model, these fluxes are parameterized with a trans-fer coefficient C H for the turbulent fluxes using the bulk Richardson number (Essery et al., 1999): where k is the von Kármán's constant (=0.4), z a is the observational height of the wind speed and the air temperature, and z 0 is the aerodynamic surface roughness length. A typical roughness length for snow of z 0 =10 −3 m is used in the model.

Modification of the snowpack model for subarctic snow
CROCUS was originally designed to simulate the physical snowpack properties in the French Alps (Brun et al., 1989(Brun et al., , 1992, and it has been successfully validated against measurements obtained there (Durand et al., 1999;Essery et al., 1999;Boone and Etchevers, 2001). Therefore, in its original form it can be considered as an excellent tool for the simulation of a warm and deep snowpack close to the melting point with regular melting events and features typical for so-called maritime snow according to the snow classification of Sturm  et al. (1995). In addition, the model has successfully been applied outside its originally planned domain of application like in the Arctic and Antarctic and in a tropical mountain region (Dang et al., 1997;Genthon et al. 2001Genthon et al. , 2007Lejeune et al., 2007). Since CROCUS is used here for a subarctic snowpack, the performance of the model was carefully examined. The subarctic snowpack (corresponding to the taiga snowpack in the classification of Sturm et al., 1995) exhibits quite different properties compared to the maritime snowpack: it is a thin and low-density snow cover dominated by depth hoar formed during the winter due to large temperature gradients (Sturm and Benson, 1997). A correct description of many properties of the snow cover depends critically on its thermal budget. The snowpack temperature and the temperature gradient in the snow govern many important processes including melting, re-freezing, and metamorphism of the snow determining the snow type and snow grain size (Brun et al., 1992). Moreover, the parameterizations of the SSA depend directly or indirectly on the thermal budget considering either the temperature profiles or the snow type (see below).
After performing initial simulations with the original CROCUS model large deviations between the calculations and the observed temperatures at the soil-snow interface became obvious. Using the ground heat fluxes based on the geographic location (see above), the calculated temperatures at the interface of the soil and the snow were too high compared to observed temperatures throughout the season. Moreover, the thermal conductivity of the soil normally decreases during the winter season because of the ongoing dehydration. Therefore, the imposed ground heat fluxes were reduced by 30% in December, 45% in January, 75% in February, and 80% in March and April. As a result, the average ground heat flux for the entire season was approximately 1.5 W m −2 . However, after adding the latent heat flux in the calculations, the simulated total heat flux at the snow-soil interface is larger. The monthly averages of the total fluxes begin with 9.7 W m −2 for November and decrease to 2.4 W m −2 for March resulting in a mean value of 6.7 W m −2 for the entire season. As a result the average heat flux due to the latent heat generated by the freezing of the soil amounts to around 5 W m −2 . The largest daily average of more than 28 W m −2 is calculated for 17 November. This large heat flux occurs right after the first period of very low air temperatures below −30 • C, while the snow cover is still relatively thin (around 15 cm). This leads to a substantial release of latent heat in the calculations caused by the extensive freezing of soil moisture. Overall, the calculated heat fluxes at the soil-snow interface are in excellent agreement with previous measurements performed in the same area. Sturm (1991) reported a range of 0 to 30 W m −2 during the winters of 1984 to 1987 at the University of Alaska Fairbanks Agricultural Experiment Station (64.82 • N, 147.87 • W) resulting in average total heat fluxes between 6 and 10 W m −2 for the three winter seasons (Sturm, 1991;Sturm and Johnson, 1991).

Implementation of SSA calculations
Two different parameterizations are implemented regarding the SSA simulations. First, the SSA calculations are based on a set of diagnostic equations using the simulated snow types and snow densities as input variables. The second parameterization relates the decrease of the SSA to the age of the snow, the snow temperature, and the temperature gradient in the snowpack using prognostic equations.
In the first case, different equations are used to calculate the SSA using snow densities for five different snow types as recommended by Domine et al. (2007b): fresh snow, recognizable particles, aged and rounded crystals, aged and faceted crystals, and depth hoar. These equations are at least partly based on measurements in Fairbanks. The different snow types are distinguished using the snow model parameters sphericity, dendricity, and snow grain size (Brun et al., 1992). These modeling parameters are employed to describe the metamorphism of the snow and they evolve depending on snowpack temperatures and temperature gradients. The dendricity value always starts with an initial value of 1 ascribed to fresh snow and decreases with time to zero. Simultaneously, the sphericity is initialized with a value of 0.5 for fresh snow and can either increase to 1 or decrease to 0 depending on the temperature conditions in the snowpack. Finally, when the dendricity reaches zero, the diameter of the snow grains is computed. Full details of this scheme were described by Brun et al. (1992). Since these values are assigned to each snow layer in the CROCUS model, the following criteria were employed to determine the snow type for each layer: fresh snow corresponds to a dendricity higher than 0.7; recognizable particles correspond to a dendricity between 0.2 and 0.7; aged, rounded snow corresponds to a dendricity smaller than 0.2 and a sphericity larger than 0.5; aged, faceted snow corresponds to a dendricity smaller than 0.2 and a sphericity smaller than 0.5; and depth hoar is recognized by a sphericity smaller than 0.3 and grain diameters larger than 1.2 mm. After determining the snow type, the corresponding Eqs. (4) to (8)  where SSA is in cm 2 g −1 and ρ snow is in g cm −3 . Simulated snow density values delivered by the CROCUS model are used in all cases. Domine et al. (2007b) proposed different equations for snow types like depth hoar and faceted crystals depending on the snowpack type. Equations (7) and (8)  the recommended equations for the taiga snowpack (Domine et al., 2007b), which corresponds to the subarctic snowpack encountered in Fairbanks. Domine et al. (2007b) further estimated that the errors of the parameterizations are between 21 and 31%. The second parameterization employs the two Eqs. (9) and (10) proposed by Taillandier et al. (2007) for conditions with either a low (<9 K m −1 ) or a strong (>20 K m −1 ) temperature gradient: where SSA(t) is again in cm 2 g −1 , t is the time in hours after deposition, SSA 0 is the initial SSA of the snow layer considered, i.e. the SSA just after precipitation, and T snow is the snow temperature in • C. Based on the comparison of predicted and experimental values presented by Taillandier et al. (2007), we estimated in average error of around 20% applying the two Eqs. (9) and (10). The two equations probably correspond to the two metamorphic regimes described by Sommerfeld and LaChapelle (1970) as equi-temperature (ET) and temperature gradient (TG) conditions. According to these equations the SSA of a given snow sample decreases continuously with time, while the rate of decrease depends on the age and the temperature of the sample, the temperature gradient and the initial SSA. Thus, we used the following procedure to implement such calculations. First, we calculate the age of the snow in hours corresponding to the elapsed time since deposition for each snow layer. Second, according to the snow temperature and the temperature gra-dient between adjacent snow layers Eqs. (9) or (10) is applied to calculate the SSA decrease SSA(t+ t) according to: where t corresponds to the snow age in hours and t to a time step of =0.25 h. Finally, the SSA decrease SSA(t+ t) is subtracted from the previous SSA value.
To distinguish between ET and TG conditions, we use a threshold of 15 K m −1 for the temperature gradient within the snowpack because it is close to the middle of the range of the boundary between ET and TG conditions of 9 to 20 K m −1 as proposed by Taillandier et al. (2007). Moreover, the same threshold is also employed in the model to allow the development of depth hoar (Brun et al., 1992). Nevertheless, the threshold remains uncertain since no experimental data are available in this range. The temperature gradient is calculated each time before applying Eqs. (9) or (10) for all snow layers taking into account the difference in the temperatures of adjacent layers and their thickness.
At long time scales of the order of several months, Eqs. (9) and (10) lead to negative values without any physical meaning. To overcome this difficulty for the Fairbanks snowpack that transforms almost entirely into depth hoar Taillandier et al. (2007) proposed to calculate the SSA of depth hoar until a value of 65 cm 2 g −1 is reached and to use a constant value of 65 cm 2 g −1 afterwards. This procedure is consistent with their observation that the SSA of depth hoar stabilized near that value (Taillandier et al., 2007). Moreover, as discussed by Taillandier et al. (2007) at this stage the proposed equations are recommended to be applied only to alpine and taiga snowpack types because the validity of the equations has been demonstrated so far only for low-and mediumdensity snow.

Model runs
In this section, we present results of a standard CROCUS model run that is characterized by the above described The Cryosphere, 4, 35-51, 2010 www.the-cryosphere.net/4/35/2010/ parameters for the simulation of the thermal budget of the snowpack including the modified ground heat fluxes. We also present results of sensitivity runs by changing the three parameters (1) ground heat flux, (2) thermal conductivity, and (3) turbulent fluxes to investigate the impact on the thermal budget of the snow. We performed simulations, in which one of the parameters was multiplied by either 0.5 or 1.5. The different model runs are briefly summarized in Table 1.
The results of the different runs regarding snowpack properties like snowpack height, temperatures, and temperature profiles are presented in the following sections. Finally, simulations of the snow density, the SSA, and the depth hoar formation in the standard run are presented and discussed.
To give a quantitative comparison between simulations and observations, we report the correlation coefficient from linear regressions, the root mean squared error (RMSE), or the bias between simulated and observed values. The error bars of the observed and simulated SSA values shown in the figures below represent the error of the measurements and the estimated errors of the SSA parameterization. They do not take into account further errors for the simulated values introduced by differences between observed and simulated values for the snowpack height, the snow density, the snow age, or the snow temperature. Figure 2 shows the time series of simulated snowpack heights. A comparison with the observed heights indicate a good agreement. It demonstrates that the model is well able to reproduce the dynamics of the snowpack height including the decrease in thickness after fresh snow fall events. The linear regression of modeled versus observed snowpack heights SH mod and SH obs results in the following regression line with a correlation coefficient of R 2 =0.940 and a RMSE of 49 mm:

Standard CROCUS run
Nevertheless, in the second half of the winter the simulated snowpack heights are generally lower than the observations. This is expressed by the negative bias of −36 mm. The largest deviation occurs on 19 March, when the simulated snowpack height is almost 100 mm less than observed. One obvious disagreement in the simulations occurs at the end of the season, when the simulated snowpack disappears completely between 5 and 13 April. Such a quick melting of the snowpack was not observed because snowpack heights were still around 150 mm on 14 April. Due to the generally good agreement between simulated and observed snowpack heights, we neglected to normalize the snowpack heights before further comparisons of snowpack parameters. This introduced an additional error in the comparisons, which is more pronounced for the second half of the season. The thermal budget of the snowpack can be described using three different parameters: the snow surface temperature and the temperature at the soil-snow interface (corresponding to the temperature at the upper and lower limits) and the temperature profile within the snowpack. Unfortunately, the temperature of the snow surface was not measured at LARS. Therefore, the comparison of the thermal budget is limited to the temperature at the soil-snow interface and to the temperatures within the snowpack as shown in Fig. 1. Figure 2 shows a comparable plot of the simulated snowpack temperatures for the standard CROCUS run. The overall pattern of the snowpack temperature is reproduced by the model with the fast temperature fluctuations at the snow surface and the damped variations at the soil-snow interface. Nevertheless, distinct differences between simulations and observations remain. These become more obvious in the bottom panel of Fig. 2, which shows differences between simulated and observed snowpack temperatures. These differences indicate that during extended periods the standard model underestimates observed snowpack temperatures in almost all layers. For example, the results indicate a negative average bias for the simulations for the entire winter season increasing from −1.2 K to −3.1 K between the heights of 15 and 37.5 cm (Table 2). Part of these deviations in the top layers of the snowpack can be attributed to errors in the measurements induced by radiative heating of the sensors as described above.
Closer to the bottom of the snowpack the agreement is much better. As a consequence of the insulating property of the snowpack, the ground heat flux constitutes the dominant parameter for the temperature at the soil-snow interface. Thus, the ground heat flux, which has been adjusted (see above) based on the temperature at the soil-snow interface, obviously leads to a good agreement between the simulated and observed temperatures at this interface and also in the lower snow layers (Fig. 2). This agreement is displayed in a negligible bias of −0.02 K between the seasonal average of the simulated and measured temperatures at the soil-snow interface in the standard run (Table 3).
Contrary to the general underestimation of the snowpack temperatures in the simulation, some periods occur when the snowpack temperatures are overestimated. This is particularly obvious for the cooling event between 5 and 9 March. During this period the temperatures between 7.5 and 30 cm are overestimated by 2.6 to 6.6 K (Fig. 2) although snowpack temperatures at these levels are normally underestimated (Table 3). Such efficient cooling events are probably due to strong convective circulation processes within the snowpack, which were extensively investigated at the same location (Sturm, 1991). This circulation strongly enhances the turbulent exchange within the snowpack and is caused by fast changes in air temperature and/or high wind speeds. Accordingly, the periods encountered are characterized by decreases in air temperatures below −20 • C (Fig. 1). However, in all cases the wind speed remained low. The maximum wind speed recorded at ∼3 m height remained be-   leading to an enhanced air flow through the porous snow (Colbeck, 1989). Since such a mechanism becomes only efficient at high wind speeds probably much larger than 5 m s −1 and since the snow surface was always smooth, windpumping can be neglected here. Sturm (1991) described possible further mechanisms like pressure fluctuations caused by gravity waves or synoptic-scale weather systems. Nevertheless, the exact mechanism for such vigorous processes at low wind speeds remains currently unknown and is, thus, not in-cluded in the snow model. As a result, efficient temperature changes throughout the snowpack caused by convection are not well reproduced in the simulations. In addition to the continuous recordings of the snowpack temperature, further individual temperature profiles were obtained manually in snow pits. These observations provided 13 additional profiles in the period between 25 November and 25 March. Therefore, a further comparison regarding the snowpack temperatures is shown in Fig. 3   Period 10/11-2/4 10/11-12/4 20/11-10/4 23/11-7/4 28/11-7/4 26/12-2/4 6/2-2/4  Figure 3 includes measurements of the temperature strings and the manual measurements performed in the snow pits. In general, the agreement between simulations and observations gets better closer to the ground. There, the deviations between simulated and observed temperature profiles are in most cases comparable to the differences of the measured profiles using the two different approaches. The RMSE between simulated and observed profiles using the strings amounts to 4.4, 8.6, and 1.4 K for the profiles shown in Fig. 3. The average of the RMSE for all 13 profiles corresponds to (3.7±2.7) K.

Sensitivity runs
The simulation of the thermal budget of the snowpack depends on three important parameters: the ground heat flux, the thermal conductivity, and the turbulent fluxes. Unfortunately, the simulation of all of these parameters exhibit considerable uncertainty. In the CROCUS model the exchange of latent and sensible heat between the snow surface and the atmosphere is calculated using the bulk Richardson number (Eq. 3). It is well known that the method using the bulk Richardson number generally calculates transfer coefficients that are too low if the small roughness lengths normally assumed for snow covers are used (e.g. Boone and Etchevers,  length (e.g. King and Anderson, 1994;Calanca, 2001;Andreas, 2002). This discrepancy is especially pronounced during stable boundary layer conditions with a smooth aerodynamic flow (e.g. Andreas, 2002), which were typical for the Fairbanks snowpack due to the low wind speeds encountered (Taillandier et at., 2006). In addition, Eq. (1) used in the CROCUS model delivers relatively high coefficients for the thermal conductivity compared to the more recent expression recommended by  for low-density snow. For example, at ρ snow =0.2 g cm −3 k eff corresponds to 0.076 W m −1 K −1 using the equation recommended by  compared to 0.16 W m −1 K −1 obtained with Eq. (1). In addition, for certain snow types only a weak relationship exists between the thermal conductivity and the density of the snow. Relevant for the study here is mainly the case of depth hoar.  reported rather constant experimental values of k eff =(0.072±0.025) W m −1 K −1 for this snow type that is considerably lower than the minimum conductivity coefficient imposed in the model.
To investigate the impact of these different parameters on the thermal budget of the snowpack, sensitivity runs were performed. In these runs, either the heat ground flux, the thermal conductivity coefficient, or the coefficient for the turbulent flux was decreased or increased by 50% (Table 1). For these six different model runs the results concerning the snowpack height, the snowpack temperatures, the temperature profiles, and the SSA calculations were analyzed. Tables 2 to 6 provide a full summary of the results of all runs.
Changing one of the three parameters did not substantially change the agreement between simulated and observed snowpack heights (Table 2). For example, the RMSE varies between 41 and 65 mm in the sensitivity runs compared to 49 mm in the standard run. The simulated melting of the snowpack at the end of the winter season occurred in all cases too early compared to the observations. Moreover, the general underestimation of the snowpack temperatures remained either almost unchanged in the sensitivity runs or increased even further (Table 3). For example, a reduced thermal conductivity leads to stronger temperature gradients in the snowpack reducing differences in the snowpack temperatures at higher levels and increasing differences in the lower layers.
The sensitivity runs further indicate that the temperature at the snow-soil interface is a very sensitive parameter. The seasonal average varies substantially after changing the ground heat flux as well as the thermal conductivity. Decreasing or increasing the ground heat flux by 50% enhances the average bias between simulated and observed temperatures to either −2.6 K or +1.3 K. Nevertheless, changing the thermal conductivity has a similar, yet opposite impact leading to an average bias of +2.0 or −1.6 K ( Table 3). As discussed above, the heat flux from the ground constitutes the most important term in the energy budget at the soil-snow interface. The thermal conductivity controls how efficient the heat from the ground is released upward into the snowpack. However, even changing the turbulent fluxes at the interface between the snow and the atmosphere influences the snow-soil temperature slightly.
Only the increase in the turbulent flux by 50% leads to a slight reduction of the deviations between simulated and observed temperatures almost throughout the snowpack. According to Eq. (3), the increase of the transfer coefficient for the turbulent flux corresponds to an increase of the roughness length from 1 to 4 mm, which is still less than the value of 10 mm applied by Boone and Etchevers (2001) to account for the effective roughness of the snow surface. Even with the increased turbulent flux, the average bias at the highest snowpack level at 450 mm still amounts to −1.6 K ( Table 3). This small improvement in the simulations is also displayed in the calculated temperature profiles. Accordingly, compared to the standard the average RMSE for all 13 profiles decreased by 0.4 K. Similar to the standard run, the cooling events in the snowpack (for example between 5 and 9 March) are also not reproduced in the sensitivity runs.
In summary, the transfer of heat between the snow and the atmosphere at Fairbanks seems to be more efficient than described by the processes taken into account in CROCUS. It is possible to reduce this limitation to some extent by increasing the turbulent fluxes. Nevertheless, a full representation of all physical process occurring in the snowpack would be needed to further refine the simulation of the thermal budget of the snowpack. This may concern the representation of convective processes together with the effect of windpumping on the airflow in the snowpack. For example, Lehning et al. (2002b) proposed a parameterization of the windpumping effect based on wind speed, snow density, pore length, and the ratio of the tortuosity to the porosity of the snow. However, such further model developments are beyond the scope of this study. Figure 4 presents a comparison of simulated and observed snow density profiles corresponding to the temperature profiles shown in Fig. 3. Observed snow densities are relatively low with maximum values around 0.25 g cm −3 , which are typical of the Fairbanks snowpack (Sturm and Benson, 1997). In most cases these low values are reasonably well reproduced by the model. However, especially in the second half of the winter season densities in the bottom layers are overestimated by the model, while densities in higher layers are underestimated. This is probably due to the fact that water vapor transport within the snowpack is neglected in CRO-CUS, while the mass transport in the snowpack at Fairbanks is likely to be efficient due to the strong temperature gradients. Moreover, convective processes that can enhance the heat transfer within the snowpack as described above probably also increase the water vapor transport. The different sensitivity runs show only negligible differences regarding the simulated snow densities. Thus, the parameters analyzed in the sensitivity runs only have a small impact on the simulations of the snow densities since the above mentioned water vapor transport is missing in the model. Since the density is directly used in the diagnostic SSA parameterizations, errors in the simulated densities lead to additional errors in the simulated SSA values. This additional error is largest in the bottom layers of the snowpack and in the second half of the season. Figure 5 presents the simulated SSA values for the entire snowpack. The values are calculated using either the snow density or the snow age as described above. The SSA exhibits highest values in the top layers of the snow that consist of the youngest snow and that have the lowest densities. Large fractions of the snowpack show SSA values between 300 cm 2 g −1 and the imposed lower limit of 65 cm 2 g −1 . Values higher than 600 cm 2 g −1 are only simulated in the top layers for periods of one to two weeks after snowfall events.

SSA
In the calculations based on the density, values continuously decrease from higher to lower layers. In contrast, in the simulations based on the snow age the lowest values are simulated for an intermediate layer. This is due to the fact, that the SSA decreases faster under TG conditions compared to ET conditions (Taillandier et al., 2007). Figure 3 indicates that the simulated temperature profiles exhibit smaller gradients closer to the ground compared to the higher levels. Therefore, the SSA decreases more quickly in the intermediate layer than in the bottom layers. Obviously, this effect is not compensated by the higher snowpack temperatures in the bottom layers. The same effect can be observed in Fig. 6 that presents simulated SSA profiles together with observations for the three snow profiles previously shown. If the SSA is calculated using the snow age, the three profiles show higher SSA values in the lowest snow layer compared to the layers above. The method based on the snow density always delivers profiles with decreasing values closer to the ground, which is in better agreement with the shape of the observed profiles. An increase of the SSA closer to the ground has only been measured on 25 March (Fig. 6), when a small increase was encountered in the lowest investigated layer. All other observed profiles exhibit SSA values decreasing with depths. Figure 7 compares simulated and observed SSA values for all 13 snow profiles studied between 25 November and 25 March. The simulated SSA values are integrated over the same depth range of the analyzed snow samples. If the depth range of the samples included more than one simulated snow layer, thickness-weighted averages were calculated using the SSA values of the concerned layers. The correlation coefficients for the linear regression between simulated and observed SSA values correspond to values of R 2 =0.789 and 0.827 using the snow density and the snow age, respectively. Using the snow density, the calculated SSA values show a bias of +97 cm 2 g −1 and a RMSE of 162 cm 2 g −1 .
They are significantly reduced to +55 and 86 cm 2 g −1 using the snow age to determine the SSA values. Note, that six out of 21 experiments used to develop Eqs. (9) and (10) were performed in the same snowpack at Fairbanks (Taillandier et al., 2007). The comparison indicates a reasonable agreement between the simulations and the observations considering that the accuracy of the SSA measurements was estimated to be 12% (Legagneux et al., 2002). Nevertheless, the simulated SSA values are generally higher than the observed values. The differences between simulated and observed SSA values slightly improve in the sensitivity runs with reduced ground heat flux, increased thermal conductivity, or increased turbulent fluxes. In the first two cases this is probably due to a better agreement between simulated and observed snowpack heights. In the last case, the differences in the thermal budget of the simulations are reduced leading to an improved simulation of the SSA values. This may indicate that an improved representation of the snowpack heights, the snowpack temperatures and, thus, the temperature profiles within the snowpack advances the SSA simulations. Nevertheless, even in the best run regarding the comparison of the SSA values, the improvements in the RMSE for both parameterizations are less than 10% compared to the standard run.
One disadvantage of the application of the parameterization based on the snow density and snow type in SSA modeling becomes more evident in Fig. 8. It shows a plot of the simulated SSA values based on the two different types of parameterization. The SSA values based on the density and snow type can be clearly grouped according to the snow type used. The values decrease for the different types such as fresh snow to recognizable particles to aged snow. Lowest values are calculated for depth hoar. Using the different snow www.the-cryosphere.net/4/35/2010/ The Cryosphere, 4, 35-51, 2010 types as a condition for the application of the Eqs. (4) to (8) in the model must lead to discontinuities because a change in snow type leads to a sharp decrease in the SSA as long as the snow density remains unchanged. Such discontinuities do not occur using the snow age approach. Moreover, the fast decrease of SSA for fresh snow seems to be better represented using the snow age parameterization. Figure 8 shows that the SSA quickly decreases with the snow age by more than 400 cm 2 g −1 , while the density based parameterization leads to decreases of less than 250 cm 2 g −1 . Only if the snow is then classified as recognizable particles the SSA decreases to values below 650 cm 2 g −1 . Similar differences are found for the other snow types as well. The faster decrease of SSA for fresh snow is also obvious in Fig. 5 showing smaller areas of SSA values above 700 cm 2 g −1 if the simulations are based on the snow age. For the subarctic snowpack investigated here, the prognostic equations for the SSA implemented in the current version of CROCUS deliver in the standard run (and in all further runs) a better agreement with the observations than the diagnostic SSA parameterizations (Table 5). The larger deviations for the SSA values based on the density in all model runs are at least partly due to the fact that the simulated snow densities also differ from observations (Fig. 4). Therefore, it can be expected that improving the simulations of the snow densities will also improve the agreement between the diagnostic SSA simulations and the observations. The difficul-ties of the diagnostic equations based on the snow density apply only to snowpack modeling, because in field experiments snow density and snow type can easily be determined with low error margins. In contrast, the snow age is not directly observable in an existing snowpack. It may be estimated if recordings of accumulation are available. The snow age can also be calculated using a snow physics model, if it simulates the development of a seasonal snowpack during the entire winter season or during the (year-round) accumulation period for permanently snow-covered regions.

Snow Area Index (SAI)
The SAI represents the total surface area of the interface between the snow crystals and the air per unit of ground surface area. It is calculated according to: where n is the total number of snow layers, and SSA(i), h(i), and ρ(i) are the SSA, the height, and the density of the snow layer i. Observed SAI values calculated from measured SSA profiles are shown in Fig. 9 (Taillandier et al., 2006  fresh snow in the second half of December. Starting early January, the snowpack height ( Fig. 1) and the SAI developed differently. While the snowpack height still increased until the beginning of April, the observed SAI remained relatively constant within a range of 1000 to 1150 m 2 m −2 . This is due to the aging of the snow leading to faceted snow types with smaller SSA values (Taillandier et al., 2006). Figure 9 also shows daily values of the simulated SAI for the entire winter season obtained from results of the standard run. Like in the observations, the simulated SAI reach high values in late November followed by a subsequent decrease in the SAI to reach a minimum around mid-December. While early in the season SAI values based on the snow age are higher compared to values based on the density, simulated values show only small differences regardless of the SSA parameterizations between late November and late December. Subsequently, the simulated SAI values start to differ again with higher values for the SAI based on the snow density. During this part of the season, experimental SAI values fall within the range of simulated SAI values using the two different parameterizations. In summary, the observations are reasonably well represented by the simulations. This is at least partly due to the fact that the overestimation of the simulated SSA values is compensated by the underestimated snowpack height compared to the observational data.
The simulations indicate further that short-term fluctuations of the SAI are probably larger than captured by the observations. In all model runs each fresh snow event leads to increases of the SAI, which are clearly recognizable in the time series. However, in the second half of the winter season these increases are not sufficient to compensate for the loss of surface area due to the formation of faceted crystals leading to a decoupling of the SAI and the snowpack height (Taillandier et al., 2006).

Depth hoar
Since the subarctic snowpack is characterized by the widespread formation of depth hoar (Sturm and Benson, 1997), we also compare the simulated and observed depth hoar formation. Figure 10 shows a time series of the depth hoar fraction of the snowpack using the snow type observations and the snowpack height of the 13 snow pits. The results demonstrate that depth hoar formation started in the second half of November. The snow pit sampled on 25 November indicated already a significant depth hoar layer. At that time, more than 20% of the snowpack consisted of complete or mixed depth hoar layers. The depth hoar fraction increased quickly and from mid-December until the end of March the observed complete or mixed depth hoar layers constituted around 60% of the entire snowpack. The same seasonal trend is obtained from the simulations. In the standard run depth hoar appears for the first time on November (Table 6) followed by a quick increase of the depth hoar layer (Fig. 10). Throughout the winter season, the calculated depth hoar fraction falls almost always into the range of observations of complete and mixed depth hoar layers. This demonstrates that the model is well able to accurately track snow metamorphism and to predict the formation of the correct snow types. This is especially satisfying since CROCUS has mainly been evaluated using data obtained for the maritime snow type of the French Alps (Brun et al., 1992;Boone and Etchevers, 2001;Durand et al., 1999)

Conclusions
We implemented two different empirical parameterizations for SSA calculations into an existing snow physics model. The comparison with field data from a site near Fairbanks, Alaska with a subarctic (or taiga) snowpack reveal a general overestimation of the SSA values compared to the observations. Nevertheless, a good agreement between SSA simulations and measurements is obtained if the SSA is calculated based on prognostic equations using the snow age. The agreement is less strong using diagnostic parameterizations based on the snow density. Moreover, the latter parameterization leads to discontinuities in the SSA values because different equations are used for different snow types. Changes in the simulated snow type due to the metamorphism of the snow cause abrupt changes in the simulated SSA. However, these disadvantages only apply to snowpack modeling, because in field experiments snow density and type are easily observable. In contrast, the snow age is more difficult to determine from observation of an existing snowpack. Since the applied SSA parameterizations are at least partly based on measurements at Fairbanks (Taillandier et al., 2007), the differences between observed and SSA values reflect deficiencies in the simulation of snow parameters like height, density, and temperature. The accurate prediction of the SSA using the diagnostic parameterization requires a reliable prediction of snow density. Since CROCUS does  (Taillandier et al., 2006). Error bars indicate an estimated average error of 16% (Taillandier et al., 2006). The simulated SAI is calculated using simulated values for the same parameters. The calculations are performed with the SSA obtained either from snow density (blue circles) or snow age (red squares) in the standard CROCUS run. Observed depth hoar fractions are calculated from stratigraphic information obtained from 13 snow pits throughout the winter season and includes the fraction consisting entirely of depth hoar (black circles) and the fraction including depth hoar mixed with other snow types (grey squares). The simulated depth hoar fraction for the standard CROCUS run is calculated using simulated snow types, grain diameter, and snowpack heights.
not describe mass transfer within the snowpack through sublimation-condensation cycles, the densities are not very well reproduced contributing significantly to the error of the SSA simulations. The analysis of the sensitivity runs indicates that the calculation of the thermal budget for the subarctic snowpack still remains challenging. The major deficiency in the calculations concerns the exchange of heat at the snow-atmosphere interfaces and within the snowpack. For a location like Fairbanks that is characterized by low wind speeds throughout the winter season the parameterization of the turbulent transfer using the method based on the Richardson number seems to be too limited. An additional heat transfer process due to convection inside the snowpack (Sturm and Johnson, 1991) is ignored in the current snowpack model and should be taken into account by further modeling studies as well as the transport of water vapor within the snowpack. However, the development of a parameterization of additional processes would need to be constrained by more measurements than those presented in this study. For example, direct observations of heat and water vapor fluxes above the snow surface may help to obtain more robust parameterizations of such fluxes. We tried to overcome this shortcoming of the model by enhancing the turbulent fluxes by 50% to increase artificially the heat exchange between the snow and the atmosphere. However, only a small part of the deviations between observed and simulated snowpack temperatures was compensated.
The results of the sensitivity runs indicate a slightly better agreement between simulated and observed SSA values if the simulated snowpack height and temperatures agrees better with observations. This applies to both parameterizations used. Therefore, it can be expected that SSA simulations improve further if the thermal budget of the snowpack as well as the snowpack height is better represented than in this study. This can be the case for a different location, where the Richardson method delivers better results for the turbulent fluxes. For example, it has been demonstrated that CROCUS is well able to simulate the snow surface temperature and temperature gradients in the snowpack at different sites in the French Alps and Pyrenees (Durand et al., 1999;Essery et al., 1999;Boone and Etchevers, 2001). In contrast, the implemented SSA parameterizations have been validated for the cold and dry subarctic snowpack, while the SSA evolution under warmer conditions may better be described differently www.the-cryosphere.net/4/35/2010/ The Cryosphere, 4, 35-51, 2010 due to different snow metamorphism processes in the presence of liquid water within the snowpack. Data on the development of the SSA in a wet and/or melting snow are currently lacking. As a result, a thorough validation of the SSA simulation for other snow types than the subarctic snowpack and for a melting snowpack is required to verify the robustness of the SSA simulations under different climatic conditions. Regarding the modeling of the sub-arctic snowpack, the implementation of additional processes like those mentioned above to improve the simulations of the thermal budget of the snowpack will certainly result in better SSA simulations. The SSA constitutes a crucial parameter for the modeling of the radiation transfer within the snowpack and the interaction between the snowpack and the atmosphere. It allows a more sophisticated formulation of the wavelength-dependent radiation transfer in the snow. The SSA can be used to compute the effective radius of the snow grains, which in turn delivers reliable information about the optical properties of the snow in the infrared Flanner and Zender, 2006). In further model developments this can be exploited to refine the current approach of the radiation transfer in the CROCUS model that is currently based on constant absorption coefficients using three different bands for the solar radiation (Brun et al., 1989). Such a refinement would introduce a feedback mechanism by which the SSA and the thermal budget of the snowpack would interact in both directions.
The SSA also delivers needed information for the successful modeling of the exchange of impurities between the snow and the atmosphere. For example, the number of volatile or semi-volatile molecules adsorbed in the snowpack can be quantified using the SSA and measured adsorption coefficients (e.g. Herbert et al., 2005;Taillandier et al., 2006, Domine et al., 2007a. Moreover, the majority of the chemical reactions in the snow probably take place at the surface or in the surface layer of the single snow crystals (Domine et al., 2008). While the surface area can directly be used to quantify the available area for heterogeneous reactions, the product of the thickness of the surface layer and the SSA can possibly be considered as a "reaction volume" available for reactions in the snow. Therefore, the implementation of the SSA calculation in the snow physics model CROCUS offers a significant step towards the development of a fully coupled physics and chemistry model of the snowpack.