A two-layer canopy model with thermal inertia for an improved snowpack energy balance below needleleaf forest (model SNOWPACK, version 3.2.1, revision 741)

A new, two-layer canopy module with thermal inertia as part of the detailed snow model SNOWPACK (version 3.2.1) is presented and evaluated. As a by-product of these new developments, an exhaustive description of the canopy module of the SNOWPACK model is provided, thereby filling a gap in the existing literature. In its current form, the two-layer canopy module is suited for evergreen needleleaf forest, with or without snow cover. It is designed to reproduce the difference in thermal response between leafy and woody canopy elements, and their impact on the underlying snowpack or ground surface energy balance. Given the number of processes resolved, the SNOWPACK model with its enhanced canopy module constitutes a sophisticated physics-based modeling chain of the continuum going from atmosphere to soil through the canopy and snow. Comparisons of modeled sub-canopy thermal radiation to stand-scale observations at an Alpine site (Alptal, Switzerland) demonstrate improvements induced by the new canopy module. Both thermal heat mass and the two-layer canopy formulation contribute to reduce the daily amplitude of the modeled canopy temperature signal, in agreement with observations. Particularly striking is the attenuation of the nighttime drop in canopy temperature, which was a key model bias. We specifically show that a single-layered canopy model is unable to produce this limited temperature drop correctly. The impact of the new parameterizations on the modeled dynamics of the sub-canopy snowpack is analyzed. The new canopy module yields consistent results but the frequent occurrence of mixed-precipitation events at Alptal prevents a conclusive assessment of model performance against snow data. The new model is also successfully tested without specific tuning against measured tree temperature and biomass heatstorage fluxes at the boreal site of Norunda (Sweden). This provides an independent assessment of its physical consistency and stresses the robustness and transferability of the chosen parameterizations. The SNOWPACK code including the new canopy module, is available under Gnu General Public License (GPL) license and upon creation of an account at https://models.slf.ch/.


Introduction
In the Northern Hemisphere, around 19 % of the annually snow-covered areas are forested (Rutter et al., 2009).As this type of ecosystem has considerable implications for the mass and energy balance of the surface snowpack (e.g., Harding and Pomeroy, 1996;Otterman et al., 1988), the proper understanding and representation of the snow-canopy interactions is crucial whenever realistic estimates of snow cover and melt dynamics in forested environments are needed.This is specifically of concern for hydrological modeling at all scales, runoff estimates from poorly gauged catchments, flood and drought forecasting, global water budget assessment, and in support of local water resources management including irrigation, provision of drinking water,as well as industrial, touristic, or hydropower applications.
Also, the snowpack insulates the underlying soil from winter cold air temperature, with implications for the ecosystem in terms of vegetation cover and dynamics (Rasmus et al., 2011;Grippa et al., 2005), litter decomposition (e.g., Saccone et al., 2013), or carbon cycling (e.g., Kelley et al., 1968).The representation of this insulation is one of the critical uncertainties of the modeling of the global soil carbon cycle and its evolution in permafrost environments (Lawrence and Slater, 2010;Gouttevin et al., 2012).The northwards migration of shrubs observed in the last decades at high latitudes (e.g., ACIA, 2005) also indicates that snow-forest interactions are to become more and more a concern for climate modeling in the context of global warming.
The insulation properties of snow depend on snow depth and snow thermal conductivity, which in the end relates to the type, characteristics, and spatial arrangement of snow crystals within the snowpack.The realistic description of these parameters can hence be a prerequisite for a reliable representation of soil thermal regime and microbiological processes.Snow stratigraphy is also of concern for specific local activities like reindeer grazing in northern countries (Tyler et al., 2010;Vikhamar-Schuler et al., 2013).At present, to the authors' knowledge, such a description is rarely provided by modeling tools for sub-canopy snowpacks (Rasmus et al., 2007;Tribbeck et al., 2006).
Several processes affect the snow cover in sub-canopy environments when compared to open sites.Snow interception by dense canopies and subsequent sublimation or melt of intercepted snow can reduce sub-canopy snow accumulation by up to 60 % (Hardy et al., 1997).Conversely, canopy shading from solar shortwave (SW) radiation can lead to longer lasting snow cover in forested environments, while enhanced long-wave (LW) emission from sunlit trees with low albedo can have the reverse effect (Sicart et al., 2004;Strasser et al., 2011;Lundquist et al., 2013).In such an environment, effects by topographical shading, solar angle, canopy structure, and understory further complicate matters.Sub-canopy snow is additionally sheltered from wind, thereby experiencing reduced turbulent fluxes.Finally, canopy debris tends to accumulate over the snow and modify its optic properties on the course of the season.
This complexity makes the understanding and prediction of the sub-canopy snow cover evolution a challenging task.As a result, generations of modelers have worked to capture how canopies affect the micro-meteorological conditions above snow, and predict the resulting evolution of subcanopy snow (e.g., Essery, 1998;Pomeroy et al., 1998;Durot 1999;Liston and Elder, 2006;Rutter et al., 2009;Strasser et al., 2011).The first focus of this modeling was usually on snow interception, interception evaporation, and (lack of) snow redistribution, as major features shaping the amount of snow beneath canopies (Essery, 1998;Pomeroy et al., 1998).Later, the sub-canopy or within-canopy micro-meteorology was increasingly refined to include a representation of temperature, radiation, and turbulent fluxes: Durot (1999) initiated a treatment of the meteorological fields provided by the French analysis SAFRAN (Vidal et al., 2010) to translate them into sub-canopy fields.This involved an homothetic transformation of the air temperature, an increase in air humidity, a formulation for turbulent fluxes in canopydampened wind conditions, and the shading from solar radiation through a Beer-Lambert or a linear law.Later, Jansson and Karlberg (2001), Yamazaki (2001), Tribbeck et al. (2004Tribbeck et al. ( , 2006)), Liston and Elder (2006), and Strasser et al. (2011) set up dedicated snow models designed for forest environments, featuring different strengths and weaknesses.For instance, the COUP model (Jansson and Karlberg, 2001) features an advanced representation of snow-canopy processes but lacks a detailed, layered snowpack and the associated physical processes.Oppositely, the SNOWCAN model (Tribbeck et al., 2004(Tribbeck et al., , 2006) ) couples a robust radiative transfer model for canopies to the detailed snowpack model SNTHERM, but their treatment of interception is coarse and experimental.
In 2004 and 2009, two significant international intercomparison exercises were compiled to compare the skills of a broad scope of snow models, ranging from land surface model snow schemes, to very sophisticated models designed for local catchments or point-scale applications (SnowMIP: Etchevers et al., 2004;SnowMIP2: Essery et al., 2008;Rutter et al., 2009).They demonstrated the increasing skill of the snow models to capture the dynamics of the sub-canopy snow cover, but highlighted some remaining challenges: misrepresentation of mid-winter melt events, lack of timetransferability of calibration at forested sites, and difficulty in capturing the maximum snow accumulation in warm environments, sometimes due to unreliable precipitation data.According to Rutter et al. (2009), the former and the latter could be in part due not only to a coarse representation of mixed-precipitation events and rain-on-snow, but also to the misrepresentation of ablation events driven by air temperature rising above 0 • C, when models diverge from observations due to their treatment of sub-canopy long-wave radiation.In a recent publication, Lundquist et al. (2013) pointed to long-wave canopy emissions as the main cause of an early sub-canopy snow melt in snow regions where mean winter temperature exceeds −1 • C and mid-winter melt events are frequent.This effect, and the importance of accounting for the thermal structure of different canopy elements, has been pinpointed before by other observation-based studies (Ellis et al., 2013;Pomeroy et al., 2009;Sicart et al., 2004).However, most snow models of the current generation fail to capture it due to an inappropriate treatment of the canopy thermal regime: most of them, like the COUP model or the SNOW-PACK model before our work (Bartelt and Lehning, 2002;Lehning et al., 2002aLehning et al., , b, 2006)), use the above-canopy air temperature as a substitute for canopy temperature.There are though a few exceptions: in their design of a land surface model dedicated to intensively cold regions, Yamazaki et al. (1992) and Yamazaki (2001) resolve a separate energy balance for two canopy layers (crown and trunks).However, they neither compare their results to canopy temperature or radiation data, nor do they assess the added value of this specific model design for the sub-canopy snow surface energy balance.In the aforementioned SNOWCAN and in Snow-Model (Liston and Elder, 2006), an observed or hypothesized canopy temperature can be used as a model input to compute the thermal emission of the canopy and its impact on sub-canopy snowmelt.However, this is not a comprehensive modeling approach, which would compute the canopy temperature by itself.In ADMUNSEN, Strasser et al. (2011) used the heuristic formulation by Durot (1999) which accounts for thermal dampening by the canopy, but does not propose a physical formulation of the canopy temperature.This approach may show some limitations in specific meteorological conditions or in discontinuous forest where trunks are exposed to direct solar radiation.
Our work here builds on the hypothesis by Rutter et al. (2009).We aim to test how the sub-canopy energy balance is improved by a physically consistent formulation of the canopy thermal structure and the distinction between woody and leafy elements.Focus is mostly on snow-covered environments, but not exclusively.For our purpose, we develop a two-layer canopy representation (leaves and wood) in the aforementioned SNOWPACK model, which proposes a very detailed physical and microphysical representation of the snowpack.Before our work, it included a simple, onelayer canopy module where radiation and precipitation interception by forest elements were represented (e.g., Musselman et al., 2012), but the equivalence between air and canopy temperature was assumed.Because SNOWPACK was initially developed for alpine environment and is still mostly used in alpine or boreal context where conifers are dominant, our new canopy module is for now only suited for needleleaf, evergreen forest.
Micro-meteorological and sub-canopy, stand-scale radiation data collected during the SnowMIP2 experiment provide a proper data set for the evaluation of our developments.We complement them with tree trunk temperature and biomass flux measurements collected in a summer boreal environment, which is an interesting test case for model transferability and robustness.
Our contribution is hence structured in the following way: 1.An exhaustive documentation of the new model and the canopy module is included for the sake of clarity and knowledge dissemination.Earlier versions of the canopy module had been only partially described in Stähli et al. (2006) and in Appendix A of Musselman et al. (2012).
2. Existing simultaneous observations of sub-canopy radiation, snow evolution, and meteorological conditions from Alptal (Switzerland) are used to validate the new model and demonstrate its robustness and improvement over simpler canopy formulations and with respect to observations.
3. Model validity and transferability is finally tested against observations of components of the canopy energy balance taken from a different coniferous environment (Norunda, Sweden).

The SNOWPACK/Alpine3D snow model
SNOWPACK is a one-dimensional, physics-based snowcover model originally dedicated to avalanche risk assessment.Driven by standard meteorological observations, the model describes the stratigraphy, snow microstructure, snow metamorphism, temperature distribution, and settlement as well as surface energy exchange and mass balance of a seasonal snow cover.It has been extensively described in Bartelt and Lehning (2002) and Lehning et al. (2002a, b).Since 2005, it has also included the effect of vegetation above and within or below the snowpack.Snowpack can be wrapped into an open-source, spatially distributed, three-dimensional model for analyzing and predicting the dynamics of snow-dominated surface processes in complex alpine topographies: Alpine3D (Lehning et al., 2006).In addition to SNOWPACK, Alpine3D includes a preprocessing and interpolation module for meteorological fields (Bavay and Egger, 2014), a module computing the spatially distributed radiations as affected by topography (Helbig et al., 2009), an optional snow transport model (Groot Zwaaftink et al., 2011), and an optional runoff model (Zappa et al., 2003;Comola et al., 2015).The interpolated or provided spatial meteorological fields drive the energy and mass balance of the surface snowpack, computed by SNOWPACK.The canopy module and its new features described hereafter can run within Alpine3D.They are included in the SNOW-PACK model version 3.2.1 (revisions from 741) which is available under GPL license and upon creation of an account at https://models.slf.ch/.

The canopy model structure
The canopy module of SNOWPACK calculates the upper boundary conditions for the snowpack or bare soil surface below the canopy.It is based on an energy balance approach in order to be consistent with the distributed radiation scheme used in Alpine3D.Interception and throughfall of precipitation, transpiration, and evaporation of intercepted snow or rain, as well as the influence of the canopy on radiative and

I. Gouttevin et al.: A two-layer canopy model with thermal inertia
turbulent heat fluxes at the snow or soil surface, are included in the model.
In its one-layer version, the model represents the vegetation canopy as a single big leaf with state variables (i) canopy temperature T can (K) and (ii) storage of intercepted water or snow I (mm).All canopy processes are then computed based on three basic input parameters: canopy height z can (m), leaf area index (LAI) or plant area index (PAI) (m 2 m −2 ), and direct throughfall fraction c f (-).PAI has more of a physical sense as non-leafy canopy elements play a role in radiative extinction and turbulent fluxes, but PAI and LAI can usually be derived from each other via a factor depending on stand characteristics; thus, the switch between both just affects parameter values in our formulations.The description here uses LAI; the direct throughfall fraction can be set to zero if LAI is provided as a stand-scale average including canopy gaps of moderate size (up to ∼ 1 m).These three model parameters intend to describe differences between forest stands without further tuning.
The consideration of the thermal inertia of the forest stand in the one-layer version with heat mass (1LHM) and the twolayer version (2LHM) imposes the use of an additional input parameter, the mean stand basal area B (m 2 m −2 ).The different parameters used by the SNOWPACK canopy module are listed in Table 1, distinguishing between the ones to be provided by users according to forest specificities, and the ones internal to the model.
The idea behind the two-layer version of the canopy module is to capture the thermal contrast between two distinct compartments of the canopy: -the upper (or outer) canopy compartment (leaves or needles) which is most directly exposed to the atmosphere; -the lower (or inner) canopy compartment (twigs, branches, trunks, some leaves), for which energy and mass fluxes have already been altered by the upper canopy compartment.
This modeling choice relies on observational data highlighting this contrast and its relevance for the sub-canopy energy balance (Pomeroy et al., 2009).With respect to the one-layer version, one state variable is added, namely, the temperature of the trunk or lower canopy compartment T trunk (K).T can is then replaced by T leaves , the temperature of the upper canopy.
The coupled water and heat balances of the canopy layer are calculated in three steps: 1. First, a preliminary mass balance is calculated including interception and throughfall of precipitation.
2. Second, the canopy temperature T can is calculated by solving the energy balance of the canopy.For this purpose, all the non-linear energy fluxes to the canopy have been linearized in terms of canopy temperature via the Taylor series.The radiation transfer and turbulent exchange of sensible and latent heat are then deduced.For the two-layer version, the energy balance of the upper canopy also includes thermal emission from the lower canopy which is similarly linearized in terms of T can via the explicit formulation of an energy balance for the lower canopy.
3. Third, the mass balance of the canopy is updated by the evaporation (or condensation) calculated in step two.
The two-layer version affects the canopy energy balance and computation of net radiation in each layer.For the sake of simplicity the one-layer canopy module is first fully described.The specificities implied by the consideration of two layers are then dealt with in the last part of this section.

Interception parameterization
The mass balance of the canopy layer includes three fluxes of water: interception of precipitation I (mm day −1 ), interception evaporation E int (mm d −1 ), and water unloading from the canopy U (mm d −1 ): where I (mm) is the interception storage.
A fraction (1 − c f ) of the precipitation P (mm day −1 ) is available for interception at each time step.The interception rate is calculated as a function of canopy storage saturation with an equation originally proposed by Merriam (1960), in the form given by Pomeroy et al. (1998): where the parameter c (-) is a model time-step dependent parameter known as the unloading coefficient.Pomeroy et al. (1998) suggested a value of c = 0.7 appropriate for hourly time steps.Canopy interception capacity I max (mm) is assumed to be proportional to leaf area index: where the parameter i LAI (mm) is either set to a constant corresponding to the interception capacity for liquid precipitation when these occur, or parameterized as a function of snow density during snowfall events, following Pomeroy et al. (1998): (4) Schmidt and Gluns (1991) reported estimates of the parameter i max (mm) for spruce (5.9) and pine (6.6).The density of the intercepted snow ρ s,int (kg m −3 ) is estimated as a function of air temperature (Lehning et al., 2002b).Different values have been reported for the interception capacity of snow, depending on forest type and climate (e.g., Koivusalo and Kokkonen, 2002;Essery et al., 2003).Most important is to recognize the large difference between solid and liquid precipitation.The phase of the intercepted water is assumed to be equal to the phase of precipitation at each time step.A mixture of liquid water and snow can therefore form the interception storage, and unloading proceeds as the interception capacity of the needles decreases with the enhanced density of the intercepted mixture or with a shift towards positive temperatures.
The partition of precipitation into snowfall and rainfall in SNOWPACK depends on available data.Usually precipitation with undistinguished phase is used, and a temperature threshold disentangles the phases with linear or logistic smoothing around the threshold (Kavetski and Kuczera, 2007).When phase information is available and mixed events occur, the interception capacity is calculated according to Eq. ( 4), but using the weighed sum of liquid water and new snow density instead of the density of snow.For rain-only or snow-only events, Eqs. ( 3) and ( 4) are respectively used without change.
Different approaches have been proposed for calculations of snow unload from the canopy: Essery et al. (2003) set the unload rate equal to a fraction (40 %) of calculated melt of intercepted snow.Koivusalo and Kokkonen (2002) assumed that all intercepted snow unloads as soon as the air temperature rises above 0 • C. We have chosen to calculate snow unload U (mm day −1 ) only when the interception storage exceeds the actual interception capacity: which happens when the interception capacity is reduced due to the precipitation of heavy, wet snow or due to an increase in air temperature.Sudden release of a large amount of snow is thus avoided since the intercepted snow density is increased gradually towards the threshold air temperature for snowfall.This is favorable for the numerical stability of the snowpack simulation.This simple parameterization also re-spects the fact that individual branches usually release snow at a time and total unloading of a whole tree is not very frequent.
Throughfall T (mm day −1 ) to the forest floor is thus equal to Evaporation of intercepted water is calculated as part of the canopy energy balance (cf.below) and added to the water balance at the end of the model time step.

Canopy energy balance
The canopy temperature is directly derived from the canopy energy balance.The one-layer canopy module with no heat mass (1LnoHM, e.g., the version used in previous modeling studies; Rutter et al., 2009;Musselman et al., 2012) relies on an assumption of stationarity, whereby net radiation of the canopy R net,can (W m −2 ) is assumed to equal the sum of sensible H can (W m −2 ) and latent LE can (W m −2 ) heat fluxes neglecting any storage or sources/sinks of heat within the canopy: In the new canopy module, one-layer version with heat mass (1LHM), the thermal inertia of trees is accounted for via a biomass storage flux BM can (W m −1 ), modifying the canopy energy balance: (8)

Radiation transfer
A radiation transfer model for a single canopy layer above a snow or bare soil surface has been adopted from Taconet et al. (1986) by Stähli et al. (2009).The model assumes a fractional absorption of radiation in the canopy layer given by the absorption factor σ f (-).A fraction of the absorbed radiation is reflected, as defined by the reflection factors for shortwave (albedo) and long-wave radiation.Radiation transmitted to the surface below is absorbed and reflected according to the corresponding reflection factors for the surface.Following these basic assumptions, and integrating n multiple reflections between the canopy layer and the underlying surface, the net shortwave radiation absorbed by the canopy layer SW net,can (W m −2 ) is given by where SW ↓ (W m −2 ) is the incoming shortwave radiation above the canopy layer, and α can (-) and α surf (-) are the albedo of the canopy and the snow/soil surface below, respectively.The first three terms on the right hand side are the incident, reflected, and transmitted downward radiation with regard to the canopy layer.The remaining three terms are the sums of incident, reflected, and transmitted upward radiation, as a result of multiple reflections between the canopy and the surface below.Equation ( 9) can be simplified to by mathematical relationships for geometric series.The same procedure can be applied for net shortwave radiation absorbed by the ground surface SW net,surf (W m −2 ) which thus can be written as The calculation of the long-wave radiation is further simplified by assuming an emissivity equal to 1, giving the following equations for net long-wave radiation absorbed by the canopy LW net,can (W m −2 ) and the ground surface LW net,surf (W m −2 ): where σ is the Stefan-Boltzmann constant 5.67 × 10 −8 W m −2 K −4 and LW ↓ is the thermal radiation from the sky.Neglecting the emissivity might overestimate the loss and gain of thermal radiation from the canopy.On the other hand, the absorption factor σ f (-) has a similar effect on the net adsorption/emittance, and it may be difficult to separate these two properties.
The net radiation to the canopy is then the sum of the LW and SW net contributions: The albedo of the canopy α can (-) is equal to where f wet (-) is the fraction of the canopy covered by intercepted water calculated as and α wet (-) and α dry (-) are the albedo of wet and dry canopy, respectively.The albedo for the wet part of the canopy can be set differently for liquid and solid interception (Table 1).
The canopy absorption factor σ f (-) is assumed to be equal for long-wave and diffuse shortwave radiation, independent of interception storage and phase, and is calculated as a function of LAI: where k LAI (-) is an extinction parameter with values normally between 0.4 and 0.8.For direct shortwave radiation, it can optionally be a function of solar elevation angle θ elev , following Chen et al. (1997): where θ elev is limited to the range [0.001 − π/2] to ensure a positive value of σ f,dir .Direct and diffuse SW radiations are in this case disentangled by the model after Erbs et al. (1982).
For the sake of completeness, the effective surface albedo, α total (-), and radiative surface temperature, T eff (K), above the canopy layer are given as and respectively.These variables have no influence on the onedimensional simulations presented here, but are used to estimate the contribution of long-wave and shortwave radiation from surrounding terrain when the SNOWPACK model is used within the distributed Alpine3D model.
Finally, the radiation fluxes calculated by the canopy module are only applied to the fraction of the surface covered by the canopy, assumed to be the complement of the direct throughfall parameter: (1 − c f ).An exception to that occurs for direct shortwave radiation which is collimated in the solar direction: when sun is not at the zenith, the sun beams are not parallel to the tree trunks and the projected surface occupied by the canopy along their trajectory is higher than (1 − c f ).This higher fraction of canopy shading (1 − c f,dir ) is derived following Gryning et al. (2001) from the mean canopy height z can (m) and an average canopy diameter D can (1 m by default): In the remaining fraction of the surface, the exchange of long-wave and shortwave radiation between the atmosphere and the ground surface is calculated without influence of the canopy.

Turbulent fluxes
The turbulent fluxes of sensible and latent heat from the canopy to the reference level of the meteorological input (above the canopy) are calculated using the bulk formulation: where ρ (kg m −3 ) and c p (J kg −1 K −1 ) are the density and heat capacity of air, T can (K) is the canopy-layer temperature, T air (K) and e air (Pa) are the air temperature and the actual vapor pressure in the air at a reference level z ref (m) above the ground surface, L (J kg −1 ) is the latent heat of vaporization of water (or sublimation when T air < 273.15 K), R a is the specific gas constant for air (J kg −1 K −1 ), and e sat [T can ] (P a ) is the saturated vapor pressure corresponding to the canopy temperature.Furthermore, the turbulent transfer coefficients for heat and vapor are expressed in terms of the aerodynamic resistances r H (s m −1 ) and r E (s m −1 ) (further described below).Latent heat flux is the sum of transpiration E tr (mm s −1 ) and evaporation of intercepted water E int (mm s −1 ).The partitioning of the components from partly wet canopies can be a delicate problem.To simplify the numerical solution of the energy balance, we have chosen to formulate an effective aerodynamic resistance r E for latent heat calculated as an average of the corresponding values for transpiration r Etr and interception evaporation r Eint , weighted by the fraction of wet canopy f wet : The total evaporation E can (m day −1 ) is calculated directly (Eq.23), and its components are derived as secondary results: The derivation of the aerodynamic resistances for transpiration and interception evaporation is given in the next section.
Transpiration is not allowed if the achieved E can is negative (condensation).In such cases, the solution of the energy balance has to be re-calculated using f wet = 1.At temperatures below the freezing point the modeled canopies do not transpire anymore.If the canopy energy balance forces, through Eq. ( 24), an evaporation that cannot be sustained by the interception storage, the latter limits the possible evaporation and the canopy energy balance is recalculated accordingly.

Aerodynamic resistances
The aerodynamic resistances for sensible and latent heat fluxes are calculated using a two-layer model adapted from Blyth et al. (1999) which for simplicity assumes logarithmic or log-linear wind profiles both above, within, and below the canopy.More elaborate models have been suggested by for instance Shuttleworth and Wallace (1985); however, the remaining uncertainties in the representation of the withincanopy turbulent exchange call for a simple approach.The aerodynamic resistance for scalars from the canopy level, defined by the displacement height d (m), to the reference level of the wind and temperature measurements z ref above the canopy, is calculated as where u * (m s −1 ) is the friction velocity: k is the Kármán constant (0.4), z 0h (m) and z 0m (m) are the canopy roughness lengths for heat and momentum, ψ m (-) and ψ h (-) are functions correcting for atmospheric stability following Högström (1996) and Beljaars and Holtslag (1991).In addition to Blyth et al. (1999), and following, e.g., Koivusalo and Kokkonen (2002), we introduce an additional parameter c h0 (W m −2 K −1 ) representing a minimum heat exchange coefficient for windless conditions.Displacement height, and canopy surface roughness length of momentum and heat are related to the canopy height through the parameters f d (-), f z0m (-), and f z0h/z0m (-) with values

I. Gouttevin et al.: A two-layer canopy model with thermal inertia
given in Table 1: In addition to the resistance between the canopy air (canopy reference level) and the reference level for meteorological measurements (above the canopy), excess resistances from the canopy surface, and from the soil/snow surface (beneath the canopy), to the canopy level are defined as There, a multiplicative increase of the resistance below the canopy f surf (-) is introduced as a function of the leaf area index: with a maximum value of 1 + r a,LAI (-).The excess surface resistance below the canopy, r surf , affects the heat and latent fluxes computed from the ground to the reference level.This resistance is corrected for atmospheric stability by applying the same stability functions as in Eqs. ( 27) and ( 28), but in this case using the temperature difference between the canopy and the snow or bare soil surface instead of the temperature difference between the canopy and the air.With the current choice of parameter values, the excess resistance for the canopy surface is almost zero, but the theoretical framework for a later use/optimization of this parameter based on observational data is set.
In the end, the total aerodynamic resistances for heat from the reference level to the canopy and the ground surface, respectively, are given by The aerodynamic resistances for sensible and latent heat from the ground surface are assumed to be equal.For evaporation from intercepted snow, the resistance from the canopy to the canopy layer can be increased with a factor f ra,snow (-) compared to rain following Lundberg et al. (1998) and Koivusalo and Kokkonen (2002): The total resistance for transpiration also takes the stomatal control into account: where the stomata resistance r stomata (-) is calculated as a function of a minimum resistance r smin (-), incoming solar shortwave radiation, vapor pressure deficit and soil water content soil as suggested by Jarvis (1976), and soil temperature T soil following Mellander et al. (2006) and Axelsson and Ågren (1976): The functions f 1 -f 4 in Eq. ( 39) all take values between 0 and 1, specifying optimal conditions for root water uptake corresponding to the response of the leaf stomata to conditions in the atmosphere and the root zone.

Biomass heat flux
Due to their thermal inertia, trees can store energy over periods of high exposure to solar radiation, and release it at night.This biomass heat flux is accounted for in the 1LHM version of the canopy module via the areal heat mass of trees HM can (J K −1 m −2 ): where T t can (K) and T t−1 can are the canopy temperature at the model t and t−1 time steps, and t (s) is the model time step.HM can is here derived from parameters commonly observed by foresters: LAI, mean stand basal area B (m 2 m −2 ) and mean canopy height (z can ).The leaf thickness e leaf (m), biomass density ρ biomass (kg m −3 ) and biomass specific heat mass C p biomass (J kg −1 K −1 ) are fixed parameters with values 10 −3 , 900 and 2800, respectively (Lindroth et al., 2010).In Eq. ( 43), the volume of woody biomass (referred to as "trunk" but comprising trunks and branches assimilated to the lower canopy layer) is calculated from mean tree basal area and height assuming a conical profile for trunks.In this study, areal heat masses will be expressed as "water equivalent areal heat masses" HM eq (kg m −2 ), e.g., as the areal mass of water yielding the same heat mass than HM (J K −1 m −2 ): where C p water = 4181 J kg −1 K −1 is the liquid water specific heat mass.

Two-layer canopy version
With respect to the one-layer canopy module, the two-layer formulation induces changes in the formulation of radiative transfer, turbulent and biomass fluxes, and in the end the energy balance of the canopy.These differences are the focus of the present section, whereby the upper canopy layer is equivalently referred to as "leaves" while the lower canopy layer is labeled "trunk".The formulation of the radiative and turbulent components of the two-layer module is illustrated in Fig. 1.

Radiative transfer
In a real forest the trunk-layer intercepts parts of the shortwave and long-wave radiation transmitted, reflected, and emitted by the uppermost canopy layer and upwelling from the soil surface.
Our model features a simplified representation of the following: -For SW radiation, only the transmitted radiation from the upper canopy (with absorption factor σ fleaves and albedo α leaves ) are intercepted or reflected by the trunk layer (with the respective factors σ ftrunk and α trunk ).Radiation undergoing multiple reflections between ground surface and upper canopy are unaffected by the trunk layer (Fig. 1).The SW flux reaching the ground and both canopy layers are expressed accordingly: Obviously, the biomass responsible for SW and LW extinction now has to be split into the two canopy layers so that the total extinction for SW is similar in both versions.Equating the first-order radiation from Eqs. ( 11) and ( 47) yields or equivalently, based on Eq. ( 17): where LAI leaves and LAI trunk are the respective portions of the total LAI attributable to the upper and lower canopies.We denote hereafter and express the leaf-layer and trunk-layer absorption factors as functions of LAI and f LAI : Similarly to the one-layer version (Eq.18), these factors can be adapted to enhance absorption of direct SW radiation based on solar elevation angle.
f LAI is an a priori undetermined parameter of our model due to the difficulty of deriving it from existing data sets for different forest types and structures.In Sect.4, we show that the calibration of the model at Alptal against this parameter yields f LAI = 0.5, which means equal contribution from the woody and leafy parts of the forest to shortwave extinctions.This value is adopted as default value in the model (see Sect. 5 for discussion).
-For LW radiation, the choice of an emissivity of 1.0 for ground and canopy suppresses multiple reflections.Thermal emission from the upper canopy layer and from the ground is attenuated by the trunk layer with the same absorption factor as for SW radiation σ ftrunk .The trunk layer then radiates thermally towards the ground and the As for the one-layer version, this radiation balance is only valid on the canopy-covered fraction of the model grid cell, which is (1 − c f ) for diffuse SW radiation and LW radiation, and (1 − c f,dir ) for direct SW.

Turbulent fluxes
Sensible heat exchange between the lower or upper canopy layer and the atmosphere is parameterized the same way as in the one-layer model version, e.g., via the resistance r H,can .We consider that latent heat exchange between canopy and atmosphere only occurs through interception evaporation and transpiration at the leaf level, e.g., via the upper canopy layer only.

Biomass heat flux
The upper and lower canopy layers are respectively attributed the HM leaves and HM trunk heat masses from Eqs. ( 42) and ( 43) which are used in the biomass heat flux parameterization Eq. ( 40) in the place of HM can .

Energy balance
An energy balance is formulated separately for each layer according to the energy balance equation with heat mass (Eq.8), where all terms are linearized as functions of T leaves and T trunk .The coupled system is then iteratively solved for both temperatures.
The values of all the model parameters as used in the SNOWPACK canopy module are listed in Table 1.

Data
The data from two field sites are used here.

Alptal site
The first data set is from the Alptal forest site (47 • 03 N, 8 • 43 E; Erlenbach sub-catchment, Switzerland; site 1012 in the Fig. 1 of Stähli et al., 2006) that served as test site for the SnowMIP intercomparison study (Rutter et al., 2009) and builds on a long tradition of snow and meteorological investigations (e.g., Stähli et al., 2006Stähli et al., , 2009)).The site features an ∼ 11 • west-orientated slope at 1185 m a.s.l. and is dominated by Norway spruce (85 %) and silver fir (15 %), with a basal area of 41 m 2 ha −1 and a maximum height of typically 25 m.The site LAI (including slope corrections and corrections for clumping) ranges from 3.4 to 4.6 with a mean value of 3.9 m 2 m −2 (Stähli et al., 2009).
At this site, the SNOWPACK model is run using meteorological data derived from observations: -downward shortwave and long-wave radiation measured on a 35 m high mast above the canopy: a heated, non-ventilated CNR1 from Kipp and Zonen ( 2002) comprising two pyranometers CM3 (for SW) and two pyrgeometers CG3 (for LW); -precipitation measured by a heated gauge placed at 25 m height on the high mast: the highest trees providing sheltering similar to a fence; -wind speed recorded by a cup anemometer at 35 m on the mast; -air temperature measured at 35 m by a ventilated thermo-hygrometer Thygan (Meteolabor) also integrating a dew point hygrometer; -relative air humidity at 35 m height, derived from the air temperature and dew point.

Validation data include
-Downward SW and LW radiation measured below the canopy (LW ↓BC , SW ↓BC ) by a second CNR1 radiation sensor, but mounted on a carrier constantly moving along a 10 m transect at 2 m altitude above ground at 1 m min −1 speed.This transect was previously shown to have a representative LAI for the stand (Stähli et al., 2009).Great care was put in the collection and preprocessing of this data set, as below-canopy SW radiation is typically close to zero.This effort is well described in Stähli et al. (2009).
As a post-treatment to this data, LW radiation was masked in cases when snow interception on the sensor was suspected.A typical case is illustrated in Fig. 2: from the evening of 19 February to midday of 21 February; the radiation measured by the heated pyrgeometer is close to the emission level of a blackbody at 0 • C (snow emissivity is around 0.98), whereas the air temperature is much colder and modeled canopy temperature closely follows the air temperature signal.The precipitation record (Fig. 2b) features almost continuous snowfall over that period.It is hence suspected that the measured radiation originates from snow at temperature -Snow depth, snow density, and snow water equivalent (SWE) that were measured below the canopy on a weekly basis, at 1 m intervals along a 30 m transect adjacent to the trajectory of the radiometer carrier.More details of the exact procedure are available in Stähli et al. (2009).We use the spatial average of the measurements to come up with stand-representative values.
Meteorological and validation data are available for four consecutive winter seasons between 2003 and 2007.

Norunda site
The second data set is from the Norunda forest site (60 • 05 N, 17 • 28 E), located in a quite level region about 30 km north of Uppsala, Sweden, at 45 m a.s.l.Since June 1994 it has been equipped with meteorological instruments which were complemented by biomass thermometers in June and July 1995.The forest stand is composed of Scots pine (61 %), Norway spruce (34 %), and birch (5 %) with a stand LAI between 4 and 5 m 2 m −2 , a mean basal area of ∼ 34.7 m 2 ha −1 , and a maximum tree height of ∼ 28 m.At this site, SNOWPACK is driven by observed meteorological variables: -Downwelling LW and SW radiation are measured by a combination of a ventilated CM21 pyranometer (Kipp and Zonen) placed at 102 m above ground at the top of a Fluxnet tower (http://fluxnet.ornl.gov/site/730)and a ventilated LXV055 net radiometer placed at 68 m on the same tower.
-Air temperature is recorded at 37 m height above ground by a copper-constantan thermocouple placed in the ventilated radiation shields.
-Air humidity is measured at 28 m by a HP100 TST probe (Robotronic).
-Wind speed is recorded at 37 m by a sonic anemometer.
-Precipitation data were unfortunately not available at the site.We therefore made use of precipitation data recorded at the Uppsala Aut WMO station (WMO no.2-462) openly provided by the Swedish Meteorological and Hydrological Institute (SMHI; http: //opendata-catalog.smhi.se/explore/).This station is 26 km away from the Norunda site and the nearest station in operation at the time of the measurements used here.
The specificity of the Norunda site lies in the continuous measurement, over a summer, of the biomass temperature at different heights and depths within the trunks and branches of the dominant tree species: pines and spruces.They were complemented by a detailed calculation of tree-level and stand-level biomass heat storage, which builds a unique data set to evaluate a physics-based canopy model with heat mass.The details of the tree temperature measurements and heat storage calculations can be found in Lindroth et al. (2010).
In the present study we make use pine trunk temperature at 1.5 m height, which has been measured close to the trunk surface (1 cm deep within the bark).Indeed, we are mostly interested in the ability of the model to reproduce the trunk surface temperature which generates the thermal emission of the trunk layer.We also provide an assessment of the canopy energy balance modeled by SNOWPACK by comparing the stand-scale modeled biomass storage flux to the one inferred from observations by Lindroth et al. (2010).

Methods: Model calibration
Three versions of the canopy module, corresponding to activation of the different features of the new developments (bi-layered canopy and heat mass; Table 2), are calibrated at Alptal in order to evaluate the model in its best-performance setup.Calibration is performed against the observed incoming long-wave and shortwave radiation below the canopy (LW ↓BC , SW ↓BC ).The former is specifically affected by our new developments.the treatment of rain-on-snow events and the parameterization of interception) could compromise a proper calibration of the canopy module.
Depending on the version, one or two model parameters are calibrated, consistently with our modeling choices: k LAI and/or f LAI (Table 2).
Canopy heat mass also affects the LW radiation downwelling to the ground surface.Heat mass is a physical property of a forest stand, and not a free parameter of the model.However, its value is difficult to measure and our model only proposes a coarse estimation of it (see Sect. 2).In each of the versions with heat mass, we therefore try to optimize its value considering it as an additional calibration parameter (versions 1LHM * and 2LHM * ; Table 2).This procedure is designed to assess the physical consistency of our formulation, by comparing its performance to results obtained with unrealistic heat mass values.
Calibration is performed by minimizing the error function CC (calibration criterion) which is the sum of the model-todata RMSE (root mean square error) and MB (mean bias) for the two observed variables LW ↓BC , SW ↓BC .

CC = |MB(LW ↓BC )| + |MB(SW ↓BC |
(56) We prefer CC to the more common Nash and Sutcliffe (1970) efficiency (NSE) because LW ↓BC and SW ↓BC exhibit a strong diurnal cycle: for such cyclic variables, even a lowperformance representation of the cycles yields a high NSE, and the NSE sensitivity to further improvements is typically low (Schaefli and Gupta, 2007).For all versions, the calibrated extinction coefficient k LAI is within the [0.4-0.8]range of expected values (Stähli et al., 2009).Both LW ↓BC and SW ↓BC are affected by k LAI but LW ↓BC is less sensitive to radiation extinction (as atmospheric LW extinction by canopy is partly compensated by canopy thermal emission in the same range of magnitudes).k LAI is therefore mostly determined by calibration against SW ↓BC and is the same for most versions, which differ only in their modeling of LW ↓BC .

Model calibration
The calibration of the f LAI parameter yields the value of 0.5.f LAI partitions the forest LAI between the uppermost and lowermost canopy layers in the 2LHM version.The value of 0.5 would have been a reasonable first choice for that partition.
The successive addition of heat mass (1LHM) and a twolayer partition of the canopy (2LHM) to the default 1LnoHM simulation improves the general model performance, as reflected in the decrease of the CC error function and its components (MB, RMSE).
In the two versions where canopy heat mass is optimized (1LHM * , 2LHM * ), optimization yields unrealistically high heat mass values (HM = 90 kg m −2 and HM = 60 kg m −2 , respectively, whereby field data indicate 30 kg m −2 ).However, while optimizing heat mass quite significantly improves the performance of the one-layer versions (from CC = 23.6 W m −2 for 1LHM to CC = 19.3W m −2 for 1LHM * ), it only marginally affects the performance of the two-layer version (from CC = 18.4 W m −2 for 2LHM to CC = 17.5 W m −2 for 2LnoHM).These are encouraging results for the two-layer canopy formulation: on the one hand, this model version shows a better performance than the onelayered canopy model, even with the physically estimated heat mass.With the one-layered version, such a performance can only be approached with an unrealistic canopy heat mass.On the other hand, the performance of 2LHM shows a considerably reduced sensitivity to the prescribed areal heat mass of the canopy, a physical parameter which can be spatially variable and hard to retrieve with precision in noninstrumented forests.
For the two snow seasons when all model versions simulate a reasonable dynamics for the Alptal snowpack (2004-2005and 2006-2007;see Sect. 4.1.3),the sensitivity of the modeled SWE to the calibration parameters k LAI and f LAI was determined (not shown).The modeled SWE is sensitive to k LAI and f LAI , but the calibrated (k LAI , f LAI ) values lead to RMSE to observed SWE close to the absolute minimum obtained when varying k LAI and f LAI over their full range (13 mm vs. 10 mm at minimum for 2006-2007).In the surrounding of the calibrated (k LAI , f LAI ) values, the modeled SWE has furthermore a reduced sensitivity to variations in k LAI and f LAI .This result enhances our confidence in the model robustness.The performance of all model versions after calibration over 2003-2004 slightly degrades over the longer 2003-2007 time period when observations are available.Especially the MB in LW ↓BC , and (to a smaller degree) in SW ↓BC , are increased over 2003-2007, questioning the transferability of our 2003-2004 calibration.We therefore calibrate the 1LnoHM, 1LHM, and 2LHM versions over the 2003-2007 period and analyze the changes in best-fit parameters and performance (Table 4).
The calibration over 2003-2007 yields a slightly different best-fit parameter value for the extinction coefficient in the 1LHM and 2LHM versions (k LAI = 0.85 vs. k LAI = 0.75 when calibrated over [2003][2004]: this enhanced radiation extinction improves the MB for SW ↓BC over the 2003-2007 period, but slightly degrades the results over 2003-2004.The overall picture is however not changed after this new calibration: -Over both periods, 2LHM performs better than 1LHM which also performs better than 1LnoHM: this is an indication of the added value of our new parameterizations.
-For all model versions, performance is better over 2003-2004 than over the full 2003-2007 period, especially for LW ↓BC .This may indicate that our model is still too simple to capture the full range of snow-forest processes.
-Over both periods, the two, slightly different calibrations yield thoroughly comparable model performances.This gives confidence in the validity of our calibration and in the possibility of calibrating the model over only 1 year of data.
In the simulations discussed in the rest of the paper, calibration over 2003-2007 is used.

Model evaluation against thermal radiation
In Fig. 3 are compared observed and modeled LW ↓BC as computed by the different model versions without heat mass optimization (1LnoHM, 1LHM, 2LHM) over 2003-2004.Similarly to the performance metrics of Table 3, this figure illustrates gradually increasing model performances from the 1LnoHM to the 2LHM model versions.
With respect to 1LnoHM, the consideration of tree heat mass in 1LHM slightly delays and reduces the canopy cooling at night and warming up in the morning: this translates into a slight delay and smoothing of the diurnal cycle of LW ↓BC , part of which originates from canopy thermal emission.More striking, however, is the attenuation of the daily amplitude of LW ↓BC induced by 2LHM, which brings the modeling results in closer agreement to observations: especially, the nighttime (18:00-06:00 LT) mean bias in LW ↓BC is considerably reduced in 2LHM with respect to other model versions, amounting to −10.8, −7.8, and −2.8 W m −2 in 1LnoHM, 1LHm and 2LHM respectively.
When only one bulk layer of canopy is considered, this layer is exposed at night to intense radiative cooling towards the sky, whose thermal emissivity is low.With two layers of canopy, only the uppermost layer experiences this uncompensated cooling.The lower layer receives thermal radiation from the upper layer which has a higher emissivity than the sky.This thermal sheltering yields higher temperature and LW emission at night from the lower canopy towards the ground surface.This mechanism proves to efficiently reproduce the daily cycles (Fig. 3a, b) and daily averages (Fig. 3c) of the thermal radiation affecting the snowpack.

Impact on the underlying snowpack
Over the four winters of interest here, a similar ranking of sub-canopy SWE modeled by 1LnoHM, 1LHM, and 2LHM is observed, with 1LHM accumulating more snow and 2LHM generally featuring the smallest SWE (except for the 2005-2006 winter; Figs.3c and 4).With respect to the thermal behaviors of the different model versions, such a result is somehow counterintuitive as 1LHM and 2LHM generally deliver greater amounts of LW radiation to the snowpack than does 1LnoHM (Fig. 3c.), hence contributing more energy to mid-winter ablation events (e.g., Fig. 3d, December to January).In 1LHM, this increased ablation is, however, compensated by a side effect of the thermal canopy mass: as a result of the high thermal mass of the bulk canopy in 1LHM, the canopy temperature and hence interception evaporation is reduced, and more snow unloads than in the two other versions, resulting in higher sub-canopy snow accumulation.In 2LHM, the high diurnal temperature variations of the upper canopy temperature combine with stronger LW radiation to the snowpack, resulting in a thinner snowpack.
Noteworthy, the model ability to represent SWE (as typically assessed by the RMSE to observations) is degraded in 1LHM and improved in 2LHM with respect to the original canopy module 1LnoHM.Hence, the LW-enhanced ablation in 2LHM (and small associated changes in interception evaporation) does not deteriorate the overall model skills.
In some specific ablation periods, 2LHM also proves to reproduce the observed snowpack dynamics better: one such event is the early February 2004 severe ablation, when high thermal exposure of the snowpack is better reproduced by 2LHM (Fig. 3c) while the concomitant ablation is also stronger in 2LHM, which matches the observations better (Fig. 3d).Similarly, the LW-enhanced ablation in 2LHM leads to a sub-canopy SWE dynamics in closer agreement with observations in the 2005 ablation phase and in early 2007 (mid-winter complete snow disappearance).
As mentioned in the methods, we do not trust the modeling of the accumulation phases, where high uncertainties in precipitation phase and interception (enhanced by the warm temperatures at Alptal) can initiate a permanent bias in the modeled snow cover.Therefore, we do not use observed SWE records to validate or evaluate our model.However, the capability of a model to better reproduce observed, relative ablation events when precipitations are absent (like in the 2005 main ablation phase) reliably means enhanced performances: our results are therefore encouraging for the overall consistency of the 2LHM canopy module.

Norunda: tree temperature and biomass storage flux
At the Norunda site, SNOWPACK is run using the Alptal calibration from 2003 to 2007, and a canopy basal area and areal heat mass derived from local data (Sect.3).The difference in latitudes (hence in solar angle), tree species (mostly Scots Pine at Norunda), and context (Alpine winter vs. boreal summer) between both sites constitutes a huge challenge and an excellent benchmark to test one desired feature of a physically based model, e.g., its transferability to different climate and ecosystem types.We here specify that SNOW-PACK includes all the necessary features to be used as a soilvegetation-atmosphere transfer (SVAT) model in the absence of a snow cover: a soil water balance, a surface and canopy energy balance, and a temperature diffusion scheme in the soil.The model has also been used as such in continuous multi-year simulations in previous studies (e.g., Bavay et al., 2013).We compare observed tree trunk temperature to modeled temperature of the bulk canopy (for 1LnoHM and 1LHM) or of the lower trunk layer (for 2LHM) over summer 1995 at Norunda (Fig. 5; Table 5).The modeled trunk-layer temperature of 2LHM shows an improved ability to reproduce the observed tree trunk temperature signal: similar to the improvements seen at Alptal, radiative loss of energy from the lower layer at night is considerably reduced with 2LHM, bringing nighttime modeled temperature in closer agreement to observed data at Norunda.Also, the reduced SW insolation received by the lower canopy layer during daytime in 2LHM prevents too high midday temperature of the trunks, an observation that 1LHM and 1LnoHM cannot reproduce.Finally, the combination of thermal sheltering of the lowermost canopy layer and its thermal inertia delays the tree trunk cooling (warming) at evening (morning) times, improving the temporal correlation with observations.Heat fluxes to canopy elements are a substantial, though not dominant, component of the canopy energy balance (Lindroth et al., 2010, their Fig. 6): they can amount to ∼ 7 % of the daily net radiation received by the canopy.To assess the consistency of the SNOWPACK canopy module, we compare the modeled canopy heat fluxes to the ones derived by Lindroth et al. (2010) from field measurements and extrapolated at the stand scale.Note that 1LnoHM, having no heat mass, does not consider any such fluxes.
Both the 1LHM and 2LHM versions overestimate the daily amplitude of biomass heat fluxes with respect to observations, with an increased bias for 1LHM (Fig. 6; Table 5).This is in line with an overestimation of the daily amplitude of canopy temperature (or of the temperature of the lower canopy layer for 2LHM) which is stronger with 1LHM (Fig. 5).Also, the model biomass heat fluxes peak ∼ 2 h earlier than the observed ones.We interpret this as an artifact of modeling the canopy with only one or two thermally homogeneous layers, whereas it is in reality a continuous medium experiencing thermal diffusion at scales smaller than our layers.In reality, the low thermal inertia of a bark surface layer provokes quick surface heating as a result of solar energy input (e.g., in the morning).This temporarily limits further heating from turbulent and radiative fluxes, until the surface heat has diffused into the trunk.In other words, heat uptake by trunks is diffusion limited.Contrarily, the bulk, thermally inert trunk layer of our model heats up to a smaller temperature because the heat flux is accommodated by the whole layer and not only by its uppermost surface: further heating by turbulent and radiative fluxes is then still possible and the heat flux towards the biomass keeps being sustained.As a result, our modeled canopy accommodates incoming energy more rapidly than a real one during the first part of the diurnal cycle.The aforementioned mechanism can also cause the accommodation of more heat energy by the modeled trunk layer than in reality: in reality, the capacity of the canopy to accommodate heat is limited by thermal diffusion within the wood.Heat uptake stops when available solar energy starts going down.At that time, the wooden medium may not have reached an homogeneous, high temperature yet (e.g., Fig. 1 from Lindroth et al., 2010).
As such, the representation of the biomass storage fluxes by 1LHM and 2LHM yield only a moderate improvement to the model: they feature a reasonable (though slightly shifted) diurnal cycle (cf.the correlation coefficients in Table 5) but their RMSE to observations is on the order of magnitude of the standard deviation of the observed biomass fluxes (Table 5, first row).
However, model performance, especially for 2LHM, is improved if the total heat-storage flux towards the biomass and canopy air space is considered (thick black line in Fig. 6; Table 5).The air heat-storage flux corresponds to the changes in latent and sensible heat stored in the within-canopy air space.Lindroth et al. (2010) provided estimates of these heat storage terms based on air temperature and humidity measurements at seven heights within the canopy air space.On a daily basis, the air heat-storage term reacts more rapidly to solar heating than the biomass heat-storage flux.The air heat-storage flux is not specifically accounted for in SNOW-PACK.However, the increased correlation coefficient and reduced RMSE, obtained when the SNOWPACK canopy heat flux is compared to the sum of estimated air and biomass heat fluxes, indicate that the canopy module produces a bulk representation of the observed fluxes.Such a result should be confirmed against further observational data sets.

Discussion
Our results show that the new features implemented in the SNOWPACK canopy module, especially the two-layer scheme, improve the representation of the radiation budget at the sub-canopy level.The importance of accounting for the canopy temperature and for contrasts between different canopy elements has often been underlined (Sicart et al., 2004;Pomeroy et al., 2009).Pomeroy et al. (2009) compared three sub-canopy thermal irradiance models based on measurements of (i) air temperature only; (ii) air, trunk and needles temperatures; and (iii) air temperature and empirical shortwave-long-wave conversion function by the canopy.The two latter formulations exhibited distinctively better per-formances than the first one in terms of mean bias and RMSE in both the uniform and discontinuous stands investigated.Our present work confirms these findings also propose a seamless physics-based canopy model to account for this effect.This has to our knowledge never been brought to the scientific literature.
Radiation can be an important driver of the springtime subcanopy snow energy balance and subsequent melt.Garvelmann et al. ( 2014) reported a contribution of about 50 % from the net long-wave radiation to the sub-canopy energy balance during two cloudy-sky rain-on-snow events closely monitored in the Black Forest, Germany.In open environments the long-wave contribution to the surface energy balance does not exceed 20 % for the same meteorological conditions, and is at times negative.According to Lundquist et al. (2013), forest regions with average December-January-February (DJF) temperatures greater than −1 • C experience 1 to 2 weeks reduction in snow-cover duration compared to adjacent open areas because of increasing long-wave radiation from the canopies.There, the sub-canopy snowpack is mostly impacted during mid-winter warm events, when the long-wave enhancement by the canopy dominates over the shadowing from shortwave radiations.This is precisely the situation observed at Alptal, where mean DJF temperature is about 1.5 • C (1975-2005), and our improved canopy module with more realistic thermal heating from the canopy yields a better reproduction of the two mid-winter ablation events from 2004 to 2005 and 2006 to 2007.In that sense our work contributes an enhanced capability to model sub-canopy melt conditions that can jointly benefit water and forest management (Lawler and Link, 2011;Ellis et al., 2013).
In line with Rutter et al. (2009), Essery et al. (2008), and many others, we agree that the mass balance of snow is most relevant when it comes to snow hydrological applications or to the assessment of snow hazards.In the forest, snow interception and subsequent sublimation of intercepted snow majorly affects the sub-canopy mass balance; reductions up to 60 % for sub-canopy snowfall have been reported as a result of these processes (Hardy et al., 1997).Radiation, as the main driver of the melt, shapes the end-of-season snow mass balance, but the latter also critically depends on peak accumulation from the accumulation phase which our developments barely touch.We see our contribution as a necessary step in a sequential, multi-directional development and validation process, whereby the careful and independent validation of each component of the snow model will gradually improve its skills; SNOWPACK, now equipped with a more reliable and sophisticated radiative transfer scheme for the canopy, diagnosing flaws originating from other processes (mixed-precipitations, rain-on-snow events, or misrepresented canopy interception), should be easier.Promising work has just been published very recently as to new ways to parameterize canopy interception in alpine forests (Moeser et al., 2015).The proposed methodology should later serve the improvement of snow models like SNOW-PACK in aspects of crucial interest for the snow mass balance.
Other physical processes could further improve SNOW-PACK, in the present version the within-canopy air humidity is equal to the above-canopy one, while Durot (1999) show a 10 to 20 % increase in within-canopy air humidity in spring, with peaks during unloading.This should impact the modeled sub-canopy turbulent and radiative fluxes.Also, the reduction of albedo as a result of canopy debris should be considered, though tests performed at Alptal did not show much change upon a specific parameterization of sub-canopy snow aging.
Further wintertime assessment of model performance in colder, controlled environments, where mixed-precipitation events are scarce but radiation play an important role, would help confirm the added value of our new canopy formulation for the representation of sub-canopy snow dynamics.Data from the SnowMIP sites could be used for that purpose provided they are combined with site knowledge and expertise, and ancillary data that help fit important model parameters to the local canopy conditions (interception, radiation extinction).
SNOWPACK has a multi-layer and detailed representation of snow and soil, which features a highly resolved modeling of energy and mass balance in thin layers including, e.g., snow metamorphism and freezing point depressions during phase change in the soil (Wever et al., 2014).This detailed and physics-based description should have a corresponding representation of canopy processes, which has not been the case in earlier versions of SNOWPACK.The more detailed model described in this contribution is therefore a consistent extension of SNOWPACK and leads to an overall more balanced representation of processes in the air-canopy-snowsoil continuum.A physics-based, integrated modeling chain featuring such level of homogeneity and detail is rare.Sivapalan et al. (2003) and Rutter et al. (2009) underlined that such process-based model (rather than calibrated parametric models) offer the best possibility to address the current hydrological and ecosystemic challenges related to snow in a manner that ensures site transferability and robustness with respect to changing climate.The new version of SNOW-PACK with the two-layer canopy module builds a sound basis for such investigations.The current two-layer formulation of the canopy is also a suitable basis for a future model adaptation to deciduous forest environments.
Our two-layer canopy model exhibits robustness in two ways: -First, it shows little sensitivity to physical parameters that are hard to assess from standard forestry metrics or for non-investigated forests.The canopy heat mass is one of such parameters, as stated in Sect. 4. The other one is the fraction of LAI attributed to the top-most ("leafy") canopy layer, as illustrated in Fig. 7.The ratio of woody to total plant area is hard to measure op- tically, especially for evergreen canopies (Weiss et al., 2004).Pomeroy et al. (2009) used a formulation somewhat similar to ours to attribute LW radiation to emission from leafy or woody elements.They conclude that, depending on the forest structure and type, the needlebranch fraction as seen from a ground observer would range from 0.6 to 0.75 of the total plant elements.Our Alptal calibration attributing 50 % of canopy LAI to the upper "leafy" layer is consistent with this model-based estimate for leaves and branches.
-Second, the model exhibits a good performance at the Norunda site, while its free parameters (k LAI and f LAI ) have been calibrated in a different forest ecosystem and climatic context at Alptal.In both forests, coniferous species are dominant and it is suspected that extrapolation of our parameterizations to deciduous forests requires further adaptation.However, our results give confidence in the possibility of using our physics-based model without prior tuning in different alpine and subarctic catchments majorly covered by conifers.
Finally, it is a quite general finding that two-layer formulations of physical continuums often bring substantial improvements over single-layer ones.The step from big-leaf soilvegetation-atmosphere transfer models to dual-source models (e.g., Blyth et al., 1999;Bewley et al., 2010) is a typical illustration of this phenomenon for the computation of the land surface energy balance.Similarly, Dai et al. (2004) improved their modeling of forest CO 2 absorption by considering different regimes for sunlit and shaded leaves.Our results here are in line with this more general observation.

Conclusion
Our new canopy model demonstrates ability to simulate the difference in the thermal regimes of the canopy leafy and woody compartments, as assessed by comparison to observed canopy temperature and thermal radiation.This is achieved via the separation of the canopy in two layers of different heat masses, radiatively interacting with each other.In comparison, a one-layered version of the canopy module always yields poorer results despite optimization attempts.The most striking improvement is a reduction of the nighttime canopy cold bias, which can only be achieved via the two-layer formulation and results from the sheltering role of the upper canopy layer.The robustness of the new canopy model is confirmed by the successful evaluation of the model without prior tuning at a boreal coniferous site.Besides, the new formulation shows a weak sensitivity to biomass areal heat mass, a forest-dependent input parameter that can be hard to estimate locally.Model evaluation against snow water equivalent data indicate that the new parameterizations do not degrade the overall model skills while improving the representation of LW-enhanced ablation events.These are an important feature of the mid-winter and spring snow dynamics in temperate alpine regions.
The improved representation of the radiative components of the sub-canopy energy balance achieved here opens the path to the tracking, understanding, and modeling of further processes relevant for the underlying snowpack such as turbulent fluxes or heat advection by rain.In the end, enhanced models and process understanding should help obtain better hydrological simulation tools for crucial purposes like climate change impact assessment.

Figure 1 .
Figure 1.Radiative and turbulent fluxes in the two-layer canopy module.Ellipses feature radiation sources, dotted lines indicate radiation absorption within the layer with the indicated absorption factor; albedos at the border between layers are underlined.For turbulent fluxes, arrows denote aerodynamic resistance.

Figure 2 .
Figure 2. Typical event when snow on sensor is suspected.(a) Observed and modeled LW ↓BC .(b) Observed precipitation record.

Figure 3 .
Figure 3. LW ↓BC and SWE as represented by the different model versions over the calibration period; (a, b) subsets of daily cycles; (c) 24 h running means over the calibration period; (d) SWE.

Figure 5 .
Figure 5.Comparison between observed pine trunk temperature at 1.5 m height, 1 cm deep into the trunk, and modeled canopy temperatures: bulk canopy temperature for 1LnoHM and 1LHM, lowermost canopy-layer temperature for 2LHM.

Figure 6 .
Figure 6.Comparison between biomass (and biomass + air) storage fluxes inferred from observations (obs) and biomass fluxes modeled by the different SNOWPACK versions (model) at Norunda.

Figure 7 .
Figure 7. Sensitivity of model performance over 2003-2007 (with k LAI = 0.85) to f LAI .The MB and RMSE are for the variables SW ↓BC (SW in the legend) and LW ↓BC (LW in the legend).

Table 1 .
Parameters used by the SNOWPACK canopy module.

Table 2 .
Model versions and their calibration/optimization parameters.

Table 3
↓BC and SW ↓BC data from the snow season2003- 2004.

Table 3 .
Model performance after calibration and optimization over2003-2004.The calibration criterion (CC) is in bold.The * denotes versions where heat mass is optimized and not physically derived.

Table 5 .
Statistics of model evaluation at Norunda; "corr" is the correlation coefficient.The mean modeled and observed biomass (and biomass + air) heat fluxes are null over a period between two equal thermal states.