Land-atmosphere interactions in sub-polar and alpine climates in the CORDEX FPS LUCAS models: I. Evaluation of the snow-albedo effect

In the Northern Hemisphere, the Sseasonal snow cover plays a major role in the climate system of the Northern Hemisphere via its effect on land surface albedo and fluxes. In climate models Tthe parameterization of these interactions between snow and -atmosphere interactions in climate models remains a source of uncertainty and biases in the representation of the local and global climate. Here, we first evaluate the ability of an ensemble of regional climate models (RCMs) coupled to different land surface models to simulate the snow-atmosphere interactions cover over Europe, in winter and spring. We use a previously defined index, the Snow Albedo Sensitivity Index (SASI), to quantify the radiative forcing due associated with snow cover anomalies the snow albedo effect. By comparing RCM-derived SASI values with SASI calculated from reanalyses and satellite retrievals, we show that an accurate simulation of snow cover is essential for correctly reproducing the observed forcing over mid-and high-latitudes in Europe. The choice of parameterizations, and with first and foremost primarily the choice of the land surface model, but also the convection scheme and the planetary boundary layer strongly influences the representation of SASI as it affects the ability of climate models to simulate snow cover accurately. The degree of agreement between the datasets differs between the accumulation and ablation periods, with the latter one presenting the greatest challenge for the RCMs. Given the dominant role of land surface processes in the simulation of snow cover during the ablation period, the results suggest that, during this time period, the choice of the land surface model is more critical for the representation of SASI than the atmospheric model during this time period.


Introduction
Snow is an important part of the climate system, as it regulatinges the temperature of the Earth's surface via its effect on surface albedo and surface fluxes.In mid-and high-latitude regions, snow is the main interface through which land interacts with the atmosphere during the cold season and the importance of snow-atmosphere interactions in modulating the energy budget at high latitudes during winter has been demonstrated (Diro and Sushama, 2018;Henderson et al., 2018;Xu and Dirmeyer, 2013a).Snow cover extent and depth can modify both surface energy and moisture budgets, triggering complex feedback mechanisms that impact both local and remote climates (Diro and Sushama, 2018).
In particular, snow can have a strong impact on climate due to its high albedo, primarily because of the contrast in the surface energy balance between snow-covered and snow-free land surfaces (Qu and Hall, 2014).Reciprocally, with climate change, rising temperatures are already altering the Earth's snow amount and occurrences, for example shortening, for example, the snow season in Eurasia (Ye and Cohen, 2013;Gobiet et al., 2014;Mioduszewski et al., 2015;Beniston et al., 2018;Matiu et al., 2020).
In this context, it is crucial to better understand snow-atmosphere processes and to the evaluate the ability of climate models to represent them.
The direct impact of snow on the atmosphere is known as the snow albedo effect (SAE; Xu and Dirmeyer, 2013b2011, 2013), where the presence of snow affects the land surface energy budget and influences the local climate, modifying near-surface air temperature.To quantify the contribution from the SAE to the snow-atmosphere coupling, Xu and Dirmeyer (2013b1) developed the Snow Albedo Sensitivity Index (SASI).This index combines incoming shortwave radiation with snow cover variability to quantify the snow-albedo coupling strength, i.e.SASI estimates the degree to which the atmosphere responds to anomalies in snow cover.Applying SASI to satellite observations, Xu and Dirmeyer (2013b1) found that the coupling between snow and albedo is particularly strong during the snowmelt period in the Northern Hemisphere.At high-latitudes, for example, the effects of snow cover on the climate is strongly related to the way vegetation cover is prescribed.Removal of boreal forests locally reduces surface air temperature and precipitation by increasing surface albedo and decreasing plant evapotranspiration (Snyder et al. 2004).The strength of the coupling between snow and the atmosphere is determined by processes involving radiative fluxes but also hydrology.Therefore, Xu and Dirmeyer (2013b) also defined the snow hydrological effect (SHE), which includes the effects is a result of soil moisture anomalies from snowmelt.Through land-atmosphere interactions, they soil moisture anomalies have a delayed impact on the atmosphere.Besides these direct and indirect effects, positive and negative snow-atmosphere feedbacks, such as the snow-albedo feedback (SAF; Qu and Hall, 2007;Fletcher et al., 2015;Thackeray et al., 2018) can amplify or damp anomalies.The SAF represents changes in surface albedo from cooling (warming) that can cause decreases (increases) in absorbed solar radiation, amplifying the initial cooling (warming).Hence, SAFIt is an important driver for regional climate change in Northern Hemisphere land areas.Here, we focus on the one-way direct impact of snow on the atmosphere through, SAE.To quantify the contribution from the SAE to the snowatmosphere coupling, Xu and Dirmeyer (2013b) developed the Snow Albedo Sensitivity Index (SASI).This index combines incoming shortwave radiation with snow cover variability to quantify the snowalbedo coupling strength, i.e.SASI estimates the degree to which the radiative forcing responds to anomalies in snow cover.Applying SASI to satellite observations, Xu and Dirmeyer (2013b) found that the coupling between snow and albedo is particularly strong during the snowmelt period in the Northern Hemisphere.At high-latitudes, for example, the effects of snow cover on the climate is strongly related to the way vegetation cover is prescribed.Removal of boreal forests locally reduces surface air temperature and precipitation by increasing surface albedo and decreasing plant evapotranspiration (Snyder et al. 2004).
While some previous studies have investigated snow-atmosphere processes in climate models for specific regions (e.g.European Alps; Magnusson et al., 2010;Diro et al. 2018;Matiu et al., 2019;Lüthi et al., 2019), the literature remains limited.Here, we build on earlier work from Xu andDirmeyer (2011, 2013a,b), investigatinge the ability of an ensemble of regional climate models (RCMs) to represent snow cover and the radiative forcing associated with from snow cover anomalies the snow albedo effect by evaluating (SASI) over Europe, including a comparison between mid-and high-latitude regions.We derive SASI using radiative fluxes and snow cover from satellites, reanalyseis and climate model outputs.WBuilding on findings by Xu andDirmeyer (2011, 2013a,b), we focus on winter and Formatted: Not Highlight spring seasons, i.e. transitioning duringfrom the accumulation andto the ablation period, when SASI is reachesing itsa maximum.While some previous studies have investigated snow-atmosphere processes in climate models for specific regions (e.g.European Alps; Magnusson et al., 2010;Diro et al. 2018;Matiu et al., 2019;Lüthi et al., 2019), the literature remains limited.Here, wWe use the RCMs outputs from the flagship pilot study Land Use and Climate Across Scale (LUCAS;Rechid et al., 2017;Breil et al., 2020;Davin et al., 2020;Reinhart et al., 2020;Sofiadis et al., 2021).LUCASIt is endorsed by the Coordinated Regional Climate Downscaling Experiment (CORDEX) of the World Climate Research Programme (WCRP) over the European domain (EURO-CORDEX, Jacob et al., 2020) and it enables us to perform a broader assessment of several RCMs within a consistent framework.Our assessment is carried out in two parts and published in companion articles.In Part I, we investigate the ability of these RCMs to represent snow cover and the SASI under present-day land cover distribution, while in Part II we explore the effects of large-scale changes in vegetation cover.In LUCAS, each RCM performed three coupled land-atmosphere experiments at the European scale: two idealized and intensive land use change experiments (GRASS and FOREST) and a control experiment (EVAL).The GRASS and FOREST experiments will be examined in the companion paper (Part II), while here, we use simulationsten models from the EVAL experiment only, where RCMs which employ their standard land use and land cover maps.
Section 2 introduces the modeling and observational datasets used in this study as well as the derivation of SASI, while Section 3 examines and discusses the ability of climate models to represent snow cover and SASI compared with satellite observations and reanalyses, focusing on the strength and timing of the signal.Further, the origin of the differences in SASI between the models are explored by evaluating potential common biases in the ensemble of simulations as well as individual model biases.
The analysis also explores the differences in SASI between mid-and high-latitude regions, opening the discussion on the impacts of different land cover for the simulation of SASI, which will be further explored in Part II.Finally, Section 4 the last sections offers some concluding remarks.

LUCAS experiments and models
Commented [1]: The same sentence is written in the first paragraph of the methods.

The LUCAS experiments
The simulations from the flagship pilot study LUCAS simulations cover the standard EURO-CORDEX domain (Jacob et al., 2014) with a horizontal grid resolution of 0.44° (around 50 km).All RCMs in LUCAS, except the RegCM model, use a rotated coordinate system, which is a cartographic projection to transform coordinates from a 3D sphere to a 2D plane (the model domain).The RegCM model except the RegCM model, which applies a Lambert conformal projection (suitable for midlatitudes) on a regular grid.Here we use outputs from the EVAL experiment, which employ standard land use and land cover maps.; the GRASS and FOREST experiments will be examined in the companion paper (part II).The time resolution at which outputs have been stored varies from one variable to another and follows the CORDEX protocol.All simulations span the period 1986-2015 (with a spin-up period ranging from one up to six years depending on the model) and take lateral and boundary conditions from the ERA-Interim reanalysis (Dee et al., 2011).More details can be found in Davin et al. (2020).

Models and configurations
We use the outputs from ten coupled surface-atmosphere RCM simulations that participated in the LUCAS project and were available at the time when we performed the analysis.The main model characteristics that are important for snow albedo coupling are summarized in Table 1, while a detailed description of the RCMs is provided by Davin et al. (2020).The model ensemble presents five different RCMs: COSMO-CLM version 5.0-clm9 (Sørland et al., 2021), WRF version 3.8.1 (Skamarock et al., 2008), RegCM versions 4.6 and 4.7 (Giorgi et al., 2012), RCA4 (Strandberg et al., 2015) and REMO (Jacob et al., 2012).These RCMs contributed with different setups and configurations as described in Table 1.For example, the same RCM is coupled to different land surface models (LSMs): COSMO-CLM is coupled to three distinct LSMs, which are CLM5.0(Lawrence et al., 2020), VEG3D (Breil and Schadler, 2017) and TERRA-ML (Schrodin and Heise, 2002).WRF is coupled with either CLM4.0 (Oleson et al., 2010) or NOAH-MP (Niu et al., 2011).Vice versa, the same LSM is combined with different versions of RCMs.The LSM CLM4.5 (Oleson et al., 2013) LSM is coupled to two distinct versions of RegCM (4.6 and 4.7) which also differ in their choice of convection schemes.There are also two ensemble members institutes withwhere the same RCM and LSM are used (WRF and Noah-MP)  Davin et al. (2020).The time resolution at which model outputs have been stored, varies from one variable to another and follows the CORDEX protocol.For the analyses in the present study, we use daily and monthly model outputs for incoming shortwave radiation and snow cover.For deriving SASI, the native grid of the models was kept, minimising data loss.The other fields were interpolated to a common 0.5°x0.5°grid using Climate Data Operators (CDO) bilinear remapping.The time resolution at which outputs have been stored varies from one variable to another and follows the CORDEX protocol.

Snow schemes across different land surface modelsRepresentation of Ssnow-buried fraction of vegetation in models
At high-latitudes, the effects of snow cover on regional climate is strongly modulated by strongly depend on vegetation cover.Removal of boreal forests locally reduces surface air temperature and precipitation by increasing the duration of the snow cover and, subsequently, In VEG3D, the snow cover fraction is internally calculated as a function of snow depth and vegetation height, and is used to update surface parameters, such as albedo.However, since fsno is not a default model output in VEG3D, the snow cover fraction has been computed for analysis purpose as a snow flag in case of a snow height above a certain threshold, producing a value that is equal to one or zero (i.e., the grid box is covered by snow or not).reproduce the effect on radiation of trees masking the ground snow, TERRA-ML applies a reduction factor for the snow albedo when vegetation (e.g., forest canopies) masks the snow.

In the ensemble, sSome
In terms of snow schemes, some LSMs contain more sophistication than others.Compared to previous CLM versions (i.e., CLM4.0 and CLM4.5),CLM5.0 used by CCLM-CLM5.0counts more snow layers (12 instead of 5), treats separately canopy intercepted snow and more realistically captures temperature and wind effects on the density of fresh snow (Lawrence et al., 2020;van Kampenhout et al., 2017).The RCA4 model system and its internal LSM, used in RCA, include sub-grid orography in the snow cover to capture inhomogeneous snow cover in mountainous areas.Noah-MP allows for 3 snow layers, depending on the total snow depth.To provide a better representation of the ground heat fluxes, the first very layer is only 0.045 m thick.Noah-MP also considers snow interception by the canopy, accounting for wind and temperature effects on snow accumulation and precipitation from the canopy, snow melting and refreezing (Niu and Yang, 2004).The ground snow cover fraction is a function of the snow depth and density and ground roughness (Niu et al., 2007)

Reanalyses and remote sensing data
Reanalysis data from ERA5-Land (Muñoz Sabater, 2019;Muñoz Sabater et al., 2021) and MERRA-2 (Gelaro et al., 2017) as well as satellite data from the Moderate Resolution Imaging Spectroradiometer (MODIS; Hall and Riggs, 2016) are used to evaluate the modelled snow distribution and radiation in the RCMs.Specifically, we use monthly data for snow cover (variable "fractional area of land snow cover" in ERA5-Land MERRA-2) and, incoming shortwave radiation from ERA5-Land and MERRA-2, and daily snow cover data from the MODIS sensors AQUA (MYD10C1) and TERRA (MOD10C1).The reanalysis data are interpolated bilinearly to the common 0.5°x0.5°grid (see Section 2.2).Reanalysis data cover the time period 1986-2015 and MODIS data the period 2003-2015.Only MODIS-AQUA data are will be displayed in the figures of the main part of the article, while data from MODIS-TERRA are will be included in the supplemental material.For MODIS data, the following processing steps are applied: Formatted: Font color: Black 1.Since heavy cloud cover prevents a correct estimation of snow cover, Ddata are masked according to the prevailing cloud cover by applying a threshold of 50% to the percent of clouds in each grid cell.since high cloud cover prevents a correct estimation of snow cover.We apply a threshold of two different thresholds (250% and 50%) to the percent of clouds in each cell.
For comparison, we also show the results when applying a threshold of 20% in Supplementary Figure S1.
2. Only data flagged as "best", "good", and "ok" are used while all other data are masked.
3. Data are conservatively remapped to the common 0.5°x0.5°grid.Conservative remapping is chosen due to the large difference in resolution between the original MODIS data (0.05°) and the target grid (0.5°), as it.It considers all grid points in the interpolation while, e.g., bilinear interpolation would only consider the neighbouring grid cells of the target grid.
4. A land-sea mask is applied to make sure that only land grid points are included in the analysis.
Only grid points with more than 50% land fraction are included.4.5.Data are averaged to monthly resolution.
The masking for MODIS data implies that single grid points can contribute differently to the average over one region.To make the models and reanalyses comparable, each grid point is weighted by the amount of available MODIS data (individually for each month of the whole time period).

Snow Albedo Sensitivity Index (SASI) and geographical scope
SASI is an index that quantifies the climate forcing due to the snow albedo effect (Xu and Dirmeyer, 2013b).It is defined as: where  is the incidentnet shortwave radiation at the surface, (  ) is the standard deviation of snow cover fraction, which represents the interannual variation of monthly-mean snow cover monthlymean values, and ∆ is the average difference between the albedo of a snow-covered surface and the albedo of a snow-free surface.∆ is a constant value of 0.4 as assumed in Xu and Dirmeyer (2013b).
SASI is given in Wm -2 and high values of SASI, such as 10 Wm -2 , indicate a strong climate forcing from the snow albedo effect (Xu and Dirmeyer, 2013b).
To better understand geographical differences in the role of snow for land-atmosphere coupling,

Europe
We start by giving an overview of exploring spatiotemporal differences in snow cover between the different datasets.In Figure 2, we first shows the geographical distribution of snow cover over Europe from January to June based on satellite observations averaged over the 2003-2015 period, and data from the ERA5-Land, MERRA2 and reanalysis, and the LUCAS models from January to June, averaged over the 1986-2015 period.Here we want to start exploring the systematic differences that can be observed between the different datasets in terms of snow cover by giving an overview of where they agree and where they differ.MODIS-AQUAFocusing first on the satellite observations, ERA5-Land, and MERRA2 , it appears that all these datasets capture follow a all show a similar Formatted: Font: 11 pt spatiotemporal cycle, albeit even if with some differences exist in terms of amplitude, e.g. higher snow cover in spring in ERA5-land compared to the MODIS-AQUAsatellite observations and MERRA-2.
Ssnow cover is high during the first months of the year when snow is accumulating (accumulation period), and then decreasing when while snow is melting (ablation period).The satellite and reanalysis datasets capture the At higher latitudes later snow melts at higher latitudes later than at mid-latitudes, showing high snow cover values during spring, while over the rest of Europe theysnow cover values are very low., as shown in Figure 2.Then, looking at Most of the models simulations, most of them also exhibit the same overall spatiotemporal cycle in snow cover, similarly to the satellite observations and reanalyses.H. Hhowever, large differences exist across models can be seen between the simulations regarding the in terms of amplitude and or pattern of snow cover, especially during the ablation period.
This figure aims at illustrating the systematic differences that can be observed between model simulations in terms of snow cover.To facilitate further investigations into inter-model differences, we can consider modelsthem byLooking at the different by atmospheric modeling groups (i.e., WRF, CCLM, RegCM, and others), highlighted shown by which are symbolized by the different colors in the labels of Figure 2., interestingwe can see that Looking at each group individually, we can see that l Large dissimilarities can appear between the members of each group members (e.g.WRFa-NoahMP, WRFb-CLM4.0and WRFc-NoahMP).Therefore, uusing the same atmospheric model does not guarantee producing a similar representation of snow cover.For example, iIn the WRF modeling group, WRFb-CLM4.0has presents much higher values in snow cover than the other members., for example, the three configurations provide very different representationsrepresentation of snow cover.This is also the case for the CCLM group, with CCLM-VEG3D showing higher snow cover values than the rest of the CCLM models.For the RegCM group, the differences between the two group members are less pronouncedobvious.Comparing across Interestingly, the comparison between the atmospherice model groups shows that RegCM models tend to have snow staying longer on the ground, with higher values in snow cover in May and June compared to the other models.Generally, the comparison between the different atmospheric model groups indicates that using the same atmospheric model does not guarantee producing a similar representation of snow cover, emphasizing that the specificities of each configuration (e.g.parameterization, land surface model) can have a large impact on the representation of snow variables.andRegCM.Indeed Unfortunately, tTh These dissimilarities aredissimilarities es in terms of the simulation of snow in climate models are not limited to snow cover, but to other variables such as snow depth, which also exhibit large inter -model differences (see.Supplemental Material S3).
SuchThese variations canvariations differences can have a large effects on the representation of the local climate through, for example, the impact on the surface energy budget.

SASI in satellite observations, reanalyses and RCMs over Europe
In Figure 32, similarly to Figure 2, we first show tFigure 3 shows tThe geographical distribution of SASI over Europe from January to June is shown in Figure 3 for based on satellite observations (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015), the ERA5-Land, MERRA-2 reanalysis, and the LUCAS models  from January to June, averaged over the 1986-2015 period.Focusing first on the satellite observations, MERRA-2 and ERA5-Land, an increase in SASI can be observed during the first months of the year when solar radiation increases and snow accumulates.The maximum is reachedThen it reaches a maximumit is accumulating (accumulation period)reaches, reaching a maximum during the ablation period, in March or April, depending on the region examined, and finallythen SASI decreasesing as when snow completely meltsstarts melting, during the (ablation period).At higher latitudes snow melts later than at mid-latitudes (see, as shown in Figure 2), causing giving rise to high SASI values during spring.The SASI maximum shows a rather sharp peak in East Baltic, while it is more spread out in East Europe, and in Scandinavia., as shown in Figure 2.Then, SASI reaches very low values in May and June when the snow has melted almost entirelyT.This is as expected, and tThise overall seasonal trend is consistent with Xu and Dirmeyer (2013b).Most of tThe models data exhibits a similar the same overall spatiotemporal cycle in SASI as the satellite observations, and ERA5-Land, and MERRA-2.However, large differences can be seen acrossbetween the models in terms of pattern and amplitude.the modeling groups and each of their members the simulations in terms of amplitude or pattern, especially during the ablation period.For the WRF modeling group in April, large differences appear over the domain, with higher SASI values over central Europe for WRFb-CLM4.0and over Scandinavia for WRFa-NoahMP and WRFc-NoahMP.This type of dissimilarities also exist in the CCLM and RegCM groups.In March over the Carpathian Mountains, for example, SASI varies between 1 Wm -2 for WRFa-NoahMP and RCA, and 10 Wm -2 for CCLM-CLM5.0and RegCMa-CLM4.5.It is also noteworthy that for almost all the models, SASI is close to zero everywhere in continental Europe in May and June except , as the snow has almost entirely melted, while in May for RegCMb-CLM4.5 and CCLM-VEG3D, which there are sstill have high values of SASI (~10 Wm -2 ).This point will be further discussed in Section 3.4.
In each atmospheric model group, the simulations show large dissimilarities in terms of SASI amplitude or pattern, especially during the ablation period.For example, WRFac-NoahMP and WRFca-NoahMP show noticeable differences in the amplitude and pattern of SASI (Fig. 32), even though they use the same LSM (Noah-MP) and atmospheric model (WRF).Their differences come from their distincts parameterizations of (planetary boundary layer and convection), affecting the simulatedion of temperature and precipitation, which in return can influence their representation of snow cover and SASI.Thisthus demonstratesdemonstrating the importance of atmospheric processes and their model representation for representing snow processes.FurthermoreThen, when WRF is in a configuration coupled with the LSM CLM4.0 (WRFb-CLM4.0)italso shows different results compared to the simulations from when it is coupled with NoahOAH-MP.IFor example, in WRFa-NoahMP snow melts about one month earlier hshows an earlier poleward migration of high SASI values compared to than in WRFb-CLM4.0,moving north about one month before WRFb-CLM4.0.This type of Such dissimilarities also exist in the CCLM and RegCM and CCLM groups., as well as .Large differences can also be observed between CCLM-CLM5.0,CCLM-TERRA, and CCLM-VEG3D, which ; they all use the same RCM but different LSMs.However, iIn contrast to the two other CCLMosmo configurations, CCLM-VEG3D uses a snow flag for snow cover (i.e., only indicatinges if snow is present or not; Section 2.3), likely explaining its different representation of SASI.This suggests that SASI is very sensitive to the model configurations of and process parameterizations in the climate model.In particular, the choice of the LSM or certain parameterizations (e.g. the convection scheme) can stronglyhighly influence the representation of the climate forcing from the snow albedo effect.The role of the LSM in this context will be investigated further belowin the coming sections of the article.

Investigating the origin of the differences in the representation of SASI
Formatted: Indent: First line: 1,27 cm

Transition between the accumulation and ablation periods
To further investigate the differences in snow albedo coupling strength between the simulations and the observation-based datasets during the accumulation and ablation periods, a time-series of SASI from January to June is presented in Figure 43 for the three sub-regions East Europe, East Baltic and Scandinavia (see Figure 1 for their extents).Before looking at the differences between the different datasets, it is interesting to SASI values are generally higher in East Europe and East Baltic (mid-latitude regions) than in Scandinavia (high-latitude region)in compare the amplitude of SASI between East Baltic and East Europe (mid-latitude regions) with Scandinavia (high-latitude region), which shows slightly higher values of SASI over the mid-latitude regions in satellite observations, MERRA-2, ERA5-Land and most of the RCMs.This confirms previous findings from Xu and Dirmeyer (2013b), which estimated higher values of SASI in mid-versus high-latitude regions in satellite observations.Part of this difference is likely due to might come from the lower values in incoming solar radiation in for Scandinavia compared to the other regions.However, even with lowerhigher SASI values at high-versus mid-latitudes, this result suggests that the radiative forcing due to the snow albedo effect is not negligible over high-latitude regions in winter and spring, highlighting .This result shows again the importance of the snow-atmosphere processes in mid-and high-latitudes in the Northern hemisphere.
ReturningThen, coming back to the comparison of the different datasets, in all three regions, the models and observations indicate a pronounced springtime peak in SASI in all three regions.As already mentioned, Tthe maximum in SASI marks occurs during the ablation period, when snow is melting the transition between the accumulation and ablation periods.The exact timing of this transition depends on the latitude of the respective region examined due to, for example, latitudinal differences in incoming solar radiation.There is also a difference in the timing of the peak between the satellite observations, MERRA-2, and ERA5-Land, in particular over Scandinavia and East Baltic, Aalthough the amplitude of the peak is very similar between the satellite observations and ERA5-Land, it is interesting to see that the timing differs between them, over Scandinavia and East Baltic.Over East Europe the peak occursit happens in March for both the satellite observations and ERA5-Land, for East Baltic in March (satellites) or April (ERA5-Land), and for Scandinavia in April (satellites) or May (ERA5-Land).The origin of these differences remains unclear, but it is not related to the difference in time period between the Commented [8]: I once got a comment that reanalysis are not observation based.If we have enough time (and I still have some time on Monday), we could define the satellites and reanalysis as "reference datasets" and use that name throughout the manuscript.
satellite observations and the reanalysis (see Figure S1) butThe origin of these differences it has not been clarified yet.This might be due to the higher elevations of these two regions compared to East Europe as complex orography is a driving factor for the spatial heterogeneity of precipitation (Grunewald et al., 2014).
The LUCAS simulations also show a pronounced peak in SASI in all regions (Fig. 3), however they do not all agree on the timing and the amplitude of the signal.This is true both within the in the different groups of modelsing groups but also but also for comparing all the simulations.For example, in the East Baltic region, some models (WRFc-NoahMP and WRFa-NoahMP) simulate a peak in March, others in April (WRFb-CLM4.0 and CCLM-CLM5.0)or even in May (RegCMb-CLM4.5 and CCLM-VEG3D).In general, RegCMb-CLM4.5 and CCLM-VEG3D tend to present the latest peak in SASI as well as the highest amplitude in the signal.On the other hand, WRFa-NoahMP tends to produce an earlier peak and lower values of SASI, especially over East Europe.These differences might be related to the way snow remainings longer on the ground (Section 3.1),snow and meltings later in the different models, and will be further explored in the next section.More generally, we see that during the accumulation period, all the datasets are in better agreement compared to the ablation period (Fig. 43).
For East Europe and East Baltic, the spread largely increases in March and for Scandinavia from April until the end of the season, when the snow is melting.
This large model spread during the ablation period is further confirmed by examining Figure 4 showing the pattern correlation between the simulations and ERA5-Land from January to June (not shown).For many models, the correlation is high at the beginning of the season but strongly decreases in March or April, when the snow starts melting.These results are in agreement with previous studies showing the difficulties of climate models to represent snow processes during the ablation period (Essery et al. 2009).Given the dominant role of land surface over atmospheric processes during the ablation period, this suggests that the choice of the LSM is more critical for the representation of the climate forcing from the snow albedo effect than the atmospheric model in spring.For simulatingcalculating snow-covered areas at different stages of ablation, a correct representation of the landscape type is important (Pomeroy et al., 1998).In this context, it is interesting that no systematic differences can be observed between the PFT-dominant versus PFT-tile models representation of the sub-grid scale surface Formatted: Strikethrough heterogeneity (Table 1); as it does not seem to affect the ability of RCMs to represent snow cover or SASI. Figure 4 The derivation of the pattern correlation (not shown) also indicatesindicates shows that the behavior of the RCMs is different between East Europe and East Baltic versus Scandinavia.Over the latter region, most RCMs differ from the reanalysis, as indicated by low correlations.Earlier studies showed that snow accumulates or melts very differently in an open region compared to a forested region (Jonas and Essery, 2014;Moeser et al., 2016).Our results suggest that RCMs represent snow processes better in open spaces like the East Baltic than in forest-covered regions like Scandinavia.The relationship between the representation of SASI and land cover will be further explored in the companion article, Part II.The mountains in Scandinavia could also be a source of biases since the resolution of the RCM simulations (0.44°) can be considered insufficient to represent the more complex topography of Scandinavia.

Inter-model differences in SASI
To better understand the origin of the differences in SASI across RCMs, we explore the relationship between SASI and its components, surface snow cover and shortwave radiation, during the accumulation and ablation periods.Figure 5 presents a comparison of the averaged monthly surface snow cover for the LUCAS simulations, the reanalyses MERRA-2, and ERA5-Land, and MODIS-AQUAas well as the satellite observations from MODIS, averaged over the our three regions of interest, from January to May.As shown in Section 3.1, differences can be observed between the reanalyses and the satellite observations as these different datasets have their own limitations or biases.For example, as each reanalysis dataset is based on a different dynamical core, each model may parameterize or resolve physical processes differently (Daloz et al. 2020).TFor example, the surface snow cover in East Baltic in March is ~0.6 for MODIS, ~0.7 for MERRA-2 and ~0.8 for ERA5-Land.It is therefore important to include several reference observation-based datasets to evaluate the ability of climate models to represent snow cover and estimate the uncertainties associated with this variable.However, it is interesting that there are no systematic differences between the PFT-dominant versus PFT-tile models representation of the sub-grid scale surface heterogeneity (Table 1); as it does not seem to affect the ability of RCMs to represent snow cover or SASI.Based on Figure 43, RegCMb-CLM4.5 and CCLM-VEG3D were identified as models with higher values in SASI during the ablation period and later peaks for all regions.Figure 5 shows that this behavior can be at least partly attributed to their representation of snow cover.During the ablation period, both they all tend to produce higher values of snow cover compared to the other models and also to keep high values later in the season, when they lie completely outside the range of the reference datasets (indicated by the black dots in Figure 5).During the abolation period both This behavior is confirmed by the black dots under these two models during the ablation period as they indicate when the models are outside the range of the reference datasets (MERRA-2, ERA5-Land and MODIS).This is particularly striking for CCLM-VEG3D.Similarly, the low SASI peaks for WRFa-NoahMP, which also occur earlier than the peaks for other models (Figure 43), might be related to the comparatively low lower values in snow cover (WRFa-NoahMP lying outside the range of the reference datasets in all months and all regions) and the small interannual snow cover variability compared to the other RCMs, particularly in East Europe (Figure 5).Again, this is confirmed by the black dots indicated under the model.The differences in snow cover are also reflected by the timing of snowmelt (reduction in snow mass (SWE)) rate of snow melting for the different RCMs (Supplemental Material; Figure S41).The models having high snow cover late in spring (RegCMb-CLM4.5 and CCLM-VEG3D) tend to have later snow melt than the other models while WRFa-NoahMP, showing reduced snow cover earlier than the other models, also tends to have an earlier snowmelt sooner.
Another component of SASI is shortwave radiation at the surface, shown in which is presented in Figure 6 for the LUCAS simulations, the reanalyses and MERRA-2 and ERA5-Land, averaged over our three regions of interest, from January to May.The comparison between the RCMs and the reanalysis shows noticeable differences for some models.Both REMO-iMOVE and WRFa-NoahMPCCLM-VEG3D exhibit very different results in terms of surface shortwave radiation compared to the datasets as shown by the black dots on the figure, showing much lower and higher values than the reference datasets, respectively.However, even with these discrepancies, they both reproduce SASI reasonably well.This seems to indicate that the differences in the representation of the forcing from the snow albedo effect are mostly driven by differences in the representation of snow cover in the models.This is confirmed by Figure 7 showing the average correlation across models between SASI and shortwave radiation (left) as well as SASI and snow cover (right) for the LUCAS models.Scandinavia and East Baltic present similar results with significant, positive correlations between SASI and snow cover for almost all months, associated with positive but not statistically significant correlations between SASI and shortwave radiation.For East Europe, the correlation between SASI and snow cover is lower and not significant in January and February but remains high and significant the rest of the time period.
In parallel, the correlation between SASI and downward shortwave radiation at the surface is negative for almost all months (but and not statistically significant).Overall, high and significant correlations often appear between SASI and snow cover for the three regions from January to June.On the other hand, the correlations between SASI and shortwave radiation are low and usually not significant.This indicates that the differences in the representation of the forcing from the snow albedo effect are mostly driven by differences in the representation of snow cover in the models.

Conclusion
Previous work has demonstrated already shown ed the difficulty for climate models to represent snow variables or processes, such as snow cover and depth (Matiu et al., 2020) or the snow-albedo feedback (SAF;Fletcher et al., 2015), buthowever the origin of the differences between the models is not clear yet.In this work, we focus on the ability of RCMs to first simulate snow cover and then the radiative forcing associated with the interannual variations in snow cover anomalies from the snow albedo effect in winter and spring over Europe and explore the origin of the differences across between the RCMs.Theis radiative forcing associated with snow cover anomalies is represented by the index SASI (Xu and Dirmeyer, 2013b), which quantifies the strength of the coupling between snow and surface net shortwave SW radiation albedo.Ten RCMs from the CORDEX Flagship Pilot Study LUCAS are compared to satellite observations and the reanalysis datasets including ERA5-Land and MERRA-2.These simulations are part of the control experiment of LUCAS.
The results show that climate models are able to reproduce well some of the SASI characteristics (e.g.existence of a peak, amplitude of the peak) compared to reanalyseis and satellite observations (Section 3.21), even if large differences appear between the RCMs, for all groups of modelsing groups.
The climate models' ability to represent SASI is highly related to their representation of snow cover (Section 3.3), which can be difficult to represent for climate models (Section 3.1Matiu et al., 2020).Our results also suggest that the models' capability highly differs between the accumulation and ablation periods.Most models have much lower agreement with reanalyses and satellite observations in the ablation period, with some exceptions (e.g.CCLM-CLM5.0over East Europe), indicating a systematic common bias regarding snow cover in spring, in turn pointing towards a bias from LSMs.This bias seems to be common to most LSMs even if they are based on different assumptions and parameterizations (see Section 2.3).It is also interesting that even though CCLM-TERRA is not as advanced in terms of snow modeling compared to the other models (e.g.Section 2.1.3),it still manages to represent SASI reasonably well over Europe.In addition, there werewas no systematic differencest between the PFT-dominant versus PFT-tile models representation of the sub-grid scale surface heterogeneity (Table 1); PFT-dominant versus PFT-tile) as it does not seem to affect the ability of RCMs to represent snow cover or SASI.
Although it is difficult to identify the origin of the bias in the RCMs, an increase in spatial resolution might improve the simulation of snow cover and therefore the representation of SASI.For example, over Scandinavia, an increase in spatial resolution would provide a better representation of the complex topography of the region as well as its forested areas, which may lead to an improved simulation of the coupling between snow and albedo.The coming phases of LUCAS, phases 2 and 3, could help answer this question as they will produce simulations at a higher spatial resolution, 12 km and convection-permitting (<3km), respectively.Taking advantage of the different configurations of the LUCAS simulations, we have also explored the role of distinct parts of the models in their ability to represent SASI.The first part of this work has already emphasized the role of the LSMs, but other components can also play an important role.WRFc-NoahMP and WRFa-NoahMP, even though using the same RCM and LSM, show noticeable differences in the amplitude and pattern of SASI.Their differences in parameterizations (planetary boundary layer and convection) are certainly affecting the way they represent SASI, highlighting the impact of such choices and the role of atmospheric processes.
Mid-and high-latitude areas are also specifically examined looking at three sub-regions:  (2013b), which estimated higher values of SASI in mid-versus high-latitude regions in satellite observations.Our results also suggest that the climate forcing due to the snow albedo effect is not negligible over high-latitude regions in winter and spring.This is important since often the landatmosphere coupling is considered weaker at higher latitudes (Xu and Dirmeyer, 2011) but it is also possible that this coupling happens through snow and is therefore underestimated.
Although it is difficult to identify the origin of the bias in the RCMs, an increase in spatial resolution might improve the simulation of snow cover and therefore the representation of SASI.For example, over Scandinavia, an increase in spatial resolution would provide a better representation of the complex topography of the region as well as its forested areas, which may lead to an improved simulation of the coupling between snow and albedo.The coming phases of LUCAS, phases 2 and 3, could help answer this question as they will produce simulations at a higher spatial resolutions (, 12 km in phase 2 and convection-permitting (<3km) in phase 3), respectively.Taking advantage of the different configurations of the LUCAS simulations, we have also explored the role of distinct parts of the models in their ability to represent SASI.The first part of this work has already emphasized the role of the LSMs, but other components can also play an important role.For example, WRFc-NoahMP and WRFa-NoahMP, even though using the same RCM and LSM, show noticeable differences in the amplitude and pattern of SASI.
the surface albedo and by decreasing plant evapotranspiration (Snyder et al. 2004).Today, tThe important role of forest albedo on winter and -spring climate in the high-latitudes is highlighted supported by bothwell acknowledged based on field campaigns, such as the Boreal Ecosystem-Atmosphere Study (BOREAS; Betts et al., 2001), and on modeling studies (e.g., Betts and Ball, 1997; Betts et al., 1996; Betts et al., 2001; Bonan, 2008; Davin and Noblet-Ducoudré, 2010; Mooney et al., 2021), which has .These studies led to the implementationing of more sophisticated snow sub-models in LSMs that account for the burial of vegetation by snow.All LSMs in the LUCAS ensemble derive the fraction of vegetation buried by snow, adopting similar approaches that account for snow depth, vegetation height and snow cover fraction.The snow cover fraction fsno depends on the snow cover (which measures the snow amount in water equivalent) accumulated at the surface over bare soil or vegetation and influences the calculation of surface albedo and fluxes.Canopy-intercepted snow does not contribute to the Commented [2]: I propose to remove this sentence since it is already in the Intro the internal LSM of in the RCA4 model (Samuelsson et al., 2015) separately calculate the snow cover fraction during snowfall and snow melting processes, accounting for sub-grid orography when snow melting occurs.In NoahMP, the snow cover fraction depends on snow depth, ground roughness length and snow density (Niu and Yang, 2004).
LSMs arecontain more sophisticatedion than others.CLM5.0 (Lawrence et al., 2020) and NoahMP (Niu et al., 2007) separately treats separately canopy-intercepted snow and more realistically captures temperature and wind effects on snow processes.In addition, LSMs differ in the number of additional layers for snow calculation: CLM5.0 uses 12 snow layers; CLM4.0,CLM4.5 and TERRA-ML (Tolle et al., 2018) use five;, three in NoahMP uses three;, and two in VEG3D and iMOVE use two ;, and one in the RCA4 uses onemodel.The iMOVE model adopts the snow parameterisation from the global climate model ECHAM4 (Roeckner, et al., 1996) and reproduces the snow albedo as a linear function of the snow surface temperature and of the forest fraction in a grid cell, with fixed maximum and minimum snow albedo at temperatures lower than -10°C and at 0°C, respectively (Kotlarskis, 2007).In the VEG3D model, the snow scheme is based on the Canadian Land Surface Scheme (CLASS) (Verseghy, 1991) and ISBA (Douville et al., 1995) and accounts for changes of surface albedo and emissivity as well as processes like compaction, destructive metamorphosis, the melting of snow, and the freezing of liquid water.The TERRA-ML LSM model is a bulk/1D LSM that applies an infinitesimal vegetation layer on top of the soil surface and has no canopy (i.e., vegetation lays flat on the surface).Therefore, the snow always stays on top of the vegetation and there is no snow under the trees.To correctly Commented [5]: I added this paragraph since one reviewer pointed on the fact that snow schemes are not fully described for all RCMs Commented [6]: Thank you Susanna.Formatted: Font: (Default) Arial, Not Bold, Font color: Violet Formatted: Font: Not Bold, Italic, Subscript Formatted: Font: (Default) Arial, Not Bold, Font color: Violet Formatted: Space Before: 12 pt Formatted: Font: Not Bold Scandinavia, East Europe and East Baltic (Section 3.2).The comparison of the three sub-regions shows the difficulties for models to simulate SASI over Scandinavia during the accumulation and ablation periods.The simulation of snow processes in a forested region is more challenging than in an open region(Jonas and Essery, 2014;Moeser et al., 2016).Thus, potentially climate models can potentially have more difficulties representing snow processes in forest-covered regions like Scandinavia compared to open-land regions like East Baltic.The relationship between the representation of SASI and landcover will be further explored in the companion article (Part II), analyzing the other LUCAS experiments (GRASS and FOREST) from LUCAS.Finally, the comparison of mid-versus high-latitude regions shows slightly higher values of SASI over the mid-latitude regions in satellite observations, ERA5-Land, MERRA-2, and most of the RCMs.This confirms previous findings from Xu and Dirmeyer

Figure 1 :
Figures and Tables

Figure 2 :
Figure 2: Spatial maps of snow cover for the satellite observations MODIS-AQUA, the reanalyses ERA5-Land and MERRA2, and the ten RCMs regional climate simulations from the EVAL experiment of LUCAS from January to June.Data show monthly averages from January to June over, averaged over the time period 1986-2015 for models and reanalyses and over 2003-2015 for satellite observations.See

Figure 32 :
Figure 32: Spatial maps of SASI (Wm -2 ) for satellite observations (MODIS-AQUA), the reanalysis ERA5-Land and the ten RCMsregional climate simulations from the EVAL experiment of LUCAS from January to June, averaged over the time period 1986-2015.

Figure 43 :
Figure 43: Time series of the spatial average of SASI for the satellite observations, the reanalysis ERA5-Land and the ten RCMsregional climate simulations from the EVAL experiment of LUCAS in Scandinavia, East Europe and East Baltic (see Figure 1 for their spatial extent).Data are averaged over the time period 1986-2015.

Figure 54 :
Figure 54: As in Figure 3 but for the pattern correlation between SASI and ERA5-Land for the LUCAS simulations.

Figure 5 :
Figure 5: Snow cover for the 10 RCMs, the reanalyses MERRA-2, ERA5-Land and MERRA-2, and the satellite observations, and MODIS-AQUA llite observations (using only data from days and pixels with less than 50% cloud cover) for January to May.The box-and-whisker-plots show the interannual variability of snow cover over 1986-2015 for models and reanalyses and over 2003-2015 for satellite observations.Bars , with the bar representing the median, boxes the interquartile range, and whiskers the minimum/maximum values.Dots indicate models lying outside the range of the reference datasetsMERRA-2, ERA5-Land, and MODIS-AQUA, and MODIS-TERRA, (i.e., the 25th (75th) model percentile is higher (lower) than the highest 75th (lowest 25th) quantile of the reference datasets).For satellite observations we only use data from days and pixels with less than 50% cloud cover.

Figure 7 :
Figure 7: Pearson correlation between SASI and shortwave radiation (left), and SASI and standard deviation of snow cover (right) calculated across RCMs for the three regions Scandinavia, East Baltic, and East Europe for the months January to June during 1986-2015.The values represent the variable (shortwave radiation or variability in snow cover) to which the inter-model variability of SASI is predominantly related to.Bold values indicate statistical significance at the 0.05 level (two-tailed pvalue).
we focus on three sub-regions over Europe, with different climate, vegetation cover, topography andor latitudes: Scandinavia [5 o E-30 o E, 55 o N-70 o N], East Europe [16 o E-30 o E, 44 o N-55 o N] and East Baltic [20 o E-40 o E, 50 o N-62 o N] (see Figure1).The first two regions, Scandinavia and East Europe correspond to regions 8 and 5 of the PRUDENCE project (Prediction of Regional scenarios and Uncertainties for Christensen and Christensen, 2007)isk and Effects;Christensen and Christensen, 2007).The three selected regions differ in terms of climate but also in terms of vegetation: vegetation in Scandinavia is mostly needle-leaved evergreen forests dominate in Scandinaviatrees while the two other regions are covered by cropland and more deciduous trees cover the other two regions.The Scandinavian region also stands out because of its geographical location coverstretching over high latitudes, where the incoming shortwave radiation is very small or zero during winter.In comparison with the plain region of the East Baltic region, which is covered by plains, the East Europe and Scandinavia regions have a more complex topography as they encompass the Carpathian and Scandinavian mountains, respectively.

Table 1 :
Summary of participating RCMs and their LSMs.

Table 1 :
Summary of participating Regional Climate Models and their Land Surface Models.

Table Formatted :
Font: Not Bold