Articles | Volume 18, issue 1
Research article
12 Jan 2024
Research article |  | 12 Jan 2024

Modeling snowpack dynamics and surface energy budget in boreal and subarctic peatlands and forests

Jari-Pekka Nousu, Matthieu Lafaysse, Giulia Mazzotti, Pertti Ala-aho, Hannu Marttila, Bertrand Cluzet, Mika Aurela, Annalea Lohila, Pasi Kolari, Aaron Boone, Mathieu Fructus, and Samuli Launiainen

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., effects of vegetation and topography, as well as limitations in the theoretical understanding of turbulent transfer in the stable boundary layer. Studies that evaluate snow, hydrology 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 different climate and snow conditions, representing the southern and northern subarctic zones, as well as the contrasting forest and peatland ecosystems typical for the boreal landscape. We tested different turbulent flux parameterizations implemented in the Crocus snowpack model. In addition, we examined common alternative approaches to conceptualize soil and vegetation, and we assessed their performance in simulating surface energy fluxes, snow conditions and soil thermal regime. Our results show that a stability correction function that increases the turbulent exchange under stable atmospheric 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 simulations 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 implications 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 (Serreze et al.2009; Screen and Simmonds2010; Boisvert and Stroeve2015; 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 (Tyler2010; Pedersen et al.2021) and plant distribution (Rasmus et al.2011; Kreyling et al.2012; Rissanen et al.2021). Recent studies show that especially cold-dwelling species have been shifting 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 consequences on the quantity and properties of snow may define the destiny of many species and human activities in the boreal region. To predict future snow conditions, environmental change, and the consequences for water resources, ecosystems and people, predictive and process-based models possess 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 models (Douville et al.1995; Niu et al.2011; Lawrence et al.2019) and have more recently become common tools for interdisciplinary 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 (hereafter albedo) and the surface heat fluxes (Cohen and Rind1991; 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. Correctly representing the snowpack is thus essential for simulating 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 Heikinheimo1999; 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 vegetation, the snow cover can exclusively determine the wintertime albedo (Aurela et al.2015). With minimal solar radiation during winter months on these open snow fields, turbulent 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 (Lackner et al.2022; Conway et al.2018). Simulation of turbulent fluxes under stable atmospheric conditions is known as one of the major sources of uncertainty in snow models (Lafaysse et al.2017; Menard et al.2021). In LSMs the turbulent fluxes are commonly computed with bulk aerodynamic approaches, where sensible and latent heat fluxes are proportional to the turbulent exchange coefficient according to the Monin–Obukhov similarity theory. These approaches typically use atmospheric stability correction functions based either on the bulk Richardson number (Martin and Lejeune1998; Lafaysse et al.2017; Clark et al.2015) or the Obukhov length scale (Jordan et al.1999). It is established that the Monin–Obukhov similarity theory does not well represent low-wind and stable atmospheric conditions above aerodynamically smooth surfaces such as snow (Conway 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 (Derbyshire1999; Andreas et al.2010). To circumvent this effect, stability correction functions have been modified to permit turbulent fluxes above critical stability thresholds (Lafaysse et al.2017) by manipulating the wind speed (Martin and Lejeune1998; Andreas et al.2010) or including a windless turbulent exchange coefficient (Jordan et al.1999). Evaluations of these modifications often rely on validation with observed surface temperatures and snow depths (e.g., for the detailed snowpack model Crocus; Martin and Lejeune1998; 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 turbulent exchange is attenuated by the canopy, and the snowpack 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 angles, the albedo dynamics of seasonally snow-covered boreal forests is complex (Malle et al.2021). The absorption 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 represent the large-scale impact of forest on the snowpack energy budget (Niu et al.2011; Lawrence et al.2019; Boone et al.2017) are still prone to large errors, due to the complexity and unresolved spatial scales of the underlying physical processes (Loranty et al.2014; Thackeray et al.2019). The forest snow model evaluations against concurrent snowpack 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 neighbor 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 conductivity from high values at the top to low values at the subsurface, 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 reference sites. In addition, many spatial snow simulations neglect peat soils or SOC altogether, and their hydrological and thermal 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 SURFEX LSM (Masson et al.2013) to describe the surface energy balance and its drivers in boreal and subarctic peatlands and forests. We evaluate the effect of alternative turbulent exchange and snowpack parameterizations and examine the skills of alternative model configurations to represent the soil–snow–vegetation interactions. The modeling framework includes flexible parameterizations for different processes within the Crocus snowpack model (Vionnet et al.2012), and its coupling to ISBA (Noilhan and Mahfouf1996; 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 simulations 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 turbulence 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 composite soil–vegetation and MEB big-leaf forest scheme to assess 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).

Figure 1(a) Study area locations inside the boreal land biome (green area, Olson et al.2001), (b) study site locations in Finland (Esri2023) and (c–f) aerial images of each site (NLSF2020).

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 lasiocarpa) 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 forest is dominated by Norway spruce (Picea abies) with some deciduous trees, mainly birch (Betula pubescens) but also aspen (Populus tremula) and pussy willow (Salix caprea). According 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 distance 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 Sphagnum 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.

Table 1General site information.

Download Print Version | Download XLSX

2.2 Models

We use components from the SURFEX LSM (Surface Externalisée, Masson et al.2013) platform. SURFEX was selected 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 Planton1989; Noilhan and Mahfouf1996) for composite 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 coupled to Crocus is used for both peatland and forest experiments (Fig. 2a, b), whereas MEB coupled to Crocus is only used for the forest experiments (Fig. 2c). In the next subsections, we briefly describe these model components. Parameterizations 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 Atmosphere) is the soil and vegetation component of SURFEX (Noilhan and Planton1989; Noilhan and Mahfouf1996). 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 models 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, W m−2) is computed as the residual of the sum of all surface/atmosphere energy fluxes:

(1) G 0 = SWD ( 1 - LSA ) + ϵ ( LWD - σ T s 4 ) + H + LE ,

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, σ is the Stefan–Boltzmann constant and Ts (K) is the surface 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:

(2) H = ρ a c ρ C H V a ( T s - T a ) ,

where the air density, the specific heat capacity, the wind speed and the air temperature are denoted with ρa (kg m−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; evaporation of intercepted water on the canopy, Ec; transpiration from the vegetation, Etr; and sublimation from bare soil ice, Si:

(3) LE = L v ( E g + E c + E tr ) + ( L f + L v ) ( S i ) ,

where Lv (Jkg−1) and Lf (Jkg−1) are the latent heat of vaporization and fusion, respectively. Total evapotranspiration (ET) is computed as

(4) ET = E g + E c + E tr = ( 1 - veg ) ρ a C H V a [ h u q sat ( T s ) - q a ] + veg ρ a C H V a h v ,

where veg is the fraction of vegetation cover, qsat(Ts) (kgkg−1) is the saturated specific humidity at the surface, qa(Ts) (kg kg−1) is the atmospheric specific humidity, hu is the dimensionless relative humidity at the ground surface related 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 details)

The turbulent exchange coefficient CH is based on the formulation of Louis (1979):

(5) C H = k 2 ln ( z u / z 0 t ) ln ( z a / z 0 t ) f ( R i ) ,

where zu (m) is the reference height of Va, za (m) is the reference 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 number (Ri)  (Louis1979).

Instead of separate treatment of the vegetation canopy and ground, ISBA considers the composite soil–vegetation energy 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. Similar to Napoly et al. (2020), we use the Ags stomatal conductance formulation derived from the coupling of photosynthetic 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 inertia, 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 Mahfouf1996; 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 surface 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 surface/litter layer in MEB are

(6) C v T v t = R nv - H v - LE v + L f ϕ v C g , 1 T g , 1 t = ( 1 - ρ sng ) ( R ng - H g - LE g ) + ρ sng ( G gn + τ n , N n SW nn ) - G g , 1 + L f ϕ g , 1 C n , 1 T n , 1 t = R nn - H n - LE n - τ n , 1 SW nn + ϵ n , 1 - G n , 1 + L f ϕ n , 1 ,

where Cv, Cg,1, and Cn,1 (Jm-2K-1) and Tv, Tg,1, and Tn,1K are the effective heat capacities and temperatures of the canopy, ground surface/litter layer, and snowpack, respectively. 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 conduction heat flux at the bottom of the uppermost soil or snow layer. Ggn is the conduction heat flux at the soil–snow interface. Rnv, Rng and Rnn (Wm−2) are net radiation, i.e., the sum of net shortwave radiation and net longwave radiation from/to the corresponding layer. The shortwave radiation 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 stability 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 hazard forecasting in the French mountain ranges for the past 3 decades (Morin et al.2020). It aims to mimic the vertical layering of snowpacks with a Lagrangian discretization system, avoiding the aggregation of snow layers with highly different physical properties. A detailed description of Crocus and its integration in SURFEX can be found in Vionnet et al. (2012).

The snowpack surface energy budget is the sum of net radiation, turbulent fluxes and advective fluxes from precipitation. Over the snow, the sensible heat flux is computed similarly as in ISBA (Eq. 2) for soil surface, while the latent heat flux (sublimation/deposition), LEs, is computed as

(7) LE s = ( L f + L v ) ρ a C H V a [ q sat ( T s q a ) ] ,

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 system in which both components are solved separately with an implicit approach considering that the state of the second system 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 Mahfouf1996; 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 vegetation 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 regardless of the snow depth (Fig. 2a). It is the common approach 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)

(8) p snv = min 1.0 , HS HS + w sw z 0 ,

(9) p sng = min 1.0 , HS HS g ,

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 prescribed in the subsequent Crocus snowpack and ISBA-DIF soil modules (Fig. 2c).

Figure 2Three 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, Ensemble System Crocus, Lafaysse et al.2017) to evaluate the impact and associated uncertainties of the different parameterizations of snow processes and turbulent exchange. In ESCROC, the main physical processes and properties of snowpack, as well as the turbulent fluxes, can be represented by several alternative options. These include density of new snow, snow metamorphism, absorption of solar radiation, 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 members (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 ensemble spread corresponds to simulated values between ensemble 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 Lejeune1998). Therefore, 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 conditions. 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 effective 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 deciduous 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 neglected in our simulations as major part of the simulation period covers time before the thinning. Vegetation types in ISBA are characterized according to ECOCLIMAP (Champeaux 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 following the standard methods of ISBA (Noilhan and Mahfouf1996; 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 uppermost 1 m, in accordance with field measurements (Väliranta and Mathijssen2021; Muhic et al.2023), while the deeper layers were assigned as mineral soil similar to the contiguous forests. Although peat profiles may be deeper, the soils 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 estimates. The thermal and water retention parameters are subsequently derived from the pedotransfer functions of ISBA (Noilhan and Mahfouf1996; Peters-Lidard et al.1998).

Champeaux et al. (2005)Champeaux et al. (2005)Aurela et al. (2015); Alekseychik et al. (2017)Kolari et al. (2022)Aurela et al. (2015); Alekseychik et al. (2017)Kolari et al. (2022)FMI (2021)Aurela et al. (2015); Mammarella et al. (2019)Alekseychik et al. (2022a)Hari et al. (2013); Alekseychik et al. (2022a)FMI (2021)Lindroos et al. (2022)Muhic et al. (2023)Väliranta and Mathijssen (2021)Lindroos et al. (2022)Muhic et al. (2023)Väliranta and Mathijssen (2021)Väliranta and Mathijssen (2021)Muhic et al. (2023)

Table 2Main model parameters for the study sites. Vegetation types BOGR and BONE correspond to boreal grass and boreal needleleaf evergreen, respectively.

Download Print Version | Download XLSX

2.4 Data

2.4.1 Model forcing

Meteorological forcing consist of hourly observations of air temperature, wind speed, precipitation rate, air humidity, downward shortwave and longwave radiation, and atmospheric pressure. The available meteorological observations from the nearest meteorological stations were obtained from Finnish Meteorological Institute (FMI) open database (FMI2021) (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 function of the cosine of the sun zenith angle, μ. More specifically, 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 typical winter conditions, 0.09 for the aerosol optical thickness, 300 DU for the ozone column and 0.854 g cm−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 meteorological 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 agreement 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 pyrgeometers, 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 measured by the eddy-covariance (EC) technique. The EC systems 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 descriptions 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 instrument 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.1ms−1 and 0.3≤ FST ≤1.0;

  • flag 0: u*0.1ms−1 and FST ≤0.3.

At S-WET, S-FOR and N-FOR, automated snow depth observations are directly used (FMI2021). On N-WET the automated snow depth measurement is at 0.7 m height, and therefore the exceeding snow depths were taken from biweekly 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 (Aurela 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 differences 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). Moreover, the influence of soil texture was explored by comparing M98-SOC to the simulation where soil was parameterized as mineral soil, similar to the contiguous forest site (referred to as M98-MIN).

On the forest sites, we examine the skills of the different 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 compare ESCROC-E2 simulations with these three configurations 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). Additional ISBA-VS simulations are performed to assess the sensitivity of the wsw parameter, using a value of 0.2 in Eq. (8). For a more detailed comparison of the simulated and observed above-canopy surface energy fluxes, we conduct deterministic simulations with both ISBA-VS and MEB, the default snow cover fraction and Crocus parameterizations (as in Fig. 2. in Lafaysse et al. (2017)). Additionally, we perform a deterministic ISBA-VS simulation with the default snow cover fraction parameterization but the vegetation fraction 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 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 demonstrate the results by using the winter season 2018–2019 as an example time series period thanks to its best coverage of energy 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 compare 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 (December, January, February), longwave radiation balance to a large extent determines Rn, particularly in northern Finland. Daily average Rn is negative down to −50Wm−2 and lower, which implies considerable radiative cooling. Towards spring the radiation budget is gradually counterbalanced by shortwave radiation. On the peatlands, a large fraction of SWD is reflected during snow cover, and daily Rn 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 vegetation. In summer, high solar elevation and the absence of the reflective snow surface cause daily Rn to be up to 200 Wm−2.

Figure 3Observed 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).


Rn is balanced mostly by H and LE and to a lesser extent snowpack–ground heat flux (Fig. 3). The residual line represents the amount of energy that would be required to close the observed energy budget (Fig. 3). It includes changes in internal energy of the snowpack and vegetation but also reflects the common energy balance closure problem in EC measurements (Mauder et al.2020) (see Sect. 4.4). The energy balance closure in snow-free conditions was typical for EC measurements, 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 downward H, 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 downward H 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 observed H often lies near the limit or even outside the simulated range at both S-WET and N-WET. Modeled wintertime 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.

Figure 4Time 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.


3.2.1 Impact of alternative turbulence (CH) parameterizations

To assess the sensitivity of snow depth simulations to alternative turbulence parameterizations and to alternative snow process options, we examined simulations where the ESCROC-E2 members are grouped according to their turbulent flux option (Fig. 5). During snow accumulation periods, the spread is small and the groups are consistently overlapping on both sites, indicating that the differences in snow accumulation and maximum snow depth are driven mostly by the uncertainty of snow process descriptions. The spread increases 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.

Figure 5Time 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.


These findings are consistent with the comparison of simulated 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 underestimated, while this bias is to a large extent corrected by using M98-SOC. Improved simulation of H and surface temperature also entail improved LWU (Fig. 6), but the modeled H fluxes still only moderately correlate (R2) with observations. 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.

Figure 6Scatterplots 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 turbulence parameterizations.


3.2.2 Radiative fluxes

We compare the simulated and observed albedo, SWU, LWU and surface temperatures with snow-free and snow conditions in Fig. 7. These experiments correspond to the deterministic M98-SOC simulation.

The modeled SWU generally matches the observations well, but the scatter increases with increasing SWD, indicating uncertainties in simulated albedo when shortwave forcing is high over the snowpack in spring. These cause a slight underestimation of simulated spring albedo, also visible in the time series especially on N-WET (Fig. 7). Moreover, simulated 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 incorrect 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.

Figure 7Scatterplots 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.


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 simulated across sites and ground conditions at least with the presented time intervals. During snow cover, the upper tail of the radiation distribution is slightly higher than simulated; however, the mean biases are generally very low. There are no other visible biases in LWU simulations and the other metrics are also very good, consistent with Fig. 4. The mean absolute errors in simulating LWU are similar for snow and bare ground (∼3 to 8 Wm−2).

3.2.3 Soil thermal regime

The effect of soil parameterization on simulated soil temperature dynamics at S-WET is shown in Fig. 8. Due to shallow water table, the soil profile remains nearly saturated throughout 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 insulating 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).

Figure 8Effect 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 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 (composite 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 parameter (wsw=5) with wsw=0.2 (used for hydrological modeling in Le Moigne et al.2020) yields slightly better snow depth dynamics for N-FOR, but the results remain unsatisfactory (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 balance) 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 observations. 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 process representations prevail in these simulations.

Figure 9Effect of alternative configurations of ISBA and MEB on the snow depth (HS). The envelopes visualize the corresponding ensemble spreads between minimum and maximum values.


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 likely caused by the lack of explicit litter and canopy layers. All model versions tend to overestimate day-to-day temperature variability.

Figure 10Effect 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.


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 energy fluxes at N-FOR. The snow cover periods are defined according to agreement between MEB simulations and observations, 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 overestimated 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 vegetation 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 coefficient (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 assigning the vegetation fraction to unity (i.e., full vegetation coverage and no soil evaporation; see Fig. S7).

Figure 11Simulated 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).


3.3.4 Evolution of albedo

Figure 12 illustrates the time evolution of modeled and observed 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, consistent with results in Sect. 3.3.3. The likely reason for winter albedo underestimation is because the models do not represent changes in albedo due to intercepted snow. The overestimation in spring is presumably due to representing the effective 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.

Figure 12Simulated 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.


3.4 Summary: surface energy budget on peatland and forest sites

Finally, to sum up the whole surface energy budget, we compare how the simulated Rn and turbulent fluxes (H+ LE) match the observations at the four sites. These deterministic simulations are conducted with simulation setups that provided the best fit to data: the deterministic simulation as M98-SOC for the peatland sites and deterministic MEB simulation for the forest sites (Fig. 13).

Despite the challenges in simulating snow depth evolution at the forest sites, the energy budget simulations are generally 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 peatland 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 winter 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 simulations, 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 capture 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. Performance of the simulated summer energy fluxes is very good (Fig. S4).

Figure 13Simulated (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.


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 prevailing stable atmospheric conditions (Table 3), the radiative 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 wintertime energy budget (Lackner et al.2022; Langer et al.2011; Westermann et al.2009).

Table 3Occurrence 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 conditions as 0Ri0.25 and strongly stable conditions as Ri>0.25.

Download Print Version | Download XLSX

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 support 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 contributed 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 turbulent fluxes over snowpack and identified the turbulent exchange parameterizations (M98-SOC and RI2-SOC) that improve the simulated surface energy fluxes and snowpack dynamics at high-latitude winter conditions (Figs. 5, 6). In stable (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 sensible 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 simulations 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 parameterizations 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 default ISBA turbulent flux parameterization (RIL), although widely used, e.g., in numerical weather prediction and general 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.4C at S-WET and −1.1C at N-WET). This finding is consistent with ESM-SnowMIP (Menard et al.2021), where the default configuration of Crocus had one of the lowest skills for surface temperature simulations (−2C mean cold bias) among the compared snow models. Even with the M98-SOC simulation we found a rather low skill of turbulent flux simulations. To summarize, our findings highlight the limitations of the Monin–Obukhov similarity theory to simulate turbulent fluxes under stable atmospheric conditions and emphasize the need for further model developments with observations in various environments.

The soil temperature simulations confirmed that it is necessary 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 thermal regime and consequent freezing–thawing processes in peatlands (Fig. 8) (Dankers et al.2011; Lawrence and Slater2008; 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 dynamics and lateral flow could further improve the soil temperature simulations on boreal and subarctic peatlands. The thermal state and ice/liquid water content also have major cascading 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 development 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 applications that first and foremost require grid-cell-averaged surface energy fluxes (i.e., numerical weather prediction). However, we agree with Napoly et al. (2020) that the snow cover fraction approach of ISBA (Fig. 2b) is essentially a compromise that attempts to retain the insulating impact of the snowpack over the soil while still simulating turbulent exchange from the vegetation. The energy exchange between the atmosphere and the soil–vegetation composite directly impacts the snowpack, and it led to strongly biased snow depth simulations, similar to Napoly et al. (2020). This fractional approach 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 surface 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 snowpack dynamics and soil temperature simulations, ignoring snow–vegetation interactions is a better compromise than using varying snow cover fraction as it is currently implemented in SURFEX. Consequently, ISBA-FS should be preferred 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, neglecting 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 coupling 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 underestimated in winter and overestimated in spring (Fig. 12). The winter albedo underestimation was most likely because intercepted snow increased the observed albedo, a process that is not accounted for in MEB (Napoly et al.2020). This assumption is based on Pomeroy and Dion (1996), while the increase in forest albedo by intercepted snow has been more recently shown (Webster and Jonas2018), and simple descriptions 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 parameterizations may fail to fully capture canopy shading, particularly at low solar elevation angles typical of high latitudes (Malle et al.2021).

In contrast to Napoly et al. (2020), we found MEB to simulate snowmelt too early, especially on N-FOR (Fig. 9). These errors are partly explained by inaccuracies in canopy radiative transfer, but they also suggest errors in simulated below-canopy surface heat fluxes. The differences in snow simulations were rather small between MEB and ISBA-FS, especially 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 decreased snow depths at the immediate vicinity of tree trunks but high snow depth between trees at N-FOR.

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 systematic uncertainties (e.g., Aubinet et al.2012). The absolute values of winter H and LE are small, and their relative uncertainty is high; compared to summertime measurements, the wintertime energy balance closure ratio is typically poorer, particularly at the northern ecosystems (Reba et al.2009; Molotch et al.2009; Launiainen2010). 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, undercatch, etc., Stuefer et al.2020), and data gaps further add up possible sources of errors. Uncertainties in model forcing 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 Essery1999; Meriö et al.2023; Liston and Sturm2002) 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. Omission of internal water vapor transfer by diffusion and/or convection 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 consequently 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. Complementary 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 conditions 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 confirmed by our data (Fig. 9). The small-scale forest structure 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; Rutter et al.2009). More realistic below- and above-canopy heat flux simulations could be achieved by more sophisticated canopy representations, including multiple layers and species (Bonan et al.2021; McGowan et al.2017; Launiainen 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 alternative to traditional LSMs (Broxton et al.2015; Mazzotti et al.2020b).

The generality of our findings should be tested by additional 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 observations of all components of surface energy balance, particularly 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 evaluate turbulent exchange parameterizations and alternative approaches 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 parameterizations (ESCROC-E2) ensures that uncertainties in snow processes (not evaluated in this study) do not affect the robustness of our main conclusions.

Peatland simulations showed that using a stability correction function that increases the turbulent exchange under stable atmospheric conditions is imperative to simulate the snowpack energy budget. This adjustment led to major improvements 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 dynamics, 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) performance was satisfactory only after an adjustment of a sensitive vegetation fraction parameter. In particular, shortwave and longwave radiation balances were simulated well by both approaches, whereas the turbulent fluxes had significantly higher uncertainty. Only the explicit vegetation model (MEB) was able to simultaneously simulate realistic surface energy budget, snow and soil conditions. The composite approaches only succeeded in simulating either the correct surface energy budget (ISBA-VS) or snow and soil conditions (ISBA-FS). Furthermore, we demonstrated that the composite approaches rely on a previously poorly documented parameterization of the snow cover fraction with high sensitivity on model outputs despite a limited physical interpretation.

With well-selected model configuration and parameterization, 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 version 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 predictions of forest snow and soil conditions necessary for hydrological applications. We expect that the future inclusion of MEB in operational systems will reconcile these applications. Our results can be used to inform the choice of model configuration for studies of subarctic and boreal region ecology, hydrology and biogeochemistry under the ongoing environmental change.

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

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
Appendix B: Model evaluation data availability and radiation forcing evaluation

The availability periods of the model evaluation data are presented in Table B1.

Table B1Model evaluation data availability for each site.

Download Print Version | Download XLSX

Figure B1Comparison of longwave radiation forcing data between ERA5 and site observations (OBS) on S-FOR.


Appendix C: Effect of snow cover fraction parameter wsw
Bengtsson et al. (2017)Courtier et al. (1991)Decharme et al. (2019)Caillaud et al. (2021)Nabat et al. (2020)Le Moigne et al. (2020)Verrelle et al. (2021)Vernay et al. (2022)Brun et al. (2013)Morin et al. (2020)

Table C1Summary 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.

Download Print Version | Download XLSX

Figure C1Effect of wsw parameter on snow depths simulated by ISBA-VS. The envelopes visualize the corresponding ensemble spreads between minimum and maximum values.


Figure C2The 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.


Code availability

SURFEX is an open-source project (, Le Moigne, 2018), but it requires registration. The full procedure and instructions are available at (Lafaysse et al., 2023), and the SURFEX version used in this study is archived at (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 (Cluzet and Nousu, 2023).

Data availability

Meteorological data, evaluation data, and SURFEX specific namelist and forcing files are available as an electronic supplement at (Nousu et al.2023) under the terms and conditions of the Creative Commons Attribution 4.0 International license.


The supplement related to this article is available online at:

Author contributions

JPN, ML, SL, PA and HM designed the research. 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 writing 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.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


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 acknowledge the GreenFeedBack project from the EU Horizon Europe 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 Aurela 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 internationalization support grant. The authors would like to thank Bertrand Decharme, Christine Delire, and Bertrand Bonan at CNRM/GAME and Marie Dumont at CNRM/CEN for their great support 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 Nationalfonds 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.


Aalto, J., Aalto, P., Keronen, P., Kolari, P., Rantala, P., Taipale, R., Kajos, M., Patokoski, J., Rinne, J., Ruuskanen, 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,, 2022. a, b

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 seasonally frozen ground on hydrological partitioning? A systematic review, Environ. Res. Lett., 16, 4,, 2021. a

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,, 2022a. a, b, c

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,, 2022b. a

Alekseychik, P. K., Korrensalo, A., Mammarella, I., Vesala, T., and Tuittila, E. S.: Relationship between aerodynamic roughness length and bulk sedge leaf area index in a mixed-species boreal mire complex, Geophys. Res. Lett., 44, 5836–5843,, 2017. a, b, c

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,, 2010. a, b

Aubinet, M., Vesala, T., and Papale, D.: Eddy covariance: a practical guide to measurement and data analysis, Springer Science & Business Media, 2012. a

Aurela, M., Riutta, T., Laurila, T., Tuovinen, J. P., Vesala, T., Tuittila, 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,, 2007. a

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. a, b, c, d, e, f, g, h

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,, 2017. a

Bengtsson, L., Andrae, U., Aspelien, T., Batrak, Y., Calvo, J., Rooy, W. d., Gleeson, E., Hansen-Sass, B., Homleid, M., Hortal, 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,, 2017. a, b

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,<3324:TROASI>2.0.CO;2, 2001. a

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 Surface Modelling, Current Climate Change Reports, 7, 45–71,, 2021. a

Boisvert, L. N. and Stroeve, J. C.: The Arctic is becoming warmer and wetter as revealed by the Atmospheric Infrared Sounder, Geophys. Res. Lett., 42, 4439–4446,, 2015. a

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,, 2021. a

Boone, A., Masson, V., Meyers, T., and Noilhan, J.: The influence of the inclusion of soil freezing on simulations by a soil-vegetation-atmosphere transfer scheme, J. Appl. Meteorol., 39, 1544–1569,<1544:TIOTIO>2.0.CO;2, 2000. a

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,, 2017. a, b, c, d, e, f, g, h, i, j

Bouchard, B., Nadeau, D. F., and Domine, F.: Comparison of snowpack structure in gaps and under the canopy in a humid boreal forest, Hydrol. Process., 36, e14681,, 2022. a

Broxton, P. D., Harpold, A. A., Biederman, J. A., Troch, P. A., Molotch, N. P., and Brooks, P. D.: Quantifying the effects of vegetation structure on snow accumulation and ablation in mixed-conifer forests, Ecohydrology, 8, 1073–1094,, 2015. a

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,, 2013. a, b, c, d

Brunet, Y.: Turbulent Flow in Plant Canopies: Historical Perspective and Overview, vol. 177, Springer Netherlands, ISBN 0123456789,, 2020. a, b

Caillaud, C., Somot, S., Alias, A., Bernard-Bouissières, I., Fumiè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,, 2021. a

Calvet, J. C., Noilhan, J., Roujean, J. L., Bessemoulin, P., Cabelguenne, M., Olioso, A., and Wigneron, J. P.: An interactive vegetation SVAT model tested against data from six contrasting sites, Agr. Forest Meteorol., 92, 73–95,, 1998. a

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,, 2013. a, b

Chadburn, S., Burke, E., Essery, R., Boike, J., Langer, M., Heikenfeld, M., Cox, P., and Friedlingstein, P.: An improved representation of physical permafrost dynamics in the JULES land-surface model, Geosci. Model Dev., 8, 1493–1508,, 2015. a

Champeaux, J. L., Masson, V., and Chauvin, F.: ECOCLIMAP: A global database of land surface parameters at 1 km resolution, Meteorol. Appl., 12, 29–32,, 2005. a, b, c

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 hydrologic and land-surface models: A review, Water Resour. Res., 47, W07539,, 2011. a

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 modeling: 1. Modeling concept, Water Resour. Res., 51, 2498–2514,, 2015. a, b

Cluzet, B. and Nousu, J.-P.: bertrandcz/CrocO_toolbox: CrocO_v1.2 (CrocO_v1.2), Zenodo [code],, 2024. 

Cohen, J. and Rind, D.: The Effect of Snow Cover on the Climate, J. Climate, 4, 689–706,<0689:TEOSCO>2.0.CO;2, 1991. a

Conway, J. P., Pomeroy, J. W., Helgason, W. D., and Kinar, N. J.: Challenges in modeling turbulent heat fluxes to snowpacks in forest clearings, J. Hydrometeorol., 19, 1599–1616,, 2018. a, b, c, d

Couet, J., Marjakangas, E. L., Santangeli, A., Kålås, J. A., Lindström, Ã., and Lehikoinen, A.: Short-lived species move uphill faster under climate change, Oecologia, 198, 877–888,, 2022. a

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. a

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,, 2011. a

Decharme, B., Boone, A., Delire, C., and Noilhan, J.: Local evaluation of the Interaction between Soil Biosphere Atmosphere soil multilayer diffusion scheme using four pedotransfer functions, J. Geophys. Res.-Atmos., 116, 1–29,, 2011. a, b

Decharme, B., Brun, E., Boone, A., Delire, C., Le Moigne, P., and Morin, S.: Impacts of snow and organic soils parameterization on northern Eurasian soil temperature profiles simulated by the ISBA land surface model, The Cryosphere, 10, 853–877,, 2016. a, b, c, d, e

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,, 2019. a, b, c, d

Derbyshire, S. H.: Boundary-layer decoupling over cold surfaces as a physical boundary-instability, Bound.-Lay. Meteorol., 90, 297–325,, 1999. a

Deschamps-Berger, C., Cluzet, B., Dumont, M., Lafaysse, M., Berthier, E., Fanise, P., and Gascoin, S.: Improving the Spatial Distribution of Snow Cover Simulations by Assimilation of Satellite Stereoscopic Imagery, Water Resour. Res., 58, e2021WR030271,, 2022. a

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,, 2018. a

Douville, H., Royer, J. F., and Mahfouf, J. F.: A new snow parameterization for the Météo-France climate model, Clim. Dynam., 12, 21–35,, 1995. a, b

Esri: ESRI Satellite (ArcGIS/World_Imagery), (last access: 19 September 2023), 2023. a

Essery, R., Rutter, N., Pomeroy, J., Baxter, R., Stahli, M., Gustafsson, D., Barr, A., Bartlett, P., and Elder, K.: SnowMIP2: An evalution of forest snow process simulation, B. Am. Meteor. Soc., 90, 1130–1135,, 2009. a, b

Eugster, W., Rouse, W. R., Pielke, R. A., Mcfadden, J. P., Baldocchi, D. D., Kittel, T. G., Chapin, F. S., Liston, G. E., Vidale, 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,, 2000. a

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,, 2013. a

FMI: Finnish Meteorological Institute past weather observations, (last access: 19 September 2023), 2021. a, b, c, d

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,, 2005. a

Gouttevin, I., Lehning, M., Jonas, T., Gustafsson, D., and Mölder, M.: A two-layer canopy model with thermal inertia for an improved snowpack energy balance below needleleaf forest (model SNOWPACK, version 3.2.1, revision 741), Geosci. Model Dev., 8, 2379–2398,, 2015. a

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,, 2023. a

Hamdi, R., Degrauwe, D., Duerinckx, A., Cedilnik, J., Costa, V., Dalkilic, T., Essaouini, K., Jerczynki, M., Kocaman, F., Kullmann, L., Mahfouf, J.-F., Meier, F., Sassi, M., Schneider, S., Váňa, F., and Termonia, P.: Evaluating the performance of SURFEXv5 as a new land surface scheme for the ALADINcy36 and ALARO-0 models, Geosci. Model Dev., 7, 23–39,, 2014. a

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 Ecology, 9789400756, 471–487,, 2013. a

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., Flemming, 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,, 2020. a

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,, 1999. a, b

Koivusalo, H. and Heikinheimo, M.: Surface energy exchange over a boreal snowpack: Comparison of two snow energy balance models, Hydrol. Process., 13, 2395–2408,<2395::AID-HYP864>3.0.CO;2-G, 1999. a

Kolari, P., Aalto, J., Levula, J., Kulmala, L., Ilvesniemi, H., and Pumpanen, J.: Hyytiälä SMEAR II site characteristics, Zenodo [data set],, 2022. a, b

Kreyling, J., Haei, M., and Laudon, H.: Absence of snow cover reduces understory plant cover and alters plant community composition in boreal forests, Oecologia, 168, 577–587,, 2012. a

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,, 2018. a

Lackner, G., Domine, F., Nadeau, D. F., Parent, A.-C., Anctil, F., Lafaysse, M., and Dumont, M.: On the energy budget of a low-Arctic snowpack, The Cryosphere, 16, 127–142,, 2022. a, b, c, d, e

Lafaysse, M., Hingray, B., Etchevers, P., Martin, E., and Obled, C.: Influence of spatial discretization, underground water storage and glacier melt on a physically-based hydrological model of the Upper Durance River basin, J. Hydrol., 403, 116–129,, 2011. a

Lafaysse, M., Cluzet, B., Dumont, M., Lejeune, Y., Vionnet, V., and Morin, S.: A multiphysical ensemble system of numerical snow modelling, The Cryosphere, 11, 1173–1198,, 2017. a, b, c, d, e, f, g, h, i, j, k

Lafaysse, M., Fructus, M., Vernay, M., Radanovics, S., Dumont, M., and Viallon-Galinier, L.: Procedure for new users of Crocus model,, (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 northern Siberia – Part 2: Winter, The Cryosphere, 5, 509–524,, 2011. a

Lapo, K., Nijssen, B., and Lundquist, J. D.: Evaluation of Turbulence Stability Schemes of Land Models for Stable Conditions, J. Geophys. Res.-Atmos., 124, 3072–3089,, 2019. a, b

Launiainen, S.: Seasonal and inter-annual variability of energy exchange above a boreal Scots pine forest, Biogeosciences, 7, 3921–3940,, 2010. a, b

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,, 2015. a

Launiainen, S., Katul, G. G., Leppä, K., Kolari, P., Aslan, T., Grönholm, 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,, 2022. a

Lawrence, D. M. and Slater, A. G.: Incorporating organic soil into a global climate model, Clim. Dynam., 30, 145–160,, 2008. a

Lawrence, D. M., Fisher, R. A., Koven, C. D., Oleson, K. W., Swenson, S. C., Bonan, G., Collier, N., Ghimire, B., van Kampenhout, 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., Hoffman, 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 Martin, 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,, 2019. a, b

Le Moigne, P.: SURFEX scientific documentation, SURFEX v8.1, Issue 3, Météo-France, Toulouse, France, (last access: 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., Leroux, D., and Rousset-Regimbeau, F.: The latest improvements with SURFEX v8.0 of the Safran–Isba–Modcou hydrometeorological model for France, Geosci. Model Dev., 13, 3925–3946,, 2020. a, b, c, d, e

Lindroos, A. J., Mäkipää, R., and Merilä, P.: Soil carbon stock changes over 21 years in intensively monitored boreal forest stands in Finland, Ecol. Indic., 144, 109551,, 2022. a, b, c

Liston, G. E. and Sturm, M.: Winter precipitation patterns in arctic Alaska determined from a blowing-snow model and snow-depth observations, J. Hydrometeorol., 3, 646–659,<0646:WPPIAA>2.0.CO;2, 2002. a

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. a

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,, 2014. a

Louis, J. F.: A parametric model of vertical eddy fluxes in the atmosphere, Bound.-Lay. Meteorol., 17, 187–202,, 1979. a, b, c, d

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,, 2013. a

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 lacking appropriate temperature sensitivities, Hydrol. Process., 35, 1–20,, 2021. a

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 Preliminary Results, American Meteorological Society, 8, 2039–2057, 1995. a, b

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,, 2021. a, b, c

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 processing steps in two contrasting ecosystems, Atmos. Meas. Tech., 9, 4915–4933,, 2016. a

Mammarella, I., Rannik, Ü., Launiainen, S., Alekseychik, P., Peltola, O., Keronen, P., Kolari, P., Laakso, H., Matilainen, T., Salminen, T., Levula, J., Pohja, T., Siivola, E., and Vesala, T.: SMEAR II Hyytiälä forest eddy covariance, Fairdata,, 2019. a, b

Martin, E. and Lejeune, Y.: Turbulent fluxes above the snow surface, Ann. Glaciol., 26, 179–183,, 1998. a, b, c, d, e, f, g

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 carbon cycling – Leading the way for future studies using integrated datasets at Pallas, Finland, Hydrol. Process., 35, 1–19,, 2021. a, b, c, d

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., Salgado, 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,, 2013. a, b

Mauder, M., Foken, T., and Cuxart, J.: Surface-Energy-Balance Closure over Land: A Review, vol. 177, Springer Netherlands, ISBN 0123456789,, 2020. a, b

Mazzotti, G., Essery, R., Moeser, C. D., and Jonas, T.: Resolving Small-Scale Forest Snow Patterns Using an Energy Balance Snow Model With a One-Layer Canopy, Water Resour. Res., 56, e2019WR026129,, 2020a. a

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 Resour. Res., 56, 1–25,, 2020b. a, b, c

Mazzotti, G., Webster, C., Essery, R., and Jonas, T.: Increasing the Physical Representation of Forest-Snow Processes in Coarse-Resolution Models: Lessons Learned From Upscaling Hyper-Resolution Simulations, Water Resour. Res., 57, 1–21,, 2021. a

McGowan, L. E., Paw U, K. T., and Dahlke, H. E.: What Does a Multilayer Canopy Model Tell Us About Our Current Understanding of Snow-Canopy Unloading?, in: AGU Fall Meeting Abstracts, New Orleans, 11–15 December 2017, vol. 2017, C43B–06,, 2017. a

Ménard, C. B., Essery, R., Barr, A., Bartlett, P., Derry, J., Dumont, M., Fierz, C., Kim, H., Kontu, A., Lejeune, Y., Marks, D., Niwano, M., Raleigh, M., Wang, L., and Wever, N.: Meteorological 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,, 2019. a

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., Hagemann, S., Haverd, V., Kim, H., Lafaysse, M., Marke, T., Nasonova, O., Nitta, T., Niwano, M., Pomeroy, J., Schädler, G., Semenov, 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,, 2021. a, b, c, d

Menberu, M. W., Marttila, H., Ronkanen, A. K., Haghighi, A. T., and Kløve, B.: Hydraulic and Physical Properties of Managed and Intact Peatlands: Application of the Van Genuchten-Mualem Models to Peat Soils, Water Resour. Res., 57, 1–22,, 2021. a, b

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,, 2023. a, b, c

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. a

Morin, S., Horton, S., Techel, F., Bavay, M., Coléou, C., Fierz, C., Gobiet, A., Hagenmuller, P., Lafaysse, M., Ližar, M., Mitterer, 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,, 2020. a, b

Morris, P. J., Davies, M. L., Baird, A. J., Balliston, N., Bourgault, 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.: Saturated Hydraulic Conductivity in Northern Peats Inferred From Other Measurements, Water Resour. Res., 58, e2022WR033181,, 2022. a, b

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,, 2023. a, b, c, d

Mustamo, P., Ronkanen, A. K., Berglund, Ã., Berglund, K., and Kløve, B.: Thermal conductivity of unfrozen and partially frozen managed peat soils, Soil Till. Res., 191, 245–255,, 2019. a

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 atmospheric circulation over the Euro-Mediterranean region, Atmos. Chem. Phys., 20, 8315–8349,, 2020. a

Napoly, A., Boone, A., Samuelsson, P., Gollvik, S., Martin, E., Seferian, R., Carrer, D., Decharme, B., and Jarlan, L.: The interactions between soil–biosphere–atmosphere (ISBA) land surface 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,, 2017. a, b, c

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,, 2020. a, b, c, d, e, f, g, h, i

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,, 2007. a

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 description and evaluation with local-scale measurements, J. Geophys. Res.-Atmos., 116, 1–19,, 2011. a, b

NLSF: National Land Survey of Finland Topographic Database, (last access: 29 July 2023), 2020. a

Noilhan, J. and Mahfouf, J. F.: The ISBA land surface parameterisation scheme, Global Planet. Change, 13, 145–159,, 1996. a, b, c, d, e, f, g, h

Noilhan, J. and Planton, S.: A simple parameterization of land surface processes for meteorological models, Mon. Weather Rev., 117, 536–549,<0536:ASPOLS>2.0.CO;2, 1989. a, b

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 Catchment Water Resources Research, Water Resour. Res., 59, 1–19,, 2022. a, b

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,, 2019. a

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 surface energy budget in boreal and subarctic peatlands and forests, Zenodo [data set],, 2023. a

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],, 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.: Terrestrial ecoregions of the world: A new map of life on Earth, BioScience, 51, 933–938,[0933:TEOTWA]2.0.CO;2, 2001. a

Park, H., Launiainen, S., Konstantinov, P. Y., Iijima, Y., and Fedorov, A. N.: Modeling the Effect of Moss Cover on Soil Temperature and Carbon Fluxes at a Tundra Site in Northeastern Siberia, J. Geophys. Res.-Biogeo., 123, 3028–3044,, 2018. a

Pedersen, S. H., Bentzen, T. W., Reinking, A. K., Liston, G. E., Elder, K., Lenart, E. A., Prichard, A. K., and Welker, J. M.: Quantifying effects of snow depth on caribou winter range selection and movement in Arctic Alaska, Movement Ecology, 9, 1–24,, 2021. a

Peters-Lidard, C. D., Blackburn, E., Liang, X., and Wood, E. F.: The effect of soil thermal conductivity parameterization on surface energy fluxes and temperatures, J. Atmos. Sci., 55, 1209–1224,<1209:TEOSTC>2.0.CO;2, 1998. a, b, c

Pomeroy, J. W. and Dion, K.: Winter radiation extinction and reflection in a boreal pine canopy: Measurements and modelling, Hydrol. Process., 10, 1591–1608,<1591::aid-hyp503>;2-8, 1996. a

Pomeroy, J. W. and Essery, R. L.: Turbulent fluxes during blowing snow: Field tests of model sublimation predictions, Hydrol. Process., 13, 2963–2975,<2963::AID-HYP11>3.0.CO;2-9, 1999. a

Raleigh, M. S., Lundquist, J. D., and Clark, M. P.: Exploring the impact of forcing error characteristics on physically based snow simulations within a global sensitivity analysis framework, Hydrol. Earth Syst. Sci., 19, 3153–3179,, 2015. a

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 the globe since 1979, Commun. Earth Environ., 3, 1–10,, 2022. a

Rasmus, S., Lundell, R., and Saarinen, T.: Interactions between snow, canopy, and vegetation in a boreal coniferous forest, Plant Ecol. Divers., 4, 55–65,, 2011. a

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,, 2009. a

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 evaluation of snowpack simulations in complex alpine terrain using satellite and in situ observations, Remote Sensing, 10, 1–32,, 2018. a

Ricchiazzi, P., Yang, S., Gautier, C., and Sowle, D.: SBDART: A Research and Teaching Software Tool for Plane-Parallel Radiative Transfer in the Earth's Atmosphere, B. Am. Meteor. Soc., 79, 2101–2114,<2101:SARATS>2.0.CO;2, 1998. a

Richter, B., Schweizer, J., Rotach, M. W., and Van Herwijnen, A.: Modeling spatially distributed snow instability at a regional scale using Alpine3D, J. Glaciol., 67, 1147–1162,, 2021. a

Rinne, J., Tuittila, E. S., Peltola, O., Li, X., Raivonen, M., Alekseychik, 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,, 2018. a

Rissanen, T., Niittynen, P., Soininen, J., and Luoto, M.: Snow information is required in subcontinental scale predictions of mountain plant distributions, Global Ecol. Biogeogr., 30, 1502–1513,, 2021. a

Rutter, N., Essery, R., Pomeroy, J., Altimir, N., Andreadis, K., Baker, I., Barr, A., Bartlett, P., Boone, A., Deng, H., Douville, 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,, 2009. a, b, c

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. a, b

Screen, J. A. and Simmonds, I.: The central role of diminishing sea ice in recent Arctic temperature amplification, Nature, 464, 1334–1337,, 2010. a

Serreze, M. C., Barrett, A. P., Stroeve, J. C., Kindig, D. N., and Holland, M. M.: The emergence of surface-based Arctic amplification, The Cryosphere, 3, 11–19,, 2009. a

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 Rosnay, 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 representation of snow in land surface schemes: Results from PILPS 2(d), J. Hydrometeorol., 2, 7–25,<0007:TROSIL>2.0.CO;2, 2001. a

Stiegler, C., Lund, M., Christensen, T. R., Mastepanov, M., and Lindroth, A.: Two years with extreme and little snowfall: effects on energy partitioning and surface energy exchange in a high-Arctic tundra ecosystem, The Cryosphere, 10, 1395–1413,, 2016. a

Stigter, E. E., Steiner, J. F., Koch, I., Saloranta, T. M., Kirkham, J. D., and Immerzeel, W. W.: Energy and mass balance dynamics of the seasonal snowpack at two high-altitude sites in the Himalaya, Cold Reg. Sci. Technol., 183, 103233,, 2021. a

Stuefer, S. L., Kane, D. L., and Dean, K. M.: Snow Water Equivalent Measurements in Remote Arctic Alaska Watersheds, Water Resour. Res., 56, 1–12,, 2020. a

Tayleur, C. M., Devictor, V., Gaüzère, P., Jonzén, N., Smith, H. G., and Lindström, Ã.: Regional variation in climate change winners and losers highlights the rapid loss of cold-dwelling species, Divers. Distrib., 22, 468–480,, 2016. a

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,, 2019. a

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,, 2006. a

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 particles during two contrasted snow seasons at Col du Lautaret (2058 m a.s.l., French Alps), The Cryosphere, 14, 4553–4579,, 2020. a

Tyler, N. J.: Climate, snow, ice, crashes, and declines in populations of reindeer and caribou (Rangifer tarandus L.), Ecol.Monogr., 80, 197–219,, 2010. a

Väliranta, M. and Mathijssen, P. J. H.: Geochemistry of Siikaneva peat core from Finland, PANGAEA,, 2021. a, b, c, d

Varhola, A., Coops, N. C., Weiler, M., and Moore, R. D.: Forest canopy effects on snow accumulation and ablation: An integrative review of empirical results, J. Hydrol., 392, 219–233,, 2010. a

Vernay, M., Lafaysse, M., Monteiro, D., Hagenmuller, P., Nheili, R., Samacoïts, R., Verfaillie, D., and Morin, S.: The S2M meteorological and snow cover reanalysis over the French mountainous areas: description and evaluation (1958–2021), Earth Syst. Sci. Data, 14, 1707–1733,, 2022. a, b, c, d, e

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,, 2021. a

Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le Moigne, P., Martin, E., and Willemet, J.-M.: The detailed snowpack scheme Crocus and its implementation in SURFEX v7.2, Geosci. Model Dev., 5, 773–791,, 2012. a, b, c, d, e

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., Fernandez, E., Madec, G., Maisonnave, E., Moine, M. P., Planton, S., Saint-Martin, D., Szopa, S., Tyteca, S., Alkama, R., Belamari, S., Braun, A., Coquart, L., and Chauvin, F.: The CNRM-CM5.1 global climate model: Description and basic evaluation, Clim. Dynam., 40, 2091–2121,, 2013. a, b

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., Franchistéguy, L., Geoffroy, O., Lévy, C., Madec, G., Meurdesoif, 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,, 2019.  a, b

Webster, C. and Jonas, T.: Influence of canopy shading and snow coverage on effective albedo in a snow-dominated evergreen needleleaf forest, Remote Sens. Environ., 214, 48–58,, 2018. a

Webster, C., Rutter, N., and Jonas, T.: Improving representation of canopy temperatures for modeling subcanopy incoming longwave radiation to the snow surface, J. Geophys. Res.-Atmos., 122, 9154–9172,, 2017. a

Westermann, S., Lüers, J., Langer, M., Piel, K., and Boike, J.: The annual surface energy budget of a high-arctic permafrost site on Svalbard, Norway, The Cryosphere, 3, 245–263,, 2009. a

Short summary
The snowpack has a major impact on the land surface energy budget. Accurate simulation of the snowpack energy budget is difficult, and studies that evaluate models against energy budget observations are rare. We compared predictions from well-known models with observations of energy budgets, snow depths and soil temperatures in Finland. Our study identified contrasting strengths and limitations for the models. These results can be used for choosing the right models depending on the use cases.