The Cryosphere, 18, 231–263, 2024
https://doi.org/10.5194/tc-18-231-2024
© Author(s) 2024. This work is distributed under
the Creative Commons Attribution 4.0 License.
Modeling snowpack dynamics and surface energy budget in boreal
and subarctic peatlands and forests
Jari-Pekka Nousu1,2,3, Matthieu Lafaysse2, Giulia Mazzotti2, Pertti Ala-aho1, Hannu Marttila1, Bertrand Cluzet4,
Mika Aurela5, Annalea Lohila5,6, Pasi Kolari6, Aaron Boone7, Mathieu Fructus2, and Samuli Launiainen3
1Water, Energy and Environmental Engineering Research Unit, University of Oulu, Oulu, Finland
2Univ. Grenoble Alpes, Université de Toulouse, Météo-France, CNRS, CNRM,
Centre d’Études de la Neige, Grenoble, France
3Bioeconomy and Environment, Natural Resources Institute Finland, Helsinki, Finland
4Snow and Atmosphere, WSL Institute for Snow and Avalanche Research SLF, Davos, Switzerland
5Climate System Research, Finnish Meteorological Institute, Helsinki, Finland
6Institute for Atmospheric and Earth System Research INAR, University of Helsinki, Helsinki, Finland
7CNRM, Université de Toulouse, Météo-France, CNRS, Toulouse, France
Correspondence: Jari-Pekka Nousu (jari-pekka.nousu@oulu.fi)
Received: 24 February 2023 – Discussion started: 8 March 2023
Revised: 25 October 2023 – Accepted: 22 November 2023 – Published: 12 January 2024
Abstract. The snowpack has a major influence on the land
surface energy budget. Accurate simulation of the snowpack
energy and radiation budget is challenging due to, e.g., ef-
fects of vegetation and topography, as well as limitations
in the theoretical understanding of turbulent transfer in the
stable boundary layer. Studies that evaluate snow, hydrol-
ogy and land surface models against detailed observations of
all surface energy balance components at high latitudes are
scarce. In this study, we compared different configurations
of the SURFEX land surface model against surface energy
flux, snow depth and soil temperature observations from four
eddy-covariance stations in Finland. The sites cover two dif-
ferent climate and snow conditions, representing the southern
and northern subarctic zones, as well as the contrasting for-
est and peatland ecosystems typical for the boreal landscape.
We tested different turbulent flux parameterizations imple-
mented in the Crocus snowpack model. In addition, we ex-
amined common alternative approaches to conceptualize soil
and vegetation, and we assessed their performance in simu-
lating surface energy fluxes, snow conditions and soil ther-
mal regime. Our results show that a stability correction func-
tion that increases the turbulent exchange under stable at-
mospheric conditions is imperative to simulate sensible heat
fluxes over the peatland snowpacks and that realistic peat
soil texture (soil organic content) parameterization greatly
improves the soil temperature simulations. For accurate sim-
ulations of surface energy fluxes, snow and soil conditions
in forests, an explicit vegetation representation is necessary.
Moreover, we demonstrate the high sensitivity of surface
fluxes to a poorly documented parameter involved in snow
cover fraction computation. Although we focused on models
within the SURFEX platform, the results have broader impli-
cations for choosing suitable turbulent flux parameterization
and model structures depending on the potential use cases for
high-latitude land surface modeling.
1 Introduction
The boreal zone, characterized by a mosaic of seasonally
snow-covered peatlands, forests and lakes, is the largest land
biome in the world. Snow conditions in the boreal zone are
rapidly changing due to climate warming, which is found to
be the strongest during the cold seasons in the Arctic (Ser-
reze et al., 2009; Screen and Simmonds, 2010; Boisvert and
Stroeve, 2015; Rantanen et al., 2022). Evidently snow has
an important role for water resources and human activities
in the cold regions, but it is also known that the snowpack
characteristics affect animal movement (Tyler, 2010; Peder-
sen et al., 2021) and plant distribution (Rasmus et al., 2011;
Published by Copernicus Publications on behalf of the European Geosciences Union.
232 J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget
Kreyling et al., 2012; Rissanen et al., 2021). Recent studies
show that especially cold-dwelling species have been shift-
ing towards higher latitudes and altitudes in search for more
suitable habitats (Tayleur et al., 2016; Couet et al., 2022).
Therefore, the rapid warming of the Arctic and its conse-
quences on the quantity and properties of snow may define
the destiny of many species and human activities in the bo-
real region. To predict future snow conditions, environmental
change, and the consequences for water resources, ecosys-
tems and people, predictive and process-based models pos-
sess great potential (Clark et al., 2015; Boone et al., 2017).
Land surface models (LSMs) have been used for decades in
numerical weather prediction and in global circulation mod-
els (Douville et al., 1995; Niu et al., 2011; Lawrence et al.,
2019) and have more recently become common tools for in-
terdisciplinary impact studies (Blyth et al., 2021).
Snowpack has a major impact on the wintertime energy
budget due to its influence on the land surface albedo (here-
after albedo) and the surface heat fluxes (Cohen and Rind,
1991; Eugster et al., 2000). The heat diffusion within the
snowpack is determined by the surface heat fluxes, internal
properties of the snowpack and soil thermal regime. Cor-
rectly representing the snowpack is thus essential for sim-
ulating energy and mass exchange between the snow surface
and the atmosphere, as well as below the snowpack (e.g.,
surface temperatures and soil freezing–thawing dynamics,
Koivusalo and Heikinheimo, 1999; Slater et al., 2001). The
snowpack energy budget is partitioned into downward and
upward shortwave and longwave radiation, turbulent fluxes
of sensible and latent heat, snowpack–ground heat flux, and
phase changes in the snow. The snowpack energy balance
and energy partitioning among the flux components vary
strongly across diurnal and seasonal timescales, as well as
between different ecosystems (Clark et al., 2011; Stiegler
et al., 2016; Stigter et al., 2021). It is essential that LSMs
are able to correctly reproduce this variability.
On the vast boreal and arctic peatlands with shallow veg-
etation, the snow cover can exclusively determine the win-
tertime albedo (Aurela et al., 2015). With minimal solar ra-
diation during winter months on these open snow fields, tur-
bulent fluxes make an important component in the energy
budget of the snowpack, as they compensate for the radiative
cooling processes and further contribute to snowmelt (Lack-
ner et al., 2022; Conway et al., 2018). Simulation of tur-
bulent fluxes under stable atmospheric conditions is known
as one of the major sources of uncertainty in snow mod-
els (Lafaysse et al., 2017; Menard et al., 2021). In LSMs
the turbulent fluxes are commonly computed with bulk aero-
dynamic approaches, where sensible and latent heat fluxes
are proportional to the turbulent exchange coefficient ac-
cording to the Monin–Obukhov similarity theory. These ap-
proaches typically use atmospheric stability correction func-
tions based either on the bulk Richardson number (Martin
and Lejeune, 1998; Lafaysse et al., 2017; Clark et al., 2015)
or the Obukhov length scale (Jordan et al., 1999). It is estab-
lished that the Monin–Obukhov similarity theory does not
well represent low-wind and stable atmospheric conditions
above aerodynamically smooth surfaces such as snow (Con-
way et al., 2018). In such conditions, the simulated surface
temperatures have been found to be unrealistically low, as
the turbulent boundary layer tends to decouple from the snow
surface (Derbyshire, 1999; Andreas et al., 2010). To circum-
vent this effect, stability correction functions have been mod-
ified to permit turbulent fluxes above critical stability thresh-
olds (Lafaysse et al., 2017) by manipulating the wind speed
(Martin and Lejeune, 1998; Andreas et al., 2010) or includ-
ing a windless turbulent exchange coefficient (Jordan et al.,
1999). Evaluations of these modifications often rely on vali-
dation with observed surface temperatures and snow depths
(e.g., for the detailed snowpack model Crocus; Martin and
Lejeune, 1998; Lafaysse et al., 2017), while comparisons
against turbulent energy flux data remain scarce (Lapo et al.,
2019; Conway et al., 2018).
The energy budgets of forest canopies and below-canopy
snowpack are different to those on open peatlands, as tur-
bulent exchange is attenuated by the canopy, and the snow-
pack energy budget and snowmelt are mostly driven by the
radiation balance (Rutter et al., 2009; Essery et al., 2009;
Varhola et al., 2010). However, due to heterogeneous canopy
structures and canopy processes (radiation transmittance,
snow interception and unloading) together with low solar an-
gles, the albedo dynamics of seasonally snow-covered bo-
real forests is complex (Malle et al., 2021). The absorp-
tion of the shortwave radiation can be highly heterogeneous,
having direct implications on canopy temperatures (Webster
et al., 2017) and on the resulting longwave radiative fluxes
between canopy, snowpack and the atmosphere (Mazzotti
et al., 2020b). Forest snow modeling has been identified as
a priority in advancing cold region climate and hydrological
models (Rutter et al., 2009; Krinner et al., 2018; Lundquist
et al., 2021). Various models that have been proposed to rep-
resent the large-scale impact of forest on the snowpack en-
ergy budget (Niu et al., 2011; Lawrence et al., 2019; Boone
et al., 2017) are still prone to large errors, due to the com-
plexity and unresolved spatial scales of the underlying phys-
ical processes (Loranty et al., 2014; Thackeray et al., 2019).
The forest snow model evaluations against concurrent snow-
pack and surface energy balance data are also surprisingly
scarce (e.g., Tribbeck et al., 2006). For instance, the explicit
forest scheme of SURFEX LSM, MEB (Boone et al., 2017),
has so far been evaluated only against data from three neigh-
bor sites in Saskatchewan, Canada (Napoly et al., 2020). This
considerably limits knowledge of the model skill to represent
snow–forest interactions in regional or global applications.
The texture and thermal properties of the underlying soil
can strongly impact the snowpack–ground heat exchange,
snowpack energy fluxes and snowpack dynamics (Decharme
et al., 2016). Peatlands have high soil organic content (SOC)
and are characterized by high porosity, shallow water table,
weak hydraulic suction, strong gradient in hydraulic conduc-
The Cryosphere, 18, 231–263, 2024 https://doi.org/10.5194/tc-18-231-2024
J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget 233
tivity from high values at the top to low values at the sub-
surface, low thermal conductivity and large heat capacity
(Decharme et al., 2016; Marttila et al., 2021; Morris et al.,
2022; Menberu et al., 2021). These properties result in a wet
soil profile resistant to temperature variations, while the drier
top peat and moss layer can also provide effective insulation
particularly during summertime (Beringer et al., 2001; Park
et al., 2018; Chadburn et al., 2015). The importance of the
soil texture is still often overlooked even in detailed snow
models. For instance, in model comparisons of the ESM-
SnowMIP project (Ménard et al., 2019), no SOC information
was used to parameterize the participating LSMs to the refer-
ence sites. In addition, many spatial snow simulations neglect
peat soils or SOC altogether, and their hydrological and ther-
mal characteristics are derived from fractions of sand, silt and
clay (Vernay et al., 2022; Brun et al., 2013; Mazzotti et al.,
2021; Richter et al., 2021)
The goal of this study is to evaluate the ability of SUR-
FEX LSM (Masson et al., 2013) to describe the surface en-
ergy balance and its drivers in boreal and subarctic peat-
lands and forests. We evaluate the effect of alternative tur-
bulent exchange and snowpack parameterizations and exam-
ine the skills of alternative model configurations to represent
the soil–snow–vegetation interactions. The modeling frame-
work includes flexible parameterizations for different pro-
cesses within the Crocus snowpack model (Vionnet et al.,
2012), and its coupling to ISBA (Noilhan and Mahfouf,
1996; Decharme et al., 2016) and MEB (Boone et al., 2017;
Napoly et al., 2017) models enables assessments of soil–
snow–vegetation interactions. We compare the model sim-
ulations against observed surface energy fluxes, snow depth,
and soil temperatures from two forest and two peatland sites
in Finland. We focus on the snow cover period but cover also
the snow-free season for reference. At the peatland sites, we
test the sensitivity of the surface heat fluxes to different tur-
bulence and snow parameterizations and assess how sensitive
soil temperature and snowpack dynamics are to SOC. At the
forest sites, we compare the simulations of the ISBA com-
posite soil–vegetation and MEB big-leaf forest scheme to as-
sess the suitability of different forest snow model structures.
2 Materials and methods
2.1 Study sites
We consider coniferous forest and peatland ecosystems in
southern and northern Finland. Both areas are located in the
boreal biome and have seasonal snow cover (Fig. 1, Table 1).
Site photos can be found in the Supplement (Fig. S8).
2.1.1 Pallas supersite (N-WET and N-FOR)
The Pallas area represents northern subarctic conditions and
is characterized by pine and spruce forests, wetlands, fells,
and lakes (Aurela et al., 2015; Lohila et al., 2015; Marttila
et al., 2021). In this study, we use data from its two eddy-
covariance (EC) flux stations.
Lompolojänkkä (northern peatland, N-WET) is a pristine
northern boreal mesotrophic sedge fen where the wetter parts
are dominated by sedges (Carex rostrata (most abundant),
Carex chordorrhiza, Carex magellanica and Carex lasio-
carpa) and the drier parts consist of shallow deciduous trees
(Betula nana and Salix lapponum). Moreover, the fen has a
fairly low coverage of shrubs, mainly Andromeda polyfolio
and Vaccinium oxycoccos. The vegetation height is shallow
(∼ 0.4 m), with the exception of isolated trees/bushes on the
drier edges of the peatland.
Kenttärova (northern forest, N-FOR) is a northern boreal
spruce forest, located on a hill-top plateau with mineral soil
approximately 60 m above Lompolojänkkä wetland. The for-
est is dominated by Norway spruce (Picea abies) with some
deciduous trees, mainly birch (Betula pubescens) but also as-
pen (Populus tremula) and pussy willow (Salix caprea). Ac-
cording to the classification by Brunet (2020), Kenttärova is
a sparse forest. Both sites and their measurements have been
described in detail by Aurela et al. (2015).
2.1.2 Hyytiälä and Siikaneva (S-WET and S-FOR)
The sites are located in southern subarctic conditions in the
Pirkanmaa region in southern Finland, at about 5 km dis-
tance from each other. Siikaneva fen (southern wetland, S-
WET) is a southern boreal oligotrophic fen dominated by
sedges (Eriophorum vaginatum, Carex rostrata and Carex
limos) and has an extensive Sphagnum cover (mainly Sphag-
num balticum, Sphagnum majus and Sphagnum papillosum).
The site has been described in detail in Aurela et al. (2007),
Alekseychik et al. (2017) and Rinne et al. (2018).
Hyytiälä (southern forest, S-FOR) is a managed boreal
Scots pine (Pinus sylvestris) dominated forest on mineral
soil, described in detail by Launiainen (2010); Launiainen
et al. (2022). According to the classification by Brunet
(2020), the site is a dense forest.
2.2 Models
We use components from the SURFEX LSM (Surface Ex-
ternalisée, Masson et al., 2013) platform. SURFEX was se-
lected as its modularity and vast range of model structures
and incorporated process parameterizations enable its use in
diverse applications. Specifically, we used ISBA (Noilhan
and Planton, 1989; Noilhan and Mahfouf, 1996) for compos-
ite soil–vegetation, MEB (Boone et al., 2017; Napoly et al.,
2017) for the canopy and Crocus (Vionnet et al., 2012), and
its ensemble/multiphysics version ESCROC (Lafaysse et al.,
2017) for the snowpack simulations. Specifically, ISBA cou-
pled to Crocus is used for both peatland and forest experi-
ments (Fig. 2a, b), whereas MEB coupled to Crocus is only
used for the forest experiments (Fig. 2c). In the next subsec-
tions, we briefly describe these model components. Param-
https://doi.org/10.5194/tc-18-231-2024 The Cryosphere, 18, 231–263, 2024
234 J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget
Figure 1. (a) Study area locations inside the boreal land biome (green area, Olson et al., 2001), (b) study site locations in Finland (Esri,
2023) and (c–f) aerial images of each site (NLSF, 2020).
Table 1. General site information.
Site Code Coordinates Ecosystem Soil type
Lompolojänkkä N-WET 67◦59.835′ N, 24◦12.546′ E mesotrophic fen peat
Siikaneva S-WET 61◦49.961′ N, 24◦11.567′ E oligotrophic fen peat
Kenttärova N-FOR 67◦59.237′ N, 24◦14.579′ E sparse spruce forest podzol
Hyytiälä S-FOR 61◦50.471′ N 24◦17.439′ E dense pine forest podzol
eterizations and different configurations of ISBA, MEB and
Crocus models are detailed in Sect. 2.3.
2.2.1 ISBA
ISBA (Interactions between the Soil, Biosphere and At-
mosphere) is the soil and vegetation component of SUR-
FEX (Noilhan and Planton, 1989; Noilhan and Mahfouf,
1996). It simulates the mass and energy fluxes in the soil–
vegetation composite, as well as the exchanges between
the soil–vegetation and the overlying atmosphere/snowpack
(Fig. 2a, b). ISBA is used for the general circulation mod-
els by Météo-France (Mahfouf et al., 1995; Douville et al.,
1995; Salas-Mélia et al., 2005; Voldoire et al., 2013, 2019)
and for numerical weather prediction in numerous countries
(e.g., Hamdi et al., 2014; Bengtsson et al., 2017).
In ISBA, the surface heat flux between the atmosphere
and the soil–vegetation composite (G0, Wm−2) is computed
as the residual of the sum of all surface/atmosphere energy
fluxes:
G0 = SWD(1−LSA)+ (LWD− σT 4s )+H +LE, (1)
where SWD (Wm−2) and LWD (Wm−2) are the incoming
shortwave and longwave radiations, respectively. The land
surface albedo is denoted as LSA,  is the surface emissivity,
The Cryosphere, 18, 231–263, 2024 https://doi.org/10.5194/tc-18-231-2024
J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget 235
σ is the Stefan–Boltzmann constant and Ts (K) is the sur-
face temperature. The sum of the radiation terms is hereafter
denoted as Rn (Wm−2).
The sensible heat flux H (Wm−2) is computed with the
bulk aerodynamics approach:
H = ρacρCHVa(Ts− Ta), (2)
where the air density, the specific heat capacity, the wind
speed and the air temperature are denoted with ρa (kgm−3),
cp (Jkg−1K−1), Va (ms−1) and Ta (K), respectively. CH is the
turbulent exchange coefficient described later and is one of
the parameters that is the focus of this study. When the soil is
not fully covered by snow, the latent heat flux LE (Wm−2) is
the sum of evaporation from the bare soil surface, Eg; evap-
oration of intercepted water on the canopy, Ec; transpiration
from the vegetation, Etr; and sublimation from bare soil ice,
Si:
LE= Lv(Eg+Ec+Etr)+ (Lf+Lv)(Si), (3)
where Lv (Jkg−1) and Lf (Jkg−1) are the latent heat of va-
porization and fusion, respectively. Total evapotranspiration
(ET) is computed as
ET= Eg+Ec+Etr = (1− veg)ρaCHVa
[huqsat(Ts)− qa] + vegρaCHVahv, (4)
where veg is the fraction of vegetation cover, qsat(Ts)
(kgkg−1) is the saturated specific humidity at the surface,
qa(Ts) (kgkg−1) is the atmospheric specific humidity, hu is
the dimensionless relative humidity at the ground surface re-
lated to the superficial soil moisture content, and hv is the
dimensionless Halstead coefficient describing the Ec and Etr
partitioning between the leaves covered and not covered by
intercepted water (see Noilhan and Mahfouf (1996) for de-
tails)
The turbulent exchange coefficient CH is based on the for-
mulation of Louis (1979):
CH =
[
k2
ln(zu/z0t ) ln(za/z0t )
]
f (Ri), (5)
where zu (m) is the reference height of Va, za (m) is the ref-
erence height of Ta and humidity, z0t (m) is the roughness
height for heat, k (–) is the von Kármán constant, and f (Ri)
(–) describes the decrease in CH as a function of increasing
atmospheric stability, represented through Richardson num-
ber (Ri) (Louis, 1979).
Instead of separate treatment of the vegetation canopy and
ground, ISBA considers the composite soil–vegetation en-
ergy budget (Fig. 2a, b). In the most detailed soil scheme
ISBA diffusion (ISBA-DIF, Boone et al., 2000; Decharme
et al., 2011), used in this study, the 1D Fourier law is used
to solve the soil heat diffusion, while a mixed-form Richards
equation is applied for the 1D soil water movements. Simi-
lar to Napoly et al. (2020), we use the A− gs stomatal con-
ductance formulation derived from the coupling of photo-
synthetic CO2 demand and stomatal function (Calvet et al.,
1998).
ISBA uses parameters such as one-sided leaf area index
(LAI, m2 m−2), vegetation height, vegetation thermal iner-
tia, albedos of soil and vegetation, fractions of sand and
clay, and SOC content to characterize the composite soil–
vegetation column. These parameters may be defined by the
user or obtained from global or regional databases (e.g.,
Faroux et al., 2013) and pedotransfer functions (Noilhan and
Mahfouf, 1996; Peters-Lidard et al., 1998). In the presence
of full snow cover, the surface energy budget is solved by
Crocus (Sect. 2.2.3). For partial snow cover, Crocus is used
to solve the snow-covered fraction while the energy balance
of the snow-free fraction is computed by ISBA, and total sur-
face energy fluxes are computed as weighted averages of the
snow and snow-free fractions (Sect. 2.3.1).
2.2.2 MEB
MEB (multi-energy balance) is a recent ISBA development
to explicitly describe vegetation and soil energy and mass
balances. It was developed initially for forests (Boone et al.,
2017; Napoly et al., 2017) and found to yield improved snow
and soil temperature simulations (Napoly et al., 2020) but
has not been evaluated for boreal and subarctic conditions.
MEB simulates surface energy budget separately for soil and
vegetation canopy. When the ground is snow covered, the
energy budget of the snowpack is also explicitly represented.
We used the MEB option, where the forest floor is covered
by a litter layer instead of the bare soil surface (Napoly et al.,
2020) (Fig. 2c). MEB describes vegetation canopy as a single
big leaf (Boone et al., 2017). The respective energy balance
equations for the canopy, the snowpack and the ground sur-
face/litter layer in MEB are
Cv
∂Tv
∂t
= Rnv−Hv −LEv +Lfφv
Cg,1
∂Tg,1
∂t
= (1− ρsng)(Rng−Hg −LEg)
+ρsng(Ggn+ τn,NnSWnn)−Gg,1+Lfφg,1
Cn,1
∂Tn,1
∂t
= Rnn−Hn−LEn− τn,1SWnn
+n,1−Gn,1+Lfφn,1
, (6)
where Cv, Cg,1, and Cn,1 (Jm−2 K−1) and Tv, Tg,1, and
Tn,1 K are the effective heat capacities and temperatures of
the canopy, ground surface/litter layer, and snowpack, re-
spectively. In these equations, the subscripts g,1 and n,1
represent the uppermost layer for the soil and the snowpack,
respectively. Gg,1 and Gn,1 are, respectively, the conduc-
tion heat flux at the bottom of the uppermost soil or snow
layer. Ggn is the conduction heat flux at the soil–snow in-
terface. Rnv, Rng and Rnn (Wm−2) are net radiation, i.e.,
the sum of net shortwave radiation and net longwave radi-
ation from/to the corresponding layer. The shortwave radia-
https://doi.org/10.5194/tc-18-231-2024 The Cryosphere, 18, 231–263, 2024
236 J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget
tion scheme used in MEB is described in Carrer et al. (2013).
Light transmission through the canopy is computed with a
sky view factor, which depends on LAI, solar angle and a
vegetation-dependent constant (see Eq. 45 in Boone et al.,
2017). H and LE flux parameterization as well as the stabil-
ity correction functions are detailed in Boone et al. (2017).
Obviously Tv, Tg,1 and Tn,1 are also involved in the radiative
and turbulent terms, providing a linear system of equations
to be solved by an implicit numerical scheme. In this study,
MEB is coupled to the ISBA-DIF soil scheme (Sect. 2.2.1)
and the snowpack model Crocus (Sect. 2.2.3). Energy fluxes
between the canopy and the ground surfaces are calculated
within MEB and prescribed as upper-boundary conditions in
the subsequent Crocus and ISBA-DIF calculations.
2.2.3 Crocus
Crocus is a 1D physically based multilayer snowpack model
(Vionnet et al., 2012). It is the most detailed snow scheme
in ISBA and has been used for operational avalanche haz-
ard forecasting in the French mountain ranges for the past
3 decades (Morin et al., 2020). It aims to mimic the verti-
cal layering of snowpacks with a Lagrangian discretization
system, avoiding the aggregation of snow layers with highly
different physical properties. A detailed description of Cro-
cus and its integration in SURFEX can be found in Vionnet
et al. (2012).
The snowpack surface energy budget is the sum of net ra-
diation, turbulent fluxes and advective fluxes from precipita-
tion. Over the snow, the sensible heat flux is computed simi-
larly as in ISBA (Eq. 2) for soil surface, while the latent heat
flux (sublimation/deposition), LEs, is computed as
LEs = (Lf +Lv)ρaCHVa[qsat(Ts−qa)], (7)
where Ts (K) is the snow surface temperature. The bottom of
the snowpack and the uppermost soil layer of ISBA are fully
coupled with a mass- and energy-conserving semi-implicit
solution. The semi-implicit solution refers to a coupled sys-
tem in which both components are solved separately with an
implicit approach considering that the state of the second sys-
tem remains constant during the solving of the first system.
The heat conduction flux Ggn at the snow–soil interface is
explicitly computed using the Fourier equation and depends
on the temperature gradient between the bottom snow layer
and the uppermost soil layer (Eq. 4 in Decharme et al., 2011).
The soil thermal conductivity and heat capacity are described
using pedotransfer functions of ISBA (Noilhan and Mahfouf,
1996; Peters-Lidard et al., 1998).
2.3 Model configurations and parameterization
2.3.1 Model configurations
We use three different configurations of ISBA, MEB and
Crocus modules (Fig. 2). The first two configurations use
the composite soil–vegetation conceptualization of ISBA
(Fig. 2a, b) and differ only in how the snow cover fraction is
represented. ISBA aggregates the properties of soil and veg-
etation depending on vegetation fraction (veg) that covers a
given grid cell.
The first configuration (referred to as ISBA-FS) assumes
the snowpack fully covers the soil–vegetation composite re-
gardless of the snow depth (Fig. 2a). It is the common ap-
proach for snow simulations over shallow vegetation or bare
soil (Vernay et al., 2022; Nousu et al., 2019), for large-scale
reanalyses (Brun et al., 2013), and for some hydrological
applications (Lafaysse et al., 2011; Revuelto et al., 2018).
In addition, most site-level evaluations of SURFEX snow
schemes rely on ISBA-FS configuration (Decharme et al.,
2016; Lafaysse et al., 2017).
In the second configuration (referred to as ISBA-VS), a
part of the soil–vegetation composite is covered by snow
while the remaining (non-snow) fraction stays in constant
contact with the atmosphere. This proportion is governed
by snow depth. The effective snow cover fraction is the
weighted average between the snow fraction of vegetation
(psnv) and snow fraction of the bare ground (psng), calculated
as (Decharme et al., 2019; Napoly et al., 2020)
psnv =min
(
1.0,
HS
HS+wswz0
)
, (8)
psng =min
(
1.0,
HS
HSg
)
, (9)
where HS (m) is the snow depth, HSg is the threshold value
for snow depth (0.01 m by default) and z0 (m) denotes the
surface roughness. The coefficient wsw is supposed to relate
to scale-dependent vegetation characteristics and is set to 5
by default in SURFEX and in numerical weather prediction
configurations (as well as in this study). However, without
clear consistency, highly different values of wsw have been
used, e.g., in climate simulations (wsw = 2 by Decharme
et al., 2019) and hydrological applications (wsw = 0.2 by
Le Moigne et al., 2020). We present a summary of the
application-specific treatment of the snow cover fraction and
wsw in the Appendix (Table C1). This summary shows that
selection of the wsw value seems arbitrary and the fractional
concept is only loosely linked to any physical relationships
between soil, vegetation and snow. Yet, it is necessary for
such a composite approach.
The third configuration (later denoted as MEB) is the
big-leaf approach where the fluxes between the canopy and
snowpack/ground are explicitly computed by MEB and pre-
scribed in the subsequent Crocus snowpack and ISBA-DIF
soil modules (Fig. 2c).
The Cryosphere, 18, 231–263, 2024 https://doi.org/10.5194/tc-18-231-2024
J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget 237
Figure 2. Three different model configurations used in the study: (a) ISBA-FS with full snow cover fraction, (b) ISBA-VS with varying
snow cover fraction and (c) MEB big-leaf approach. Considered energy fluxes between domains are represented with arrows.
2.3.2 ESCROC parameterizations for snow processes
and turbulent exchange
We use the multiphysics version of Crocus (ESCROC, En-
semble System Crocus, Lafaysse et al., 2017) to evaluate
the impact and associated uncertainties of the different pa-
rameterizations of snow processes and turbulent exchange.
In ESCROC, the main physical processes and properties
of snowpack, as well as the turbulent fluxes, can be repre-
sented by several alternative options. These include density
of new snow, snow metamorphism, absorption of solar ra-
diation, turbulent fluxes, thermal conductivity, liquid water
holding capacity, snow compaction and surface heat capacity
(Eqs. 1–17 in Lafaysse et al., 2017). Lafaysse et al. (2017)
showed that an optimized standard subensemble of 35 mem-
bers (E2 subensemble) is sufficient to provide a spread of
the appropriate magnitude compared to model errors. We use
this subensemble (hereafter ESCROC-E2) similar to recent
studies quantifying the model uncertainty (e.g., Deschamps-
Berger et al., 2022; Tuzet et al., 2020). The presented ensem-
ble spread corresponds to simulated values between ensem-
ble minimum and maximum.
In Crocus, the default turbulent exchange parameterization
(Eq. 5) has been found to underestimate the turbulent fluxes
under stable conditions (Martin and Lejeune, 1998). There-
fore, different stability dependencies of the CH have been
implemented in ESCROC-E2. They differ mainly in the Ri
thresholds below which CH is assigned a constant value to
enable turbulent heat and mass transport under stable con-
ditions. As shown by Fig. 4 in Lafaysse et al. (2017), these
parameterizations are the (a) classical Louis (1979) formula
(later referred to as RIL) with threshold at Ri = 0.2; (b) RIL
with threshold at Ri = 0.1 (RI1); (c) RIL with threshold at
Ri = 0.026 (RI2); and (d) modified formulation with effec-
tive roughness length for heat (10−3 m), with minimum wind
speed (0.3 ms−1), and with threshold at Ri = 0.026 (M98)
by Martin and Lejeune (1998). The RIL parameterization is
widely used in SURFEX applications (e.g., Decharme et al.,
2019; Le Moigne et al., 2020), RI2 is applied in operational
snow modeling in the Alpine area (Vernay et al., 2022) and
M98 was recently used in the Canadian Arctic by Lackner
et al. (2022). However, evaluations of the different Crocus
turbulent flux parameterizations against surface flux data are
still lacking. MEB uses a different stability correction term
(Boone et al., 2017) and applies only the RIL option for the
stable conditions.
2.3.3 Site parameters
The parameterization of ISBA and MEB for the study sites
is given in Table 2. Summer LAI and vegetation height were
obtained from literature, while winter LAI (and monthly LAI
cycle) was estimated according to the proportion of decidu-
ous and coniferous vegetation on each site. The LAI of S-
FOR refers to conditions before forest thinning in early 2020.
The thinning, resulting in ca. 35 % reduction in LAI, was ne-
glected in our simulations as major part of the simulation
period covers time before the thinning. Vegetation types in
ISBA are characterized according to ECOCLIMAP (Cham-
peaux et al., 2005); the forest sites in this study classify as
boreal needleleaf evergreen, while the peatland sites are best
represented as boreal grass. Additional parameters based on
LAI, vegetation height and vegetation type are computed fol-
lowing the standard methods of ISBA (Noilhan and Mahfouf,
1996; Carrer et al., 2013).
Soil texture (sand and clay fractions) for the forest sites is
based on in situ measurements. The peat soils at S-WET and
N-WET were parameterized as fully organic for the upper-
most 1 m, in accordance with field measurements (Väliranta
and Mathijssen, 2021; Muhic et al., 2023), while the deeper
layers were assigned as mineral soil similar to the contigu-
ous forests. Although peat profiles may be deeper, the soils
https://doi.org/10.5194/tc-18-231-2024 The Cryosphere, 18, 231–263, 2024
238 J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget
below the damping depth of annual temperature fluctuations
(ca. 1.1 m for saturated peat soil with porosity ca. 90 %) are
assumed not to have significant impact on surface energy flux
dynamics. The SOC values for mineral soils of N-FOR and
S-FOR were taken from Lindroos et al. (2022). The rest of
the parameters presented in Table 2 were assigned as esti-
mates. The thermal and water retention parameters are sub-
sequently derived from the pedotransfer functions of ISBA
(Noilhan and Mahfouf, 1996; Peters-Lidard et al., 1998).
2.4 Data
2.4.1 Model forcing
Meteorological forcing consist of hourly observations of
air temperature, wind speed, precipitation rate, air humid-
ity, downward shortwave and longwave radiation, and atmo-
spheric pressure. The available meteorological observations
from the nearest meteorological stations were obtained from
Finnish Meteorological Institute (FMI) open database (FMI,
2021) (station IDs: N-WET 778135, N-FOR 101317, S-FOR
101987). Meteorological observations at the S-WET site
come from the SMEAR database (Alekseychik et al., 2022a).
At S-WET and S-FOR the shortwave and longwave radiation
were obtained from the SMEAR database, while at N-WET
and N-FOR data from FMI stations were used. The diffuse-
to-total-shortwave-radiation ratio, r , was estimated as a func-
tion of the cosine of the sun zenith angle, µ. More specif-
ically, a third-degree polynomial fit between r and µ was
obtained using the atmospheric model SBDART (Ricchiazzi
et al., 1998) to simulate diffuse and total solar radiation in
clear-sky conditions. The atmospheric profile was set to typ-
ical winter conditions, 0.09 for the aerosol optical thickness,
300 DU for the ozone column and 0.854 gcm−2 for the water
vapor column.
The data gaps in meteorological observations were first
filled by the contiguous sites (e.g., N-FOR for N-WET and
vice versa) and the remaining gaps by other nearby meteo-
rological stations (IDs: N-WET/N-FOR 101932, S-WET/S-
FOR 101520). The missing radiation observations were first
filled by the contiguous sites, and the remaining gaps by
ERA5 reanalysis data (Hersbach et al., 2020). Only less
than 10 h of ERA5 data was used for N-WET, N-FOR and
S-WET. However, S-FOR radiation observations contained
more gaps, specifically LWD in 2008–2012. A good agree-
ment between site observations and ERA5 estimates of LWD
is shown in the Appendix (Fig. B1). Furthermore, the fraction
of snow to total precipitation Pice/(Pice+Pliq) is assumed to
linearly decrease from 1 to 0 at air temperatures between 0
and 1 ◦C.
2.4.2 Model evaluation data
We use surface energy flux observations, snow depth and
soil temperatures in model evaluation. The availability period
of each variable is given in Table B1. On all sites, upward
shortwave radiation (SWU) and upward longwave radiation
(LWU) were measured using pyranometers and pyrgeome-
ters, while ground heat flux (G) was measured using soil heat
flux plates between 5 and 10 cm depths.
The sensible heat (H ) and latent heat (LE) were mea-
sured by the eddy-covariance (EC) technique. The EC sys-
tems consist of USA-1 (METEK) three-axis and Gill HS-50
sonic anemometers as well as closed-path LI-7000 and LI-
7200 (LI-COR, Inc.) CO2/H2O analyzers. The detailed de-
scriptions of the instrumentation, footprint analysis and the
procedures for obtaining the turbulent heat fluxes from raw
eddy-covariance data are detailed in the original data and site
publications (N-WET/N-FOR, Aurela et al., 2015; and S-
WET/S-FOR, Mammarella et al., 2016, 2019; Alekseychik
et al., 2022b). In short, H and LE were screened for instru-
ment failure and data outliers, and data were quality flagged
according to friction velocity (u∗) and flux stationarity (FST)
criteria (Foken et al., 2005) as follows
– flag 2: all data (after screening of instrument failures
and outliers);
– flag 1: u∗ ≥ 0.1 ms−1 and 0.3≤ FST ≤ 1.0;
– flag 0: u∗ ≥ 0.1 ms−1 and FST ≤ 0.3.
At S-WET, S-FOR and N-FOR, automated snow depth
observations are directly used (FMI, 2021). On N-WET the
automated snow depth measurement is at 0.7 m height, and
therefore the exceeding snow depths were taken from bi-
weekly manual measurements. To account for the spatial
variability of snow depth in the forests, manual snow depth
measurements from a snow course in close proximity of
the automated measurements were used (Aalto et al., 2022;
Marttila et al., 2021). Each site has different configuration
of soil temperature sensors. At N-FOR and N-WET stations,
soil temperatures are measured at 5 and 20 cm depths (Au-
rela et al., 2015). Soil temperatures at S-FOR and S-WET
are measured at depths of 0, 5, 10, 30, 50 and 75 cm (Aalto
et al., 2022).
2.5 Model experiments
On the peatland sites, we evaluate the skill of ISBA-FS
(Sect. 2.3.1) and effect of ESCROC-E2 parameterizations
(Sect. 2.3.2) on surface heat fluxes over snowpack and bare
ground. The simulations are further used to assess the dif-
ferences in snow depth simulations between ESCROC-E2
turbulent exchange options. For a more detailed evaluation
of two contrasting turbulent exchange options, we conducted
deterministic ISBA-FS simulations with site parameters and
default ESCROC-E2 snow parameterizations but different
treatments of turbulent exchange, RIL and M98 options (see
Sect. 2.3.2; referred to as RIL-SOC and M98-SOC). More-
over, the influence of soil texture was explored by comparing
M98-SOC to the simulation where soil was parameterized as
The Cryosphere, 18, 231–263, 2024 https://doi.org/10.5194/tc-18-231-2024
J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget 239
Table 2. Main model parameters for the study sites. Vegetation types BOGR and BONE correspond to boreal grass and boreal needleleaf
evergreen, respectively.
Parameter N-WET S-WET N-FOR S-FOR Source
Vegetation type BOGR BOGR BONE BONE ECOCLIMAP: Champeaux et al. (2005)
Vegetation fraction (only 0.95 0.95 0.95 0.95 ECOCLIMAP: Champeaux et al. (2005)
with ISBA-VS) (–)
Vegetation height (m) 0.4 0.25 13 15 Aurela et al. (2015); Alekseychik et al. (2017);
Kolari et al. (2022)
LAImax (m2 m−2) 1.3 0.6 2.1 3.0 Aurela et al. (2015); Alekseychik et al. (2017);
Kolari et al. (2022)
LAImin (m2 m−2) 0.3 0.1 1.9 2.4 assigned
Vegetation albedo (NIR/VIS) (–) 0.136 0.187 0.145 0.145 assigned
Soil albedo (NIR/VIS) (–) 0.136 0.187 0.145 0.145 assigned
Ta measurement height (m) 2 2 2 2 FMI (2021)
Wind measurement height (m) 13 3 23 16.8 Aurela et al. (2015); Mammarella et al. (2019);
Alekseychik et al. (2022a)
Elevation (m) 270 162 347 181 Hari et al. (2013); Alekseychik et al. (2022a);
FMI (2021)
Clay (%) (below 1 m at peatlands) 9 7 9 7 measurements
Sand (%) (below 1 m at peatlands) 76 65 76 65 measurements
SOC (0–30 cm) (kgm−2) 93.5 93.5 3.0 3.5 Lindroos et al. (2022);
Muhic et al. (2023);
Väliranta and Mathijssen (2021)
SOC (30–70 cm) (kgm−2) 93.5 93.5 1.75 0.75 Lindroos et al. (2022);
Muhic et al. (2023);
Väliranta and Mathijssen (2021)
SOC (70–100 cm) (kgm−2) 93.5 93.5 0 0 Väliranta and Mathijssen (2021);
Muhic et al. (2023)
Start of simulation (yyyy–mm) 2013–09 2016–09 2013–09 2008–09 –
End of simulation (yyyy–mm) 2021–07 2021–07 2021–07 2021–07 –
mineral soil, similar to the contiguous forest site (referred to
as M98-MIN).
On the forest sites, we examine the skills of the differ-
ent alternatives to represent the energy and mass budgets
of soil and vegetation (ISBA-VS, ISBA-FS, and MEB in
Sect. 2.3.1), as well as their implications on snow depth,
soil temperature and surface energy fluxes. First, we com-
pare ESCROC-E2 simulations with these three configura-
tions focusing on the snow depth and soil temperature. These
ISBA-VS simulations are conducted with the default snow
cover fraction parameterization (wsw = 5 in Eq. 8). Addi-
tional ISBA-VS simulations are performed to assess the sen-
sitivity of the wsw parameter, using a value of 0.2 in Eq. (8).
For a more detailed comparison of the simulated and ob-
served above-canopy surface energy fluxes, we conduct de-
terministic simulations with both ISBA-VS and MEB, the de-
fault snow cover fraction and Crocus parameterizations (as
in Fig. 2. in Lafaysse et al. (2017)). Additionally, we per-
form a deterministic ISBA-VS simulation with the default
snow cover fraction parameterization but the vegetation frac-
tion set to unity.
Model simulation periods for each site are in Table 2. For
each site, the model initial state was obtained by a spin-up
simulation from the start date (Table 2) to September 2020.
A total of 361 ensemble and deterministic simulations were
conducted.
2.6 Model evaluation metrics
Time series of daily averages are used to represent the results,
whereas mean absolute error (MAE), mean bias error (MBE)
and coefficient of determination (R2) are used in quantitative
model–data comparison. To detect possible biases in model
https://doi.org/10.5194/tc-18-231-2024 The Cryosphere, 18, 231–263, 2024
240 J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget
simulations, we use scatterplots and quantile–quantile plots
of sorted observations against sorted simulations. The sign
convention is so that the surface energy fluxes are presented
relative to the surface (i.e., negative flux means that surface is
losing energy). In time series, the turbulent flux observations
include all EC data (Sect. 2.4.2, quality flag≤ 2). We demon-
strate the results by using the winter season 2018–2019 as an
example time series period thanks to its best coverage of en-
ergy flux data (least gaps) and typical snow conditions on all
sites. For the scatter and quantile–quantile plots, only flux
data with quality flag ≤ 1 are used, and the results computed
as aggregated 6 h means include periods where observations
are available (referred to later as evaluation period). We com-
pare snow and snow-free conditions by grouping the results
into time windows where models and observations agree of
the ground conditions (snow or snow-free).
3 Results
3.1 Observed energy balance at peatland and forest
sites
The energy budget at high latitudes has a strong seasonal
variation driven by solar radiation (Fig. 3). In winter (De-
cember, January, February), longwave radiation balance to a
large extent determines Rn, particularly in northern Finland.
Daily average Rn is negative down to −50 Wm−2 and lower,
which implies considerable radiative cooling. Towards spring
the radiation budget is gradually counterbalanced by short-
wave radiation. On the peatlands, a large fraction of SWD is
reflected during snow cover, and dailyRn turns positive in the
late melting season (Fig. 3a, c). At the forest sites, the timing
of Rn becoming positive is less sensitive to the presence of
snow, as a large proportion of SWD is absorbed by vegeta-
tion. In summer, high solar elevation and the absence of the
reflective snow surface cause dailyRn to be up to 200 Wm−2.
Rn is balanced mostly by H and LE and to a lesser extent
snowpack–ground heat flux (Fig. 3). The residual line repre-
sents the amount of energy that would be required to close the
observed energy budget (Fig. 3). It includes changes in inter-
nal energy of the snowpack and vegetation but also reflects
the common energy balance closure problem in EC measure-
ments (Mauder et al., 2020) (see Sect. 4.4). The energy bal-
ance closure in snow-free conditions was typical for EC mea-
surements, ranging from 0.81 to 0.99 (Mauder et al., 2020).
The lack of snowpack heat flux and/or temperature profile
measurements did not enable assessing the closure during
snow cover periods. In winter, LE and G are small and the
radiative cooling is counterbalanced mostly by downwardH ,
corresponding to warming of the snowpack and/or vegetation
and cooling of the ambient air. Rn during winter falls lower
(more negative) on the northern sites, and thus also down-
wardH becomes stronger (daily average up to 50 Wm−2). In
summer, both H and LE are negative (upward), heating the
atmosphere, while downward G drives the warming of soil
profile. At all sites, LE increases along the growing season
and peaks approximately in July. In autumn, the turbulent
fluxes decrease in response to reduced Rn.
3.2 Peatland simulations
The sensitivity of surface heat fluxes and snow depth to
different ESCROC-E2 model parameterizations is shown in
Fig. 4. The spread corresponds to the difference between the
minimum and maximum of the ensemble. Notably, H has
relatively high spread, especially on N-WET, and the ob-
served H often lies near the limit or even outside the sim-
ulated range at both S-WET and N-WET. Modeled winter-
time LE is low and, as for H , the observed values are near
the limit or outside the simulated range, especially in spring.
LWU has strong day-to-day variation well captured by the
model, and the spread is rather small relative to the total flux.
3.2.1 Impact of alternative turbulence (CH)
parameterizations
To assess the sensitivity of snow depth simulations to al-
ternative turbulence parameterizations and to alternative
snow process options, we examined simulations where the
ESCROC-E2 members are grouped according to their turbu-
lent flux option (Fig. 5). During snow accumulation periods,
the spread is small and the groups are consistently overlap-
ping on both sites, indicating that the differences in snow ac-
cumulation and maximum snow depth are driven mostly by
the uncertainty of snow process descriptions. The spread in-
creases during and after snowmelt events, indicating higher
importance of turbulent fluxes on snowmelt dynamics. While
it is difficult to identify a group that fits observed snow depths
best, the winter melt event in 2018–2019 on N-WET is only
captured by the M98 and RI2 parameterizations.
These findings are consistent with the comparison of sim-
ulated H and LWU by the two deterministic runs (RIL-SOC
and M98-SOC) against observations (Fig. 6). With the RIL-
SOC parameterization, the magnitude of H is largely under-
estimated, while this bias is to a large extent corrected by
using M98-SOC. Improved simulation ofH and surface tem-
perature also entail improved LWU (Fig. 6), but the modeled
H fluxes still only moderately correlate (R2) with observa-
tions. In terms of LE, the simulations are not improved by
the M98-SOC (see Fig. S1), possibly due to low magnitude
and high relative uncertainty of wintertime LE over snow.
3.2.2 Radiative fluxes
We compare the simulated and observed albedo, SWU, LWU
and surface temperatures with snow-free and snow condi-
tions in Fig. 7. These experiments correspond to the deter-
ministic M98-SOC simulation.
The modeled SWU generally matches the observations
well, but the scatter increases with increasing SWD, indicat-
The Cryosphere, 18, 231–263, 2024 https://doi.org/10.5194/tc-18-231-2024
J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget 241
Figure 3. Observed daily averaged radiation budget (a, c, e, g) and surface energy budget (b, d, f, h) of hydrological year of 2018–2019.
Colored stacks represent the observed fluxes relative to the surface as shown in legends (i.e., incoming fluxes are positive and outgoing fluxes
negative). The dashed line in the energy budget plot corresponds to the residual after the sum of each energy component, whereas the dashed
line in the radiation plot shows the net radiation (Rn). Note the different scale in the left and right columns. Ground heat flux (G) is missing
on N-WET. The observed evolution of the snow depth (HS) is shown in gray (not to scale).
ing uncertainties in simulated albedo when shortwave forcing
is high over the snowpack in spring. These cause a slight un-
derestimation of simulated spring albedo, also visible in the
time series especially on N-WET (Fig. 7). Moreover, simu-
lated albedo tends to be overestimated during shallow snow
depth both in spring and autumn. This is because the ISBA-
FS approach assumes snow to completely cover the ground
regardless of the snow depth, while in reality the fractional
snow cover can lower the albedo. In contrast in May 2019,
the underestimation of albedo on N-WET is due to an incor-
rect timing of snowmelt (too early snow disappearance in the
simulation). The mean absolute errors in simulating SWU are
small and of similar magnitude (from ∼ 4 to 9 Wm−2) both
for snow and snow-free conditions.
Warmer surface temperatures during the snow-free season
result in higher LWU compared to winter and spring (Fig. 7).
The surface temperatures and LWU are generally well simu-
lated across sites and ground conditions at least with the pre-
sented time intervals. During snow cover, the upper tail of the
radiation distribution is slightly higher than simulated; how-
ever, the mean biases are generally very low. There are no
other visible biases in LWU simulations and the other met-
rics are also very good, consistent with Fig. 4. The mean ab-
solute errors in simulating LWU are similar for snow and
bare ground (∼ 3 to 8 Wm−2).
https://doi.org/10.5194/tc-18-231-2024 The Cryosphere, 18, 231–263, 2024
242 J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget
Figure 4. Time series of daily averaged surface heat flux spread simulated by ISBA-FS with ESCROC-E2 35 ensemble members against
corresponding observed values during the 2018–2019 snow season. H , LE and LWU correspond to sensible heat, latent heat and upward
longwave radiation fluxes, respectively. The observed and simulated evolution of snow depth (HS) are shown in gray.
Figure 5. Time series of snow depths (HS) simulated by ISBA-FS ESCROC-E2. The 35 ensemble members are grouped by their turbulent
flux parameterization, and the spread of each group is presented in colored ranges. Observed snow depths are presented in black dots and
dashed lines.
The Cryosphere, 18, 231–263, 2024 https://doi.org/10.5194/tc-18-231-2024
J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget 243
Figure 6. Scatterplots and quantile–quantile plots of sensible heat
flux (H ) and upward longwave radiation (LWU) during snow cover
for the evaluation period with the RIL-SOC and M98-SOC turbu-
lence parameterizations.
3.2.3 Soil thermal regime
The effect of soil parameterization on simulated soil temper-
ature dynamics at S-WET is shown in Fig. 8. Due to shallow
water table, the soil profile remains nearly saturated through-
out the year. As the porosity and field capacity in the M98-
SOC parameterization are much higher than in the M98-
MIN, the former also has significantly higher heat capacity
and smaller thermal diffusivity. This means soil temperature
variations are attenuated in M98-SOC compared to M98-
MIN, and this attenuation becomes increasingly important in
deeper soil layers (Fig. 8). The results show that including a
realistic soil profile (SOC) greatly improves the peatland soil
temperature simulations at depths 50–70 cm but only slightly
close to the surface (0–10 cm) (see Fig. S2 for comparisons
of more soil depths). On both sites, the simulated surface soil
temperature variations in summer are greater than observed.
This is presumably because ISBA does not include the insu-
lating moss/litter layer on top of the peat soil, as well as due
to water table dynamics, potentially affected by lateral flows
not accounted for. Due to the weak influence on the surface
soil temperatures, the soil parameterization (M98-SOC vs.
M98-MIN) does not significantly affect the simulated snow
depth (Fig. S2).
3.3 Forest simulations
3.3.1 Impact of vegetation representation on snow
depth
The three different vegetation representations (Sect. 2.3.1)
have a highly contrasting effect on the forest energy budget,
snowpack and soil temperature simulations. In general, the
snowpack simulations for the forest sites are poorer than for
the peatland sites; however, the observed snow depths also
vary considerably within the forests (see OBS in Fig. 9 and
Sect. 4.4).
The simulated snow depth with the ISBA-VS (compos-
ite soil–vegetation and varying snow cover fraction, Fig. 2b)
does not agree with the observations; the model version
heavily overestimates accumulation on N-FOR in 2021 and
predicts extremely rapid, strong and too early melt events
(Fig. 9). Replacing the default snow cover fraction param-
eter (wsw = 5) with wsw = 0.2 (used for hydrological mod-
eling in Le Moigne et al., 2020) yields slightly better snow
depth dynamics for N-FOR, but the results remain unsatis-
factory (Fig. C1 in the Appendix). The different sensitivity
of the wsw parameter for S-FOR and N-FOR simulations is
explored via soil temperature simulations in Fig. C2; with
the default snow cover fraction parameter, particular warm
events on N-FOR heat up the soil, causing the snowpack
to melt, while simulation with wsw = 0.2 manages to retain
freezing soil temperatures.
MEB (explicit canopy, ground and snowpack energy bal-
ance) simulates the snow accumulation periods at N-FOR
very well, but peak snow is reached too early, and maximum
snow depths are underestimated. This is due to the combined
impact of overestimated compaction and too early start and
progression of the snowmelt. The role of both processes was
evident from the comparison of modeled and observed snow
water equivalent (see Fig. S6).
ISBA-FS performs better during the snow accumulation
period, with simulated snow depths very close to observa-
tions. However, the ablation of snow is too rapid, and the
final melt-out dates are close to those simulated by MEB.
On S-FOR, MEB captures snow accumulation (including
peak snow depths), melt dynamics and final melt-out dates
rather well. ISBA-FS predictions are generally close to MEB.
As MEB only considers one option for turbulent exchange
(RIL), the spread of the ensemble is smaller than for the
ISBA configurations (Fig. 9). The uncertainties of other snow
processes accounted for in ESCROC-E2 are not sufficient to
explain the discrepancies between simulated and observed
snow depths, suggesting that uncertainties in the canopy pro-
cess representations prevail in these simulations.
https://doi.org/10.5194/tc-18-231-2024 The Cryosphere, 18, 231–263, 2024
244 J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget
Figure 7. Scatterplots and quantile–quantile plots of modeled against observed upward shortwave radiation (SWU) (a, b) and upward
longwave radiation (LWU) (c, d) on peatland sites with snow cover (w/ snow) and without snow cover (w/o snow) as well as the time evolution
of 5 d rolling means of albedos (LSA) and surface temperature (Ts) as simulated and observed from September 2018 to September 2019 (e,
f). The evolution of the snow depth (HS) is not to scale.
Figure 8. Effect of soil parameterization on soil temperatures (a–c) during a hydrological year at the S-WET site. M98-MIN refers to mineral
soil and M98-SOC to peat soil. The observed soil temperatures are compared to the closest model layer; the depths of measurements and
simulations are presented in each panel.
3.3.2 Impact of vegetation representation on soil
temperature
Similar to the snow depth, soil temperature predictions by
ISBA-VS are erroneous, with drastically underestimated
temperatures and unrealistic dynamics (Fig. 10) (see more
soil depths in Fig. S3). While MEB and ISBA-FS provided
very similar snow depth, the soil temperatures simulated by
MEB agree better with the observations, although there is a
cold bias in autumn and a warm bias in summer (Fig. 10). On
N-FOR the warm bias in winter by MEB may be important
for determining the soil frost regime. Interestingly, ISBA-FS
seems to capture the winter soil temperatures better on N-
FOR, but this may be due to the larger cold bias in autumn
The Cryosphere, 18, 231–263, 2024 https://doi.org/10.5194/tc-18-231-2024
J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget 245
Figure 9. Effect of alternative configurations of ISBA and MEB on the snow depth (HS). The envelopes visualize the corresponding ensemble
spreads between minimum and maximum values.
likely caused by the lack of explicit litter and canopy layers.
All model versions tend to overestimate day-to-day tempera-
ture variability.
3.3.3 Impact of vegetation representation on surface
energy fluxes
Figure 11 compares the deterministic simulations (Sect. 2.5)
by MEB and ISBA-VS against observed above-canopy en-
ergy fluxes at N-FOR. The snow cover periods are defined ac-
cording to agreement between MEB simulations and obser-
vations, and thus the ISBA-VS simulations are often snow-
free as seen in Fig. 9.
MEB is superior to ISBA-VS in simulating energy fluxes.
SWU simulations with snow cover are clearly improved by
MEB, but the spread remains relatively large and albedo is
underestimated when incoming radiation is small and over-
estimated when incoming shortwave radiation is higher. The
time evolution of albedo on N-FOR and S-FOR is presented
in Sect. 3.3.4. The LWU is very well simulated by both model
configurations. Turbulent fluxes are clearly better simulated
by MEB, but the performance metrics of turbulent fluxes are
worse than for radiative fluxes. ISBA-VS uses the vegeta-
tion fraction parameter to scale the partitioning of latent heat
flux between vegetation and soil (Eq. 4). However, because
the same roughness length and thus turbulent exchange co-
efficient (CH) is used for both soil and vegetation, the soil
evaporation and snow sublimation are likely overestimated
and result in clearly wrong partitioning between H and LE
(Fig. 11c, d and g, h). In the case of N-FOR, especially the
summer energy fluxes were majorly improved by simply as-
signing the vegetation fraction to unity (i.e., full vegetation
coverage and no soil evaporation; see Fig. S7).
3.3.4 Evolution of albedo
Figure 12 illustrates the time evolution of modeled and ob-
served albedo and the shortwave components in 2018–2019
on both forest sites. Compared to the measurements, the
modeled early and midwinter albedos are underestimated,
while the spring albedos are slightly overestimated, consis-
tent with results in Sect. 3.3.3. The likely reason for winter
albedo underestimation is because the models do not repre-
sent changes in albedo due to intercepted snow. The overesti-
mation in spring is presumably due to representing the effec-
tive albedo of snow and forest canopy with only bulk canopy
parameters, as well as effect of spring needle and litter fall
decreasing snow albedo. Moreover, the simulated albedo is
dominated by the vegetation albedo parameter, and thus it is
not highly sensitive to snowpack albedo dynamics.
3.4 Summary: surface energy budget on peatland and
forest sites
Finally, to sum up the whole surface energy budget, we com-
pare how the simulated Rn and turbulent fluxes (H+ LE)
match the observations at the four sites. These determin-
istic simulations are conducted with simulation setups that
provided the best fit to data: the deterministic simulation as
https://doi.org/10.5194/tc-18-231-2024 The Cryosphere, 18, 231–263, 2024
246 J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget
Figure 10. Effect of alternative configurations of ISBA and MEB on soil temperatures. The envelopes visualize the corresponding ensemble
spreads. The observed evolution of the snow depth (HS) is not to scale. The observed soil temperatures are compared to the closest model
layer; the depths of measurements and simulations are presented in each panel.
Figure 11. Simulated against observed upward shortwave and longwave radiation (SWU and LWU, columns 1 and 2) and turbulent fluxes
(H and LE, columns 3 and 4) on the N-FOR site for the full evaluation period. Ground conditions are presented as (i) with snow cover (w/
snow, row 1) and (ii) without snow cover (w/o snow, row 2).
M98-SOC for the peatland sites and deterministic MEB sim-
ulation for the forest sites (Fig. 13).
Despite the challenges in simulating snow depth evolution
at the forest sites, the energy budget simulations are gener-
ally better than on the peatlands. Due to the challenges to
accurately simulate albedo and surface temperatures on the
open sites, the simulated Rn is considerably worse on peat-
land sites than on forests (see Fig. 7). In particular, the high
Rn periods, representing the spring conditions, are biased on
peatland sites, while the negative Rn period (i.e., the win-
ter conditions) are simulated rather well. The challenges in
describing forest wintertime albedo and thus SWU (as in
Figs. 12 and 11) do not significantly bias the Rn simula-
tions, as in wintertime the shortwave radiation balance has a
small role compared to the longwave radiation balance. The
results propose that canopy temperature, which particularly
in dense forests (e.g., S-FOR) has a central role for upward
longwave radiation, must be adequately simulated by MEB.
When it comes to the turbulent fluxes, the simulations cap-
ture the main seasonal patterns. However, there are still high
uncertainties (scatter) both on peatland and forest sites. The
relative uncertainties in simulated and observed energy fluxes
are significantly greater in winter than in summer. Perfor-
mance of the simulated summer energy fluxes is very good
(Fig. S4).
The Cryosphere, 18, 231–263, 2024 https://doi.org/10.5194/tc-18-231-2024
J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget 247
Figure 12. Simulated and observed albedo (LSA) and downward and upward shortwave radiation (SWD and SWU) in 2018–2019. LSA is
presented as 5 d rolling means. The observed evolution of the snow depth (HS) is not to scale.
Figure 13. Simulated (MOD) against observed (OBS) daily surface energy budget during winter 2018–2019. The left column shows net
radiation (Rn) and the right column presents the sum of turbulent fluxes (H+ LE). The scatterplots represent the full simulation periods
when snow cover was present. The observed evolution of the snow depth (HS) is not to scale.
https://doi.org/10.5194/tc-18-231-2024 The Cryosphere, 18, 231–263, 2024
248 J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget
Table 3. Occurrence of different turbulence regimes at S-WET and
N-WET. The regimes are defined based on the bulk Richardson
number (Ri): unstable conditions as Ri < 0, weakly stable condi-
tions as 0≤ Ri ≤ 0.25 and strongly stable conditions as Ri > 0.25.
Turbulence regimes
Site Surface Unstable Weakly stable Stable
[%] [%] [%]
N-WET all 35.1 15.1 49.8
N-WET snow 13.3 16.5 70.2
N-WET ground 63.2 15.0 21.8
S-WET all 59.7 33.7 6.6
S-WET snow 26.5 54.6 18.8
S-WET ground 78.9 20.5 0.6
4 Discussion
4.1 Insights on energy flux partitioning in boreal
environments
The energy budget observed with this novel dataset over
boreal and subarctic peatlands showed that despite the pre-
vailing stable atmospheric conditions (Table 3), the radia-
tive cooling was mostly counterbalanced by sensible heat
flux (H ) (Fig. 3). During the snow season, the dominating
regimes were strongly stable at N-WET (70.2 %) and weakly
stable at S-WET (54.6 %). Despite the stronger stability at
N-WET, we observed higher H fluxes compared to S-WET
and considerably higher H than Lackner et al. (2022) at the
Canadian site dominated by weakly stable conditions. Thus,
we presume that greater radiative cooling leads to a stronger
near-surface air temperature gradient and larger downward
sensible heat flux.
The small role of the ground heat flux (G) at the studied
peatlands is due to the high heat capacity of the peat soil
and its large water storage which progressively freezes from
the top, keeping the temperatures in the soil–snow interface
nearly constant at a minimum of 0 ◦C (Fig. 8). Other studies
in tundra environments of the Canadian Arctic, Siberia and
Svalbard have reported larger contributions of G to the win-
tertime energy budget (Lackner et al., 2022; Langer et al.,
2011; Westermann et al., 2009).
We observed a shorter snowmelt period in the peatlands
compared to adjacent forest sites (see Fig. S5), in line with
Lundquist et al. (2013), who proposed that increased canopy
shading (delaying melt) outweighs the impact of longwave
radiation enhancement (accelerating melt). Our datasets sup-
port this (Fig. S5); however, the forest sites tended to also
accumulate more snow than the peatland sites (wind erosion
is presumably higher on peatlands), which may have con-
tributed to longer snow duration in the forest.
4.2 Implications for simulating snow and energy
balance at peatland sites
Our results highlighted the uncertainties in modeling turbu-
lent fluxes over snowpack and identified the turbulent ex-
change parameterizations (M98-SOC and RI2-SOC) that im-
prove the simulated surface energy fluxes and snowpack dy-
namics at high-latitude winter conditions (Figs. 5, 6). In sta-
ble (winter) conditions, the uncertainties in turbulent fluxes
are in line with Menard et al. (2021), Conway et al. (2018),
and Lapo et al. (2019) and larger than in unstable (summer)
conditions (Fig. S4). Moreover, the turbulent fluxes of sensi-
ble and latent heat have greater uncertainty than the radiation
balance components (Fig. 13). The ESCROC-E2 simulations
showed that the turbulent exchange parameterizations also
have an impact on snowmelt simulations, in line with sim-
ulations at Col de Porte, France, and ESM-SnowMIP sites
(Menard et al., 2021). In contrast, Lackner et al. (2022) found
only small differences between the Crocus turbulence param-
eterizations in their study in the Canadian Arctic, most likely
due to smaller sensible heat fluxes and less frequent strongly
stable conditions.
Improved surface temperature simulations by the M98-
SOC (absolute biases decreased by 0.3 ◦C at S-WET and
0.4 ◦C at N-WET compared to RIL-SOC) provide support
to Martin and Lejeune (1998) and Gouttevin et al. (2023),
who adjusted the turbulent fluxes under stable conditions to
reproduce surface and air temperature observations. The de-
fault ISBA turbulent flux parameterization (RIL), although
widely used, e.g., in numerical weather prediction and gen-
eral circulation models (Mahfouf et al., 1995; Salas-Mélia
et al., 2005; Voldoire et al., 2013, 2019), provided the poorest
fit with the observed surface heat fluxes and produced a cold
bias in snow surface temperature (−0.4 ◦C at S-WET and
−1.1 ◦C at N-WET). This finding is consistent with ESM-
SnowMIP (Menard et al., 2021), where the default configu-
ration of Crocus had one of the lowest skills for surface tem-
perature simulations (−2 ◦C mean cold bias) among the com-
pared snow models. Even with the M98-SOC simulation we
found a rather low skill of turbulent flux simulations. To sum-
marize, our findings highlight the limitations of the Monin–
Obukhov similarity theory to simulate turbulent fluxes un-
der stable atmospheric conditions and emphasize the need
for further model developments with observations in various
environments.
The soil temperature simulations confirmed that it is nec-
essary to realistically describe the organic peat soil hydraulic
and thermal properties (Menberu et al., 2021; Morris et al.,
2022; Mustamo et al., 2019) to accurately simulate soil ther-
mal regime and consequent freezing–thawing processes in
peatlands (Fig. 8) (Dankers et al., 2011; Lawrence and Slater,
2008; Nicolsky et al., 2007). This is line with Decharme et al.
(2016), who implemented SOC parameterization in ISBA
and showed improved soil temperature simulations across
northern Eurasia. The implementation of water table dynam-
The Cryosphere, 18, 231–263, 2024 https://doi.org/10.5194/tc-18-231-2024
J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget 249
ics and lateral flow could further improve the soil tempera-
ture simulations on boreal and subarctic peatlands. The ther-
mal state and ice/liquid water content also have major cas-
cading effects on runoff generation during snowmelt (Ala-
Aho et al., 2021). Moreover, the interactions between low
vegetation and snow are likely improved by using explicit
vegetation (MEB in SURFEX). However, as MEB has never
been applied on snow-covered low vegetation, additional de-
velopment and evaluation beyond the scope of this study
would have been required.
4.3 Implications for simulations at forest sites
4.3.1 ISBA-VS
We found the ISBA-VS to be largely biased towards correctly
simulated surface energy fluxes at the expense of poor soil
temperature and snow depth simulations, as a major part of
the composite was always directly coupled to the atmosphere
(Figs. 9, 10). ISBA-VS with correct tuning (e.g., setting the
veg fraction to unity on N-FOR) may be an imperfect but
sufficient compromise to simulate snow-free forests in appli-
cations that first and foremost require grid-cell-averaged sur-
face energy fluxes (i.e., numerical weather prediction). How-
ever, we agree with Napoly et al. (2020) that the snow cover
fraction approach of ISBA (Fig. 2b) is essentially a compro-
mise that attempts to retain the insulating impact of the snow-
pack over the soil while still simulating turbulent exchange
from the vegetation. The energy exchange between the at-
mosphere and the soil–vegetation composite directly impacts
the snowpack, and it led to strongly biased snow depth sim-
ulations, similar to Napoly et al. (2020). This fractional ap-
proach and high sensitivity of empirically based parameters
(e.g., veg fraction) highlights the uncertainties of ISBA-VS
to provide accurate year-round lower-boundary conditions of
boreal forests for numerical weather prediction and general
circulation model applications. Also considering the very
low skill obtained in snow depth and soil temperatures for
this configuration, its use in hydrological applications or sur-
face offline reanalyses (Le Moigne et al., 2020) is highly
questionable. Nevertheless, local-scale evaluation might not
directly translate to large-scale spatial simulations, as further
discussed in Sect. 4.4.
4.3.2 ISBA-FS
We found that snow and soil simulations in forests were
strongly improved when the snow cover fraction was set to
unity (ISBA-FS), allowing the snowpack to fully insulate
the soil, similarly to the simulations at the peatland sites
(Fig. 9). These results suggest that if the focus is on snow-
pack dynamics and soil temperature simulations, ignoring
snow–vegetation interactions is a better compromise than
using varying snow cover fraction as it is currently imple-
mented in SURFEX. Consequently, ISBA-FS should be pre-
ferred to ISBA-VS in surface reanalyses as in Brun et al.
(2013) and Vernay et al. (2022). However, ISBA-FS reaches
its conceptual limits when forest energy balance, snow and
soil state variables are all of interest. For instance, neglect-
ing snow interception and subsequent canopy snow losses
may cause large errors in simulated snow water equivalent
in dense forests, and unrealistic contribution of the canopy
evapotranspiration may be expected if reanalyses are further
used for hydrological modeling. Obviously, highly biased
surface energy fluxes would also be expected for any cou-
pling with an atmospheric model.
4.3.3 MEB
MEB appears as a better compromise than ISBA-VS and
ISBA-FS for modeling forest energy exchanges; the snow
depth and soil temperature simulations were highly improved
compared to ISBA-VS, while surface–atmosphere energy
fluxes are obviously much more realistic than with ISBA-FS.
Significant improvements in energy flux simulations were
also obtained compared to ISBA-VS (Fig. 11).
Nevertheless, two systematic biases affecting upward
shortwave radiation were identified: albedo was underesti-
mated in winter and overestimated in spring (Fig. 12). The
winter albedo underestimation was most likely because in-
tercepted snow increased the observed albedo, a process that
is not accounted for in MEB (Napoly et al., 2020). This as-
sumption is based on Pomeroy and Dion (1996), while the
increase in forest albedo by intercepted snow has been more
recently shown (Webster and Jonas, 2018), and simple de-
scriptions can already be found in some forest snow models
(Mazzotti et al., 2020a). Although our results propose that the
intercepted snow has a clear impact on the albedo, its impact
on Rn was small. The spring albedo bias is in line with Malle
et al. (2021), who found albedo at sparse boreal forests to be
overestimated by the LSM CLM5. Big-leaf canopy parame-
terizations may fail to fully capture canopy shading, partic-
ularly at low solar elevation angles typical of high latitudes
(Malle et al., 2021).
In contrast to Napoly et al. (2020), we found MEB to simu-
late snowmelt too early, especially on N-FOR (Fig. 9). These
errors are partly explained by inaccuracies in canopy radia-
tive transfer, but they also suggest errors in simulated below-
canopy surface heat fluxes. The differences in snow simu-
lations were rather small between MEB and ISBA-FS, es-
pecially at N-FOR, suggesting that sparse canopies did not
majorly alter simulated snow accumulation and ablation, at
least when compared to snow depth observations between
the trees. Consistently, Meriö et al. (2023) have shown de-
creased snow depths at the immediate vicinity of tree trunks
but high snow depth between trees at N-FOR.
https://doi.org/10.5194/tc-18-231-2024 The Cryosphere, 18, 231–263, 2024
250 J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget
4.4 Limitations and outlook
Our eddy-covariance fluxes are among the longest datasets
ever used for the evaluation of turbulent flux simulations over
snow. The EC data, however, contain both random and sys-
tematic uncertainties (e.g., Aubinet et al., 2012). The abso-
lute values of winter H and LE are small, and their rela-
tive uncertainty is high; compared to summertime measure-
ments, the wintertime energy balance closure ratio is typi-
cally poorer, particularly at the northern ecosystems (Reba
et al., 2009; Molotch et al., 2009; Launiainen, 2010). As our
analysis used numerous site years from multiple sites, and we
used established quality criteria for filtering the EC fluxes,
we expect that uncertainties in flux data do not significantly
affect the study results. Moreover, the conclusions regarding
the validity of each model version were not affected by the
selected flux quality flag (Sect. 2.4.2). Intrinsic uncertainties
in meteorological forcing are known to exist, especially in
northern conditions (instrument freezing, snow blocking, un-
dercatch, etc., Stuefer et al., 2020), and data gaps further add
up possible sources of errors. Uncertainties in model forc-
ing can affect model–data comparisons, especially during the
gap-filled periods (Raleigh et al., 2015).
Potentially important snow processes on subarctic sites are
still absent in Crocus, including wind-induced snow transport
and internal water vapor transfer due to a large temperature
gradient in the snowpack. Wind can move and change the
properties of snow (Pomeroy and Essery, 1999; Meriö et al.,
2023; Liston and Sturm, 2002) and is especially noticeable
on open peatlands. In Crocus, wind modifies the properties
of falling snow (Vionnet et al., 2012) but without any lateral
transport or modifications of the mass. Although we achieved
satisfactory model performance without accounting for this
process, Meriö et al. (2023) showed notable wind transport
in transition zones between open peatland and forest at the
N-WET site. Although the spatial scale of wind transport
prevents an explicit simulation of this process in large-scale
LSMs, improved parameterizations of the wind impact on
snow properties should be considered in the future. Omis-
sion of internal water vapor transfer by diffusion and/or con-
vection has been suspected to be responsible for errors in
simulated snow properties (density, microstructure) in Arctic
snowpack (Barrere et al., 2017; Domine et al., 2018) and con-
sequently in thermal conductivity and soil thermal regime. A
realistic implementation of water vapor transfer within the
snowpack is lacking in most state-of-the-art LSMs. Comple-
mentary observations and model developments are required
to understand if the simulated snow properties are affected by
this kind of error. Furthermore, the spring and autumn con-
ditions on the peatlands are particularly difficult to correctly
simulate; in addition to the snow cover, also, e.g., ponding of
liquid water and refreezing of the ponds are not uncommon
(Noor et al., 2022) and can alter the albedo.
In forests, the spatial heterogeneity of snow cover can be
high, as demonstrated by numerous studies (Marttila et al.,
2021; Mazzotti et al., 2020b; Noor et al., 2022) and con-
firmed by our data (Fig. 9). The small-scale forest struc-
ture has an important role in the evolution of the snow cover
and may affect the representativeness of point measurements
(Bouchard et al., 2022). Consequently, the comparison of
point observations and models intended for forest stand and
larger scales (such as the big-leaf approach of MEB) can
suffer from this scale mismatch (Essery et al., 2009; Rut-
ter et al., 2009). More realistic below- and above-canopy
heat flux simulations could be achieved by more sophisti-
cated canopy representations, including multiple layers and
species (Bonan et al., 2021; McGowan et al., 2017; Launi-
ainen et al., 2015; Gouttevin et al., 2015). For site-level or
limited-area modeling, high-resolution models that explicitly
resolve tree-scale canopy structure are a promising alterna-
tive to traditional LSMs (Broxton et al., 2015; Mazzotti et al.,
2020b).
The generality of our findings should be tested by addi-
tional snow model and LSM evaluation studies, extended to
more contrasting climates and a wider range of ecosystem
types. For this purpose, model evaluation datasets should be
complemented with more boreal and Arctic sites and obser-
vations of all components of surface energy balance, partic-
ularly turbulent fluxes.
5 Conclusions
We used eddy-covariance-based energy flux data, radiation
balance, and snow depth and soil temperature measurements
in two boreal and subarctic peatlands and forests to evalu-
ate turbulent exchange parameterizations and alternative ap-
proaches to represent the soil and vegetation continuum in
land surface models. We used the SURFEX platform, but
our findings are largely transferable to other model systems.
Our evaluation with the ensemble snowpack parameteriza-
tions (ESCROC-E2) ensures that uncertainties in snow pro-
cesses (not evaluated in this study) do not affect the robust-
ness of our main conclusions.
Peatland simulations showed that using a stability cor-
rection function that increases the turbulent exchange under
stable atmospheric conditions is imperative to simulate the
snowpack energy budget. This adjustment led to major im-
provements under stable conditions during snow cover, but
the model performance still remained lower than in snow-
free conditions. Furthermore, correct hydraulic and thermal
parameterization of organic peat soils was found necessary
to reproduce the observed soil thermal regime in peatlands.
The findings have direct implications for modeling snow dy-
namics, peatland hydrology and permafrost dynamics.
The surface energy budgets of forest sites were well
simulated by the explicit big-leaf approach (MEB), while
the composite soil–vegetation approach (ISBA-VS) perfor-
mance was satisfactory only after an adjustment of a sensi-
tive vegetation fraction parameter. In particular, shortwave
The Cryosphere, 18, 231–263, 2024 https://doi.org/10.5194/tc-18-231-2024
J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget 251
and longwave radiation balances were simulated well by
both approaches, whereas the turbulent fluxes had signifi-
cantly higher uncertainty. Only the explicit vegetation model
(MEB) was able to simultaneously simulate realistic surface
energy budget, snow and soil conditions. The composite ap-
proaches only succeeded in simulating either the correct sur-
face energy budget (ISBA-VS) or snow and soil conditions
(ISBA-FS). Furthermore, we demonstrated that the compos-
ite approaches rely on a previously poorly documented pa-
rameterization of the snow cover fraction with high sensitiv-
ity on model outputs despite a limited physical interpretation.
With well-selected model configuration and parameteriza-
tion, the SURFEX model platform can realistically simulate
surface energy fluxes and snow and soil conditions in the
subarctic and boreal peatlands and forests. The common ver-
sion of ISBA (ISBA-VS) can provide rather realistic lower-
boundary conditions for numerical weather prediction and
global circulation models, at the expense of non-realistic pre-
dictions of forest snow and soil conditions necessary for hy-
drological applications. We expect that the future inclusion
of MEB in operational systems will reconcile these applica-
tions. Our results can be used to inform the choice of model
configuration for studies of subarctic and boreal region ecol-
ogy, hydrology and biogeochemistry under the ongoing en-
vironmental change.
https://doi.org/10.5194/tc-18-231-2024 The Cryosphere, 18, 231–263, 2024
252 J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget
Appendix A: Tables of abbreviations, acronyms and
mathematical symbols
A1 Table of acronyms and abbreviations
Acronyms and abbreviations Definition
LSM Land surface model
LSA Land surface albedo
LAI Leaf area index
LWD Downward longwave radiation
LWU Upward longwave radiation
SWD Downward shortwave radiation
SWU Upward shortwave radiation
H Sensible heat
LE Latent heat
Rn Net radiation
G Snowpack–ground heat flux
SOC Soil organic content
EC Eddy covariance
HS Snow depth
RIL Classical Louis (1979) formula for the turbulent exchange coefficient
RI1 RIL with threshold at Ri = 0.1
RI2 RIL with threshold at Ri = 0.026
M98 Martin and Lejeune (1998) formula for CH
BOGR Boreal grass
BONE Boreal needleleaf evergreen
FMI Finnish Meteorological Institute
SMEAR Station for Measuring Forest Ecosystem–Atmosphere Relations
FST Flux stationarity
The Cryosphere, 18, 231–263, 2024 https://doi.org/10.5194/tc-18-231-2024
J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget 253
A2 Table of mathematical symbols
Symbol Definition
 surface emissivity
σ Stefan–Boltzmann constant
Ts surface temperature
ρa air density
ρw liquid water density
ρ snow density
cp specific heat capacity
Va wind speed
Ta air temperature
CH turbulent exchange coefficient
ET total evapotranspiration
Eg evaporation from the bare soil surface
Ec evaporation of intercepted water on the canopy
Etr transpiration from the vegetation
Si sublimation from bare soil ice
Lv latent heat of vaporization
Lf latent heat of fusion
veg fraction of vegetation cover
qsat(Ts) saturated specific humidity at the surface
qa(Ts) atmospheric specific humidity
hu dimensionless relative humidity at the ground surface related to the superficial soil moisture content
hv dimensionless Halstead coefficient describing the Ec and Etr partitioning
between the leaves covered and not covered by intercepted water
zu reference height of the wind speed
za reference height of the air temperature and humidity
z0t roughness height for heat
k von Kármán constant
Ri bulk Richardson number
Cv effective heat capacity of the canopy
Cg,1 effective heat capacity of the ground surface/litter layer (uppermost layer)
Cn,1 effective heat capacity of the snowpack (uppermost layer)
Tv temperature of the canopy
Tg,1 temperature of the ground surface/litter layer (uppermost layer)
Tn,1 temperature of the snowpack (uppermost layer)
Rnv net radiation of the canopy
Rng net radiation of the ground surface/litter layer
Rnn net radiation of the snowpack
k snow effective thermal conductivity
LEs latent heat flux of the snowpack
Gn heat conduction at the snow–soil interface
psn effective snow cover fraction
psnv effective snow cover fraction of the vegetation
psng effective snow cover fraction of the soil
HSg threshold value for height of the snow
wsw coefficient relating to vegetation characteristics
E2 ESCROC optimized standard subensemble
P precipitation rate
Pice snowfall rate
Pliq rainfall rate
r diffuse-to-total-shortwave-radiation ratio
µ cosine of the sun zenith angle
u∗ friction velocity
https://doi.org/10.5194/tc-18-231-2024 The Cryosphere, 18, 231–263, 2024
254 J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget
Appendix B: Model evaluation data availability and
radiation forcing evaluation
The availability periods of the model evaluation data are pre-
sented in Table B1.
Table B1. Model evaluation data availability for each site.
Variable N-WET S-WET N-FOR S-FOR
Snow depth 2017-11–2021-05 2016-09–2021-07 2013-09–2020-09 2008-09–2021-07
Soil temperature 2013-09–2019-12 2017-06–2021-07 2016-09–2020-09 2008-09–2021-07
Upward LW flux 2017-07–2021-06 2016-09–2021-07 2013-09–2021-07 2008-09–2021-07
Upward SW flux 2017-07–2021-06 2016-09–2021-07 2013-09–2021-07 2008-09–2021-07
Sensible heat flux 2013-09–2021-06 2016-09–2020-12 2013-09–2021-07 2008-09–2021-07
Latent heat flux 2013-09–2021-06 2016-09–2020-12 2013-09–2021-07 2008-09–2021-07
Ground heat flux 2013-09–2017-07 2016-09–2021-07 2013-09–2021-02 2008-09–2021-07
Figure B1. Comparison of longwave radiation forcing data between ERA5 and site observations (OBS) on S-FOR.
Appendix C: Effect of snow cover fraction parameter
wsw
Table C1. Summary of snow cover fraction and values of wsw used in different SURFEX/ISBA applications. The parameter wsw is rarely
documented, and hence these application-specific values were obtained through communications with the authors.
Application Name Domain Resolution Snow fraction wsw Reference
Numerical weather prediction AROME, ARPEGE Europe (many) 1.3–10 km varying 5 Bengtsson et al. (2017)
Courtier et al. (1991)
Global climate modeling CNRM-CM6 Global 100 km varying 2 Decharme et al. (2019)
Regional climate modeling CNRM-AROME European Alps 2.5 km varying 1 Caillaud et al. (2021)
Regional climate modeling CNRM-ALADIN Europe, North Africa 12 km varying 2 Nabat et al. (2020)
Hydrological modeling SIM2 France 8 km varying 0.2 Le Moigne et al. (2020)
Regional reanalysis CERRA-Land Europe 5.5 km varying 0.1 Verrelle et al. (2021)
Snow cover reanalysis S2M French Alps massif-scale full (1) – Vernay et al. (2022)
Snow cover reanalysis ERA-Interim–Crocus Northern Eurasia 80 km full (1) – Brun et al. (2013)
Avalanche hazard forecasting S2M French Alps massif-scale full (1) – Morin et al. (2020)
The Cryosphere, 18, 231–263, 2024 https://doi.org/10.5194/tc-18-231-2024
J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget 255
Figure C1. Effect of wsw parameter on snow depths simulated by ISBA-VS. The envelopes visualize the corresponding ensemble spreads
between minimum and maximum values.
Figure C2. The sensitivity of the wsw parameter on 6 h surface soil temperature simulations and snow water equivalent (SWE). Simulations
are represented by one member of the ESCROC-E2 ISBA-VS.
https://doi.org/10.5194/tc-18-231-2024 The Cryosphere, 18, 231–263, 2024
256 J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget
Code availability. SURFEX is an open-source project (https:
//www.umr-cnrm.fr/surfex/IMG/pdf/surfex_scidoc_v8.1.pdf, Le
Moigne, 2018), but it requires registration. The full procedure
and instructions are available at https://opensource.umr-cnrm.fr/
projects/snowtools_git/wiki/Procedure_for_new_users (Lafaysse
et al., 2023), and the SURFEX version used in this study is
archived at https://doi.org/10.5281/zenodo.10473208 (Nousu
et al., 2023). The SURFEX version used in this work is also
available in Git (tagged as boreal_ecosystems). The ESCROC
version developed for external SURFEX users is available at
https://doi.org/10.5281/zenodo.10474747 (Cluzet and Nousu,
2023).
Data availability. Meteorological data, evaluation data, and SUR-
FEX specific namelist and forcing files are available as an
electronic supplement at https://doi.org/10.5281/zenodo.8252267
(Nousu et al., 2023) under the terms and conditions of the Creative
Commons Attribution 4.0 International license.
Supplement. The supplement related to this article is available on-
line at: https://doi.org/10.5194/tc-18-231-2024-supplement.
Author contributions. JPN, ML, SL, PA and HM designed the re-
search. JPN led the study and was in charge of the experiments.
JPN, GM and ML performed the model evaluations with scientific
contributions of SL, HM and PA. JPN was responsible for writ-
ing the article, with significant contributions from GM, ML and SL
and input from all the authors. BC and MF supported the ensemble
simulations and other technical aspects. MA, AL, PK, PA and HM
provided surface energy data and ancillary data.
Competing interests. The contact author has declared that none of
the authors has any competing interests.
Disclaimer. Publisher’s note: Copernicus Publications remains
neutral with regard to jurisdictional claims made in the text, pub-
lished maps, institutional affiliations, or any other geographical rep-
resentation in this paper. While Copernicus Publications makes ev-
ery effort to include appropriate place names, the final responsibility
lies with the authors.
Acknowledgements. This work was funded by the Research
Council of Finland (RCF) (ArcI Profi 4). Giulia Mazzotti
was funded by the Swiss National Science Foundation (grant
no. P500PN_202741). Samuli Launiainen and Jari-Pekka Nousu ac-
knowledge the GreenFeedBack project from the EU Horizon Eu-
rope Framework Programme for Research and Innovation (grant
no. 101056921). Samuli Launiainen acknowledges the support of
RCF (LS-HYDRO, 356138). Pertti Ala-aho was funded by the
RCF Research Fellow grant (no. 347348). Pasi Kolari, Mika Au-
rela and Annalea Lohila acknowledge the support of the ACCC
Flagship funded by the RCF (no. 337549 and no. 337552) and
ICOS Finland by the University of Helsinki and the Ministry of
Transport and Communication. Mika Aurela and Annalea Lohila
also acknowledge the RCF (UPFORMET, grant no. 308511). Jari-
Pekka Nousu thanks Riitta ja Jorma Takasen säätiö for the in-
ternationalization support grant. The authors would like to thank
Bertrand Decharme, Christine Delire, and Bertrand Bonan at CN-
RM/GAME and Marie Dumont at CNRM/CEN for their great sup-
port during this work. CNRM/CEN is part of Labex OSUG@2020.
Financial support. This research has been supported by the
Academy of Finland (ArcI Profi 4), the Schweizerischer Nation-
alfonds zur Förderung der Wissenschaftlichen Forschung (grant
no. P500PN_202741), Horizon 2020 (grant no. 101056921), RCF
(LS-HYDRO, 356138 and UPFORMET, grant no. 308511) and the
Academy of Finland (grant nos. 347348, 337549, and 337552).
Review statement. This paper was edited by Melody Sandells and
reviewed by two anonymous referees.
References
Aalto, J., Aalto, P., Keronen, P., Kolari, P., Rantala, P.,
Taipale, R., Kajos, M., Patokoski, J., Rinne, J., Ruuska-
nen, T., Leskinen, M., Laakso, H., Levula, J., Pohja, T.,
Siivola, E., Kulmala, M., and Ylivinkka, I.: SMEAR II
Hyytiälä forest meteorology, greenhouse gases, air quality
and soil, Fairdata, https://doi.org/10.23729/62f7ad2c-7fe0-4f66-
b0a4-8d57c80524ec, 2022.
Ala-Aho, P., Autio, A., Bhattacharjee, J., Isokangas, E., Kujala, K.,
Marttila, H., Menberu, M., Meriö, L. J., Postila, H., Rauhala,
A., Ronkanen, A. K., Rossi, P. M., Saari, M., Haghighi, A. T.,
and Klove, B.: What conditions favor the influence of season-
ally frozen ground on hydrological partitioning? A systematic
review, Environ. Res. Lett., 16, 4, https://doi.org/10.1088/1748-
9326/abe82c, 2021.
Alekseychik, P., Kolari, P., Rinne, J., Haapanala, S., Laakso,
H., Taipale, R., Matilainen, T., Salminen, T., Levula, J., and
Tuittila, E.: SMEAR II Siikaneva 1 wetland meteorology
and soil, Fairdata, https://doi.org/10.23729/7d205559-3ef9-4f34-
8e08-ea24316f50c8, 2022a.
Alekseychik, P., Peltola, O., Li, X., Aurela, M., Hatakka, J., Pihlatie,
M., Rinne, J., Haapanala, S., Laakso, H., Taipale, R., Matilainen,
T., Salminen, T., and Levula, J.: SMEAR II Siikaneva 1 wetland
eddy covariance, Fairdata, https://doi.org/10.23729/f6455f02-
905b-4bf7-a870-743bd3788bf6, 2022b.
Alekseychik, P. K., Korrensalo, A., Mammarella, I., Vesala, T.,
and Tuittila, E. S.: Relationship between aerodynamic rough-
ness length and bulk sedge leaf area index in a mixed-species
boreal mire complex, Geophys. Res. Lett., 44, 5836–5843,
https://doi.org/10.1002/2017GL073884, 2017.
Andreas, E. L., Horst, T. W., Grachev, A. A., Persson, P.
O. G., Fairall, C. W., Guest, P. S., and Jordan, R. E.:
Parametrizing turbulent exchange over summer sea ice and the
marginal ice zone, Q. J. Roy. Meteor. Soc., 136, 927–943,
https://doi.org/10.1002/qj.618, 2010.
The Cryosphere, 18, 231–263, 2024 https://doi.org/10.5194/tc-18-231-2024
J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget 257
Aubinet, M., Vesala, T., and Papale, D.: Eddy covariance: a practi-
cal guide to measurement and data analysis, Springer Science &
Business Media, 2012.
Aurela, M., Riutta, T., Laurila, T., Tuovinen, J. P., Vesala, T., Tuit-
tila, E. S., Rinne, J., Haapanala, S., and Laine, J.: CO2 exchange
of a sedge fen in southern Finland – The impact of a drought
period, Tellus B, 59, 826–837, https://doi.org/10.1111/j.1600-
0889.2007.00309.x, 2007.
Aurela, M., Lohila, A., Tuovinen, J. P., Hatakka, J., Penttilä, T., and
Laurila, T.: Carbon dioxide and energy flux measurements in four
northern-boreal ecosystems at Pallas, Boreal Environ. Res., 20,
455–473, 2015.
Barrere, M., Domine, F., Decharme, B., Morin, S., Vionnet, V.,
and Lafaysse, M.: Evaluating the performance of coupled snow–
soil models in SURFEXv8 to simulate the permafrost thermal
regime at a high Arctic site, Geosci. Model Dev., 10, 3461–3479,
https://doi.org/10.5194/gmd-10-3461-2017, 2017.
Bengtsson, L., Andrae, U., Aspelien, T., Batrak, Y., Calvo, J.,
Rooy, W. d., Gleeson, E., Hansen-Sass, B., Homleid, M., Hor-
tal, M., Ivarsson, K. I., Lenderink, G., Niemelä, S., Nielsen,
K. P., Onvlee, J., Rontu, L., Samuelsson, P., Muñoz, D. S.,
Subias, A., Tijm, S., Toll, V., Yang, X., and Køltzow, M. Ã.: The
HARMONIE-AROME model configuration in the ALADIN-
HIRLAM NWP system, Mon. Weather Rev., 145, 1919–1935,
https://doi.org/10.1175/MWR-D-16-0417.1, 2017.
Beringer, J., Lynch, A. H., Chapin, F. S., Mack, M., and
Bonan, G. B.: The representation of arctic soils in
the land surface model: The importance of mosses, J.
Climate, 14, 3324–3335, https://doi.org/10.1175/1520-
0442(2001)014<3324:TROASI>2.0.CO;2, 2001.
Blyth, E. M., Arora, V. K., Clark, D. B., Dadson, S. J., De Kauwe,
M. G., Lawrence, D. M., Melton, J. R., Pongratz, J., Turton,
R. H., Yoshimura, K., and Yuan, H.: Advances in Land Sur-
face Modelling, Current Climate Change Reports, 7, 45–71,
https://doi.org/10.1007/s40641-021-00171-5, 2021.
Boisvert, L. N. and Stroeve, J. C.: The Arctic is becom-
ing warmer and wetter as revealed by the Atmospheric
Infrared Sounder, Geophys. Res. Lett., 42, 4439–4446,
https://doi.org/10.1002/2015GL063775, 2015.
Bonan, G. B., Patton, E. G., Finnigan, J. J., Baldocchi, D. D., and
Harman, I. N.: Moving beyond the incorrect but useful paradigm:
reevaluating big-leaf and multilayer plant canopies to model
biosphere-atmosphere fluxes – a review, Agr. Forest Meteorol.,
306, 108435, https://doi.org/10.1016/j.agrformet.2021.108435,
2021.
Boone, A., Masson, V., Meyers, T., and Noilhan, J.: The
influence of the inclusion of soil freezing on simula-
tions by a soil-vegetation-atmosphere transfer scheme, J.
Appl. Meteorol., 39, 1544–1569, https://doi.org/10.1175/1520-
0450(2000)039<1544:TIOTIO>2.0.CO;2, 2000.
Boone, A., Samuelsson, P., Gollvik, S., Napoly, A., Jarlan, L., Brun,
E., and Decharme, B.: The interactions between soil–biosphere–
atmosphere land surface model with a multi-energy balance
(ISBA-MEB) option in SURFEXv8 – Part 1: Model description,
Geosci. Model Dev., 10, 843–872, https://doi.org/10.5194/gmd-
10-843-2017, 2017.
Bouchard, B., Nadeau, D. F., and Domine, F.: Compari-
son of snowpack structure in gaps and under the canopy
in a humid boreal forest, Hydrol. Process., 36, e14681,
https://doi.org/10.1002/hyp.14681, 2022.
Broxton, P. D., Harpold, A. A., Biederman, J. A., Troch, P. A.,
Molotch, N. P., and Brooks, P. D.: Quantifying the ef-
fects of vegetation structure on snow accumulation and ab-
lation in mixed-conifer forests, Ecohydrology, 8, 1073–1094,
https://doi.org/10.1002/eco.1565, 2015.
Brun, E., Vionnet, V., Boone, A., Decharme, B., Peings, Y., Valette,
R., Karbou, F., and Morin, S.: Simulation of northern Eurasian
local snow depth, mass, and density using a detailed snowpack
model and meteorological reanalyses, J. Hydrometeorol., 14,
203–219, https://doi.org/10.1175/JHM-D-12-012.1, 2013.
Brunet, Y.: Turbulent Flow in Plant Canopies: Historical Per-
spective and Overview, vol. 177, Springer Netherlands,
ISBN 0123456789, https://doi.org/10.1007/s10546-020-00560-
7, 2020.
Caillaud, C., Somot, S., Alias, A., Bernard-Bouissières, I., Fu-
mière, Q., Laurantin, O., Seity, Y., and Ducrocq, V.: Modelling
Mediterranean heavy precipitation events at climate scale: an
object-oriented evaluation of the CNRM-AROME convection-
permitting regional climate model, vol. 56, Springer Berlin
Heidelberg, ISBN 0123456789, https://doi.org/10.1007/s00382-
020-05558-y, 2021.
Calvet, J. C., Noilhan, J., Roujean, J. L., Bessemoulin, P.,
Cabelguenne, M., Olioso, A., and Wigneron, J. P.: An in-
teractive vegetation SVAT model tested against data from
six contrasting sites, Agr. Forest Meteorol., 92, 73–95,
https://doi.org/10.1016/S0168-1923(98)00091-4, 1998.
Carrer, D., Roujean, J. L., Lafont, S., Calvet, J. C., Boone,
A., Decharme, B., Delire, C., and Gastellu-Etchegorry, J. P.:
A canopy radiative transfer scheme with explicit FAPAR
for the interactive vegetation model ISBA-A-gs: Impact
on carbon fluxes, J. Geophys. Res.-Biogeo., 118, 888–903,
https://doi.org/10.1002/jgrg.20070, 2013.
Chadburn, S., Burke, E., Essery, R., Boike, J., Langer, M.,
Heikenfeld, M., Cox, P., and Friedlingstein, P.: An im-
proved representation of physical permafrost dynamics in the
JULES land-surface model, Geosci. Model Dev., 8, 1493–1508,
https://doi.org/10.5194/gmd-8-1493-2015, 2015.
Champeaux, J. L., Masson, V., and Chauvin, F.: ECO-
CLIMAP: A global database of land surface param-
eters at 1 km resolution, Meteorol. Appl., 12, 29–32,
https://doi.org/10.1017/S1350482705001519, 2005.
Clark, M. P., Hendrikx, J., Slater, A. G., Kavetski, D., Anderson,
B., Cullen, N. J., Kerr, T., Örn Hreinsson, E., and Woods, R. A.:
Representing spatial variability of snow water equivalent in hy-
drologic and land-surface models: A review, Water Resour. Res.,
47, W07539, https://doi.org/10.1029/2011WR010745, 2011.
Clark, M. P., Nijssen, B., Lundquist, J. D., Kavetski, D., Rupp,
D. E., Woods, R. A., Freer, J. E., Gutmann, E. D., Wood, A. W.,
Brekke, L. D., Arnold, J. R., Gochis, D. J., and Rasmussen,
R. M.: A unified approach for process-based hydrologic mod-
eling: 1. Modeling concept, Water Resour. Res., 51, 2498–2514,
https://doi.org/10.1002/2015WR017200, 2015.
Cluzet, B. and Nousu, J.-P.: bertrandcz/CrocO_toolbox:
CrocO_v1.2 (CrocO_v1.2), Zenodo [code],
https://doi.org/10.5281/zenodo.10474747, 2024.
https://doi.org/10.5194/tc-18-231-2024 The Cryosphere, 18, 231–263, 2024
258 J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget
Cohen, J. and Rind, D.: The Effect of Snow Cover on the
Climate, J. Climate, 4, 689–706, https://doi.org/10.1175/1520-
0442(1991)004<0689:TEOSCO>2.0.CO;2, 1991.
Conway, J. P., Pomeroy, J. W., Helgason, W. D., and Kinar,
N. J.: Challenges in modeling turbulent heat fluxes to snow-
packs in forest clearings, J. Hydrometeorol., 19, 1599–1616,
https://doi.org/10.1175/JHM-D-18-0050.1, 2018.
Couet, J., Marjakangas, E. L., Santangeli, A., Kålås, J. A., Lind-
ström, Ã., and Lehikoinen, A.: Short-lived species move up-
hill faster under climate change, Oecologia, 198, 877–888,
https://doi.org/10.1007/s00442-021-05094-4, 2022.
Courtier, P., Freydier, C., Geleyn, J.-F., Rabier, F., and Rochas, M.:
The Arpege project at Météo-France, in: Seminar on Numerical
Methods in Atmospheric Models, 9–13 September 1991, Vol. II,
ECMWF, 1991.
Dankers, R., Burke, E. J., and Price, J.: Simulation of permafrost
and seasonal thaw depth in the JULES land surface scheme, The
Cryosphere, 5, 773–790, https://doi.org/10.5194/tc-5-773-2011,
2011.
Decharme, B., Boone, A., Delire, C., and Noilhan, J.: Lo-
cal evaluation of the Interaction between Soil Biosphere At-
mosphere soil multilayer diffusion scheme using four pe-
dotransfer functions, J. Geophys. Res.-Atmos., 116, 1–29,
https://doi.org/10.1029/2011JD016002, 2011.
Decharme, B., Brun, E., Boone, A., Delire, C., Le Moigne, P.,
and Morin, S.: Impacts of snow and organic soils parameteri-
zation on northern Eurasian soil temperature profiles simulated
by the ISBA land surface model, The Cryosphere, 10, 853–877,
https://doi.org/10.5194/tc-10-853-2016, 2016.
Decharme, B., Delire, C., Minvielle, M., Colin, J., Vergnes, J. P.,
Alias, A., Saint-Martin, D., Séférian, R., Sénési, S., and Voldoire,
A.: Recent Changes in the ISBA-CTRIP Land Surface System
for Use in the CNRM-CM6 Climate Model and in Global Off-
Line Hydrological Applications, J. Adv. Model. Earth Sy., 11,
1207–1252, https://doi.org/10.1029/2018MS001545, 2019.
Derbyshire, S. H.: Boundary-layer decoupling over cold surfaces as
a physical boundary-instability, Bound.-Lay. Meteorol., 90, 297–
325, https://doi.org/10.1023/A:1001710014316, 1999.
Deschamps-Berger, C., Cluzet, B., Dumont, M., Lafaysse, M.,
Berthier, E., Fanise, P., and Gascoin, S.: Improving the Spa-
tial Distribution of Snow Cover Simulations by Assimila-
tion of Satellite Stereoscopic Imagery, Water Resour. Res.,
58, e2021WR030271, https://doi.org/10.1029/2021WR030271,
2022.
Domine, F., Belke-Brea, M., Sarrazin, D., Arnaud, L., Barrere,
M., and Poirier, M.: Soil moisture, wind speed and depth hoar
formation in the Arctic snowpack, J. Glaciol., 64, 990–1002,
https://doi.org/10.1017/jog.2018.89, 2018.
Douville, H., Royer, J. F., and Mahfouf, J. F.: A new snow param-
eterization for the Météo-France climate model, Clim. Dynam.,
12, 21–35, https://doi.org/10.1007/BF00208760, 1995.
Esri: ESRI Satellite (ArcGIS/World_Imagery), https://www.arcgis.
com/home/item.html?id=10df2279f9684e4a9f6a7f08febac2a9
(last access: 19 September 2023), 2023.
Essery, R., Rutter, N., Pomeroy, J., Baxter, R., Stahli, M., Gustafs-
son, D., Barr, A., Bartlett, P., and Elder, K.: SnowMIP2: An eva-
lution of forest snow process simulation, B. Am. Meteor. Soc.,
90, 1130–1135, https://doi.org/10.1175/2009BAMS2629.1,
2009.
Eugster, W., Rouse, W. R., Pielke, R. A., Mcfadden, J. P., Bal-
docchi, D. D., Kittel, T. G., Chapin, F. S., Liston, G. E., Vi-
dale, P. L., Vaganov, E., and Chambers, S.: Land-atmosphere
energy exchange in Arctic tundra and boreal forest: Available
data and feedbacks to climate, Glob. Change Biol., 6, 84–115,
https://doi.org/10.1046/j.1365-2486.2000.06015.x, 2000.
Faroux, S., Kaptué Tchuenté, A. T., Roujean, J.-L., Masson, V.,
Martin, E., and Le Moigne, P.: ECOCLIMAP-II/Europe: a
twofold database of ecosystems and surface parameters at 1 km
resolution based on satellite information for use in land surface,
meteorological and climate models, Geosci. Model Dev., 6, 563–
582, https://doi.org/10.5194/gmd-6-563-2013, 2013.
FMI: Finnish Meteorological Institute past weather observa-
tions, https://en.ilmatieteenlaitos.fi/download-observations (last
access: 19 September 2023), 2021.
Foken, T., Göockede, M., Mauder, M., Mahrt, L., Amiro, B., and
Munger, W.: Post-Field Data Quality Control, in: Handbook of
Micrometeorology: A Guide for Surface Flux Measurement and
Analysis, edited by: Lee, X., Massman, W., and Law, B., 181–
208, Springer Netherlands, Dordrecht, ISBN 978-1-4020-2265-
4, https://doi.org/10.1007/1-4020-2265-4_9, 2005.
Gouttevin, I., Lehning, M., Jonas, T., Gustafsson, D., and Mölder,
M.: A two-layer canopy model with thermal inertia for an im-
proved snowpack energy balance below needleleaf forest (model
SNOWPACK, version 3.2.1, revision 741), Geosci. Model Dev.,
8, 2379–2398, https://doi.org/10.5194/gmd-8-2379-2015, 2015.
Gouttevin, I., Vionnet, V., Seity, Y., Boone, A., Lafaysse, M., Deliot,
Y., and Merzisen, H.: To the Origin of a Wintertime Screen-Level
Temperature Bias at High Altitude in a Kilometric NWP Model,
J. Hydrometeorol., 24, 53–71, https://doi.org/10.1175/JHM-D-
21-0200.1, 2023.
Hamdi, R., Degrauwe, D., Duerinckx, A., Cedilnik, J., Costa, V.,
Dalkilic, T., Essaouini, K., Jerczynki, M., Kocaman, F., Kull-
mann, L., Mahfouf, J.-F., Meier, F., Sassi, M., Schneider, S.,
Vánˇa, F., and Termonia, P.: Evaluating the performance of
SURFEXv5 as a new land surface scheme for the ALAD-
INcy36 and ALARO-0 models, Geosci. Model Dev., 7, 23–39,
https://doi.org/10.5194/gmd-7-23-2014, 2014.
Hari, P., Nikinmaa, E., Pohja, T., Siivola, E., Bäck, J., Vesala, T.,
and Kulmala, M.: Station for measuring ecosystem-atmosphere
relations: SMEAR, Physical and Physiological Forest Ecol-
ogy, 9789400756, 471–487, https://doi.org/10.1007/978-94-007-
5603-8_9, 2013.
Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A.,
Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers,
D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo,
G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara,
G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flem-
ming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L.,
Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S.,
Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P.,
Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J. N.: The
ERA5 global reanalysis, Q. J. Roy. Meteor. Soc., 146, 1999–
2049, https://doi.org/10.1002/qj.3803, 2020.
Jordan, R. E., Andreas, E. L., and Makshtas, A. P.: Heat budget of
snow-covered sea ice at North Pole 4, J. Geophys. Res.-Oceans,
104, 7785–7806, https://doi.org/10.1029/1999jc900011, 1999.
Koivusalo, H. and Heikinheimo, M.: Surface energy
exchange over a boreal snowpack: Comparison of
The Cryosphere, 18, 231–263, 2024 https://doi.org/10.5194/tc-18-231-2024
J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget 259
two snow energy balance models, Hydrol. Process.,
13, 2395–2408, https://doi.org/10.1002/(SICI)1099-
1085(199910)13:14/15<2395::AID-HYP864>3.0.CO;2-G,
1999.
Kolari, P., Aalto, J., Levula, J., Kulmala, L., Ilvesniemi, H., and
Pumpanen, J.: Hyytiälä SMEAR II site characteristics, Zenodo
[data set], https://doi.org/10.5281/zenodo.5909681, 2022.
Kreyling, J., Haei, M., and Laudon, H.: Absence of snow
cover reduces understory plant cover and alters plant commu-
nity composition in boreal forests, Oecologia, 168, 577–587,
https://doi.org/10.1007/s00442-011-2092-z, 2012.
Krinner, G., Derksen, C., Essery, R., Flanner, M., Hagemann, S.,
Clark, M., Hall, A., Rott, H., Brutel-Vuilmet, C., Kim, H., Mé-
nard, C. B., Mudryk, L., Thackeray, C., Wang, L., Arduini, G.,
Balsamo, G., Bartlett, P., Boike, J., Boone, A., Chéruy, F., Colin,
J., Cuntz, M., Dai, Y., Decharme, B., Derry, J., Ducharne, A.,
Dutra, E., Fang, X., Fierz, C., Ghattas, J., Gusev, Y., Haverd, V.,
Kontu, A., Lafaysse, M., Law, R., Lawrence, D., Li, W., Marke,
T., Marks, D., Ménégoz, M., Nasonova, O., Nitta, T., Niwano,
M., Pomeroy, J., Raleigh, M. S., Schaedler, G., Semenov, V.,
Smirnova, T. G., Stacke, T., Strasser, U., Svenson, S., Turkov,
D., Wang, T., Wever, N., Yuan, H., Zhou, W., and Zhu, D.:
ESM-SnowMIP: assessing snow models and quantifying snow-
related climate feedbacks, Geosci. Model Dev., 11, 5027–5049,
https://doi.org/10.5194/gmd-11-5027-2018, 2018.
Lackner, G., Domine, F., Nadeau, D. F., Parent, A.-C., Anc-
til, F., Lafaysse, M., and Dumont, M.: On the energy bud-
get of a low-Arctic snowpack, The Cryosphere, 16, 127–142,
https://doi.org/10.5194/tc-16-127-2022, 2022.
Lafaysse, M., Hingray, B., Etchevers, P., Martin, E., and Obled,
C.: Influence of spatial discretization, underground water stor-
age and glacier melt on a physically-based hydrological model
of the Upper Durance River basin, J. Hydrol., 403, 116–129,
https://doi.org/10.1016/j.jhydrol.2011.03.046, 2011.
Lafaysse, M., Cluzet, B., Dumont, M., Lejeune, Y., Vionnet,
V., and Morin, S.: A multiphysical ensemble system of nu-
merical snow modelling, The Cryosphere, 11, 1173–1198,
https://doi.org/10.5194/tc-11-1173-2017, 2017.
Lafaysse, M., Fructus, M., Vernay, M., Radanovics, S., Dumont,
M., and Viallon-Galinier, L.: Procedure for new users of Crocus
model, https://opensource.umr-cnrm.fr/projects/snowtools_git/
wiki/Procedure_for_new_users, (last access: 26 January 2023),
2023.
Langer, M., Westermann, S., Muster, S., Piel, K., and Boike, J.:
The surface energy balance of a polygonal tundra site in north-
ern Siberia – Part 2: Winter, The Cryosphere, 5, 509–524,
https://doi.org/10.5194/tc-5-509-2011, 2011.
Lapo, K., Nijssen, B., and Lundquist, J. D.: Evaluation of
Turbulence Stability Schemes of Land Models for Sta-
ble Conditions, J. Geophys. Res.-Atmos., 124, 3072–3089,
https://doi.org/10.1029/2018JD028970, 2019.
Launiainen, S.: Seasonal and inter-annual variability of energy ex-
change above a boreal Scots pine forest, Biogeosciences, 7,
3921–3940, https://doi.org/10.5194/bg-7-3921-2010, 2010.
Launiainen, S., Katul, G. G., Lauren, A., and Kolari, P.:
Coupling boreal forest CO2, H2O and energy flows by
a vertically structured forest canopy – Soil model with
separate bryophyte layer, Ecol. Model., 312, 385–405,
https://doi.org/10.1016/j.ecolmodel.2015.06.007, 2015.
Launiainen, S., Katul, G. G., Leppä, K., Kolari, P., Aslan, T., Grön-
holm, T., Korhonen, L., Mammarella, I., and Vesala, T.: Does
growing atmospheric CO2 explain increasing carbon sink in a
boreal coniferous forest?, Glob. Change Biol., 28, 2910–2929,
https://doi.org/10.1111/gcb.16117, 2022.
Lawrence, D. M. and Slater, A. G.: Incorporating organic soil
into a global climate model, Clim. Dynam., 30, 145–160,
https://doi.org/10.1007/s00382-007-0278-1, 2008.
Lawrence, D. M., Fisher, R. A., Koven, C. D., Oleson, K. W.,
Swenson, S. C., Bonan, G., Collier, N., Ghimire, B., van Kam-
penhout, L., Kennedy, D., Kluzek, E., Lawrence, P. J., Li, F.,
Li, H., Lombardozzi, D., Riley, W. J., Sacks, W. J., Shi, M.,
Vertenstein, M., Wieder, W. R., Xu, C., Ali, A. A., Badger,
A. M., Bisht, G., van den Broeke, M., Brunke, M. A., Burns,
S. P., Buzan, J., Clark, M., Craig, A., Dahlin, K., Drewniak,
B., Fisher, J. B., Flanner, M., Fox, A. M., Gentine, P., Hoff-
man, F., Keppel-Aleks, G., Knox, R., Kumar, S., Lenaerts, J.,
Leung, L. R., Lipscomb, W. H., Lu, Y., Pandey, A., Pelletier,
J. D., Perket, J., Randerson, J. T., Ricciuto, D. M., Sanderson,
B. M., Slater, A., Subin, Z. M., Tang, J., Thomas, R. Q., Val Mar-
tin, M., and Zeng, X.: The Community Land Model Version
5: Description of New Features, Benchmarking, and Impact of
Forcing Uncertainty, J. Adv. Model. Earth Sy., 11, 4245–4287,
https://doi.org/10.1029/2018MS001583, 2019.
Le Moigne, P.: SURFEX scientific documentation, SURFEX
v8.1, Issue 3, Météo-France, Toulouse, France, https://www.
umr-cnrm.fr/surfex/IMG/pdf/surfex_scidoc_v8.1.pdf (last ac-
cess: 26 January 2023), 2018.
Le Moigne, P., Besson, F., Martin, E., Boé, J., Boone, A., Decharme,
B., Etchevers, P., Faroux, S., Habets, F., Lafaysse, M., Ler-
oux, D., and Rousset-Regimbeau, F.: The latest improvements
with SURFEX v8.0 of the Safran–Isba–Modcou hydrometeoro-
logical model for France, Geosci. Model Dev., 13, 3925–3946,
https://doi.org/10.5194/gmd-13-3925-2020, 2020.
Lindroos, A. J., Mäkipää, R., and Merilä, P.: Soil carbon
stock changes over 21 years in intensively monitored bo-
real forest stands in Finland, Ecol. Indic., 144, 109551,
https://doi.org/10.1016/j.ecolind.2022.109551, 2022.
Liston, G. E. and Sturm, M.: Winter precipitation pat-
terns in arctic Alaska determined from a blowing-
snow model and snow-depth observations, J. Hy-
drometeorol., 3, 646–659, https://doi.org/10.1175/1525-
7541(2002)003<0646:WPPIAA>2.0.CO;2, 2002.
Lohila, A., Penttilä, T., Jortikka, S., Aalto, T., Anttila, P., Asmi, E.,
Aurela, M., Hatakka, J., Hellén, H., Henttonen, H., Hänninen, P.,
Kilkki, J., Kyllönen, K., Laurila, T., Lepistö, A., Lihavainen, H.,
Makkonen, U., Paatero, J., Rask, M., Sutinen, R., Tuovinen, J. P.,
Vuorenmaa, J., and Viisanen, Y.: Preface to the special issue on
integrated research of atmosphere, ecosystems and environment
at Pallas, Boreal Environ. Res., 20, 431–454, 2015.
Loranty, M. M., Berner, L. T., Goetz, S. J., Jin, Y., and Randerson,
J. T.: Vegetation controls on northern high latitude snow-albedo
feedback: Observations and CMIP5 model simulations, Glob.
Change Biol., 20, 594–606, https://doi.org/10.1111/gcb.12391,
2014.
Louis, J. F.: A parametric model of vertical eddy fluxes
in the atmosphere, Bound.-Lay. Meteorol., 17, 187–202,
https://doi.org/10.1007/BF00117978, 1979.
https://doi.org/10.5194/tc-18-231-2024 The Cryosphere, 18, 231–263, 2024
260 J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget
Lundquist, J. D., Dickerson-Lange, S. E., Lutz, J. A., and Cristea,
N. C.: Lower forest density enhances snow retention in regions
with warmer winters: A global framework developed from plot-
scale observations and modeling, Water Resour. Res., 49, 6356–
6370, https://doi.org/10.1002/wrcr.20504, 2013.
Lundquist, J. D., Dickerson-Lange, S., Gutmann, E., Jonas, T.,
Lumbrazo, C., and Reynolds, D.: Snow interception modelling:
Isolated observations have led to many land surface models lack-
ing appropriate temperature sensitivities, Hydrol. Process., 35,
1–20, https://doi.org/10.1002/hyp.14274, 2021.
Mahfouf, A. J., Manzi, A. O., Noilhan, J., Giordani, H., and Déqué,
M.: The Land Surface Scheme ISBA within the Météo-France
Climate Model ARPEGE. Part I: Implementation and Prelimi-
nary Results, American Meteorological Society, 8, 2039–2057,
1995.
Malle, J., Rutter, N., Webster, C., Mazzotti, G., Wake, L., and
Jonas, T.: Effect of Forest Canopy Structure on Wintertime
Land Surface Albedo: Evaluating CLM5 Simulations With
In-Situ Measurements, J. Geophys. Res.-Atmos., 126, 1–15,
https://doi.org/10.1029/2020JD034118, 2021.
Mammarella, I., Peltola, O., Nordbo, A., Järvi, L., and Rannik, Ü.:
Quantifying the uncertainty of eddy covariance fluxes due to the
use of different software packages and combinations of process-
ing steps in two contrasting ecosystems, Atmos. Meas. Tech., 9,
4915–4933, https://doi.org/10.5194/amt-9-4915-2016, 2016.
Mammarella, I., Rannik, Ü., Launiainen, S., Alekseychik, P.,
Peltola, O., Keronen, P., Kolari, P., Laakso, H., Mati-
lainen, T., Salminen, T., Levula, J., Pohja, T., Siivola, E.,
and Vesala, T.: SMEAR II Hyytiälä forest eddy covari-
ance, Fairdata, https://doi.org/10.23729/40f64739-11d1-4e5f-
8dc2-da931512c91c, 2019.
Martin, E. and Lejeune, Y.: Turbulent fluxes above the snow surface,
Ann. Glaciol., 26, 179–183, https://doi.org/10.3189/1998aog26-
1-179-183, 1998.
Marttila, H., Lohila, A., Ala-Aho, P., Noor, K., Welker, J. M.,
Croghan, D., Mustonen, K., Meriö, L., Autio, A., Muhic, F.,
Bailey, H., Aurela, M., Vuorenmaa, J., Penttilä, T., Hyöky, V.,
Klein, E., Kuzmin, A., Korpelainen, P., Kumpula, T., Rauhala,
A., and Kløve, B.: Subarctic catchment water storage and car-
bon cycling – Leading the way for future studies using inte-
grated datasets at Pallas, Finland, Hydrol. Process., 35, 1–19,
https://doi.org/10.1002/hyp.14350, 2021.
Masson, V., Le Moigne, P., Martin, E., Faroux, S., Alias, A.,
Alkama, R., Belamari, S., Barbu, A., Boone, A., Bouyssel, F.,
Brousseau, P., Brun, E., Calvet, J.-C., Carrer, D., Decharme, B.,
Delire, C., Donier, S., Essaouini, K., Gibelin, A.-L., Giordani, H.,
Habets, F., Jidane, M., Kerdraon, G., Kourzeneva, E., Lafaysse,
M., Lafont, S., Lebeaupin Brossier, C., Lemonsu, A., Mahfouf,
J.-F., Marguinaud, P., Mokhtari, M., Morin, S., Pigeon, G., Sal-
gado, R., Seity, Y., Taillefer, F., Tanguy, G., Tulet, P., Vincendon,
B., Vionnet, V., and Voldoire, A.: The SURFEXv7.2 land and
ocean surface platform for coupled or offline simulation of earth
surface variables and fluxes, Geosci. Model Dev., 6, 929–960,
https://doi.org/10.5194/gmd-6-929-2013, 2013.
Mauder, M., Foken, T., and Cuxart, J.: Surface-Energy-Balance
Closure over Land: A Review, vol. 177, Springer Netherlands,
ISBN 0123456789, https://doi.org/10.1007/s10546-020-00529-
6, 2020.
Mazzotti, G., Essery, R., Moeser, C. D., and Jonas, T.: Resolv-
ing Small-Scale Forest Snow Patterns Using an Energy Balance
Snow Model With a One-Layer Canopy, Water Resour. Res.,
56, e2019WR026129, https://doi.org/10.1029/2019WR026129,
2020a.
Mazzotti, G., Essery, R., Webster, C., Malle, J., and Jonas, T.:
Process-Level Evaluation of a Hyper-Resolution Forest Snow
Model Using Distributed Multisensor Observations, Water Re-
sour. Res., 56, 1–25, https://doi.org/10.1029/2020WR027572,
2020b.
Mazzotti, G., Webster, C., Essery, R., and Jonas, T.: Increas-
ing the Physical Representation of Forest-Snow Processes in
Coarse-Resolution Models: Lessons Learned From Upscaling
Hyper-Resolution Simulations, Water Resour. Res., 57, 1–21,
https://doi.org/10.1029/2020WR029064, 2021.
McGowan, L. E., Paw U, K. T., and Dahlke, H. E.: What Does
a Multilayer Canopy Model Tell Us About Our Current Un-
derstanding of Snow-Canopy Unloading?, in: AGU Fall Meet-
ing Abstracts, New Orleans, 11–15 December 2017, vol. 2017,
C43B–06, https://doi.org/10.13140/RG.2.2.33388.62085, 2017.
Ménard, C. B., Essery, R., Barr, A., Bartlett, P., Derry, J., Dumont,
M., Fierz, C., Kim, H., Kontu, A., Lejeune, Y., Marks, D., Ni-
wano, M., Raleigh, M., Wang, L., and Wever, N.: Meteorolog-
ical and evaluation datasets for snow modelling at 10 reference
sites: description of in situ and bias-corrected reanalysis data,
Earth Syst. Sci. Data, 11, 865–880, https://doi.org/10.5194/essd-
11-865-2019, 2019.
Menard, C. B., Essery, R., Krinner, G., Arduini, G., Bartlett, P.,
Boone, A., Brutel-Vuilmet, C., Burke, E., Cuntz, M., Dai, Y.,
Decharme, B., Dutra, E., Fang, X., Fierz, C., Gusev, Y., Hage-
mann, S., Haverd, V., Kim, H., Lafaysse, M., Marke, T., Na-
sonova, O., Nitta, T., Niwano, M., Pomeroy, J., Schädler, G., Se-
menov, V. A., Smirnova, T., Strasser, U., Swenson, S., Turkov,
D., Wever, N., and Yuan, H.: Scientific and human errors in a
snow model intercomparison, B. Am. Meteor. Soc., 102, E61–
E79, https://doi.org/10.1175/BAMS-D-19-0329.1, 2021.
Menberu, M. W., Marttila, H., Ronkanen, A. K., Haghighi, A. T.,
and Kløve, B.: Hydraulic and Physical Properties of Man-
aged and Intact Peatlands: Application of the Van Genuchten-
Mualem Models to Peat Soils, Water Resour. Res., 57, 1–22,
https://doi.org/10.1029/2020WR028624, 2021.
Meriö, L.-J., Rauhala, A., Ala-aho, P., Kuzmin, A., Korpelainen,
P., Kumpula, T., Kløve, B., and Marttila, H.: Measuring
the spatiotemporal variability in snow depth in subarctic
environments using UASs – Part 2: Snow processes and
snow–canopy interactions, The Cryosphere, 17, 4363–4380,
https://doi.org/10.5194/tc-17-4363-2023, 2023.
Molotch, N. P., Blanken, P. D., Williams, M. W., Turnipseed, A. A.,
Monson, R. K., and Margulis, S. A.: Estimating sublimation of
intercepted and sub-canopy snow using eddy covariance systems,
Hydrol. Process., 21, 1567–1575, 2009.
Morin, S., Horton, S., Techel, F., Bavay, M., Coléou, C., Fierz,
C., Gobiet, A., Hagenmuller, P., Lafaysse, M., Ližar, M., Mit-
terer, C., Monti, F., Müller, K., Olefs, M., Snook, J. S.,
van Herwijnen, A., and Vionnet, V.: Application of physical
snowpack models in support of operational avalanche hazard
forecasting: A status report on current implementations and
prospects for the future, Cold Reg. Sci. Technol., 170, 102910,
https://doi.org/10.1016/j.coldregions.2019.102910, 2020.
The Cryosphere, 18, 231–263, 2024 https://doi.org/10.5194/tc-18-231-2024
J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget 261
Morris, P. J., Davies, M. L., Baird, A. J., Balliston, N., Bour-
gault, M., Clymo, R. S., Fewster, R. E., Furukawa, A. K.,
Holden, J., Kessel, E., Ketcheson, S. J., Kløve, B., Larocque,
M., Marttila, H., Menberu, M. W., Moore, P. A., Price, J. S.,
Ronkanen, A., Rosa, E., Strack, M., Surridge, B. W. J.,
Waddington, J. M., Whittington, P., and Wilkinson, S. L.: Sat-
urated Hydraulic Conductivity in Northern Peats Inferred From
Other Measurements, Water Resour. Res., 58, e2022WR033181,
https://doi.org/10.1029/2022wr033181, 2022.
Muhic, F., Ala-Aho, P., Noor, K., Welker, J. M., Klöve, B., and
Marttila, H.: Flushing or mixing? Stable water isotopes reveal
differences in arctic forest and peatland soil water seasonality,
Hydrol. Process., 37, 1–22, https://doi.org/10.1002/hyp.14811,
2023.
Mustamo, P., Ronkanen, A. K., Berglund, Ã., Berglund, K.,
and Kløve, B.: Thermal conductivity of unfrozen and par-
tially frozen managed peat soils, Soil Till. Res., 191, 245–255,
https://doi.org/10.1016/j.still.2019.02.017, 2019.
Nabat, P., Somot, S., Cassou, C., Mallet, M., Michou, M.,
Bouniol, D., Decharme, B., Drugé, T., Roehrig, R., and Saint-
Martin, D.: Modulation of radiative aerosols effects by atmo-
spheric circulation over the Euro-Mediterranean region, Atmos.
Chem. Phys., 20, 8315–8349, https://doi.org/10.5194/acp-20-
8315-2020, 2020.
Napoly, A., Boone, A., Samuelsson, P., Gollvik, S., Martin, E.,
Seferian, R., Carrer, D., Decharme, B., and Jarlan, L.: The in-
teractions between soil–biosphere–atmosphere (ISBA) land sur-
face model multi-energy balance (MEB) option in SURFEXv8 –
Part 2: Introduction of a litter formulation and model evaluation
for local-scale forest sites, Geosci. Model Dev., 10, 1621–1644,
https://doi.org/10.5194/gmd-10-1621-2017, 2017.
Napoly, A., Boone, A., and Welfringer, T.: ISBA-MEB (SURFEX
v8.1): model snow evaluation for local-scale forest sites, Geosci.
Model Dev., 13, 6523–6545, https://doi.org/10.5194/gmd-13-
6523-2020, 2020.
Nicolsky, D. J., Romanovsky, V. E., Alexeev, V. A., and Lawrence,
D. M.: Improved modeling of permafrost dynamics in a
GCM land-surface scheme, Geophys. Res. Lett., 34, 2–6,
https://doi.org/10.1029/2007GL029525, 2007.
Niu, G. Y., Yang, Z. L., Mitchell, K. E., Chen, F., Ek, M. B., Barlage,
M., Kumar, A., Manning, K., Niyogi, D., Rosero, E., Tewari,
M., and Xia, Y.: The community Noah land surface model with
multiparameterization options (Noah-MP): 1. Model descrip-
tion and evaluation with local-scale measurements, J. Geophys.
Res.-Atmos., 116, 1–19, https://doi.org/10.1029/2010JD015139,
2011.
NLSF: National Land Survey of Finland Topographic
Database, https://asiointi.maanmittauslaitos.fi/karttapaikka/
tiedostopalvelu?lang=en (last access: 29 July 2023), 2020.
Noilhan, J. and Mahfouf, J. F.: The ISBA land surface pa-
rameterisation scheme, Global Planet. Change, 13, 145–159,
https://doi.org/10.1016/0921-8181(95)00043-7, 1996.
Noilhan, J. and Planton, S.: A simple parameterization of
land surface processes for meteorological models, Mon.
Weather Rev., 117, 536–549, https://doi.org/10.1175/1520-
0493(1989)117<0536:ASPOLS>2.0.CO;2, 1989.
Noor, K., Marttila, H., Klöve, B., Welker, J. M., and Ala-
aho, P.: The Spatiotemporal Variability of Snowpack and
Snowmelt Water O and H Isotopes in a Subarctic Catch-
ment Water Resources Research, Water Resour. Res., 59, 1–19,
https://doi.org/10.1029/2022WR033101, 2022.
Nousu, J.-P., Lafaysse, M., Vernay, M., Bellier, J., Evin, G., and
Joly, B.: Statistical post-processing of ensemble forecasts of the
height of new snow, Nonlin. Processes Geophys., 26, 339–357,
https://doi.org/10.5194/npg-26-339-2019, 2019.
Nousu, J.-P., Lafaysse, M., Mazzotti, G., Ala-aho, P., Marttila, H.,
Cluzet, B., Aurela, M., Lohila, A., Kolari, P., Boone, A., Fructus,
M., and Launiainen, S.: Modelling snowpack dynamics and sur-
face energy budget in boreal and subarctic peatlands and forests,
Zenodo [data set], https://doi.org/10.5281/zenodo.8252267,
2023.
Nousu, J.-P., Lafaysse, M., Mazzotti, G., Ala-aho, P., Marttila,
H., Cluzet, B., Aurela, M., Lohila, A., Kolari, P., Boone, A.,
Fructus, M., and Launiainen, S.: SURFEX software supplement
to ”Modeling snowpack dynamics and surface energy budget
in boreal and subarctic peatlands and forests”, Zenodo [code],
https://doi.org/10.5281/zenodo.10473209, 2024.
Olson, D. M., Dinerstein, E., Wikramanayake, E. D., Burgess,
N. D., Powell, G. V., Underwood, E. C., D’Amico, J. A.,
Itoua, I., Strand, H. E., Morrison, J. C., Loucks, C. J.,
Allnutt, T. F., Ricketts, T. H., Kura, Y., Lamoreux, J. F.,
Wettengel, W. W., Hedao, P., and Kassem, K. R.: Ter-
restrial ecoregions of the world: A new map of life on
Earth, BioScience, 51, 933–938, https://doi.org/10.1641/0006-
3568(2001)051[0933:TEOTWA]2.0.CO;2, 2001.
Park, H., Launiainen, S., Konstantinov, P. Y., Iijima, Y., and Fe-
dorov, A. N.: Modeling the Effect of Moss Cover on Soil
Temperature and Carbon Fluxes at a Tundra Site in North-
eastern Siberia, J. Geophys. Res.-Biogeo., 123, 3028–3044,
https://doi.org/10.1029/2018JG004491, 2018.
Pedersen, S. H., Bentzen, T. W., Reinking, A. K., Liston, G. E., El-
der, K., Lenart, E. A., Prichard, A. K., and Welker, J. M.: Quan-
tifying effects of snow depth on caribou winter range selection
and movement in Arctic Alaska, Movement Ecology, 9, 1–24,
https://doi.org/10.1186/s40462-021-00276-4, 2021.
Peters-Lidard, C. D., Blackburn, E., Liang, X., and Wood,
E. F.: The effect of soil thermal conductivity parame-
terization on surface energy fluxes and temperatures, J.
Atmos. Sci., 55, 1209–1224, https://doi.org/10.1175/1520-
0469(1998)055<1209:TEOSTC>2.0.CO;2, 1998.
Pomeroy, J. W. and Dion, K.: Winter radiation ex-
tinction and reflection in a boreal pine canopy:
Measurements and modelling, Hydrol. Process.,
10, 1591–1608, https://doi.org/10.1002/(sici)1099-
1085(199612)10:12<1591::aid-hyp503>3.0.co;2-8, 1996.
Pomeroy, J. W. and Essery, R. L.: Turbulent fluxes during blow-
ing snow: Field tests of model sublimation predictions, Hydrol.
Process., 13, 2963–2975, https://doi.org/10.1002/(SICI)1099-
1085(19991230)13:18<2963::AID-HYP11>3.0.CO;2-9, 1999.
Raleigh, M. S., Lundquist, J. D., and Clark, M. P.: Explor-
ing the impact of forcing error characteristics on physi-
cally based snow simulations within a global sensitivity anal-
ysis framework, Hydrol. Earth Syst. Sci., 19, 3153–3179,
https://doi.org/10.5194/hess-19-3153-2015, 2015.
Rantanen, M., Karpechko, A. Y., Lipponen, A., Nordling, K.,
Hyvärinen, O., Ruosteenoja, K., Vihma, T., and Laaksonen,
A.: The Arctic has warmed nearly four times faster than
https://doi.org/10.5194/tc-18-231-2024 The Cryosphere, 18, 231–263, 2024
262 J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget
the globe since 1979, Commun. Earth Environ., 3, 1–10,
https://doi.org/10.1038/s43247-022-00498-3, 2022.
Rasmus, S., Lundell, R., and Saarinen, T.: Interac-
tions between snow, canopy, and vegetation in a bo-
real coniferous forest, Plant Ecol. Divers., 4, 55–65,
https://doi.org/10.1080/17550874.2011.558126, 2011.
Reba, M. L., Link, T. E., Marks, D., and Pomeroy, J.: An assessment
of corrections for eddy covariance measured turbulent fluxes
over snow in mountain environments, Water Resour. Res., 46,
W00D38, https://doi.org/10.1029/2008WR007045, 2009.
Revuelto, J., Lecourt, G., Lafaysse, M., Zin, I., Charrois, L.,
Vionnet, V., Dumont, M., Rabatel, A., Six, D., Condom, T.,
Morin, S., Viani, A., and Sirguey, P.: Multi-criteria evalua-
tion of snowpack simulations in complex alpine terrain using
satellite and in situ observations, Remote Sensing, 10, 1–32,
https://doi.org/10.3390/rs10081171, 2018.
Ricchiazzi, P., Yang, S., Gautier, C., and Sowle, D.: SB-
DART: A Research and Teaching Software Tool for Plane-
Parallel Radiative Transfer in the Earth’s Atmosphere, B.
Am. Meteor. Soc., 79, 2101–2114, https://doi.org/10.1175/1520-
0477(1998)079<2101:SARATS>2.0.CO;2, 1998.
Richter, B., Schweizer, J., Rotach, M. W., and Van Herwij-
nen, A.: Modeling spatially distributed snow instability at a
regional scale using Alpine3D, J. Glaciol., 67, 1147–1162,
https://doi.org/10.1017/jog.2021.61, 2021.
Rinne, J., Tuittila, E. S., Peltola, O., Li, X., Raivonen, M., Aleksey-
chik, P., Haapanala, S., Pihlatie, M., Aurela, M., Mammarella, I.,
and Vesala, T.: Temporal Variation of Ecosystem Scale Methane
Emission From a Boreal Fen in Relation to Temperature, Water
Table Position, and Carbon Dioxide Fluxes, Global Biogeochem.
Cy., 32, 1087–1106, https://doi.org/10.1029/2017GB005747,
2018.
Rissanen, T., Niittynen, P., Soininen, J., and Luoto, M.: Snow infor-
mation is required in subcontinental scale predictions of moun-
tain plant distributions, Global Ecol. Biogeogr., 30, 1502–1513,
https://doi.org/10.1111/geb.13315, 2021.
Rutter, N., Essery, R., Pomeroy, J., Altimir, N., Andreadis, K.,
Baker, I., Barr, A., Bartlett, P., Boone, A., Deng, H., Dou-
ville, H., Dutra, E., Elder, K., Ellis, C., Feng, X., Gelfan,
A., Goodbody, A., Gusev, Y., Gustafsson, D., Hellström, R.,
Hirabayashi, Y., Hirota, T., Jonas, T., Koren, V., Kuragina, A.,
Lettenmaier, D., Li, W. P., Luce, C., Martin, E., Nasonova,
O., Pumpanen, J., Pyles, R. D., Samuelsson, P., Sandells, M.,
Schädler, G., Shmakin, A., Smirnova, T. G., Stähli, M., Stöckli,
R., Strasser, U., Su, H., Suzuki, K., Takata, K., Tanaka, K.,
Thompson, E., Vesala, T., Viterbo, P., Wiltshire, A., Xia, K.,
Xue, Y., and Yamazaki, T.: Evaluation of forest snow processes
models (SnowMIP2), J. Geophys. Res.-Atmos., 114, D06111,
https://doi.org/10.1029/2008JD011063, 2009.
Salas-Mélia, D., Chauvin, F., Déqué, M., Douville, H., Gueremy,
J., Marquet, P., Planton, S., Royer, J., and Tyteca, S.: Description
and Validation of the CNRM-CM3 Global Coupled Model, Clim.
Dynam., 103, 2005.
Screen, J. A. and Simmonds, I.: The central role of diminishing
sea ice in recent Arctic temperature amplification, Nature, 464,
1334–1337, https://doi.org/10.1038/nature09051, 2010.
Serreze, M. C., Barrett, A. P., Stroeve, J. C., Kindig, D. N., and
Holland, M. M.: The emergence of surface-based Arctic ampli-
fication, The Cryosphere, 3, 11–19, https://doi.org/10.5194/tc-3-
11-2009, 2009.
Slater, A. G., Schlosser, C. A., Desborough, C. E., Pitman, A. J.,
Henderson-Sellers, A., Robock, A., Vinnikov, K. Y., Mitchell,
K., Boone, A., Braden, H., Chen, F., Cox, P. M., De Ros-
nay, P., Dickinson, R. E., Dai, Y. J., Duan, Q., Entin, J.,
Etchevers, P., Gedney, N., Gusev, Y. M., Habets, F., Kim, J.,
Koren, V., Kowalczyk, E. A., Nasonova, O. N., Noilhan, J.,
Schaake, S., Shmakin, A. B., Smirnova, T. G., Verseghy, D.,
Wetzel, P., Xue, Y., Yang, Z. L., and Zeng, Q.: The represen-
tation of snow in land surface schemes: Results from PILPS
2(d), J. Hydrometeorol., 2, 7–25, https://doi.org/10.1175/1525-
7541(2001)002<0007:TROSIL>2.0.CO;2, 2001.
Stiegler, C., Lund, M., Christensen, T. R., Mastepanov, M., and
Lindroth, A.: Two years with extreme and little snowfall: ef-
fects on energy partitioning and surface energy exchange in a
high-Arctic tundra ecosystem, The Cryosphere, 10, 1395–1413,
https://doi.org/10.5194/tc-10-1395-2016, 2016.
Stigter, E. E., Steiner, J. F., Koch, I., Saloranta, T. M., Kirkham,
J. D., and Immerzeel, W. W.: Energy and mass balance dy-
namics of the seasonal snowpack at two high-altitude sites
in the Himalaya, Cold Reg. Sci. Technol., 183, 103233,
https://doi.org/10.1016/j.coldregions.2021.103233, 2021.
Stuefer, S. L., Kane, D. L., and Dean, K. M.: Snow Water Equiva-
lent Measurements in Remote Arctic Alaska Watersheds, Water
Resour. Res., 56, 1–12, https://doi.org/10.1029/2019WR025621,
2020.
Tayleur, C. M., Devictor, V., Gaüzère, P., Jonzén, N., Smith, H. G.,
and Lindström, Ã.: Regional variation in climate change win-
ners and losers highlights the rapid loss of cold-dwelling species,
Divers. Distrib., 22, 468–480, https://doi.org/10.1111/ddi.12412,
2016.
Thackeray, C. W., Derksen, C., Fletcher, C. G., and Hall,
A.: Snow and Climate: Feedbacks, Drivers, and Indices
of Change, Current Climate Change Reports, 5, 322–333,
https://doi.org/10.1007/s40641-019-00143-w, 2019.
Tribbeck, M. J., Gurney, R. J., and Morris, E. M.: The radiative
effect of a Fir Canopy on a snowpack, J. Hydrometeorol., 7, 880–
895, https://doi.org/10.1175/JHM528.1, 2006.
Tuzet, F., Dumont, M., Picard, G., Lamare, M., Voisin, D., Nabat,
P., Lafaysse, M., Larue, F., Revuelto, J., and Arnaud, L.:
Quantification of the radiative impact of light-absorbing par-
ticles during two contrasted snow seasons at Col du Lautaret
(2058 m a.s.l., French Alps), The Cryosphere, 14, 4553–4579,
https://doi.org/10.5194/tc-14-4553-2020, 2020.
Tyler, N. J.: Climate, snow, ice, crashes, and declines in populations
of reindeer and caribou (Rangifer tarandus L.), Ecol.Monogr., 80,
197–219, https://doi.org/10.1890/09-1070.1, 2010.
Väliranta, M. and Mathijssen, P. J. H.: Geochemistry
of Siikaneva peat core from Finland, PANGAEA,
https://doi.org/10.1594/PANGAEA.927689, 2021.
Varhola, A., Coops, N. C., Weiler, M., and Moore, R. D.: For-
est canopy effects on snow accumulation and ablation: An in-
tegrative review of empirical results, J. Hydrol., 392, 219–233,
https://doi.org/10.1016/j.jhydrol.2010.08.009, 2010.
Vernay, M., Lafaysse, M., Monteiro, D., Hagenmuller, P., Nheili,
R., Samacoïts, R., Verfaillie, D., and Morin, S.: The S2M mete-
orological and snow cover reanalysis over the French mountain-
ous areas: description and evaluation (1958–2021), Earth Syst.
The Cryosphere, 18, 231–263, 2024 https://doi.org/10.5194/tc-18-231-2024
J.-P. Nousu et al.: Modeling snowpack dynamics and surface energy budget 263
Sci. Data, 14, 1707–1733, https://doi.org/10.5194/essd-14-1707-
2022, 2022.
Verrelle, A., Glinton, M., Bazile, E., and Le Moigne, P.: CERRA-
Land : A new land surface reanalysis at 5.5 km resolution over
Europe, EMS Annual Meeting 2021, online, 6–10 Sep 2021,
EMS2021-492, https://doi.org/10.5194/ems2021-492, 2021.
Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le
Moigne, P., Martin, E., and Willemet, J.-M.: The detailed snow-
pack scheme Crocus and its implementation in SURFEX v7.2,
Geosci. Model Dev., 5, 773–791, https://doi.org/10.5194/gmd-5-
773-2012, 2012.
Voldoire, A., Sanchez-Gomez, E., Salas y Mélia, D., Decharme,
B., Cassou, C., Sénési, S., Valcke, S., Beau, I., Alias, A.,
Chevallier, M., Déqué, M., Deshayes, J., Douville, H., Fernan-
dez, E., Madec, G., Maisonnave, E., Moine, M. P., Planton,
S., Saint-Martin, D., Szopa, S., Tyteca, S., Alkama, R., Bela-
mari, S., Braun, A., Coquart, L., and Chauvin, F.: The CNRM-
CM5.1 global climate model: Description and basic evaluation,
Clim. Dynam., 40, 2091–2121, https://doi.org/10.1007/s00382-
011-1259-y, 2013.
Voldoire, A., Saint-Martin, D., Sénési, S., Decharme, B., Alias, A.,
Chevallier, M., Colin, J., Guérémy, J. F., Michou, M., Moine,
M. P., Nabat, P., Roehrig, R., Salas y Mélia, D., Séférian,
R., Valcke, S., Beau, I., Belamari, S., Berthet, S., Cassou,
C., Cattiaux, J., Deshayes, J., Douville, H., Ethé, C., Fran-
chistéguy, L., Geoffroy, O., Lévy, C., Madec, G., Meurdes-
oif, Y., Msadek, R., Ribes, A., Sanchez-Gomez, E., Terray, L.,
and Waldman, R.: Evaluation of CMIP6 DECK Experiments
With CNRM-CM6-1, J. Adv. Model. Earth Sy., 11, 2177–2213,
https://doi.org/10.1029/2019MS001683, 2019.
Webster, C. and Jonas, T.: Influence of canopy shading and
snow coverage on effective albedo in a snow-dominated ev-
ergreen needleleaf forest, Remote Sens. Environ., 214, 48–58,
https://doi.org/10.1016/j.rse.2018.05.023, 2018.
Webster, C., Rutter, N., and Jonas, T.: Improving representation of
canopy temperatures for modeling subcanopy incoming long-
wave radiation to the snow surface, J. Geophys. Res.-Atmos.,
122, 9154–9172, https://doi.org/10.1002/2017JD026581, 2017.
Westermann, S., Lüers, J., Langer, M., Piel, K., and Boike,
J.: The annual surface energy budget of a high-arctic per-
mafrost site on Svalbard, Norway, The Cryosphere, 3, 245–263,
https://doi.org/10.5194/tc-3-245-2009, 2009.
https://doi.org/10.5194/tc-18-231-2024 The Cryosphere, 18, 231–263, 2024