Improvements in the land surface configuration to better simulate seasonal snow cover in the European Alps with the CNRM-AROME (cycle 46) convection-permitting regional climate model

Snow cover modeling remains a major challenge in climate and numerical weather prediction (NWP) models even in recent versions of high-resolution coupled surface–atmosphere (i.e., at kilometer scale) regional models. Evaluation of recent climate simulations, carried out as part of the WCRP-CORDEX Flagship Pilot Study on Convection (FPSCONV) with the CNRM-AROME convection-permitting regional climate model at 2.5 km horizontal resolution, has highlighted significant snow cover biases, severely limiting its potential in mountain regions. These biases, which are also found in AROME numerical weather prediction (NWP) model results, have multiple causes, involving atmospheric processes and their influence on input data to the land surface models in addition to deficiencies of the land surface model itself. Here we present improved configurations of the SURFEX-ISBA land surface model used in CNRM-AROME. We thoroughly evaluated these configurations on their ability to represent seasonal snow cover across the European Alps. Our evaluation was based on coupled simulations spanning the winters of 2018–2019 and 2019–2020, which were compared against remote sensing data and in situ observations. More specifically, the study tests the influence of various changes in the land surface configuration, such as the use of multi-layer soil and snow schemes, the division of the energy balance calculation by surface type within a grid cell (multiple patches), and new physiographic databases and parameter adjustments. Our findings indicate that using only more detailed individual components in the surface model did not improve the representation of snow cover due to limitations in the approach used to account for partial snow cover within a grid cell. These limitations are addressed in further configurations that highlight the importance, even at kilometer resolution, of taking into account the main subgrid surface heterogeneities and improving representations of interactions between fractional snow cover and vegetation. Ultimately, we introduce a land surface configuration, which substantially improves the representation of seasonal snow cover in the European Alps in coupled CNRM-AROME simulations. This holds promising potential for the use of such model configurations in climate simulations and numerical weather prediction both for AROME and other high-resolution climate models.

1 Introduction

Accurate modeling of the interactions between land surfaces and the atmosphere in mountainous regions is crucial for numerical weather prediction (NWP) and climate projections. Applications range from short-term forecasts for weather-dependent human activities (risk management, hydropower production, tourism, and traffic management) to long-term studies of the impacts of climate change on the various components of a mountain range. In these regions, an appropriate representation of the seasonal snow cover is crucial as its presence strongly affects the evolution of the surface and near-surface conditions by the modification of the albedo, the roughness of the terrain, and its insulating properties for the underlying soil.

Snow models used in NWP and climate models have been widely evaluated and tested in standalone configurations, most often at the point scale and driven by observations or reanalysis of near-surface atmospheric variables (Decharme et al., 2016; Menard et al., 2020). However, when used in coupled surface–atmosphere models, they can produce significantly different results due to the combined effect of errors arising from atmospheric modeling (Raleigh et al., 2015; Lapo et al., 2015), the use of subgrid parameterizations to account for surface heterogeneity within a grid cell of a discretized domain, and their own deficiencies. Testing snow models in standalone “offline” configurations is not sufficient, and tests in coupled configurations are required. This is particularly challenging in mountainous areas.

Indeed, modeling atmospheric and surface fluxes as well as snow cover in mountain regions is challenging in many aspects. The complex topography induces a number of atmospheric phenomena on a wide range of spatio-temporal scales (e.g., Föhn effect, convection phenomena, preferential deposition of snowfall, temperature inversions, and wind-induced snow transport), which have a major impact on surface weather conditions. In addition, the strong heterogeneities of the surface characteristics (elevation, surface type, and aspect) generate high variability in near-surface conditions at sub-kilometer scales, affecting all surface components, including the snowpack.

Most regional coupled atmosphere–surface climate models (RCMs) exhibit deviations with respect to observational references, which can be particularly substantial in mountainous areas. In the European Alps, numerous studies have evaluated the EURO-CORDEX regional climate simulation ensembles at a horizontal resolution from 12 to 50 km and identified strong biases in near-surface precipitation and temperature indicators (Kotlarski et al., 2014; Smiatek et al., 2016; Vorkauf et al., 2021), showing generally excessive precipitation and temperature values that are too low. In general, snow cover (depth, mass, and duration) is overestimated (Terzago et al., 2017; Matiu et al., 2020). One potential approach to mitigate these issues is to develop and apply kilometer-scale modeling frameworks, such as those considered in the WCRP-CORDEX Flagship Pilot Study (FPS) on Convection (Coppola et al., 2020; Ban et al., 2021; Pichelli et al., 2021). In addition to their capacity to explicitly resolve deep convection and thereby enhance the representation of precipitation extremes (Caillaud et al., 2021), models operating at the kilometer scale make it possible to better represent the topography of mountain areas and the heterogeneities that characterize the surface through higher resolution, holding great potential for mountain regions.

At Météo-France, the limited-area non-hydrostatic model AROME (Applications de la Recherche à l'Opérationnel à Méso-Echelle; Seity et al., 2011; Brousseau et al., 2016) has been used operationally for NWP since 2008, initially at 2.5 km horizontal resolution and 1.3 km since 2015 and used for climate studies, referred to as CNRM-AROME, since 2014 (Déqué et al., 2016; Fumière et al., 2019; Caillaud et al., 2021) at 2.5 km horizontal resolution. Simulation results of these models exhibit a number of issues that limit their use and relevance in mountain regions. Indeed, in a recent study comparing 30 years of past climate simulations carried out with CNRM-AROME with the S2M reanalysis (Vernay et al., 2022) over the French Alps, we highlighted a negative temperature difference on the order of 2 to 3 °C, with a maximum in winter at high elevations, and an excess amount of precipitation, particularly at high elevations (Monteiro et al., 2022).

In these climate simulations, we were also able to identify substantial snow cover biases, such as an excessive snow accumulation at intermediate and high elevations, with an overestimated snow cover extent and duration (Monteiro and Morin, 2023) and unrealistic snow accumulation in some grid cells, reaching several hundreds of meters after 30 years of simulation. A near-surface temperature bias has also been identified and analyzed in the NWP version of AROME (Vionnet et al., 2016; Arnould et al., 2021; Gouttevin et al., 2023). These dismissed issues related to the horizontal resolution of the model but rather pointed towards multiple other factors – namely, the underestimation of the cloud cover, also identified by Lucas-Picher et al. (2023) in CNRM-AROME climate simulations; the underestimation of turbulent mixing under stable conditions; and strongly underestimated sub-surface soil temperatures used to diagnose the near-surface air temperature.

There are certainly multiple reasons for the widespread overestimation of snow amount and duration in AROME model results. Monteiro et al. (2022) identified several possible factors:

• biased atmospheric forcings, such as an overestimation of snowfall and an underestimation of melting due to excessively cold near-surface temperatures and errors in downward radiation fluxes,

• the use of an overly simplified surface configuration (one-layer snow model and force–restore soil scheme),

• the lack of glacier dynamics and snow redistribution processes leading to the creation of “snow towers” in some high-elevation grid cells (Freudiger et al., 2017).

The land surface configuration used in the current version of the CNRM-AROME model (also used in the current version of AROME used for NWP applications) does not provide an adequate representation of snow cover dynamics over the French and European Alps, with potential effects on other surface variables of interest, such as the near-surface air temperature (Monteiro et al., 2022).

In this study, we investigate the representation of snow cover for a set of surface model configurations already implemented within the land surface model SURFEX-ISBA (Noilhan and Mahfouf, 1996; Masson et al., 2013) but not yet evaluated in a coupled surface–atmosphere context at high resolution, such as CNRM-AROME, especially in mountainous regions.

In this context, we document the advantages and limitations of using different levels of complexity in the representation of the snowpack and the soil: from a single-layer parameterization for snow (Douville et al., 1995) and a force–restore scheme for soil (Noilhan and Mahfouf, 1996) to explicit multi-layer modules for both snow (Boone and Etchevers, 2001) and soil (Boone et al., 2000; Decharme et al., 2011). In addition to improving the individual components of the model (soil and snow schemes), we test the use of multiple patches (i.e., a “tiling approach”) to divide the energy balance by surface types, which remains required even for kilometer-scale modeling systems, and address limitations of subgrid parameterizations, such as the partial snow cover fraction approach when only one soil column is used for both covered and uncovered snow parts. As these approaches are common for many land surface models (LSMs) used in coupled systems – e.g., HTESSEL (Balsamo et al., 2009), NOAH-MP (Niu et al., 2011), CLM5 (Lawrence et al., 2019), JULES (Best et al., 2011) – our study may provide information on the necessary content of surface configurations to correctly represent snow cover in mountainous regions in a high-resolution coupled surface–atmosphere context. Moreover, the factors proposed in the study to explain the erroneous representation of the snowpack in CNRM-AROME are strongly suspected to contribute to the shortcomings in the representation of seasonal snow cover documented in coarser-resolution coupled simulations using the SURFEX-ISBA LSM, such as CNRM-ALADIN (Termonia et al., 2018) in the Alps (Monteiro and Morin, 2023) and CNRM-CM6 in high-latitude boreal forests (Decharme et al., 2019).

The results of the experiments introduced in the present study are analyzed and compared to different sets of observational data, enabling us to assess the impact of the modifications in complementary ways:

• comparisons of snow depth values on a large set of in situ measurements collected and presented in Matiu et al. (2021a), enabling a quantitative analysis on a broad spatial scale;

• comparisons with MODIS snow durations, providing near-exhaustive spatial coverage of time-aggregated information on snowpack conditions.

In the end, we introduce a SURFEX-ISBA configuration that is relevant for coupled surface–atmosphere modeling and allows for a significant improvement in the representation of mountain snow cover.

2 Materials and methods

2.1 CNRM-AROME model

In this study, simulations are carried out using the CNRM-AROME climate model, which is the convection-permitting regional climate model (CP-RCM) used at CNRM, which includes the surface model SURFEX (Masson et al., 2013) coupled with the AROME atmospheric model. CNRM-AROME is directly based on the non-hydrostatic limited-area model AROME that has been used for NWP at Météo-France since 2008 (Seity et al., 2011; Brousseau et al., 2016). An alternative version of the AROME model, referred to as HARMONIE-AROME (Bengtsson et al., 2017), is used in NWP applications by several European meteorological services and also used for climate studies by the HARMONIE Climate community (Belušić et al., 2020; Lind et al., 2020).

In this study, the CNRM-AROME model is based on NWP AROME cycle 46t1 in operational use at Météo-France since 2022, operated for climate simulations at a horizontal resolution of 2.5 km, with 60 vertical levels. The time step of the model is 60 s. This version has a lot in common with cycle 41t1, used for the CNRM-AROME climate simulations carried out as part of FPS convection of CORDEX (Coppola et al., 2020; Pichelli et al., 2021). Detailed information about its atmospheric and surface configuration can be found in Termonia et al. (2018) and Caillaud et al. (2021). The main difference between cycles 41t1 and 46t1 relevant to our study is the use of a more recent version of SURFEX (version 8.0).

2.2 SURFEX: the surface platform

For this study, the surface modeling is ensured by the surface platform SURFEX v8.0 (Masson et al., 2013). Within SURFEX, the estimation of energy and mass fluxes of each grid cell is carried out by specific modules depending on the type of surface environments called tiles. Four distinct such environments are accounted for in SURFEX:

• nature tile – “natural” continental surfaces (i.e., including bare soil, rocky ground, permanent snow, glaciers, and natural and cultivated vegetation) using the ISBA land surface model (LSM) (Noilhan and Planton, 1989; Noilhan and Mahfouf, 1996);

• town tile – urban environments using the town energy balance (TEB) module (Masson, 2000);

• lake tile – continental waterbodies such as lakes and rivers using the Charnock formulation (Charnock, 1955);

• sea tile – seas and ocean using version 6 of ECUME (Belamari and Pirani, 2007).

NATURE land surfaces modeling is carried out by the LSM ISBA, representing the evolution of soil and vegetation biophysical variables, including the snowpack, either parameterized or explicitly represented.

2.2.1 ISBA LSM: main principles and identified weaknesses/flaws for snow representation

Three main different land surface configurations are described and analyzed in this study. Despite their differences, the calculation of the surface energy balance and the parameterization of the snow fraction are identical and play a major role in the seasonal evolution of the snowpack.

Surface energy balance

The surface energy balance is computed for a surface layer with a fixed depth of 0.01 m, which is a composite representation of the soil–vegetation system (soil–vegetation–snow in the case of the approach using the D95 single-layer snow parameterization; Douville et al., 1995).

A single surface temperature, T_s, is calculated for each grid cell, whose evolution depends on the surface heat flux into the composite layer, G; the heat flux between the surface and the soil, F_surface-soil, for which the formulation depends on the soil scheme used; and the heat flux between the surface and the snowpack, F_surface-snow, in the case of the use of an explicit snow model. Thus, the time evolution of the surface temperature, T_s, is expressed as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}{T}_{\mathrm{s}}}{\mathrm{d}t}}=C_t\times G-F_{\text{surface-soil}}-F_{\text{surface-snow}},\end{array}$$

with G (W m⁻²) being the surface heat flux between the atmosphere and the soil–vegetation composite layer given by

$$\begin{array}{}\text{(2)}& G=R_n-H-\mathit{LE},\end{array}$$

resulting from the evolution of the radiation balance, R_n, and the sensible, H, and latent, LE, heat fluxes, weighted by C_t, a composite coefficient accounting for the heat capacity of the surface layer, whose formulation depends on the soil and snow schemes used.

The radiation balance is the cumulated difference (W m⁻²) between the incoming solar radiation, SWd, and the infrared atmospherical radiation, LWd, and the reflected shortwave radiation, SWu, and emitted longwave radiation, LWu, expressed as follows:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& R_n=\mathrm{SWd}+{\mathit{\epsilon }}_{\mathrm{s}}\mathrm{LWd}\mathrm{-}\mathrm{SWu}-\mathrm{LWu},\\ & \mathrm{SWu}=\mathrm{SWd}\times {\mathit{\alpha }}_{\mathrm{s}},\\ & \mathrm{LWu}={\mathit{\epsilon }}_{\mathrm{s}}\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}},\end{array}\end{array}$$

with α_s and ε_s, respectively, being the surface albedo and emissivity and σ the Stefan–Boltzmann constant.

The turbulent fluxes are computed by means of the bulk aerodynamic formulae defined by Louis (1979) and modified by Mascart et al. (1995) to account for different roughness length values for heat and momentum.

The patch approach

In order to take into account the heterogeneity of the land surface within the NATURE tile of each model grid cell, ISBA offers the possibility of splitting the calculation of energy balances by surface types. A total of 19 surface types (called patches) are available, dividing natural surfaces into soil and vegetation categories with distinct physical characteristics. The nomenclature and categorization of the 19 patches are taken from the ECOCLIMAP physiographic database (Faroux et al., 2013) and correspond to the plant functional types (PFTs) of ECOCLIMAP. In this study, we use ECOCLIMAP version 1 (Masson et al., 2003).

The number of patches is set by the user, with a number ranging from 1 to 19. When fewer than 19 patches are used, the physical characteristics of multiple land surfaces are aggregated by grouping them by categories and weighted by their respective fractions within the cell while following the aggregation laws defined by Noilhan et al. (1995) and Noilhan et al. (1997) (e.g., logarithmic for the roughness length; linear for the albedo, the leaf area index (LAI), and the vegetation fraction; and inverse for the stomatal resistance).

For a given grid cell, the atmospheric fluxes received are thus identical for all tiles and patches, but a specific energy and mass balance is calculated for each of the patches. There are no energy and mass exchanges between the soil–snow columns of the different patches. The fluxes for each of the patches are then aggregated by weighting the relative fraction of each type of surface within the grid cell, enabling the estimation of average fluxes for all the natural surface types in the grid cell, which are provided to the atmospheric model or used as diagnostics for each grid point.

Parameterization of the snow cover fraction

The presence of snow on the ground has a major impact on the surface mass and energy balance in several ways. As the snow cover extends, the albedo of the surface increases, its roughness decreases, and the snowpack insulates the underlying ground from heat and mass exchanges with the atmosphere. The way in which the fraction of the grid cell covered by snow is calculated and influences the computation of the energy balance is therefore critical and is represented in widely different ways in different land surface models (Essery et al., 2013; Menard et al., 2020; Lalande et al., 2023).

In ISBA, for each patch, the snow fraction is calculated differently depending on whether vegetation is present or not, and the fraction is used for energy balance calculations.

In the absence of vegetation, even a small amount of snow covers the entire surface. This is represented in ISBA by the fact that the snow cover fraction in non-vegetated areas (P_sng) reaches 1 as soon as the snow water equivalent, W_s, exceeds a threshold value set to $W_{{\mathrm{s}}_{\mathrm{crit}}}=\mathrm{10}$ kg m⁻² (Eq. 4). Thus, the snow fraction over ground, P_sng, is expressed as follows:

$$\begin{array}{}\text{(4)}& P_{\mathrm{sng}}=min\left(\mathrm{1},{\displaystyle \frac{W_{\mathrm{s}}}{W_{{\mathrm{s}}_{\mathrm{crit}}}}}\right).\end{array}$$

In the presence of vegetation, the calculation of the snow cover fraction (P_snv) is done based on the snow depth (also referred to as the height of the snow) value, H_s, and takes into account the height of the vegetation through the roughness length, z₀. W_sn is a scaling factor, modulating the weight of vegetation height in the calculation of the snow cover fraction (Eq. 5). P_snv is formulated as follows:

$$\begin{array}{}\text{(5)}& P_{\mathrm{snv}}={\displaystyle \frac{H_{\mathrm{s}}}{\left(H_{\mathrm{s}}+W_{\mathrm{sn}}{z}_{\mathrm{0}}\right)}}.\end{array}$$

The total snow cover fraction, P_sn is the sum of the snow fractions for each patch weighted by their respective fraction (Eq. 6):

$$\begin{array}{}\text{(6)}& P_{\mathrm{sn}}=(\mathrm{1}-\mathrm{Veg})\times P_{\mathrm{sng}}+\mathrm{Veg}\times P_{\mathrm{snv}}.\end{array}$$

Snow-related prognostic variables are defined for each patch. Integrated diagnostics for the NATURE tile of each grid cell are computed as the weighted average using the patch fractions.

2.3 Land surface configurations

The goal of the study is to describe and evaluate new land surface configurations in order to improve the representation of seasonal snow cover in the European Alps and address some of the issues identified in Monteiro et al. (2022). Consequently, the atmospheric configurations and initialization of all experiments are similar, and we explore the impacts of changes in surface configuration mostly on the simulated snowpack. Also note that part of the content of the configurations tested here were already used in the latest version of the general circulation model (GCM) CNRM-CM6 (Decharme et al., 2019) and the RCM CNRM-ALADIN63 (Nabat et al., 2020) but had not been used in coupled model simulations using AROME.

For all configurations tested in this study, including the configuration referred to as the default one, we activate the option described in Decharme et al. (2016) that limits snow accumulation above a certain snow depth threshold (see Decharme et al., 2016). Its value is set to the default value of 33.3 m. This option is activated in all experiments and avoids the formation of the problem of snow towers identified in Monteiro et al. (2022).

Figure 1 illustrates the main characteristics of the three main surface configurations used in this study. The configurations are described in detail in Sect. 2.3.1, 2.3.2, and 2.3.3, and Table 1 summarizes their main model components.

(Boone et al., 1999)(Boone et al., 2000; Decharme et al., 2011)(Boone et al., 2000; Decharme et al., 2011)(Douville et al., 1995)(Boone and Etchevers, 2001)(Boone and Etchevers, 2001)

Table 1 Main model components for each configuration.

Figure 1 Schematic illustration of the main physical processes, flux exchanges, and prognostic variables for the three main configurations documented and tested in this study. D95-3L configuration (framed in orange), ES-DIF (framed in blue), and ES-DIF-OPT (framed in green), with the different modifications displayed and schematically illustrated. H and LE: sensible and latent heat fluxes, respectively. ROS: rain on snow. F_x: fluxes from component x. W_xn and WI_xn: respectively, the liquid and ice content of the nth layer of component x, with s for surface, g for ground, and sn for snow. W_sn: snow water equivalent. α_sn: snow albedo. ρ_sn: snow density. H_sn: snow enthalpy. T_sn: snow temperature. WI_sn: liquid content of the snow.

2.3.1 D95-3L: one-layer snow parameterization and force–restore approach for the soil

This surface configuration is the default one as it is currently in use for NWP version of AROME and CNRM-AROME for climate studies (Caillaud et al., 2021; Lucas-Picher et al., 2023; Monteiro et al., 2022). It is schematically described in Fig. 1. The evolution of biophysical soil variables is ensured by the “3-L” soil model, using a force–restore approach. Heat exchanges in the ground (temperature evolution) are represented using two layers (Noilhan and Mahfouf, 1996), and three layers are used for the evolution of hydrological variables (Boone et al., 1999). In this configuration, the snowpack is parameterized as a single layer with homogeneous physical properties, referred to as the D95 parameterization (Douville et al., 1995) and mixed with the soil–vegetation composite surface layer. Consequently, no specific energy balance is solved for the snowpack, which is taken into account by modifying the properties of the composite surface layer. The main equations governing the evolution of the surface components are given below:

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{d}{T}_{\mathrm{s}}}{\mathrm{d}t}}=C_t\cdot G-F_{\text{surface-soil}},\text{(8)}& {\displaystyle }{C}_t={\displaystyle \frac{\mathrm{1}}{\frac{(\mathrm{1}-\mathrm{Veg})\left(\mathrm{1}-P_{\mathrm{sng}}\right)}{C_{\mathrm{g}}}+\frac{\mathrm{Veg}(\mathrm{1}-P_{\mathrm{snv}})}{C_{\mathrm{v}}}+\frac{P_{\mathrm{sn}}}{C_{\mathrm{s}}}}},\text{(9)}& {\displaystyle }{F}_{\text{surface-soil}}={\displaystyle \frac{\mathrm{2}\mathit{\pi }}{\mathit{\tau }}}\left(T_{\mathrm{s}}-T_{\mathrm{g}\mathrm{2}}\right),\end{array}$$

with T_g2 being the temperature of the deep soil layer (which evolves through a relaxation term towards T_s) and C_g, C_v, and C_s, respectively, the heat capacity of the ground, vegetation, and snow.

Three prognostic variables characterize the snowpack:

• The first one is the density, an exponentially decreasing function, forced to 100 kg m⁻³ for fresh snow, limited to 300 kg m⁻³ for aged snow. The density of the entire snowpack is updated during snowfall by a weighted average of the previously present layer and that of the snow newly fallen to the ground.

• The second one is albedo, whose evolution can follow two functions, forced to 0.85 for fresh snow and limited to 0.5 for old snow, and linearly decreasing in the absence of melting and exponentially decreasing in the presence of melting (i.e., to account for wet metamorphism).

• Thirdly, the snow water equivalent (total mass) of the snowpack results from a mass balance calculation depending on snowfall, snow sublimation/evaporation, and melting.

In this configuration, the snow layer has no prognostic temperature of its own but is included in the composite soil–vegetation–snow surface layer from which the melting temperature, i.e., the temperature value used to compute the melt intensity, is derived. In the presence of vegetation, the snowmelt intensity is calculated based on a hybrid diagnostic temperature, a weighted average between the surface soil temperature and the deep soil layer temperature, with a value closer to the deep soil layer temperature as the proportion of vegetation increases (see Eq. 10).

$$\begin{array}{}\text{(10)}& T_{\mathrm{melt}}=(\mathrm{1}-\mathrm{Veg})T_{\mathrm{s}}+\mathrm{Veg}{T}_{\mathrm{p}},\end{array}$$

with T_melt being the melting temperature, T_s the instantaneous surface temperature, T_p the daily mean surface temperature, and Veg the fraction of vegetation within the grid cell.

This approach was developed to prevent unrealistic snowpack melting. Indeed, using the instantaneous value of the surface, the temperature representative of the soil–vegetation–snow system tends to be too high during daytime (i.e., due to the mixed albedo between snow and vegetation), leading to spurious snowmelt computations (Douville et al., 1995). As shown later, this approach has strong consequences for the modeling of snow conditions in forested environments.

2.3.2 ES-DIF: multi-layer snow scheme and multi-layer soil scheme

This approach uses intermediate complexity schemes for soil and snow in a multi-layer manner, allowing the resolution of specific energy balances for the soil–vegetation system and for snow as well as a more detailed representation of the processes within them. These are the schemes currently used in recent versions of the CNRM-CM6 global model (Voldoire et al., 2019; Decharme et al., 2019), the CNRM-ALADIN regional model (Nabat et al., 2020), and the most recent version of the HARMONIE Climate AROME regional climate model (Belušić et al., 2020; Lind et al., 2020). However, note that only one patch is used herein for the NATURE tile, which is not the way the configuration is implemented for the coupled systems CNRM-CM6 and CNRM-ALADIN using 12 patches and HARMONIE Climate using two patches. This configuration is illustrated in Fig. 1.

Heat and mass exchanges within the soil are computed using the ISBA diffusion scheme (ISBA-DIF; Boone et al., 2000; Decharme et al., 2011), with 14 layers from the surface to 12 m, representing explicitly heat exchanges within the different soil layers through the resolution of a 1D Fourier law. In this scheme, a single surface temperature, T_s, is calculated for the soil–vegetation system, whose evolution depends on the surface energy balance (G; Eq. 2); the surface–soil heat flux, F_surface-soil, with the second soil layer on the non-snow-covered part; and the surface–snow heat flux at the surface–snow interface on the snow-covered part.

The main equations are provided below:

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\mathrm{d}{T}_{\mathrm{s}}}{\mathrm{d}t}}=C_t\times G-F_{\text{surface-soil}}-F_{\text{surface-snow}},\text{(12)}& {\displaystyle }{C}_t={\displaystyle \frac{\mathrm{1}}{\frac{(\mathrm{1}-\mathrm{Veg})(\mathrm{1}-P_{\mathrm{sng}})}{C_{\mathrm{g}}}+\frac{\mathrm{Veg}\left(\mathrm{1}-P_{\mathrm{snv}}\right)}{C_{\mathrm{v}}}}},\text{(13)}& {\displaystyle }{F}_{\text{surface-soil}}=C_{\mathrm{g}}{\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_{\mathrm{1}}}{\mathrm{\Delta }{\stackrel{\mathrm{‾}}{z}}_{\mathrm{1}}}}\left(T_{\mathrm{s}}-T_{\mathrm{g}\mathrm{2}}\right),\end{array}$$

with $\stackrel{\mathrm{‾}}{{\mathit{\lambda }}_{\mathrm{1}}}$ (W m⁻¹ K⁻¹) being the inverse-weighted arithmetic mean of the soil thermal conductivity at the interface between the surface layer and the underlying soil layer and $\mathrm{\Delta }\stackrel{\mathrm{‾}}{z_{\mathrm{1}}}$ the thickness (m) between the two consecutive layer mid-points. F_surface-snow is a heat conduction term between the lowermost snow layer and the soil surface layer.

The snowpack evolution is carried out by the Explicit Snow scheme (ISBA-ES) (Boone and Etchevers, 2001; Decharme et al., 2016) using up to 12 snow layers for which a specific energy balance is solved unlike in D95, which computes a single energy balance for the composite soil–vegetation surface layer.

Three prognostic variables are used to describe the state of each layer at each time step:

• heat content (i.e., temperature and water/ice content), which defines the energy required to melt the layer and thus combines the information of snow temperature and liquid water content at melting point;

• density, which evolves under the effect of parameterized compaction and metamorphism (Brun et al., 1989), wind-induced densification of near-surface snow layers, and fresh snowfall (whose density is a function of air temperature and wind at the time of fall);

• the thickness of each layer, ranging from a few millimeters to several tens of centimeters, defined to be the finest close to the ground–snow and atmosphere–snow interfaces (see Decharme et al., 2016, for more details).

One additional prognostic variable for the surface layer is the albedo. As stated by Boone and Etchevers (2001), snow albedo follows a linear decrease rate for dry snow (Baker et al., 1990) and an exponential decrease rate to model the wet metamorphism (Verseghy, 1991).

The mass balance of the snowpack is expressed as the sum of precipitation on snow (solid and liquid since each layer can have a liquid water content), evaporation, and sublimation as well as a term describing the flow of water out of the snowpack at its base.

ISBA-ES includes a number of parameterizations that reproduce the effects of physical processes affecting the evolution of the snowpack:

• compaction, metamorphism, and wind-induced densification (Brun et al., 1989);

• transmission of incident solar flux through the snow layers (Brun et al., 1992);

• water percolation between snow layers;

• refreezing and melting of water contained in snow layers.

The temperature of all snow layers is computed simultaneously, following the system of equations below:

$$\begin{array}{}\text{(14)}& C_{{\mathrm{s}}_i}{D}_i{\displaystyle \frac{\mathrm{d}{T}_{{\mathrm{sn}}_i}}{\mathrm{d}t}}={\mathrm{Gs}}_{i-\mathrm{1}}-{\mathrm{Gs}}_i+{\mathrm{Rs}}_{i-\mathrm{1}}-{\mathrm{Rs}}_i-{\mathrm{Ss}}_i,\end{array}$$

with, for each layer i, Ss_i representing the heat sink/source linked to water phase changes, Rs_i the incident solar radiation transmitted (decreasing exponentially with distance to the snow surface), Gs_i the layer energy balance, and Gs_i−1 the energy balance of the layer above. For the layers below the surface, the Gs_i term corresponds to thermal diffusion in snow, while for the uppermost layer, the energy balance includes the following terms:

$$\begin{array}{}\text{(15)}& {\mathrm{Gs}}_{\mathrm{0}}=R_{\mathrm{ns}}-H-\mathit{LE}-C_{\mathrm{w}}{P}_{\mathrm{sn}}\left(P-P_{\mathrm{s}}\right)\left(T_{\mathrm{f}}-T_{\mathrm{r}}\right),\end{array}$$

with R_ns being the snow surface radiation balance and $H\left(T_{{\mathrm{sn}}_{\mathrm{1}}}\right)$ and $\mathit{LE}\left(T_{{\mathrm{sn}}_{\mathrm{1}}}\right)$ the turbulent fluxes above snow (calculated according to Louis, 1979, formulae) and a latent heat source term related to the fall of liquid precipitation in the snowpack, where C_w is the heat capacity of water; P and P_s the total and solid precipitation, respectively; T_r the temperature of the rain; and T_f the fusion temperature. Any excess heating of snow temperature above the freezing point is converted into energy available for melting. Then, the liquid water percolation follows a bucket scheme based on a liquid water retention capacity and accounting for possible refreezing in colder layers.

2.3.3 ES-DIF-OPT: multi-layer snow scheme and multi-layer soil scheme, including optimal modifications

ES-DIF-OPT stands for optimized ES-DIF configuration. Starting from the second configuration (ES-DIF), we add a series of modifications (see Fig. 1) concerning the use of multiple patches, changes in some parameterizations, input physiographic databases, and calculation of heat and mass exchanges.

3-PATCHS

The 3-PATCHS modification consists of activating three patches for energy and mass balance calculations in contrast to the D95-3L and ES-DIF configurations, which only use one patch (see Sect. 2.2.3 and illustration in Fig. 1c for further details). In the most recent version of the HARMONIE Climate AROME model (2.5 km horizontal resolution) (Belušić et al., 2020), two patches have been activated in order to distinguish between forested and open-land areas, while 12 patches are used in the latest versions of CNRM-CM6 (150 km horizontal resolution) (Decharme et al., 2019) and CNRM-ALADIN (12.5 km horizontal resolution) (Nabat et al., 2020). When three patches are used, the categories are grouped into “uncovered surface” (e.g., permanent snow, rock, and bare soil), “low vegetation” (shrubs, grass, and crops) and “high vegetation” (various types of trees), allowing for a clear distinction between vegetated and non-vegetated surface types. While the number of patches used can be as high as 19, here we activate three patches as a compromise in order to avoid increasing the computational cost and storage burden of the land surface modeling within the CNRM-AROME modeling system too much.

GFLUX

The GFLUX modification consists of reducing the heat flux between the soil surface and the base of the snow. It is designed to reduce the unrealistic soil–snow conduction heat flux due to the unrealistic assumption of an identical soil physical state between the fractions of the patch that are covered by snow and those that are not. For this configuration, the thermal conductivity of the interface between the lowermost snow layer and the uppermost soil layer, calculated as the harmonic average of the conductivity of each layer, is capped at 5 % of its value below a snow fraction of 75 %, increasing linearly and reaching its base value when the snow fraction reaches 100 %.

The thermal conductivity of the soil–snow interface Csng is thus computed as

$$\begin{array}{}\text{(16)}& \mathrm{Csng}=\mathrm{Csng}\times max(\mathrm{0.05},\mathrm{3.8}{P}_{\mathrm{sn}}-\mathrm{2.8}).\end{array}$$

WSN-1

The modification WSN-1 consists of adjusting the parameter governing the estimate of the snow cover fraction on vegetation. In this case, the value of W_sn in the formula for snow fraction on vegetation (see Eq. 5) is lowered from 5 to 1. As illustrated in Fig. A2 in Appendix A, this modification increases the sensitivity of snow fraction to snow depth, allowing it to reach higher values even with moderate amounts of snow. The motivation for this modification is similar to that explained for the GFLUX modification but achieved through the reduction in the range of snow depth values with intermediate snow-covered fraction values.

SG-LSOC

The modification SG-LSOC refers to the use of the SoilGrids v2.0 database (Poggio et al., 2021) for soil textures (proportion of sand and clay) and the activation of the soil organic carbon parameterization effect (Decharme et al., 2016) on soil heat and mass exchanges. The use of SoilGrids v2.0 is motivated by its better estimate of the soil organic carbon stock than the HWSD database (Batjes, 2016) over France and boreal regions (Tifafi et al., 2018).

Experimental design

All these modifications have been defined and tested iteratively with the aim of improving the seasonal dynamic of snow cover in CNRM-AROME simulations over the European Alps. Only one major modification was done per experiment, with the aim of moving towards a more physical configuration and/or resolving remaining issues without overtuning because the land surface model is not the only cause of errors in regional climate modeling. Modifications that consistently reduced discrepancies with the observations used as a reference were retained for the following experiment, reaching an optimum configuration for simulating snowpack in the European Alps. For clarity and brevity, only the three configurations described above are shown in the main figures of the article. Further results with intermediate configurations are provided in the Appendix.

2.4 Geographical domain and simulations setup

The CNRM-AROME simulations were performed at 2.5 km horizontal resolution over a domain that covers the alpine ridge, shown in Fig. 2, namely the ALP-3 domain, that is the domain used for the CORDEX FPS Convection (Coppola et al., 2020; Ban et al., 2021; Pichelli et al., 2021; Caillaud et al., 2021; Monteiro et al., 2022). The black highlighted contour of the map defines the mountainous region of the Alps in which all the analyses were performed and which corresponds to the boundaries of the Alpine Convention domain (Convention, 2020). Simulations were driven by atmospheric fields directly from the ERA5 (Hersbach et al., 2020) reanalysis at 50 km horizontal resolution each hour, thanks to the increasing resolution of global reanalyses. This is the first time using the CNRM-AROME climate model without the need for an intermediate RCM to downscale ERA5 fields. The simulations cover the 3-year time period from 1 January 2018 to 31 December 2020, for which we had the largest number of available observations. The CNRM-AROME atmosphere was initialized using ERA5 fields interpolated on the ALP-3 domain, and the surface fields were initialized by realizing one time step with these interpolated atmospherical fields, with both happening on 1 January 2018. A dedicated study was carried out to analyze the impact of the absence of spin-up on the simulation results; see Sect. 2.5. Note that the simulation results were evaluated over the seasons 2018–2019 and 2019–2020, i.e., from September to August of the following years. The first 8 months (and last 4 months) of the coupled simulations were therefore not used for the evaluation.

Figure 2 Simulation setup displaying the whole domain of simulation and the contour of the Alpine Convention outline of the Alps within which the evaluation was carried out, with the orography of CNRM-AROME at 2.5 km horizontal resolution.

2.5 Impact of the initialization (spin-up) approach on snowpack simulations

The high computational cost of coupled surface–atmosphere simulations at 2.5 km horizontal resolution over the whole Alpine ridge prevented us to perform a fully fledged coupled spin-up for all the configurations. We nevertheless assessed the impact of the reduced spin-up procedure that we implemented, which is potentially insufficiently long to obtain a balanced ground state. Indeed, albeit soil heat and water content are known to have a relaxation time ranging from a few years to a decade (Christensen, 1999; Cosgrove et al., 2003), performing our experiments in a transient regime for the soil could affect the results of our study. As reiterated by Jerez et al. (2020), the impact of the spin-up is largely related to the goal of the study (i.e., variables of interest, magnitude of the investigated changes in a comparison, etc.). In our case, we find that the order of magnitude of the changes in the surface configurations we performed is substantially larger than the impact of an unbalanced deep soil heat and water content over snowpack simulations. Appendix B provides comparisons of model runs using the ES-DIF-OPT configuration using either the default initialization procedure (described above) or an initialization obtained from 13 years of standalone (offline) simulations of the surface model from 1 January 2006 to 1 January 2018 driven by near-surface atmospheric fields from the CNRM-AROME coupled model run that is itself driven by the ERA-Interim–ALADIN model pair (Caillaud et al., 2021; Monteiro et al., 2022). Figure B1 in Appendix B confirms that at the initialization date (i.e., 1 January 2018), the default initialization strongly underestimates snow amounts and provides significantly too wet and warm soil conditions over most of the European Alps compared to the result of multiple years of offline simulations. Nonetheless, the snow depth time series at four sites Fig. B2 in Appendix B show that, after the first 6 months, the effect of the spin-up is negligible in comparison, with respect to the objectives of our study.

2.6 Observational references

Various observational datasets taken as reference are used to analyze different aspects of the impact of the choice of the land surface model configurations on simulated snow and atmospheric surface variables and are described in the following subsections.

2.6.1 In situ snow depth observations

A set of daily in situ snow depth observations is employed to perform an extensive point-scale evaluation of the simulated snow depth values over the 2018–2019 season (i.e., from 1 September 2018 to 31 August 2019) as it was the season for which the largest number of observations was available based on existing consolidated datasets described in Matiu et al. (2021a). Figure 3 shows the location of the in situ measurements and their distribution with respect to elevation, spanning the whole Alpine ridge over elevations ranging from 0 to 3000 m.

The observational time series used in this study are quality-checked. The greatest part of the dataset was gathered and described by Matiu et al. (2021a), to which we added Austrian stations from the Hydrographic Central Office of Austria (HZB) and GeoSphere AT (i.e., the TAWES and SNOWPAT datasets) and Swiss stations (i.e., the IMIS datasets; Measurement and IMIS, 2023) from the WSL Institute for Snow and Avalanche Research (SLF).

From this large set of in situ snow depth measurements (i.e., 1005 stations), 266 stations were selected according to the following criteria:

• The AROME grid point that includes it has less than 150 m difference with the station elevation.

• The AROME grid point that includes it is filled by less than 75 % of the high-vegetation cover type.

These criteria limit some of the representativeness issues of a point-scale comparison between a local in situ station and a model grid point representative of a square of 2.5 × 2.5 km in mountainous regions.

Figure 3 Location of the in situ snow depth observations with their associated number per bin of 300 m width elevation band. The dotted black line shows the elevation distribution using a reference digital elevation model at 100 m horizontal resolution (E.E.A., 2016).

2.6.2 Satellite (MODIS) snow cover duration

A large-scale evaluation of the snow cover duration (SCD, defined as the longest consecutive period with snow on the ground based on hydrological years from September to August) was performed using the MODIS/Terra daily normalized difference snow index (NDSI) fields at 500 m for the 2018–2019 and 2019–2020 seasons. These data from the MODIS/Terra sensor have been treated by the National Snow and Ice Data Center (NSIDC) (Hall and Riggs, 2020) and correspond to a daily gap-filled product using an algorithm described in Hall et al. (2010). In this study, MODIS NDSI data were regridded to match the CNRM-AROME horizontal resolutions of 2.5 km using a first-order conservative method.

Figure 4 shows the SCD over the 2018–2019 and 2019–2020 seasons from MODIS over our area of interest regridded at 2.5 km. The MODIS SCD is calculated upon the MODIS NDSI by converting it to a series of binary snow cover maps (absence or presence of snow) using a threshold value of NDSI >0.2. This threshold corresponds to a snow cover fraction of approximately 30 % (Salomonson and Appel, 2004). In this study, the CNRM-AROME SCD was computed using snow depth values with a threshold set to 1 cm, motivated by the minimization of error metrics as described in Monteiro and Morin (2023).

2.7 Point-scale comparison, elevation bands analyses, and used statistics

2.7.1 Point-scale comparison

Section 3.1 and Appendix D introduce point-scale comparisons between individual station measurements and the corresponding CNRM-AROME grid cell. It means that each station is compared to the CNRM-AROME grid cell representative of a 2.5 × 2.5 km square including the station location based on its geographical coordinates.

2.7.2 Elevation bands analyses

Section 3.1 and 3.2 introduce analyses performed using an elevation-based categorization. Here, we used 300 m width elevation bands, meaning that for a given elevation band at median elevation z, all stations or grid points with an elevation ranging between z ± 150 m are gathered and used. This choice is a trade-off between the heterogeneity within an elevation band and the inclusion of a maximum of grid points or observations within.

For clarity and brevity, only results for four elevation bands are presented in the main article figures, representing distinct environments: 900 m ± 150 m for the valleys and low-elevation hills, 1500 m ± 150 m and 2100 m ± 150 m for the intermediate elevation, and 2700 m ± 150 m for the high-mountain conditions. The results for the other elevation bands, not shown, are consistent with the main patterns observed across the analyzed elevation bands.

Figure 4 (a) MODIS snow cover duration (SCD) over two seasons (2018–2019 and 2019–2020) within the contour of the Alpine Convention, regridded at 2.5 km horizontal resolution over the CNRM-AROME mesh grid. (b) Box plot representing the spatial distribution of the snow cover duration values for two seasons (2018–2019 and 2019–2020) for multiple elevation bands in the European Alps.

2.7.3 Surface type analyses

The evaluation of the snow cover duration using MODIS remote sensing data is complemented by a categorical analysis by surface type. Figure 5 shows the location and the elevational distribution of points per prevailing surface type (i.e., meaning that the surface type represents more than 75 % of the cover; otherwise it is classified as mixed) for the CNRM-AROME mesh grid. The classification per vegetation type is based on the ECOCLIMAP land-use database, from which the 19 vegetation types have been gathered into three categories: no vegetation (i.e., bare ground, rock, and permanent snow and ice), high vegetation (i.e., grouping all types of high trees), and low vegetation (i.e., crops, grasslands, and shrubs).

2.7.4 Statistics

The error metrics used are defined as follows:

• mean error (ME), i.e., $\mathrm{ME}=\frac{{\sum }_{i=\mathrm{1}}^N\left(x_i-y_i\right)}{N}$;

• correlation (Pearson linear correlation), i.e., $R=r_{xy}=\frac{\sum x_i{y}_i-N\stackrel{\mathrm{‾}}{x}\stackrel{\mathrm{‾}}{y}}{\sqrt{\left(\sum x_i^{\mathrm{2}}-N{\stackrel{\mathrm{‾}}{x}}^{\mathrm{2}}\right)}\sqrt{\left(\sum y_i^{\mathrm{2}}-N{\stackrel{\mathrm{‾}}{y}}^{\mathrm{2}}\right)}}$.

Here, x_i and y_i are data x and y at time i, $\stackrel{\mathrm{‾}}{x}$ and $\stackrel{\mathrm{‾}}{y}$ the mean of x and y, and N the sample size.

3 Results

The presentation of the results is first performed by comparing simulation results with a large sample of in situ snow depth measurements covering the European Alps during the 2018–2019 season. We then evaluate the simulation results in terms of snow cover duration compared to remotely sensed (MODIS) observations for the two winters, 2018–2019 and 2019–2020.

3.1 Point-scale evaluation of snow depth values over the 2018–2019 winter

Figure 6 shows the multi-station mean time series of the height of snow of the three land surface configurations and the reference observations for multiple elevations bands of 300 m width ranging from 900 m ± 150 m to 2700 m ± 150 m over the 2018-2019 winter. For each of the elevation bands and configurations the mean errors (MEs) and the Pearson correlation (R²) calculated over multi-station mean time series are displayed.

The D95-3L configuration simulates snowpack evolution with similar behaviors for all elevation bands. During the accumulation period and until the observed annual snow depth maximum, the simulated multi-station average values remain close to the observations, with highly correlated variations. Nevertheless, after the observed annual snow depth maximum, deviations from the observed multi-station average values widen, mainly due to the underestimation of melt events (i.e., their frequency and amplitude) in the simulations compared to observations, resulting in less correlated snow depth time series. This leads to a significantly delayed and higher annual snow depth maximum (i.e., up to 2 months at intermediate and high elevations) and, to a lesser extent, a delayed end of snow season (i.e., from a few days to a month, partly compensated by faster melting) in the simulation. As stated in introduction, these overestimations of the annual snow depth maximum in its value and its timing are in line with previous studies working with CNRM-AROME climate simulations (Monteiro et al., 2022; Lucas-Picher et al., 2023; Monteiro and Morin, 2023). Overall, the simulated multi-station average values show high Pearson correlation scores (R²) for the whole snow season, from 0.81 to 0.97, but mean errors (MEs) can be large, from 2 cm at 900 m to 30 cm at 2700 m (i.e., around 10 %–15 % of the annual snow depth maximum).

Figure 5 Location of the prevailing surface type (i.e., if the surface type represents more than 75 % of the cover; otherwise, is it classified as mixed) for the CNRM-AROME mesh grid within the contour of the Alpine Convention. The horizontal bar plot on the right shows the number of grid points per bin of 300 m width elevation band classified by prevailing surface type.

The snow depth values simulated using the ES-DIF configuration in Fig. 6 exhibit large differences with the D95-3L (i.e., the default configuration). During the accumulation period and until the date of the annual snow depth maximum, the simulated multi-station average values follow similar variations to the observations and slightly underestimating the amount of snow at 2100 m and below, slightly overestimating it above. Compared to the D95-3L simulated snow depth, melt events appeared to be captured better (i.e., negative variations in the snow depth are better correlated) all along the snow season but often overestimated, notably below and at 2100 m. As a consequence, even if the timing of the simulated annual snow depth maximum often matches with the observation, its value is strongly underestimated at these elevations (from 20 cm at 900 m to 50 cm at 1500 m and 2100 m), and the time of snow disappearance is too early, from 15 days to a month. The ES-DIF R² value against in situ snow depth observations is not systematically improved and the ME scores are degraded at all elevations except at 2700 m. Indeed, the R² values only significantly increase(i.e., by more than 0.1) at 1500 m compared to the default configuration, and the ME ranges from 2 cm in the D95-3L configuration to −3 cm in the ES-DIF configuration at 900 m, from 14 to −15 cm at 1500 m, and from 19 to −21 cm at 2100 m. This shows that simply using a more complex soil or snow scheme does not warrant improved results compared to a coarser snow or soil model.

The ES-DIF-OPT configuration provides the best estimation of the snow depth values against in situ observations. All along the snow season, as shown in Fig. 6, its variations are almost identical to the ES-DIF snowpack simulation and thus similar to the observed multi-station average value until the annual snow depth maximum. Its similar variations in snow depth with ES-DIF are coherent as both configurations share the same physical basis, but the specific features of the ES-DIF-OPT option seem to attenuate the sensitivity to the melt events, leading to simulated snow depth values closer to the observations during the main melt period after the annual snow depth maximum. Apart from these improvements over the two other configurations, the melt-out date at intermediate elevations (i.e., 1500 and 2100 m) is still too early in the simulations as a result of an overestimated melt notably during springtime. Overall, at all elevations, the R² and the ME scores are improved compared to the other configurations tested. Correlation scores reach very high values, with the lowest values of 0.94 at 2100 m and above 0.98 elsewhere. ME values are also strongly reduced compared to other configurations tested, with values changing from −3 to −0.8 cm at 900 m, from −15 to −5 cm at 1500 m, from −20 to −8 cm at 2100 m, and from −8 to −5 cm at 2700 m between ES-DIF and ES-DIF-OPT, respectively.

Figure 6 Multi-station average time series of snow depth values for the 2018–2019 winter for the four elevation bands centered at 900 m ± 150 m, 1500 m ± 150 m, 2100 m ± 150 m, and 2700 m ± 150 m above sea level. Colored continuous lines correspond to the simulated multi-station average time series for each of the configurations: D95-3L in orange, ES-DIF in blue, and ES-DIF-OPT in green. Black circular markers correspond to the multi-station mean time series of the in situ measurements, with the inter-station standard deviation represented by gray-shaded areas. For each elevation band, the number of stations used to compute the mean and the standard deviation is displayed in blue font. At each elevation band and for all configurations, the correlation (R²) and the mean error (ME) computed using the multi-station time series between the simulated and the in situ measurements are displayed.

3.2 Snow cover duration evaluation using MODIS remote sensing data

In this section, we compare the simulated snow cover duration against MODIS remote sensing data as a reference. The snow cover duration of the two seasons used (2018–2019 and 2019–2020) is averaged for the analysis. Individual seasons show similar differences between each simulation and the reference. Figure 7 shows the differences in terms of the mean snow cover duration (SCD) over two seasons (i.e., 2018–2019 and 2019–2020) between each experiment and the MODIS SCD in the European Alps.

Differences obtained using the D95-3L configuration (Fig. 7a) indicated a widespread overestimation compared to the MODIS SCD. Apart from a few small areas of underestimated SCD evenly distributed over the Alpine ridge, no specific region seems to show more marked differences than others. Conversely, the SCD of the ES-DIF configuration (Fig. 7a) is largely underestimated, with an amplitude that appears to depend on elevation. Indeed, stronger negative values of the difference between simulated and observed SCD (Δ SCD) are found on the outer edge and in the northeastern part of the European Alps, while a few patch of slightly positive values of Δ SCD are located along the ridge. The last configuration, ES-DIF-OPT, provides the best match with observed SCD values, reducing the magnitude of the differences with MODIS SCD at each location of the study area compared to the D95-3L and ES-DIF configurations. The northeastern part of the Alps concentrates more areas with an underestimation of the SCD, while the rest of the Alps presents a widespread slight overestimation.

Figure 7 Snow cover duration differences between the simulation results using the different configurations and MODIS observations in the European Alps. Note that MODIS products that initially have 500 m horizontal resolution have been regridded over CNRM-AROME horizontal resolution (2.5 km) grid using a first-order conservative method. (a) Map of the average differences (mean error) of the snow cover duration (SCD) over two seasons (2018–2019 and 2019–2020) for each configuration compared to MODIS SCD. (b) Box plot representing the spatial distribution of the average differences (mean error) of the SCD over two seasons (2018–2019 and 2019–2020) compared to MODIS SCD for each dataset for the six elevation bands: 900 m ± 150 m, 1500 m ± 150 m, 1800 m ± 150 m, 2100 m ± 150 m, 2400 m ± 150 m, and 2700 m ± 150 m above sea level. Each column corresponds to the values classified by the prevailing type of surface (see Sect. 2.6.2 for details).

The elevational distribution of the differences is represented by box plots in Fig. 7b, categorized by prevailing type of surface based on the ECOCLIMAP I land surface classification used as a physiographic database in the CNRM-AROME simulations (see Sect. 2.2.3). It informs us further on the specific behavior of the simulated behavior of the snow cover regarding the elevation and the surface type.

Looking at the “all” category (i.e., gathering all grid points regardless of their surface type) in Fig. 7b quantitatively confirms what was found on the map (Fig. 7a). The results from the D95-3L and ES-DIF experiments exhibit distributions of SCD differences strongly biased towards an overestimation and an underestimation, respectively, for all elevation bands. The median SCD differences in D95-3L range from +5 to +40 d and the ES-DIF from −5 to −75 d, with both larger at 1500 m elevations than above and below (i.e., in terms of median values and larger in terms of variance). The ES-DIF-OPT configuration is at the center of the other two configurations, showing zero-centered distributions of SCD differences, with greatly reduced variance for most elevation bands. Note that the mixed surface type (i.e., all grid cells with less than 75 % prevalence in each category) shows similar distributions of differences to the all category that gathers all grid points.

While an analysis by elevation alone would lead us to interpreting that greater variations can be found at intermediate elevations, categorization by dominant surface type (i.e., no vegetation, high vegetation, and low vegetation) brings more contrasting results, supporting the hypothesis that the simulation results mainly depend on surface types.

It is on the prevailing no-vegetation surface type that the simulated SCD values show the smallest differences with the MODIS reference values as well as between the different experiments, meaning that the changing surface configuration has only a marginal effect on the simulation of the SCD for this surface type. Indeed, the median Δ SCD values of all experiments combined ranges from −15 to +10 d at most, with this value increasing slightly with elevations for each experiment and only slight improvements in the scores provided by the ES-DIF-OPT configuration.

The highest Δ SCD and the most contrasted behaviors between experiments are found in the high-vegetation surface type. On this type of surface, the D95-3L experiment strongly overestimates the SCD, with median Δ SCD values of +10 d at 900 m, increasing to +50 d at 1500 m and 2100 m. The ES-DIF experiment shows the opposite behavior, with a median Δ SCD of −5 d at 900 m, increasing from −60 to −70 d at 1500 and 2100 m, respectively. For the high-vegetation category, the ES-DIF-OPT configuration brings substantial improvements in scores, with a median difference in the SCD with MODIS of −3 d at 900 m, −1 d at 1500 m, and +10 d at 2100 m.

For the low-vegetation surface type, the simulated SCD values using the D95-3L configuration are overestimated above 900 m compared to observations, with increasing differences with MODIS SCD values, reaching +60 d at high elevations. The ES-DIF Δ SCD values are negative at all elevations but exhibit higher discrepancies at intermediate elevations (i.e., at 1500 and 2100 m), with median values between −30 and −60 d. Again, the ES-DIF-OPT configuration shows the smallest differences, with median values contained in the −10 to +15 d range.

Overall, the analysis of the differences between simulated and observed snow cover duration values (Δ SCD) demonstrates a clear added value of the ES-DIF-OPT configuration, reducing discrepancies across all surface types and elevations. Indeed, it often provides zero-centered median values of the differences as well as a smaller standard deviation of the differences than the other experiments. Analyses by elevation band show that differences are often larger in terms of median and standard deviation at intermediate altitudes (i.e., 1500 and 2100 m), which may be linked to partial or intermittent snow conditions, more sensitive to atmospheric and ground physical state and therefore more difficult to model adequately. A closer look at each type of surface also shows that the main differences between our experiments lie in the presence of vegetation and are higher for high-vegetation than for low-vegetation surface types.

4 Discussion

In this study, we analyzed the results of various surface configurations in the CNRM-AROME high-resolution regional climate model on snowpack simulations in the European Alps, tested through coupled model simulations at 2.5 km horizontal resolution, driven by the ERA5 large-scale reanalysis.

Various reference datasets and indicators were used to evaluate multiple aspects of the snowpack simulations, including the snow cover duration using remote sensing data from MODIS, a multivariate analysis at four well-instrumented sites (including air temperature and radiation balance terms), and a comparison of the snow depth on a large set of in situ stations covering the Alpine ridge.

These comparisons allowed us to gain insight into the challenges of snowpack simulation within the different AROME surface configurations and ultimately document an optimized land surface configuration. In the subsequent sections, we examine the various causes for the successful and unsuccessful modifications that we tested and remaining problems and limitations and finally propose perspectives to further improve snow simulation in the European Alps and beyond using CNRM-AROME or other regional climate models.

4.1 D95-3L: an overly simplistic configuration that fails at reproducing the seasonal snowpack evolution

The default surface configuration – namely, D95-3L in this study – is based on a force–restore approach for the soil exchanges and a single-layer parameterization for snow. This surface configuration, used in recent climate study frameworks such as the CORDEX FPS on Convection (Coppola et al., 2020; Pichelli et al., 2021; Ban et al., 2021; Caillaud et al., 2021) and Météo-France's numerical weather prediction system using AROME (Seity et al., 2011; Brousseau et al., 2016) exhibits a series of issues over mountainous regions, such as a cold bias at high elevations (Vionnet et al., 2016; Monteiro et al., 2022; Arnould et al., 2021; Gouttevin et al., 2023) and a generalized overestimation of the amount of snow (Monteiro et al., 2022; Monteiro and Morin, 2023).

These problems have been reproduced in our experiments using 2 years of regional simulation driven by ERA5, and our multiple comparisons allow us to characterize them further.

In Sect. 3.1 and in particular in Fig. 6, we demonstrated on a large set of in situ snow depth observations that this configuration is unable to provide a satisfactory seasonal evolution of the snowpack at all elevations, from the lowest studied at 900 m to the highest at 2700 m. While the first accumulations are generally well correlated with observations, the start of the melt period, from late winter at low elevations to late spring at high elevations, marks the beginning of strong divergences from observations. From this point onwards, the magnitude of melting events is often severely underestimated or even completely missed, leading to a delay and overestimation of the snow depth annual maximum and then of the end of the snow season, which can last for up to a month.

Section 3.2 confirms the underestimated magnitude of snowmelt at a larger spatial scale, displaying a generalized overestimation of the duration of snow cover at all elevations on the map and box plots Fig. 7. Nevertheless, the categorization of differences by surface type (i.e., no vegetation, low vegetation, and high vegetation) adds nuance to this analysis. The overestimation appears to be particularly linked to the presence of vegetation, as the no-vegetation category shows close to zero-centered differences compared with the reference, while the strongest overestimation cases are found for the high-vegetation category above 1500 m and the low-vegetation category above 2100 m.

Based on these pieces of evidence, multiple statements can be formulated to explain the widespread overestimation of snow in the simulation using the D95-3L configuration. The overestimation of winter snowfall at high elevation already reported in past studies (Monteiro et al., 2022; Lucas-Picher et al., 2023; Monteiro and Morin, 2023) may contribute to providing overestimated snow accumulation at the highest elevation bands. Nonetheless, a significant part of this overestimation is likely to be attributed to the design of the configuration itself and, more specifically, how melt is computed in the presence of vegetation, explaining its propensity to make larger errors over these surface types. As stated in Sect. 2.3.1, in the presence of vegetation, the calculation of the melting temperature becomes composite between the surface temperature and the deep soil layer temperature (see Eq. 10 in Sect. 2.3.1). The latter leads to a decoupling between melting intensity and the sub-daily oscillations of the energy balance, which unfortunately results in an underestimation of the snowmelt in many cases.

4.2 ES-DIF: an intermediate complexity surface configuration holding conceptual issues in coupled surface–atmosphere model if only one patch is used

In order to solve some of the issues of the original, simplified D95-3L configuration, we experimented a more detailed and physically based land surface configuration using the multi-layer soil scheme ISBA-DIF (Boone et al., 2000) together with the explicit multi-layer snow scheme ES (Boone and Etchevers, 2001) using only one patch for the NATURE tile in SURFEX.

However, the majority of the results from our study indicate no improvement and, in some cases, a degradation of the simulation of snow depth. Against the multi-station mean snow depth time series in Fig. 6 in Sect. 3.1, the ES-DIF configuration displayed degraded scores of ME compared to D95-3L except at the highest elevations studied of 2700 m. As opposed to the D95-3L experiment, the simulation using the ES-DIF configuration underestimated the snow depth value from the first melt event after the date of the snow depth annual maximum to the end of the snow season, often anticipated to be from 15 days to up to a month. Albeit this configuration generally underestimates the amount of snow, we note that it simulates variations in the snow depth (see, e.g., Fig. 6) that are much better correlated with the observations than the simulation using the D95-3L configuration, reflected in a higher R² score at all elevations excepted 2100 m. We attribute it to its explicit treatment of the snowpack and an enriched description of the physical processes within the snowpack (i.e., liquid water retention, phase change, and compaction).

Section 3.2 highlights the underestimation of the snow cover with this configuration by displaying a spatially generalized underestimation of the snow cover duration compared to MODIS except at the highest elevation, as demonstrated by the near-zero-centered differences at 2700 m in the box plots in Fig. 7b. As for the D95-3L configuration, the box plots show that the differences are enhanced in the presence of vegetation and at intermediate elevations, while the lowest differences are found at high elevations for the no-vegetation surface type.

As for the overestimation of the snowpack simulated by the D95-3L configuration discussed in Sect. 4.1, several reasons can be invoked to explain the underestimation of the amount of snow and the exaggerated snowmelt intensity at intermediate elevations and in the presence of vegetation in the ES-DIF configuration. Apart from the biased atmospheric conditions, the effects of which are discussed further in a separate subsection, we suspect conceptual choices made in ISBA for computing the surface energy balance in the presence of snow on the ground brought about numerous drawbacks on snowpack modeling, exacerbated in the case when only one patch is used.

Figure 8a illustrates in a simplified way the heat exchanges between vegetation, soil surface, and atmosphere in the presence of snow on the ground in a forested area and Fig. 8b how it is modeled using the ES-DIF configuration.

Figure 8 Schematic illustration of the heat exchanges between the different components of a forested area in the presence of snow for two cases: a simplified case to represent the reality and the way the ES-DIF configuration represents it.

What we would expect to observe in forested areas as well as in open areas is that even a small amount of snow covers most of the surface even if some of it is intercepted by trees, and tree trunks and branches remain uncovered. Consequently, the soil surface would be isolated from the canopy and the atmosphere and would interact mainly with the overlying snow and underlying soil by thermal diffusion and latent heat of phase change. Just above the soil surface, the snowpack exchanges energy with a part of the atmosphere strongly influenced by the presence of the canopy, reducing momentum and overshadowing part of the incoming shortwave but emitting longwave radiation and latent heat flux through evapotranspiration. In this case, the air temperature above the canopy, T_air, would be influenced by the surface components mainly through turbulent fluxes, strongly modified by the canopy roughness.

In ISBA, however, the increase in snow amount on the ground is accounted for by progressively covering the surface and modifying surface variables such as the roughness length, z₀, and albedo, α_s, modulating the amplitude of the radiations absorbed and the turbulent fluxes (Boone and Etchevers, 2001). However, in forested areas, as the vegetation is included in a composite soil–vegetation surface layer, this means that snow gradually “replaces” vegetation.

To ensure that the effects of the vegetation on surface variables such as lowering albedo and increasing roughness length are adequately represented, notably to preserve turbulent fluxes to the atmosphere, the parameterization of the snow fraction over vegetation (see Eq. 5) has been developed in such a way that part of the high-vegetation surface remains uncovered even when a large amount of snow is present.

Unfortunately, since a single surface temperature is used for both snow-covered and snow-free surfaces, this approach inhibits most of the insulating effect of snow cover on the underlying ground.

Napoly et al. (2020) and Nousu et al. (2024) have documented the effects of this approach on the simulated snowpack in forests and examined its effects on the radiation balance and turbulent fluxes using a standalone simulation of the surface forced by observed atmospheric variables. The authors found that the model significantly underestimated the snow extent and depth at several forested sites. The main reason is an unrealistic coupling between the soil surface and the atmosphere in the case of a partially snow-covered surface, which leads to a strong overestimation of the diurnal amplitude of the soil heat flux. During daytime, a considerable amount of energy can be absorbed over the snow-free area (mainly through incoming shortwave radiation), overestimating the warming of the surface layer below the snow-covered areas, heating the snowpack from below and likely causing excessive melting at its base. They also note that this approach tends to overestimate the latent heat fluxes, mainly due to the overestimation of soil evaporation and that the strong coupling also leads to unrealistic surface cooling during nighttime, highlighted by an average cold bias of soil temperatures on the order of −5 °C.

Initially intended as a compromise between the simulation of turbulent fluxes in the presence of vegetation and the insulating effect of the snowpack, this way of treating snow cover in vegetated areas turns out to be largely unbalanced to the detriment of soil temperature and snowpack simulation.

In our study, we suspect these feedbacks to constitute the leading mechanisms causing a large part of the underestimation of snow cover in the ES-DIF configuration, particularly in the presence of vegetation and at intermediate elevations, where we often find partially snow-covered surfaces. It is important to note that, in addition to the heating feedback from snow-free surfaces to snow-covered surfaces, other factors not directly tested in this study are likely to exacerbate this effect. Indeed, due to their effects on reducing thermal conductivity, some physical processes, such as air trapping at the soil–snow interface due to the presence of a litter layer and/or low vegetation as well as a poor representation of organic matter content in the upper soil layers may contribute to and exacerbate basal snowpack melting.

4.3 ES-DIF-OPT: an optimized configuration towards addressing conceptual issues in snow representation in the ISBA LSM

The ES-DIF-OPT simulation uses identical soil and snow schemes as the ES-DIF setup, to which a number of modifications have been added. Their primary aim is the reduction in excessive snowmelt in the ES-DIF simulations, discussed in Sect. 4.2.

Section 3.1 demonstrates that the ES-DIF-OPT configuration provides the best seasonal evolution of the snow depth at all elevations (see Fig. 6), systematically increasing R² and decreasing ME values compared to the other two configurations. Its snow depth variations during accumulation and melt periods are similar to the ES-DIF simulation, but the frequency and magnitude of melting events are in much better agreement with observations.

Comparisons with MODIS snow cover duration (Sect. 3.2) show a clear improvement in the seasonality of the snow cover compared to the other configurations. The differences are smaller for all elevation bands and surface types, with a greater reduction for vegetated surface types and at intermediate elevations, as illustrated in Fig. 7, by medians of differences close to zero and distributions that show a reduced variance compared to the other configurations.

Figures C1 and C2 in Appendix C provide a more thorough visualization of the distinct impacts of each modification and clarify the source of the biases in the snowpack simulation identified in the ES-DIF configuration.

The 3-PATCHS modification, as described in Sect. 2.2.3, splits the calculation of energy and mass balances at the subgrid scale, performing an independent calculation for each patch and summing the fluxes obtained rather than averaging the surface variables of each surface type and performing a single calculation. Its effects are significant on the simulation of the snowpack, as shown in Fig. C1, by reducing the ME and increasing the R², especially at low and intermediate elevations. We note that its impacts are limited above 1500 m and negligible at 2700 m (i.e., where most of the grid cells are devoid of vegetation; see Fig. 2.6.2). The box plots of the snow cover duration differences in Fig. C2 show that the 3-PATCHS modification has limited impact in the case of a prevailing (i.e., cover more than 75 % of the grid cell) no-vegetation surface type, while the reductions in the differences are high for the low-vegetation and high-vegetation surface types and, obviously, for the mixed surface type at low and intermediate elevations, where the proportion of vegetation is high. Thus, the significant improvements brought about by this modification show that, in many cases, the aggregated characteristics of the surfaces, when using a single patch, cause exaggerated melt. It is likely that in these cases, the resulting aggregated surface characteristics present sufficiently low surface albedo and high roughness length to trigger the undesirable mechanism described in detail in Sect. 4.2. This leads to undue melting at either the base of the snowpack through an overestimation of the heating of the soil surface under snow-covered surface or its surface through the overestimation of turbulent fluxes. This hypothesis is in line with the greater impact of this modification seen in the presence of high vegetation, where the aggregation of surface characteristics produces a higher roughness length and a lower albedo, and at intermediate altitudes, where the surface temperature is often near the freezing point.

The GFLUX modification consists of lowering the intensity of heat exchange between the surface and the overlying snow layer in the case of partially snow-covered grid cells (see Sect. 2.3.3). This is a pragmatic approach, albeit not grounded on physical principles, to decrease exaggerated melt due to excessive melting at the base of the snowpack in the case of a partially snow-covered surface as described in Sect. 4.2. Figures C1 and C2 demonstrate significant improvements at the same elevation bands and surface types as the 3-PATCHS modification (i.e., intermediate elevation bands and in the presence of vegetation). The effectiveness of the modification in these areas supports the hypothesis of an overestimation of the ground heat flux, inducing basal melt in partially snow-covered surfaces, when this modification is not implemented.

Compared to the 3-PATCHS and GFLUX, the WSN-1 modification shows only slight improvements, probably limited by the chosen value, which may still be too high to substantially increase the sensitivity of the snow fraction parameterization. However, the chosen value is a first attempt as a compromise to avoid the over-reduction in the turbulent fluxes in the near-surface atmosphere (see Sect. 4.2 for further details).

Overall, the ES-DIF-OPT configuration outperformed the other two in every aspect of the snow simulations investigated in this study. We suggest considering this configuration as a basis for future simulations using the CNRM-AROME model. Although the study concentrates on evaluating seasonal snow cover in mountainous areas, a clear improvement in the representation of snow events in lowlands is also expected. Snow events in lowland areas, which are typically less intense than those in mountainous locations, are unlikely to cover the entire surface of the grid cells within the ISBA model. Consequently, these events would most likely be underestimated (snow depth and snow cover duration) using the ES-DIF configurations with only one patch, leading to unrealistically early melting, which is greatly reduced using our optimized configuration. Specific investigations are required to assess these expected improvements explicitly.

4.4 Perspectives regarding error compensations and the effects of land surface–atmosphere coupling

The present work has only addressed a fraction of the sources of model errors, focusing on those related to the surface scheme. The Results section and, more specifically, Fig. 9 provide evidence of numerous other potential sources of errors that can have major implications for the simulation of the snowpack. This justifies that we did not attempt to achieve a “perfect” match with observations through modifications of the land surface scheme, both because observations are also affected by uncertainties and because other sources of errors, in particular in the atmospheric part of the CNRM-AROME model, certainly play a role in the overall performance of the model used.

Figure 9 (a) Panel of time series the snow depth at four well-instrumented stations – i.e., Davos (DAV), Torgnon (TOR), Weissfluhjoch (WJF), and Col du Lac Blanc (CLB) – for the period from 1 October 2019 to 30 September 2020. The solid colored lines correspond to the time series simulated for each land surface configuration at the grid cell, including the site location, while the black dots indicate the observed time series. For each graph and each experiment, the simulated snow fraction values are indicated in the color-shaded areas (the y axis on the right refers to the simulated snow fraction). (b) Mean error (ME) values for each of the configurations for the four well-instrumented stations. The ME values are calculated over daily values for the 2019–2020 season for days during which snow is present (i.e., snow depth > 1 cm) in both observed and simulated time series. “Count” refers to the number of days needed to compute the score.

Figure 9a shows the time series for the period from 1 October 2019 to 30 September 2020 of the simulated and observed snow depth at four well-instrumented stations for which a detailed description can be found in Fig. D1 and Table D1 in Appendix D. Figure 9b displays the ME of multiple variables, calculated as the differences between each configuration and the observed time series over the common snow-covered days (snow depth of > 1 cm), with the aim of highlighting the interactions between atmospheric forcings and other near-surface variables that could alter snowpack simulations.

The site of Torgnon and Col du Lac Blanc in Fig. 9a, both located above 2000 m above sea level (see Table D1 in Appendix D), show an overestimated amount of snow all along the season, with discrepancies progressively widening during the accumulation period. During this period (early in the season), only few melt events occur (in both observations and simulations) and appear to be well captured by the different model configurations, as shown by simulated snow depth variations highly correlated to observations. These findings would therefore point towards an overestimation of the snowfall amount by the CNRM-AROME model at these high-elevation sites. This overestimation of snowfall is consistent with previous studies that used the CNRM-AROME model, such as Monteiro et al. (2022), Monteiro and Morin (2023), and Lucas-Picher et al. (2023), which reported a widespread overestimation of winter precipitation at high elevations over the Alpine ridge. This overestimation is also confirmed by Haddjeri et al. (2023), who used the NWP AROME precipitation fields to force standalone SURFEX simulations, resulting in an overestimated amount of snow over multiple areas the French Alps. Nonetheless, this overestimation may not be systematic, as shown at the Weissfluhjoch site located at 2500 m, presenting snow depth values during the accumulation period that are close to the observations, and at 2100 m and 2700 m compared to the multi-station mean in Fig. 6 in Sect. 3.1.

The simulated irradiance values (incoming longwave and shortwave) at each of the well-instrumented sites display systematic biases (see Fig. 9b) – namely, an overestimation of the shortwave radiation and an underestimation of the longwave radiation. These can be substantial during the snow season, reaching −35 W m⁻² for the incoming longwave radiation and +50 W m⁻² for the incoming shortwave radiation. They corroborate previous studies that have documented these biases in the NWP version of AROME (Vionnet et al., 2016; Quéno et al., 2020; Gouttevin et al., 2023) and the CNRM-AROME climate simulation (Lucas-Picher et al., 2023), all of them attributing it to the underestimation of the cloud cover over mountainous regions. In our study, the positive incoming shortwave biases may play a key role in the exaggerated melt frequency and magnitude in the ES-DIF configuration at low and intermediate elevations. Indeed, in cases of partially snow-covered surfaces, incoming shortwave biases may strongly favor the warming of the soil surface and enhance the feedback leading to undue melt, as described in detail in Sect. 4.2.

Conversely, the negative biases of the incoming longwave radiation in Fig. 9b are likely to contribute to the underestimated melting rate at some high-elevation sites such as the Col du Lac Blanc site as Lapo et al. (2015) and Quéno et al. (2020) demonstrated that a deficit of incoming longwave radiation often leads to a large underestimation of melt intensity. These biases are also prone to lead to the over-cooling of the surface during nighttime, increasing the near-surface stability of the atmosphere and triggering the self-sustaining stability feedback energy loss described by Lapo et al. (2015) in which the increasing stability inhibits turbulent exchanges, thus accelerating the surface cooling. As reported by Gouttevin et al. (2023), the biased incoming longwave explains a large fraction of the surface temperature biases in AROME NWP simulations at Col du Lac Blanc and contributes to the cold bias in near-surface air temperature as it is determined from surface temperature. These biases and associated feedback are likely to contribute to the cold biases observed at other high-elevation sites, such as Torgnon and Weissfluhjoch, and we can even expect it to be relatively widespread across the Alpine ridge, where a generalized winter and spring cold bias at high elevations has been documented in several studies (Monteiro et al., 2022; Lucas-Picher et al., 2023; Monteiro and Morin, 2023).

In the end, in addition to the surface modeling errors, our experiments also corroborate two substantial error sources in terms of irradiance values. The literature suggests that the biases have conflicting impacts on the snowpack and are location-specific (Lapo et al., 2015; Quéno et al., 2020). The location factors that significantly influence the effects of the biases include exposure, elevation, and climate type. Therefore, disentangling the impacts of individual biases in coupled surface–atmosphere simulations is almost impossible.

5 Conclusions

In this study, we investigated 3-year-long simulation results using three main surface configurations of the coupled surface–atmosphere convection-permitting regional model CNRM-AROME over the European Alps at 2.5 km horizontal resolution. It is the first case where detailed investigations using CNRM-AROME as a regional climate model are performed using a land surface configuration strongly deviating from the land surface configuration used by AROME for NWP applications.

By leveraging different datasets used as a reference, we explore multiple aspects of the simulation of the seasonal snow cover, such as an extensive analyses of the snow depth time series over 2018–2019 and a spatially exhaustive comparison of the snow cover duration using MODIS remote sensing data over 2018–2020.

Based on this analysis, we further documented the issues of the current land surface configuration used in CNRM-AROME for climate studies (Coppola et al., 2020; Caillaud et al., 2021) and numerical weather prediction (Seity et al., 2011; Brousseau et al., 2016) (i.e., D95-3L). We shed light on the potentials and limitations of an enriched surface configuration using intermediate complexity and multi-layer soil (Boone et al., 2000; Decharme et al., 2011) and snow (Boone and Etchevers, 2001) schemes (i.e., ES-DIF). Ultimately, we introduced an optimized land surface configurations based on the ES-DIF configuration (i.e., ES-DIF-OPT).

We confirmed the documented issues of the D95-3L default configuration (Monteiro et al., 2022; Lucas-Picher et al., 2023; Monteiro and Morin, 2023) – namely, a spatially widespread overestimation of the amount of snow and delay of the end of the snow season of up to a month and a half. Using a categorical analysis of the snow cover duration by surface type and a comprehensive comparison of the energy balance at some punctual sites reveal wider discrepancies in vegetated areas and a clear underestimation of melt during most of the snow season. These issues were mainly attributed to the over-simplicity of the snow scheme, including the snowpack in a soil–vegetation–snow composite layer.

We demonstrated that the multi-layer soil- and snow-scheme configuration ES-DIF failed to reproduce the seasonality of the snow cover in the European Alps if only one patch is used. Although the many additional physical processes (compared to the D95-3L configuration) make it possible for this configuration to capture well most of the variations in the snow depth during the snow season, the simulation presents a widespread underestimation of the duration of the snow cover below 2500 m, particularly at intermediate elevations and in the presence of vegetation, resulting from an exaggerated melt. We discussed the origin of this issue, already reported in standalone surface simulations by Napoly et al. (2020) and Nousu et al. (2024) and attributed to conceptual choices in the ISBA LSM with respect to the snow cover fraction. This behavior indeed results from an underestimated snow cover fraction in vegetated areas, leading to an over-warming of the soil surface below the snowpack provoking undue melt at its base.

This issue is a major topic as it appears in many modeling systems using similar configurations. Indeed, it is identified as the cause of a significant underestimation of the snow cover in boreal regions in the CNRM-CM6 GCM (Decharme et al., 2019) results and likely explains the early melt in the latest simulations of the CNRM-ALADIN RCM, also using a similar land surface configuration, as shown by Monteiro and Morin (2023).

Finally, we introduced the ES-DIF-OPT configuration, mostly based on existing options in SURFEX but not activated hitherto, which provides the best estimation of the seasonality of the snow cover and daily evolution of the height of snow over a large sample of observations in the European Alps. We find that the reduction in the ground heat flux and the splitting of the energy balance for three surface type categories constitute the major contributions to lowering the errors. Their effectiveness confirms the hypotheses put forward to explain the exaggerated melting in the ES-DIF simulation and underlines the importance, even at kilometer resolution, of taking into account the main subgrid heterogeneities concerning the surface type in mountainous terrain. In the analysis of the preferred configuration, ES-DIF-OPT, this allows for the simulation of snow cover to be satisfactory compared with the references used in many cases.

At the end, the ES-DIF-OPT configuration consists of an adjusted configuration, which limits the shortcomings of the ES-DIF configuration, offering an interesting alternative that we recommend for future simulations using the CNRM-AROME model and other modeling systems using similar configurations, pending that a more physical solution can be implemented in coupled model configurations. Indeed, the ES-DIF-OPT configuration still fails to adequately capture the evolution of snow beneath the forest and around the mean snow line elevation, where partially snow-covered surfaces are often observed. These results advocate for the use of explicit vegetation modules instead of composite soil-vegetation approach, enabling the current snow-covered fraction parameterization to be redefined as the below-canopy snow coverage, reducing the excessive sensitivity of simulation results to this very uncertain parameterization (Nousu et al., 2024). Such an approach is already implemented in most LSMs (e.g., ECLand Boussetta et al., 2021, NOAH-MP Niu et al., 2011, CLM5 Lawrence et al., 2019, JULES Best et al., 2011), and recently, developments have been made to implemented it in ISBA through the MEB (multi-energy balance) module (Boone et al., 2017; Napoly et al., 2017). Unfortunately, this advanced version of the land surface scheme SURFEX cannot yet be activated in CNRM-AROME. The developments performed and evaluated in this work demonstrate the benefit in bridging the gap between currently used land surface configurations in AROME (for climate and NWP applications) and existing state-of-the-art configurations in SURFEX that require further work in order to be applicable in coupled model experiments and applications.

Appendix A: Sensitivity of the snow cover fraction parameterization over vegetation

Figure A1 (a) Map of the mean surface roughness length values over the two winter seasons (November to April for the 2018–2019 and 2019–2020 periods) using the D95-3L configuration. (b) Box plot representing the spatial distribution of the mean roughness length values over the two winter seasons (November to April for the 2018–2019 and 2019–2020 periods) using the D95-3L configuration, classified by the prevailing type of surface (see Sect. 2.7.3 for details). It is noteworthy that although the roughness length values are displayed for the D95-3L configuration and vary with the amount of snow on the grid cell, they are very similar for all the configurations tested.

Figure A2 Snow cover fraction over vegetation as a function of snow depth for multiple combinations of values for roughness length z0 and scaling factor W_sn. The black curves represent the sensitivity of the snow cover fraction parameterization using a value of W_sn=5 as set values for the D95-3L and ES-DIF configurations, while the orange curves use a value of W_sn=1 for the ES-DIF-OPT configuration.

Appendix B: Impact of the initialization (spin-up) approach on snowpack simulations

Figure B1 (a) The maps of the snow depth, the total water content of soil, and the soil temperature at 2 m fields of our experiments at the initialization date (1 January 2018) using the default initialization procedures. (b) The maps of the differences between the default initialization procedures and the initialization resulting from 13 years of offline spin-up for the snow depth, the total water content of soil, and the soil temperature at 2 m fields at the initialization date (1 January 2018).

Figure B2 Time series of the snow depth at four sites over the period from 1 January 2018 to 31 December 2020 using the ES-DIF-OPT configuration with either the default initialization field (continuous dark-green line) or the initialization field resulting from 13 years of offline spin-up (continuous light-red line). Shaded areas correspond to periods that were not evaluated in the study.

Figure B3 Time series of the total soil water content at four sites over the period from 1 January 2018 to 31 December 2020 using the ES-DIF-OPT configuration with either the default initialization field (continuous dark-green line) or the initialization field resulting from 13 years of offline spin-up (continuous light-red line). Shaded areas correspond to periods that were not evaluated in the study.

Figure B4 Time series of the soil temperature at 2 m below surface at four sites over the period from 1 January 2018 to 31 December 2020 using the ES-DIF-OPT configuration with either the default initialization field (continuous dark-green line) or the initialization field resulting from 13 years of offline spin-up (continuous light-red line). Shaded areas correspond to periods that were not evaluated in the study.

Appendix C: All tested experiments

C1 Large-scale evaluation of snow depth over the 2018–2019 winter

Figure C1 Multi-station mean time series of the height of snow for the 2018–2019 winter for four elevation bands of 300 m (900 m ± 150 m, 1500 m ± 150 m, 2100 m ± 150 m, and 2700 m ± 150 m). Continuous colored lines correspond to the simulated multi-station mean time series for each of the configurations. Black lines with circle markers correspond to the multi-station mean time series of the in situ measurements, with the inter-station standard deviation represented in gray-shaded areas. For each elevation band, the number of stations used to construct the mean and standard deviation is displayed in blue font. At each elevation band and for all configurations, the correlation (R²) and the mean error (ME) computed using the multi-mean time series between the simulated and the in situ measurements are displayed.

C2 Snow cover duration evaluation using MODIS remote sensing data

Figure C2 Snow cover duration differences between the different configurations and MODIS observations in the European Alps. Note that MODIS products that initially had 500 m horizontal resolution have been regridded over the CNRM-AROME horizontal resolution (2.5 km) grid using a first-order conservative method. (a) Map of the average differences (mean error) of the snow cover duration (SCD) over two seasons (2018–2019 and 2019–2020) for each configuration compared to MODIS SCD. (b) Box plot representing the spatial distribution of the average differences (mean error) of the SCD over two seasons (2018–2019 and 2019–2020) compared to MODIS SCD for each dataset for six elevation bands of 300 m (900 m ± 150 m, 1500 m ± 150 m, 1800 m ± 150 m, 2100 m ± 150 m, 2400 m ± 150 m, and 2700 m ± 150 m). Each column corresponds to the values classified by prevailing type of vegetation (see Sect. 3.2 for details).

Appendix D: Multivariate comparison at four well-instrumented sites

Location and characteristics of the station and the corresponding CNRM-AROME grid point

Figure D1 Location of the well-instrumented in situ stations together with the digital elevation model at 1 km horizontal resolution.

Table D1 Main characteristics of the well-instrumented in situ stations. Besides their longitude and latitude, the elevation and surface type are given for the stations itself and for the corresponding CNRM-AROME grid cell including it and used for the comparison.