The S2M meteorological and snow cover reanalysis over the French mountainous areas: description and evaluation (1958 - 2021)

. This work introduces the S2M (SAFRAN - SURFEX/ISBA-Crocus - MEPRA) meteorological and snow cover reanalysis in the French Alps, Pyrenees and Corsica, spanning the time period from 1958 to 2021. The simulations are made over elementary areas, referred to as massifs, designed to represent the main drivers of the spatial variability observed in mountain ranges (elevation, slope and aspect). The meteorological reanalysis is performed by the SAFRAN system, which combines information from numerical weather prediction models (ERA-40 reanalysis from 1958 to 2002, ARPEGE from 2002 to 2021) 5 and the best possible set of available in-situ meteorological observations. SAFRAN outputs are used to drive the Crocus detailed snow cover model, which is part of the land surface scheme SURFEX/ISBA. This model chain provides simulations of the evolution of the snow cover, underlying ground, and the associated avalanche hazard using the MEPRA model. This contribution describes and discusses the main climatological characteristics (climatology, variability and trends), and the main limitations of this dataset. We provide a short overview of the scientiﬁc applications using this reanalysis in various scien- 10 tiﬁc ﬁelds related to meteorological conditions and the snow cover in mountain areas. An evaluation of the skill of S2M is also displayed, in particular through comparison to 665 independent

This simulation system was later expanded to covering all of mainland France and Corsica for hydrological monitoring and forecasting purpose (SAFRAN-France -ISBA -MODCOU, SIM) (Vidal et al., 2010;Le Moigne et al., 2020), and provided inspiration for a European-scale analysis system (Soci et al., 2016;Morin et al., 2021). However, in this article we only focus on the "original" model chain addressing mountain regions of France. 60 This paper introduces the latest version of the SCM reanalysis, now referred to as SAFRAN -SURFEX/ISBA-Crocus -MEPRA (S2M). This new version differs from the previous one (Durand et al., 2009a, b) by its temporal extent (15 more years, now spanning 1958-2021), its extension to Corsica and the Pyrenees in addition to the French Alps, and an update of the observations and models involved. For example, Crocus is now fully embedded as a snow cover model of the ISBA land surface model within the SURFEX interface Masson et al., 2013). However, the major innovation is that this new dataset is now freely available (see section 3 for the dataset description and section 7 for access informations) for scientific applications using the AERIS portal (http://dx.doi.org/10.25326/37#v2020.2). The first part of this paper describes the input data and model chain used to produce this reanalysis as well as the simulated variables and details on the practical access to the reanalysis dataset. The last section gives an overview of the possible uses of this dataset, with an emphasis on its three main dimensions (spatial, temporal and altitudinal). It also presents an objective assessment of the limitations of this dataset in 70 terms of climate trends, which have never been published until now in spite of its large use in the French scientific community for climate applications. Last, an evaluation of the reanalysis performance is given in terms of total snow depth by comparison to independent observations. Although the S2M dataset is provided and available from 1958 to 2021, all the results of this study are based on the period 1958-2020 of the dataset because the last year was not available at the time it has been carried out. the SURFEX/ISBA-Crocus snow cover model (including MEPRA), which is driven by atmospheric fields from the 80 SAFRAN reanalysis.

Geometry of the S2M reanalysis
The concept of semi-distributed modelling (i.e. Hydrological Response Units) has been widely used since Beven and Kirkby (1979) to discretize a hydrological catchment following the main drivers of the spatial variability of the key processes with an optimal numerical cost (e.g. MacDonald et al. (2010), Fiddes and Gruber (2012), Ajami et al. (2016), Meng et al. (2018)). In 85 mountainous environments, elevation, aspect and slope are known to be the main drivers of the spatial variability of snow energy balance, due to the strong dependence of temperature and radiation on these topographic features. Therefore, topographic classes are the cheapest solution of numerical discretization to represent this variability and is a common modelling choice in alpine hydrology (Lafaysse et al., 2011;Tarasova et al., 2016;Garavaglia et al., 2017), consistently with the pioneer introduction of Snow Cover Units by Seidel et al. (1983) and Ehrler et al. (1997) for snow remote sensing. This concept was also 90 chosen by Durand et al. (1999) for the operational application of SAFRAN-Crocus modelling in support of avalanche hazard forecasting. Consistently, the S2M reanalysis results from simulations performed over elementary areas specifically designed to represent the main drivers of the spatial variability in mountain ranges called "massifs" (shapefiles of the different massifs are included in the dataset and a glimpse of the geographic division of the three areas can be seen on Figure 2). A massif is a conceptual object corresponding to a mountainous area (of about 1000 km 2 on average) over which the meteorological 95 conditions are considered homogeneous at a given elevation. This hypothesis simplifies the representation of a complex topography by covering the different elevations and aspects of a given massif with a minimum number of representative computation points. An example of this simplification of a real topographic massif is provided in Figure 3. The S2M reanalysis uses a 300 m  Table 8 in Appendix A for more informations. vertical resolution. The upper (resp. lower) elevation for each massif is defined as the 300 m multiple immediately above the highest point (resp. below the lowest point) of a 50 m digital elevation model from the French National Geographic Institute 100 (IGN) in the considered massif. For each elevation, two different simulations are provided: one simulation only contains flat terrain, and for the other one 8 aspects (N, NE, E, SE, S, SW, W and NW) with two different slope angles (20 • and 40 • ) for each aspect. However, users must be aware that a significant small-scale spatial variability of snow cover remain unresolved by this approach (e.g. Weber et al. (2020)).

105
In addition to simulations provided on the massif geometry, simulations are also performed over a set of 665 observation stations. These sites have been selected to be higher than 600 m a.s.l and to provide snow depth measurements. The selected sites cover the three domains of the reanalysis with 435 observations sites in the Alps, 208 in the Pyrenees and 22 in Corsica.
Each site is characterized by its altitude, slope and aspect. A shapefile containing all the stations informations is provided with the dataset. These simulations result from an interpolation between the nearest elevation bands of the corresponding massif 110 and a projection of the direct solar radiation according to the slope and aspect of the station, the time during the day and the information about solar masks from surrounding topography. The main purpose of these local simulations is to compare snow depth simulations to observations in order to assess the performance of the model (see section 4.4) with as few artifacts as pos- sible due to topographical discrepancies between observations and the simulation configuration. However, it must be kept in mind that the spatial scale of SAFRAN analyses (e.g. precipitation) remains at the massif scale even in this local configuration.

Input data
The S2M reanalysis is only fed by input data to the SAFRAN atmospheric analysis system as described in section 2.4.

120
The S2M reanalysis uses two different NWP model output as a guess for the SAFRAN analysis : -The ERA-40 reanalysis (Uppala et al., 2005) between 1958 and 2002, which is based on a uniform data assimilation system (but variable in-situ and satellite network density) over the whole period.
-Operational forecasts of the French global NWP model ARPEGE from 2002 to 2021, which evolved over this time period with on average one major evolution per year.

125
The use of the ERA-40 reanalysis instead of the most recent ERA-interim or ERA-5 is inherited from previous work (Durand et al., 2009a, b). It is planned to use ERA-5 for future updates of this reanalysis. Outputs from the above-mentioned NWP systems are used as a preliminary guess at a 6-hour time resolution of the main variables driving the evolution of the snow cover (see Table 1). This guess contains informations both at the surface and at different heights above the surface.

130
The guess for precipitation since 1958 to 1 August 2017 is obtained using the AURELHY analysis method (Bénichou and Le Breton, 1987) providing a 24-hour climatological accumulation for each massif, depending on the weather type (without any use of the precipitation fields from the NWP model). Since 1 August 2017, 24-hour cumulated ARPEGE precipitation fields are used as precipitation guess and the chronology of precipitation is taken from 6-hour ARPEGE precipitation. This change is due to the improvements in simulated precipitation and their availability at a 6-hour time resolution, which increases the 135 consistency of the analysis and meaningfulness of the computation of the phase of precipitation. However, it is not possible to extend this approach back in time, using the current NWP model output, since 6-hour resolution precipitation data in ARPEGE is only available since 2017.

140
The most influential observations analysed by SAFRAN are surface observations from various networks. This includes manual observations from the Météo-France climatological network and the dedicated snow observation network ("réseau nivométéorologique" in French) resulting from a collaboration between Météo-France and mountain stakeholders (in particular Domaines Skiables de France, Association Nationale des Maires de Stations de Montagne, Association Nationale des Directeurs de Pistes et de la Sécurité de Stations de Sports d'Hiver). The latter network has been progressively implemented since 145 the 1970s, with observations relevant to the mountain snow cover. These observations are a key part of the analysis system since they often are the only available information for a given area. They include a large range of meteorological and snow cover variables including past weather conditions, information on the rain-snow elevation or 24h height of new snow one or two times per day, which are used to ensure a consistent analysis and check automatic observations. Observations for Andorra and the Spanish side of the Pyrenees are provided by means of international collaborations for the exchange of snow and

155
The temporal evolution of the number of daily 2 m surface temperature and 24-hour precipitation observations available in the massifs and effectively used in the assimilation process is shown in Figure 4. Available observations (dashes) only refer to stations within the boundaries of the massifs while the analysis model can use observations coming from more distant low elevation areas (see Figure 2 for the extraction domains of observations). This difference explains that the daily average number of assimilated observations can be higher than the number of available observations within the simulation domain.
160 Figure 4 shows a generally growing number of available observations in the mountains over most of the period 1958-2020 in the three areas. This explains the increasing number of assimilated observations and suggests that the share of mountain observations in the reanalysis also rose. The stabilisation (for 2 m-temperature) or decrease (for precipitation) of the number of observations for the five last years of the period is due to a cost reduction in some manual measurements and automatic networks. The rapid increase of 2 m-temperature observations ( Figure 4) starting from the beginning of the 1990s can be explained 165 by the development of automatic networks providing hourly observations in mountain areas, thus adding many more observations than daily precipitation observations. It is an important source of temporal heterogeneity in the dataset (see section 5.2).
These figures also show that the three areas have different observation network densities, in particular there are significantly more observations in the French Alps than in the Pyrenees although these two areas have a similar number of massifs.   For this new version of the S2M reanalysis, a more complete set of observations than the previous one (Durand et al., 2009a) was used, using all observations available in the Météo-France database.

185
SAFRAN also uses two additional sources of information: -Radiosondes provide vertical profiles used to correct the upper part of the profiles coming from the guess. However very few radiosondes observations are available daily and they are often launched from distant areas.
-Satellite observations are used since 1991 to analyse the cloudiness based on 1 km observed cloud structures which can especially help detect total cloudiness and clouds in valleys.

Short description of SAFRAN
SAFRAN (Durand et al., 1993) is an atmospheric analysis system, which provides the main meteorological variables necessary to drive a land surface model (and in particular a snow cover model) at an hourly time step. Each SAFRAN analysis covers a 24-hour period from 6:00 on day D-1 to 6:00 on day D. It combines the gridded meteorological guess from ERA40 or ARPEGE NWP models (see section 2.2.1) providing vertical profiles of air temperature, humidity, wind speed and direction 195 every 6 hours. An estimation of the 24-hour cumulative precipitation over the analysis period is obtained either from a climatology or directly from the NWP simulated precipitation (see section 2.2.1).
The first step of the SAFRAN analysis of all variables except precipitation is to compute a meteorological guess on its massif geometry (see section 2.1) every 6 hours using the vertical profiles from the NWP model. That guess is corrected by the assim-200 ilation of a first set of surface observations (of 2 m-temperature, 2 m-humidity and 10 m-wind) available every 6 hours using an optimal interpolation. The weight of each assimilated observation is based on the horizontal distance between the observation and analysis points. These 6-hour values are then interpolated at an hourly time step. The daily evolution of 2 m-temperature is particularly important and the implemented scheme is based on Martin and Mainguy (1988). A first estimate of the maximum daily 2 m-temperature is performed using the analysis at 12:00. Then the 2 m-temperature profiles are temporally interpolated 205 using a diurnal adjustment depending on the other meteorological variables. If hourly observations are available, a variational assimilation of these observations produces the final hourly simulation for the main relevant atmospheric variables affecting surface processes (i.e., 2 m air temperature, 10 m wind speed, 2 m air humidity, cloudiness, long-wave incoming radiation, and direct and scattered solar radiation).

210
The precipitation analysis consists in an optimal interpolation between the guess and the different available observations of 24 h cumulative precipitation within the massif followed by a temporal distribution of hourly precipitation according to the NWP-model chronology when available (see section 2.2.1) or simulated hourly relative humidity. The phase of hourly precipitation is determined depending on the simulated 0 • C isotherm elevation for each aspect and potential past weather condition observations as well as an estimation of the daily fraction of solid precipitation mainly based on the manual observation net-215 work ("réseau nivo-météorologique").
Meteorological variables for each elevation of each massif are analysed using a maximum number of : Some final adjustments are made at the end of the analysis, to ensure the physical consistency, such as an adjustment of the elevation of the rain/snow limit with respect to 2 m-temperature.

225
The SURFEX/ISBA-Crocus (hereafter Crocus) model represents the snowpack as a stratified medium depicted by a dynamical number of numerical layers up to 50. The prognostic variables for each layer are snow mass, density, enthalpy (i.e. temperature and liquid water content), age, and complementary variables for snow microstructure (specific surface area and sphericity).
Evolution equations rely on the solving of the diffusion heat equation in this stratified medium with Neumann boundary condition. Phase changes (melting and refreezing) are computed assuming a decoupling with heat diffusion at the model internal

Short description of MEPRA
MEPRA is an expert model designed to estimate the avalanche hazard from the snowpack stratigraphy simulated by Crocus, 240 from mechanical diagnosis and expert rules (Giraud et al., 2002). It has been fully implemented in the SURFEX platform, and its outputs are computed and made available with the other diagnostic variables. The general concept in MEPRA is to compare the shear strength to the shear stress in each snow layer. The shear strength is parametrised as a function of density and microstructure variables. For natural release, only the weight of overlying layers is taken into account in the shear stress.
An additional load is added in order to compute accidental triggering. Expert rules are defined to compute a hazard index 245 from these mechanical stability indicators, both for natural release and accidental triggering. The system is only applied on 40 degrees slopes, at a internal time step of 3 hours. The S2M reanalysis spans the time period from 1st August 1958 at 6:00 UTC until 1st August 2021 at 6:00 UTC. The atmospheric (FORCING.nc) and snow cover (PRO.nc) datasets are each stored in 63 annual NetCDF files for each mountainous area (French Alps, Pyrenees and Corsica). Two kinds of simulation are available: one on flat terrain only, and one taking slopes into account by projecting the incoming radiation variables according to the slope value and its aspect.

Metadata
The S2M dataset consists of two shapefiles containing all informations concerning the geometry of the simulation (for both massifs and stations geometries) and NetCDF files containing the simulated data. Simulated variables have two main dimensions : a temporal one ("time", the number of days from the previous 1st of August) and the total number of simulation points ("Number_of_points"). Geometry variables that can be used to retrieve specific sub-data from the full yearly files are sum-260 marised in Table 2. For example to get snow depth simulation over south heading slope, 40 • steep in massif number 1 at 2400 m a.s.l elevation, select the following characteristics for variable "DSN_T_ISBA" : A practical code in python is provided to get the simulated snow depth evolution over south heading slope, 40 • steep in massif number 1 at 2400 m a.s.l elevation, as well as for one specific date. It also show how to plot the Alps massifs from the corresponding shapefile.

Meteorological variables 270
Meteorological fields are provided at an hourly time step. The list of SAFRAN output variables is summarised in Table 3. The 0 • C isotherm elevation and the rain snow limit are are computed using the whole vertical profile simulated by SAFRAN for each massif.

Snow cover and soil variables 275
Files containing snow cover and soil variables include data at a daily time resolution (with state variables provided at 6:00 UTC). The list of the output snow cover and soil diagnostics is given in

Impact of the change of precipitation guess
The impact of the temporal heterogeneity introduced by the different precipitation guess has been evaluated over one single season (2017-2018) by comparing a simulation made with precipitation guess based on the AURELHY analysis method to the reference one (with precipitation guess from ARPEGE). The comparison of the simulated snow depth mean deviation and root mean square deviations (RMSD) for these 2 configurations didn't show any significant impact on the performance of the 295 system. Furthermore, the simulated annual precipitation amount with the precipitation guess based on the AURELHY analysis method seems to be slightly higher (about 3% on average) than those simulated with precipitation guess from ARPEGE : for the season 2017-2018 the average total precipitation over the 665 stations of the simulation with the guess based on the AURELHY analysis method and from ARPEGE respectively is 1507 mm (resp. 1464 mm) with accumulation ranging from 728 mm (resp. 675 mm) up to 3125 mm (resp. 3242 mm).  Figure 7).
This visualisation highlights the strong inter-annual variability of the winter snow cover, with fluctuations that can be greater than 100% of the mean value over the reference period. It also points the spatial variability for one given year with anomalies that can be positive for some massifs and negative for other massifs.
330 Figure 9 explores the vertical dimension of the dataset, as well as the simulated trends for the different seasons of the year. It compares the mean of the different variables over the Alps for each season at different elevations over the periods 1960-1990 and 1990-2020 (left), as well as the difference between the two periods (right). Only elevations between 900 m a.s.l and 3000 m a.s.l are considered because the lack of input observations outside this range of elevations (see Figure 5) reduces the quality of 335 the reanalysis.

340
-In fall and winter the simulation has a slight negative trend of about -0.1 • C per decade. Beaumet et al. (2021) compared the 2 m-temperature trends of the S2M reanalysis to that of the ones simulated by a model run of the MAR regional climate model driven by the ERA-20C reanalysis and SPAZM reanalysis (Gottardi et al., 2012) and showed that the S2M reanalysis simulates smaller trends both in winter and at high and low elevations in summer. Further 345 analysis of these simulated trends and a comparison to observations is presented on section 4.3.3.
Concerning precipitation, Figures 9 (c, d) shows that the mean of total precipitation in S2M over the Alps in the last three increase of the fraction of solid precipitation (Figure 9, f). The comparatively low amplitude of the increase of the simulated snow depth (only few centimetres, Figure 9, h)) may be explained by an averaging effect and the very short lifespan of snow on ground at this season. software (Mestre et al., 2013). However it is important to consider that the corresponding raw daily observations are assimilated in the S2M reanalysis described in this paper and therefore do not constitute a fully independent evaluation dataset. Consequently, another reanalysis was performed after removing these observations from the assimilation process for an independent evaluation. A third reanalysis was made without the assimilation of any 2 m-temperature observation to identify the impact of the assimilation of these observations on the simulation performance. For the evaluation of 2 m-temperature simulations, it is necessary to take into account the bias due to the hourly resolution of SAFRAN simulations whereas the observed minimum and maximum 2 m-temperatures often stands in-between two hourly 2 m-temperatures. These biases can be estimated by comparing observation series of hourly 2 m-temperatures and hourly minimum and maximum 2 m-temperatures using the raw hourly observations of the evaluation sites (available since the mid 1990s 375 until August 1st 2020). These data show that the minimum value of hourly 2 m-temperatures is on average about 0.2 • C higher than the absolute daily minimum and the maximum value of hourly 2 m-temperatures is on average about 0.4 • C lower than the absolute daily maximum. The 0.2 • C gap between the minimum and maximum 2 m-temperatures shifts highlights that the daily 2 m-temperature evolution is not a perfect sinusoid and presents a day/night asymmetry (Martin and Mainguy, 1988) : a typical daily 2 m-temperature evolution curve for a clear sky day exhibits a narrow maximum and flat minimum. Thus,

Evaluation of minimum and maximum 2 m-temperatures and precipitation
The box plots in Figure 10 show the variability of the scores among the different sites within each elevation range and the 395 notches indicate the confidence interval of the median obtained by a bootstrap sampling of the considered stations. Root mean square deviation of maximum temperature (°C) Reference reanalysis Renalysis without evaluation observations Reanalysis with no temperature observation (b) RMSD of minimum 2 m-temperature 1960-1990 1990-2012   (see Figure 4). This temporal heterogeneity is even more significant since an important part of these new observations are hourly observations. These hourly observations are crucial to accurately simulate the diurnal cycle of 2 m air temperature. Thus the assimilation of more observations over time tends to improve the simulation of the daily variations of 2 m air temperature and to bring it closer to observations for both maximum and minimum temperatures. However, the reduction over time of the warm bias in minimum 2 m-temperature is stronger than the reduction over time of the cold bias in maximum 2 m-temperatures.

420
The evaluation of precipitation exhibits no strong systematic bias for any simulation or period (Figure 11, b), but a stronger dispersion of the mean deviation for the independent simulation. This is confirmed by the fact that the RMSD (Figure 11, a) is by far higher for the independent simulation than for the two other simulations despite an improvement over time. As expected, removing all 2 m-temperature observations from the assimilation process does not affect much the precipitation analysis (only 425 minor effects but a major impact exists on the precipitation phase) but removing few precipitation observations (blue boxplots) has a stronger negative impact due to the small number of available precipitation observations (see Figure 4).

Trends of minimum and maximum 2 m-temperatures and precipitation
Here, climatological trends are defined by the difference between the mean of a variable over two 30-year long periods (e.g. -2020 and 1960-1990).  The magnitude of this warm bias on maximum 2 m-temperature only partially balances the cold bias on minimum 2 mtemperature, resulting in an artificial cooling of the simulation relatively to observations in terms of mean air temperature at 2 m. This pattern is even more pronounced when considering only the winter season (see Figure 16 of Appendix B), which 445 explains the negative 2 m-temperature trend simulated in winter as noticed in section 4.2. Besides, Figure 12 (b) shows that the simulated trends of maximum air temperature at 2 m in summer significantly increases with elevation up to about 1800 m a.s.l.

1990
However this elevation dependency of the simulated trend of maximum 2 m-temperature is not visible in station observations.
The dispersion of the observed trends of precipitation in Figure 12 (c) indicates a high spatial variability of these trends. The fact that the simulated trend of the mean precipitation over the Alps stands in the middle of the scatter plot suggests a good 450 agreement between simulated and observed trends of precipitation both in summer (no simulated trend) and winter (slightly negative simulated trend).
We computed two scores for each observation station : mean deviation and RMSD, by taking into account data from October 1st to June 30th of each year.  This evaluation shows that the S2M reanalysis is able to simulate a variable that is not assimilated and cumulates errors from both the meteorological analysis and the snow cover model with no systematic bias and moderate deviations to observations.
Another evaluation of the performance of the simulation over time is presented in Figure 15. It compares the mean devia- Eastern Pyrenees, and 20 Corsica). See Figure 2 and Table 8 for more details.
observations assimilation to their corresponding observations for the last four decades. While the reference S2M reanalysis shows no systematic bias for the four decades, removing 2 m-temperature observations from the SAFRAN analysis introduces a constant negative bias of around 10 cm on average. In addition, the RMSD for the reference reanalysis is always lower by a few centimetres than the simulation with no 2 m-temperature observation (see Figure 15) despite a mean simulated snow 495 depth systematically higher by about 20%. For both simulations, the RMSD tends to decreases over time except for the last decade where the mean simulated snow depth is significantly higher. This consolidates the results of section 4.3.2, confirming that the reanalysis described in this paper provides an optimal simulation at all time by considering all available information at that time. The downside of this is that the simulated trends are not fully representative of the climatological trends due to the temporal heterogeneity of available observations. 500

Discussion
In this study, we introduce the open access S2M reanalysis, which provides meteorological and snow cover data for the French Alps, Pyrenees and Corsica for the time period 1958-2021. This dataset enables a large number of research and operational applications, but several limitations need to be considered. Here we review the main strengths and weaknesses of this dataset, 505 and their consequences.

Main assets and known uses
The S2M reanalysis makes it possible to take into account meteorological and snow cover conditions for any time between 1958 and 2021 in the French high elevation mountain regions. It supersedes the SCM reanalysis, which has been developed by Durand et al. (2009a) in the early 2000s. SCM and then S2M have been used to address a wide range of applications in 510 various scientific domains. The S2M dataset was used as a reference for evaluating snow cover simulations driven by the NWP model AROME (Queno et al., 2016;Vionnet et al., 2019). SCM and S2M reanalyses have been used in a number of studies addressing glacier mass balance in the French Alps (Gerbaux et al., 2005;Réveillet et al., 2018;Bolibar et al., 2020;Peyaud et al., 2020) or hydrological simulation of alpine mountainous catchments (Lafaysse et al., 2011). Pellarin et al. (2016) used S2M soil thermal state and the overlying snow cover to investigate the potential of L-band satellite measurements to improve 515 soil moisture retrievals. Although climate trends may be questionable, the S2M dataset provides a robust climatological baseline for mountain regions. It has been used as a reference for adjusting climate change projections with statistical downscaling techniques (Lafaysse et al., 2014;Verfaillie et al., 2017). This methods leads to the provision of meteorological forcing driving files corresponding to future climate time series, on the same geometry and data format as S2M reanalysis forcing files, enabling homogeneous post-processing including running SURFEX/ISBA-Crocus for natural (Verfaillie et al., 2018) and man- dation for innovative developments towards the assimilation of remotely sensed and in situ snow cover observations in the simulations (Viallon-Galinier et al., 2020;Deschamps-Berger et al., 2020;Cluzet et al., 2021).

Temporal heterogeneity
A number of studies rely on S2M for the analysis of meteorological and snow cover trends at climatic scale Lopez-Moreno et al., 2020). For these applications, the main caveat in using S2M lies in the temporal heterogeneity of the dataset, in particular the strong changes in 2 m-temperature observations across the mountain regions considered during the full 530 time period of the reanalysis. The S2M reanalysis for the period 1991-2020 is produced using significantly more observations than the simulation over the period 1961-1990 ( Figure 4) and we established in section 4.3.3 that this has a significant impact on the simulated 2 m-temperature trends. This is superimposed on changes in meteorological guess in 2002 (from ERA-40 to ARPEGE) and changes intrinsic to ERA-40 and ARPEGE: assimilation of satellite observations in ERA-40 is known to be responsible for temporal breaks (Sturaro, 2003;Sterl, 2004). Over the time period from 2002 to 2021, the ARPEGE analysis was 535 even more affected by various changes in its physical parametrisation as well as its changes of horizontal resolution over time.
Similarly to many other reanalysis systems, the original purpose of the S2M reanalysis was not to provide a system dedicated to the analysis of climate trends, but rather the best available estimate of meteorological and snow cover conditions for every day within the covered time period, hence heterogeneity was allowed in the input data to the S2M reanalysis. This choice is corroborated by the various temporal scores presented in this study showing a clear improvement of the simulation over time 540 correlated to the growing number of available observations. The downside of this temporal refinement is the introduction of artificial biases in various simulated variables that can either compensate or amplify actual climatological trends as presented in this study for 2 m-temperature in winter. Thus, trend analysis using S2M must be considered with high caution, especially for air temperature at 2 m, but also probably for other sensitive variables such as snow depth (see Figure 4 of Verfaillie et al. (2018) at Col de Porte) although very few observations series allow an accurate characterisation of long term trends.

545
A further evaluation of the simulated snow cover duration values is planned in the coming years, for example products from Hüsler et al. (2014) as evaluation data.

Future updates
The S2M dataset is intended to be updated every year to extend its time coverage period and take into account evolutions of 550 the various components of the model chain. This may impact, in the future, the simulations presented in this article. Replacing ERA-40/ARPEGE by ERA-5 is under preparation and is expected to reduce temporal heterogeneities -although changes in observation data are likely to remain the dominant source of heterogeneity. It is also planned to expand the S2M reanalysis to lower-lying mountain ranges Vosges, Jura and Massif Central.

555
This study introduces and describes the latest release of the meteorological and snow cover reanalysis S2M covering the 63-year period from 1958 to 2021 for the French mountainous areas (French Alps, Pyrenees and Corsica). It includes the description of the different models and data used to produce this dataset, a comprehensive list of the parameters forming the dataset itself in its specific geometry and the technical access to the database. An evaluation of the simulation quality is provided by comparison 560 to in-situ observations of snow depth as well as an overview of known and potential uses of this dataset, and a series of caveats associated with the use of this dataset. Yearly updates of this reanalysis will extend the period in the future and could lead to significant updates of the current S2M data by including model or input data modifications.

Data access 565
The S2M dataset is freely available on the AERIS data center on the following https://doi.org/10.25326/37. It is required that scientific publications using the S2M dataset mention in the acknowledgements: "The S2M data are provided by Météo-France -CNRS, CNRM, Centre d'Études de la Neige, through AERIS" and cite this article as reference and refer to the dataset as Vernay et al. (2022).
To access the data five fields are required : CNRM/CEN is part of LabEx OSUG@2020. 9 Appendix B (a) RMSD of maximum 2 m-temperature in winter 1960-1990 1990-2012 Time period Root mean square deviation of minimum temperature (°C) Reference reanalysis Renalysis without evaluation observations Reanalysis with no temperature observation (c) Mean deviation of maximum 2 m-temperature in winter 1960-1990 1990-2012 Time period  Root mean square deviation of minimum temperature (°C) Reference reanalysis Renalysis without evaluation observations Reanalysis with no temperature observation (c) Mean deviation of maximum 2 m-temperature in summer 1960-1990 1990-2012