A factorial snowpack model ( FSM 1 . 0 )

A model for the coupled mass and energy balances of snow on the ground requires representations of absorption of solar radiation by snow, heat conduction in snow, compaction of snow, transfer of heat to snow from the air and retention and refreezing of meltwater in snow. Many such models exist, but it has proven hard to relate their relative performances to the complexity of their process representations. This paper describes the systematic development of an open-source snowpack model with two levels of representation for each of the five processes mentioned above, allowing factorial experimental designs with 32 different model configurations. The model is demonstrated using driving and evaluation data recorded over one winter at an alpine site.


Introduction
Snow on the ground reflects solar radiation, limits surface temperatures, insulates the ground and stores water.These properties have important influences on the meteorology, hydrology and ecology of seasonally snow-covered regions; therefore, representations of snowpacks have to be included in meteorological, hydrological and ecological models.There are many surface mass and energy balance models that include snowpack processes, varying in complexity from simple modifications of land surface characteristics in global climate models (e.g.Cox et al., 1999) to multi-layer snow physics models used in regional avalanche forecasting (e.g.Bartelt and Lehning, 2002;Vionnet et al., 2012).Many studies have compared snowpack model predictions with observations (e.g.Douville et al., 1995;Dutra et al., 2010;Schmucki et al., 2014) and a few have compared multiple models in attempts to understand how differences in model structure and parametrizations determine differences in model performance (e.g.Slater et al., 2001;Etchevers et al., 2004;Essery et al., 2009).Understanding why models of coupled processes with large parameter spaces differ, however, is extremely difficult.Although useful insights have been gained, snowpack model comparisons have generally failed to find clear relationships between model complexity and performance and have failed to find a best model.
Several recent commentaries have discussed how models and data can be better used to develop understanding of complex environmental systems.Larsen et al. (2014) reviewed exploratory and "appropriately minimalist" modelling with simple representations of multiple processes alternately switched on or off in factorial experimental designs to investigate causality in geomorphological systems.Mendoza et al. (2015) argued that relaxing constraints on the choices of parametrizations and parameter values in complex process-based models can increase their agility.Gupta and Nearing (2014) advocated a systems theory approach to model building for hydrological systems, which focusses on process modelling without imposing rigid parametrizations.An approach of this kind has been put into practice by Clark et al. (2015) in developing a common framework within which multiple model representations of hydrological processes can be systematically evaluated.
For snowpacks, Essery et al. (2013) presented a model with a rather ad hoc selection of alternative process parametrizations forming a large ensemble of 1701 model configurations, each configuration having between 9 and 32 parameters.The ensemble was run for four winters at an alpine site in France.No configuration was found to give the best simulations of snow mass on the ground for every winter, but a group of configurations incorporating prognostic equations for snow albedo, density and liquid water content were found to have the best overall performances, and alternative parametrizations of fresh snow density and thermal conductivity were found to have relatively little effect.The robustness of these results in the face of parameter uncertainty was not fully investigated due to the size of the ensemble, and some of the parametrization options were incompatible with each other.This paper now describes the much more systematic development of a new snowpack model called a factorial snowpack model (FSM) with five parametrizations that can be turned on or off independently, giving an ensemble of 32 possible configurations with similar spread but much faster run times than the model of Essery et al. (2013).
The parametrizations used are all simple, and none of them are entirely new; similar parametrizations can be found in the CLASS (Verseghy, 1991), CLM (Oleson et al., 2010), HTESSEL (Dutra et al., 2010), ISBA (Douville et al., 1995;Boone and Etchevers, 2001), JULES (Best et al., 2011), MOSES (Cox et al., 1999) and ORCHIDEE (Wang et al., 2013) land surface models, and more complex parametrizations of the same processes can be found in the Crocus (Vionnet et al., 2012), SNOWPACK (Bartelt and Lehning, 2002) and SNTHERM (Jordan, 1991) snow physics models.The intention of the model development here is to allow for investigations of how snowpack process parametrizations work in combination.This could also provide a framework for evaluations of new process parametrizations within a complete snowpack model.Following a detailed description of the model in the next section, an ensemble of simulations is compared with observations and the influence of each process on the results is determined.

Model building
In the Gupta and Nearing (2014) programme, model building involves construction of a conceptual model for the system of interest, decomposition of the system into subsystems representing its spatial organization, parametrization of the processes linking the subsystems to form a closed set of equations, and specification of computational methods and approximations for solving the equations.By following this formal procedure, assumptions introduced at each stage can be clearly identified.

Conceptual model
A conceptual model can be illustrated by a system diagram, which identifies the boundaries of a system, the fluxes across the boundaries, the state variables of the system and the conservation principles linking the fluxes to changes in the state variables.For the snowpack model in this paper, the control volume is a column of snow of 1 m 2 surface area and height h shown by the system diagram in Fig. 1.The state variables are the ice mass I , the liquid water mass W , the density ρ = h −1 (I +W ) and the internal energy U of the column (internal energy and liquid water mass multiplied by the latent heat of fusion can be combined in a single heat content state variable, but they are kept separate here for clarity).Horizontal homogeneity is assumed; therefore, only vertical fluxes are considered; this limits the model to applications over large areas or at sheltered sites where divergence of horizontal heat and mass fluxes can be neglected.It is assumed that there is no vegetation protruding above the snow, but the influences of vegetation could be represented by adjusting the near-surface driving data accordingly (Hellström, 2000).
Mass is added at the snow surface by precipitation at rate P r and removed or added by vapour flux E to the atmosphere and runoff R b at the base of the snowpack; vapour fluxes between the snow and the ground are neglected.Changes in the combined ice and water mass of the column are given by a conservation equation and further constrained by the conditions I, W ≥ 0. Internal energy change is driven by heat fluxes G s and G b at the surface and the base of the snowpack, respectively, in a conservation equation (2)

Model architecture
A subsystem diagram shows the architecture used to represent the internal structure of the system.Gupta and Nearing (2014) referred to this as a "directed graph" because it can be viewed as a collection of nodes (state variables) joined by directional links (fluxes).The snowpack model is discretized by dividing the snow column into layers as shown in Fig. 2 to represent gradients in the state variables; density and liquid water content can be expected to vary vertically due to compaction of snow and drainage of water, and energy gradients are set up by heating and cooling at the surface.The number and thicknesses ( h) of the layers depends on the depth of the snow: one layer is used if the depth is less than 0.2 m, two layers with a 0.1 m top layer if the depth is between 0.2 and 0.5 m, and three layers with thicknesses 0.1, 0.2 and h−0.3 m at the base if the depth exceeds 0.5 m.Numerical subscripts are added to the state variables to identify their values in the layers.The arrows in Fig. 2 define conventions for the directions in which fluxes are taken to be positive.Precipitation has been divided into rainfall and snowfall in Fig. 2. The solid mass fluxes at the surface are snowfall (or deposition of wind-blown snow) S f and sublimation E. Solid mass fluxes between layers are included because redistribution of mass is required by the discretization when the snow depth changes.The liquid mass fluxes into and out of the snow column are rainfall R f and melt M at the surface and runoff R b at the base of the snowpack; evaporation of liquid water in the snow is neglected.The mass conservation equations become for ice and for liquid water, where R i and S i are liquid and solid mass fluxes at the base of layer i and F i is the rate of freezing for water in the layer.
The surface heat flux is divided into radiative, turbulent and melt components in a surface energy balance equation where L f and L s are the latent heats of fusion and sublimation for water (physical constants and quantities that are assumed to be constant in the model are listed in Table 1).R n is the net radiation absorbed by the surface, H is the turbulent sensible heat flux from the surface to the atmosphere and L s E is the turbulent latent heat flux from the surface to the atmosphere.Advection of heat by precipitation is neglected.
For incoming shortwave and long-wave radiation fluxes SW ↓ and LW ↓ , the net radiation absorbed by a surface with albedo α and Kelvin temperature T s is Penetration of shortwave radiation in snow is neglected and the thermal emissivity of snow is assumed to be equal to 1.The energy conservation equations for the model layers are The internal energy and the temperature of layer i are related by where is the areal heat capacity of the layer.

System parametrizations
Parametrizations are required for calculation of the fluxes in Eqs.
(3) to (7).Snow model parametrizations and parameter values are reviewed in Essery et al. (2013).The five parametrizations that can be switched on or off in FSM are three prognostic equations for the albedo, density and liquid water content of snow, and two diagnostic equations for the dependence of thermal conductivity on snow density and the dependence of turbulent fluxes on atmospheric stability.The parameters that are introduced at this stage in the model development are listed in Table 2. Defaults are set, but, following the recommendation of Mendoza et al. (2015), all of the parameters are adjustable.

Albedo
If the prognostic parametrization for snow albedo α s is switched on, decreasing albedo as snow ages and increasing albedo as fresh snow falls are parametrized by where the timescale τ α has different values τ cold and τ melt for cold and melting snow, respectively, as shown in Fig. 3a.If Although separate energy balances are not calculated for snow and snow-free ground, the effective albedo of patchy snow cover is represented by calculating the albedo in Eq. ( 6) as a weighted average where α g is the measured albedo of snow-free ground and the snow cover fraction as a function of snow depth is Snow of depth equal to parameter h f thus covers 76 % of the ground and depth 2h f covers 96 %.This rapid establishment of unform snow cover is most appropriate for simulations at level sites and on small spatial scales (Niu and Yang, 2007).

Heat conduction
Advection of heat by water movement in snow is neglected and conduction of heat is parametrized by where λ is the thermal conductivity of snow and z is depth below the snow surface.For the discretized model, the surface heat flux is and the heat flux at the base of layer i is where Thermal conductivity is calculated as if the conductivity parametrization is switched on (Fig. 3b) and set to a constant value λ 0 if the parametrization is switched off.Heat flux G b at the base of the snowpack is calculated using the temperature of the upper layer in a fourlayer soil model with heat capacities and thermal conductivities depending on liquid and frozen soil moisture contents as in Cox et al. (1999).

Snow compaction
Snow density is set to a constant value ρ 0 if the prognostic density parametrization is switched off.If the parametrization is switched on, the density of fresh snow is given by parameter ρ f and the rate of density increase is parametrized by where the maximum density ρ max that is approached has different values ρ cold and ρ melt for cold and melting snow, respectively, as shown in Fig. 3c.Higher densities can be reached if liquid water freezes in the snow, but increased compaction at depth due to the overburden of snow is neglected.

Turbulent fluxes
Calculations of turbulent fluxes are driven with measurements of air temperature T a and specific humidity Q a at height z T Q and wind speed U a at height z U .Fluxes of sensible heat and water vapour between the snow surface and the atmosphere are parametrized as and where P s is the surface air pressure, ρ a = P s /(R air T a ) is the density of air and Q sat is the specific humidity at saturation with respect to ice.The transfer coefficient in Eqs. ( 20) and ( 21) is where is the surface momentum roughness length, z 0g is the roughness length for snow-free ground and z 0h = 0.1z 0 is the roughness length for heat and moisture transfer.The balance between shear production and buoyant suppression of turbulence in the atmospheric surface layer is characterized by a bulk Richardson number If the adjustment of turbulent fluxes for atmospheric stability is switched on, the stability factor in Eq. ( 22) is from Louis et al. (1982), as shown in Fig. 3d.The stability factor is set to 1 if the stability adjustment is switched off.

Liquid water
A very simple bucket storage parametrization is used for liquid water in snow.The porosity of a layer with ice mass I i and thickness h i is and the maximum liquid water mass that can be held in the layer at 0 • C is Water in excess of this capacity drains to the layer below, and so on to the base of the snowpack.Water in a layer at a temperature below 0 • C will freeze, releasing latent heat.If the prognostic liquid water parametrization is switched off, rain and meltwater at the snow surface drain immediately to the base of the snowpack and W i = 0 for all layers.

Computational model
The conservation equations and the parametrizations form a set of simultaneous non-linear equations that cannot be solved analytically.Instead, the equations are linearized and used sequentially to update the model state variables over time steps of length δt.
First, the snow albedo is updated and the thermal conductivity and the atmospheric stability factor are diagnosed if the relevant parametrizations are switched on.If α (n) s is the snow albedo at the beginning of time step n, Eq. ( 10) is integrated to give albedo at the end of the time step, where and Next, the surface energy balance equation is solved to find the time step increment in the surface temperature while keeping the driving variables and the exchange coefficient constant.Writing substituting in Eqs. ( 6), ( 15), ( 20) and ( 21) and linearizing in δT s gives and where from the Clausius-Clapeyron equation.Equation ( 5) then gives The surface temperature increment is first calculated assuming no melt (M = 0).If this gives a surface temperature above 0 • C, the snow is melting and the temperature increment is recalculated assuming that all of the snow melts (M = I /δt).If this gives a surface temperature below 0 • C, the snow only partially melts during the time step; the fluxes are then recalculated with T s = T m and the melt rate is diagnosed from Next, the temperatures of the snow layers are updated while keeping their and thicknesses constant.Writing the time derivatives in Eq. ( 7) are approximated to first order by An implicit scheme ) is unconditionally stable and can be written as for a three-layer model.This is a tridiagonal matrix equation, which can be generalized to any number of model layers and solved by standard methods.Flux coupling between the bottom snow layer and the top soil layer is calculated explicitly.Numerical stability is maintained even for vanishingly thin snow by always calculating G s as the heat flux between the surface skin and a level h 1 /2 below the surface.Next, ice is removed from the surface layer by melting and sublimation, and liquid water is added by melting and rain.If the liquid water in a layer exceeds W max , the excess is moved to the layer below or runs off at the base of the snowpack.If liquid water enters a layer with temperature T i < T m , then an amount freezes and the layer temperature is increased by an amount If heat conduction increases the temperature of a layer above T m , an amount of liquid water is produced by melting and the layer temperature is reset to T m .Integrating Eq. ( 20) to update the density of each layer gives New snow is added as a layer with the same temperature as the existing surface layer and density ρ f or ρ 0 , depending on the density parametrization.Finally, the thicknesses of the layers are reset according to the rules in Sect.2.2.Temperatures and masses are assumed to be uniform within each layer and are averaged when layers are combined so that total ice mass I = I i , total liquid water mass W = W i and total internal energy U = C i T i are conserved.

Example results
The well-instrumented and well-maintained Météo-France experimental site at Col de Porte (45.30 • N, 5.77 • E; 1325 m elevation) provides quality-controlled data for driving and evaluating snowpack models (Morin et al., 2012).As an example of model performance, Fig. 4 compares observations and simulations of snow mass, snow depth, albedo, runoff at the base of the snowpack, surface temperature and soil temperature for the winter of 2005-2006 at Col de Porte.The 32 configurations of the model produce ensembles of simulations that encompass the observations.The ensemble spreads are particularly wide for snow mass and snow depth simulations during the melt period.The same winter was simulated with the large 1701-member ensemble in Essery et al. (2013); this produced a similar spread to the 32-member FSM ensemble but with delayed snowmelt due to the use of separate energy balances for snow and snow-free ground.
The influence of a particular process on simulations of a particular variable can be measured by taking differences between simulations averaged over all model configurations that have a process parametrization switched on and averaged over all that have it switched off, as shown in Fig. 5. Switching the prognostic albedo parametrization on gives higher albedos after snowfall, lower albedos for cold aged snow and higher albedos for melting snow, which delays snowmelt and increases snow mass and depth.Later-lying snow then gives a period in April with lower surface and soil temperatures.Switching on the thermal conductivity parametrization only makes a difference if the prognostic snow density parametrization is also switched on, and if so, the variable thermal conductivity is lower than the fixed value with the parametrization switched off until the snow density reaches 300 kg m −3 ; the measured bulk density of the snowpack at   Col de Porte did not reach that level until mid-February in 2006.With lower thermal conductivity, the diurnal range of the modelled snow-surface temperature is increased and reaches the melting point more often, increasing mid-winter melting and decreasing the snow mass.Soil temperatures are increased in winter because of the lower thermal conductivity, and surface and soil temperatures are slightly increased in spring because the surface becomes snow free and can warm to above 0 • C sooner.Not surprisingly, switching on the prognostic density parametrization has a large influence on snow depths, which are increased early in the winter and after snowfall but decreased when the density of partially melted snow exceeds the fixed density used if the parametrization is switched off.Winter soil temperatures are kept higher by the deeper snow, but the snow melts earlier in spring because of the interaction between the density and conductivity parametrizations.Switching on the atmospheric stability adjustment delays snowmelt by limiting the transfer of heat to the snow when the air temperature is higher than the snow-surface temperature.Surface and soil temperatures are decreased throughout the winter, and are also strongly reduced in the snow-free periods beyond the scope of this discussion.Switching on the prognostic liquid water parametrization has the largest and earliest impact on snow mass because it prevents the runoff of surface meltwater from mid-winter melt events.Runoff from notable events in late December and mid-February is suppressed, but runoff is increased in April because the snow melts later.Surface and soil temperatures are increased in winter and decreased in spring.

Conclusions
A snowpack model that can be run in 32 different configurations of varying complexity by switching five process parametrizations on or off independently has been presented.The model performance was demonstrated using driving meteorological data over one winter at Col de Porte.Running the model with every possible combination of parametrizations revealed rich behaviour, with some parametrizations having different behaviours at different times of year or depending on the selection of other parametrizations.All of the processes were found to have important influences on model outputs, and all are subjects of current research; for examples, see Dang et al. (2015)

Figure 1 .
Figure 1.System diagram for a snow column of height h with ice mass I , liquid water mass W and internal energy U .Arrows show mass and energy fluxes at the top and bottom of the column.

Figure 2 .
Figure 2. Subsystem diagram for a three-layer snowpack model with linked conservation equations for liquid water, ice and internal energy.The F , G, M, R and S fluxes represent freezing, heat conduction, melt, drainage of liquid water and redistribution of ice between model layers, respectively.

Figure 3 .
Figure 3. (a) Albedo decay as functions of time for cold snow (dashed line) and melting snow (solid line).(b) Thermal conductivity of snow as a function of density.(c) Snow density as functions of time for cold snow (dashed line) and melting snow (solid line).(d) Atmospheric stability factor as a function of bulk Richardson number.Dotted lines show constant values used when parametrizations are switched off.

Figure 4 .
Figure 4. Observations (black lines) and simulations (grey lines) of snow mass, snow depth, albedo, runoff, surface temperature and soil at Col de Porte in 2005-2006.Dotted lines show the envelopes of snow mass and depth simulations from the 1701-member ensemble described in Essery et al. (2013) for comparison.

Figure 5 .
Figure5.Differences between averages of simulations with prognostic snow albedo (grey lines), thermal conductivity (red), prognostic snow density (black), atmospheric stability (green) and prognostic liquid water (blue) parametrizations switched on or off.
on snow albedo, Calonne et al. (2014) on heat transfer in snow, Morris and Wingham (2014) on snow compaction, Reba et al. (2014) on snow-atmosphere interactions and Wever et al. (2014) on meltwater runoff from snow.A paper evaluating the model configurations with and without calibration of parameters for multiple years at multiple sites is in preparation.

Table 1 .
Physical and model constants.