Past changes in natural and managed snow reliability of French Alps ski resorts from 1961 to 2019

Snow reliability is a key climatic impact driver for the ski tourism industry, although there are only a few studies addressing past changes in snow reliability in ski resorts accounting for snow management practices (grooming and snowmaking, in particular). This study provides an assessment of past changes in natural and managed snow cover reliability from 1961 to 2019 in the French Alps. In particular, we used snowmaking investment figures to infer the evolution of snowmaking coverage at the ski resort scale for 16 ski resorts in the French Alps, which we used together with a detailed snow cover modelling system driven by a local atmospheric reanalysis. We find different benefits of snow management to reduce the variability and long-term decrease in snow cover reliability because of the heterogeneity of the snowmaking deployment trajectories across ski resorts. The frequency of challenging conditions for ski resort operation over the 1991–2019 period increased in November and February to April compared to the 30-year reference period 1961–1990. In general, snowmaking had a positive impact on snow reliability, especially in December to January. While for the highest-elevation ski resorts, snowmaking improved snow reliability for the core of the winter season, it did not counterbalance the decreasing trend in snow cover reliability for lower-elevation ski resorts and in the spring.

Abstract. Snow reliability is a key climatic impact driver for the ski tourism industry, although there are only a few studies addressing past changes in snow reliability in ski resorts accounting for snow management practices (grooming and snowmaking, in particular). This study provides an assessment of past changes in natural and managed snow cover reliability from 1961 to 2019 in the French Alps. In particular, we used snowmaking investment figures to infer the evolution of snowmaking coverage at the ski resort scale for 16 ski resorts in the French Alps, which we used together with a detailed snow cover modelling system driven by a local atmospheric reanalysis. We find different benefits of snow management to reduce the variability and long-term decrease in snow cover reliability because of the heterogeneity of the snowmaking deployment trajectories across ski resorts. The frequency of challenging conditions for ski resort operation over the 1991-2019 period increased in November and February to April compared to the 30-year reference period . In general, snowmaking had a positive impact on snow reliability, especially in December to January. While for the highest-elevation ski resorts, snowmaking improved snow reliability for the core of the winter season, it did not counterbalance the decreasing trend in snow cover reliability for lower-elevation ski resorts and in the spring.

Introduction
Ski tourism is a major socio-economic component of mountainous regions for many countries around the world (Vanat, 2020). Due to its reliance on the snow cover, ski tourism has long been identified to be particularly vulnerable to climate change impacts, such as decreases in snow cover duration at low elevation in mountain regions (Martin et al., 1994;Beniston et al., 1995;Abegg et al., 2017). In order to reduce the impact of natural snow cover variability and its long-term decreasing trend due to climate change, snowmaking has emerged and is nowadays routinely used in almost all ski resorts in developed countries . Recent progress has been made in the literature to better account for snowmaking in future projections of snow reliability in ski resorts in Europe at the pan-European scale , in separate European countries (Austria -e.g. Marke et al., 2015; France -e.g. Spandre et al., 2019;Norway -e.g. Scott et al., 2020a) and in North America (Scott et al., 2020b). However, while several studies have documented strong reductions in snow cover amount, depth and duration in many mountain regions of the world over the past decades (Mote et al., 2018;Klein et al., 2016;Marty et al., 2017;Matiu et al., 2021;Hock et al., 2019), explicit assessments of the impact of climate change on ski resort operations, based on past observations, have remained limited (Beaudin and Huang, 2014;Hamilton et al., 2003;Rutty et al., 2017;Steiger, 2011). In fact, the emergence in past decades and popularization in major news media of several studies addressing future climate change risks to ski tourism, together with the recurrence of snow-scarce winters in Europe (Durand et al., 2009) and North America (Cooper et al., 2016), have somehow led to the broadly accepted consideration that climate change has indeed already been having a strong impact on ski resort operations in the past and at present (Knowles and Scott, 2020). Based on interviews with ski resort managers and ski tourists, several studies have documented the fact that many stakeholders from the ski tourism industry have perceived a change in meteorological and snow conditions, which they did not always attribute to climate change (Trawöger, 2014;Bicknell and McManus, 2006;Morrison and Pickering, 2013). In fact, a substantial body of evidence on climate change impacts on ski tourism has been inferred primarily from future climate change projections rather than from the analysis of past changes in snow conditions. In a way, the link between climate change and snow cover reliability appears to be so direct and obvious (less snow in a warmer world) that it has been somehow neglected in previous studies addressing past changes, with only a few exceptions. Alongside this knowledge gap, the efficiency of snowmaking based on past observations has seldom been assessed quantitatively: this also requires in-depth analysis of ski resorts' operating conditions for past seasons in a way that allows disentangling the contribution of snowmaking from actual snow conditions.
Previous studies have led to the development of a sophisticated modelling system to simulate snow conditions in French ski resorts for the past and future, taking explicitly into account grooming and snowmaking, based on the meteorological reanalysis system SAFRAN and the snow cover model Crocus and accounting for several key features of individual ski resorts (Spandre, 2016;Spandre et al., 2019). This model chain takes as input, for each ski resort, the fraction of ski slopes covered with snow guns. However, in past studies the time evolution of snowmaking was applied uniformly to all ski resorts (Spandre et al., 2019) in the absence of resortlevel data describing the past time evolution of snowmaking fractional coverage in ski resorts. A proper assessment of the impact of the interannual variability and long-term climate change on ski resorts' operating conditions requires taking into account not only observed meteorological conditions driving the time variations in snow conditions in ski resorts but also the time evolution of their individual snowmaking capacity. In the French Alps, snowmaking emerged in the 1990s, but not all ski resorts have followed the same development pathway in terms of snowmaking, and there is no consolidated database describing this past evolution (Berard-Chenu et al., 2020). For a given change in a climatic impact driver operating at the scale of an entire mountain range, various snowmaking equipment pathways in different ski resorts are thus anticipated to lead to different impacts.
In this study, we quantitatively assessed the changes in operating conditions of ski resorts, with and without snowmaking, taking into account the variability in snowmaking development dynamics across ski resorts. To do so, we analysed snowmaking investment figures for 16 ski resorts in the French Alps from 1997 to 2018 and developed an original method to infer the time evolution of the snowmaking coverage for each ski resort. These data were used to feed the SAFRAN-Crocus model chain, making it possible to perform simulations with and without snowmaking for the past decades. Indeed, there is no long-term record of snow conditions in ski resorts that could be used for this purpose, and there is no way to reconstruct what would be the snow cover situation without snowmaking once it has been implemented in actual ski resorts. Numerical simulations make it possible to compare the state of the snow cover with and without snowmaking, everything else being equal. Section 2 introduces the methodology used for this study; Sect. 3 provides the results; Sect. 4 discusses them, and Sect. 5 concludes.

Material and methods
In the absence of a consolidated database providing the time evolution of the snowmaking coverage in each ski resort, this study uses investment figures related to snowmaking infrastructure in order to infer the time evolution of snowmaking coverage in individual ski resorts in the French Alps. To the best of our knowledge, this approach is unprecedented in the scientific community. Such estimates are then used to process snow cover simulation produced using the SAFRAN-Crocus model chain, with and without snowmaking, spanning the time period from 1961 to 2019, enabling us to quantify past changes in snow reliability at the ski resort scale, explicitly accounting for snowmaking equipment dynamics in individual ski resorts.

Snowmaking investments dataset
Investment figures come from a snowmaking investment dataset for 100 French Alps ski resorts from 1997 to 2018 (Berard-Chenu et al., 2020). The dataset originates from the professional journal Montagne Leaders, which manages a yearly investment survey completed in a declarative manner by ski lift operators. Although collected by a non-scientific and unofficial organization, these data hold significant value. For instance Falk and Vanat (2016) and Berard-Chenu et al. (2020) used this data source for academic research purposes. The investment survey covers five types of investment: new ski lift, ski lift maintenance (i.e. servicing and replacement of parts on existing ski lifts and compliance investments), snowmaking, ticketing and ski slope remodelling. Each type of investment is broken down into subgroups. Snowmaking investments are decomposed into six types: equipment and concepts, water reservoirs, electricity, civil engineering, buildings, and other. In our investment dataset we do not consider all snowmaking investments but only those referring to "equipment and concepts". Indeed, the subgroup of equipment and concepts is the most relevant proxy to capture the growth rate of the snowmaking coverage of ski resorts as it only accounts for the purchase of snow guns. The other snowmaking investments can refer to an increase in the snowmaking instantaneous production capacity, i.e. an electrical power increase or an expansion of the water storage capacity. These investments may increase the snow-making production efficiency without any growth of the ski slope areas equipped with snowmaking facilities. The investment dataset only contains ski resorts which have invested at least once in snowmaking equipment and concepts over the 1997-2018 period. Investments are in constant prices: we derived constant price investments from current price investments with the GDP deflator, a widely used inflation adjustment index.
The literature points to a broad range of figures on cost investment for snowmaking facilities. Table 1 shows several estimates related to snowmaking investment costs per unit ski slope surface area, ranging from EUR 28 000 to EUR 240 000 per hectare. Furger (2002) and RMS (Remontées Mécaniques Suisses), the professional association of the Swiss ski resort operators (RMS, 2011) mentioned a cost investment in Switzerland but expressed it per kilometre of slope equipped and not per hectare as in the figures in Table 1. Furger (2002) estimated the investment costs of snowmaking to be CHF 1 million per kilometre of slope equipped based on a sample of 10 ski resorts located in the canton of Vaud. In France, the latest estimate from Domaines Skiables de France dates back to 2018, in an internal training document intended for snowmakers, with a unit cost ranging between EUR 130 000 and EUR 150 000 per hectare.

Snowmaking fractional coverage reference dataset
Atout France, which is the national organization responsible for monitoring the tourism sector, has not published estimations of snowmaking facilities since 2009. These estimates were the only reference for professionals (Badré et al., 2009) until the study of Spandre et al. (2015). The latter assessed that 32 % of ski slope areas in the French Alps were equipped with snowmaking facilities in 2014, based on a survey involving ski resort managers. They also predicted that this proportion was likely to reach 43 % by 2020. DSF (Domaines Skiables de France), the professional association of the French ski resort operators, produces a snowmaking coverage average rate based on information collected from its members on an annual basis. In its most recent release in 2020, it indicated that 37 % of ski slope areas in France were equipped with snowmaking facilities. A comparison between the DSF snowmaking facilities rate in 2020 and the estimation of Spandre et al. (2015) is difficult since the latter only considers the French Alps and DSF does not provide any information regarding the representativeness of its sample and how the aggregated value was derived for all ski resorts in France.
In contrast to the national situation lacking quantitative estimates of snowmaking coverage for individual ski resorts, such information is available for a subregion of the northern French Alps Département de la Savoie) which corresponds to the NUTS-3 level in the European Union nomenclature of territorial units. Indeed, the Direction Départementale des Territoires de la Savoie, an administration in charge of the water policy, has initiated an observatory of snowmaking. Each year this public body collects snowmaking data from ski resorts in Savoie, among which is the snowmaking fractional coverage. Among the 32 ski resorts equipped with snowmaking facilities and registered in the Direction Départementale des Territoires de la Savoie's database, we removed 3 ski resorts with missing or questionable data. We have a sample of 29 ski resorts after data cleaning. For each ski resort, information is available for their total ski slope area and the area equipped with snow production facilities for the years 2010, 2011, 2015, 2016, 2017 and 2018.

Relationship between snowmaking investments and snowmaking fractional coverage
We developed a model analysing the relationship between the snowmaking surface area coverage with snowmaking investments. Among various possibilities, we retained a model relating them in a linear way (an exponential model was also tested, leading to insignificant differences for ski resorts with sufficiently numerous observation data).
The linear model to be estimated is where t denotes the year. The left-hand variable SMA denotes the ski slope area equipped with snowmaking facilities and I the cumulated snowmaking investment. α represents the snowmaking unit surface area equipment cost, and ε is an error term. We estimated the area equipped with snowmaking facilities based on the area equipped and the snowmaking investment increase both from the previous year. Among the 29 ski resorts equipped with snowmaking equipment, 13 showed no increase in terms of their snowmaking coverage despite the fact they made snowmaking investments over the study period. We removed these 13 ski resorts from our sampling since an increase in the snowmaking coverage rate over the 6 reference years is required to estimate the α coefficient in our modelling. Thus our final sample for the evolution of the individual snowmaking coverage encompasses 16 ski resorts.

SAFRAN-Crocus model chain, ski resorts geospatial modelling and snow reliability indicator
We used the SAFRAN-Crocus model chain to simulate snow cover characteristics from 1961 to 2019. SAFRAN is a meteorological reanalysis system, combining in situ observations and numerical weather prediction model output to provide an estimate of the meteorological conditions as a function of elevation (by 300 m elevation steps), for "massifs", i.e. mountain areas assumed to be meteorologically homogenous. This system makes it possible to reconstruct the time evolution of the meteorological conditions in the mountain regions of France since the late 1950s (Durand et al., 2009;Vernay et al., 2019Vernay et al., , 2022 and has been used as a basis for several pre- vious studies on the snow cover reliability of ski resorts in the French Alps and Pyrenees (e.g. Spandre et al., 2019, and references therein). The hourly resolution SAFRAN meteorological conditions are used as input to the detailed, multilayer snow cover model Crocus, resolving natural processes occurring in the snow cover and at its interfaces with the underlying ground and overlying atmosphere (Vionnet et al., 2012). In its Resort version , the Crocus model is also equipped with dedicated options for representing snow grooming and snowmaking, in terms of both their physical characteristics and their typical management operation timing and snowmaking production rate (Spandre et al., , 2019Hanzer et al., 2020;Morin et al., 2021). Crocus model runs were carried out for each relevant SAFRAN massif and all elevations, not only on flat terrain but also for eight main orientations (N, NE, E, SE, S, SW, W, NW) and slope angles of 10, 20, 30 and 40 • . Simulation outputs were then aggregated for individual ski resorts at a daily scale used to calculate a reliability index for different periods across the ski season. For each ski resort, a ski resort gravitational envelope is computed based on a digital elevation model (DEM) and the location of ski lifts (François et al., 2014Spandre et al., 2019). Depending on the snowmaking coverage, simulations carried out with and without snowmaking are combined, as described in Spandre et al. (2019). The location of ski slopes equipped with snowmaking is determined based on their position within the ski resort (elevation, distance to major ski lift, distance to main housing infrastructure), following Spandre et al. (2016). This entire procedure enables the generation of, for a given snowmaking coverage, a resort-level snow cover indicator, quantifying the fraction of the ski resort surface area where the amount of snow exceeds a given threshold (here, 100 kg m −2 snow water equivalent, which corresponds to 25 cm of snow with a density of 400 kg m −3 ) . Daily values of this indicator are averaged for each month from November to April, as well as for the Christmas time period from 20 December to 5 January. This approach is comparable to the "snow reliable skiing terrain" approach developed and used by Steiger and Stötter (2013) and Scott et al. (2019). We also computed the combined reliability index, corresponding to the weighted  (2019) continuously, whose values can be computed depending on the snowmaking fractional coverage value, with a regular 15 % step between 0 % and 90 % where 0 % corresponds to grooming only. The time-variable snowmaking coverage where considered a linear interpolation of the values obtained fractional coverage by steps of 15 %. In addition, based on the gravitational envelope computed for each ski resort, we compute the fraction of the ski resort above 2000 m elevation. This elevation threshold roughly corresponds to the divide between low-elevation and high-elevation mountain areas, most relevant in the context of elevation-dependent natural snow cover trends (Matiu et al., 2021) in the European Alps. It enables us to compare ski resorts' geographical characteristics, not only in terms of their minimum-maximum or mean elevation but also using a single indicator representative of their elevational coverage.

Statistical analysis of snow reliability indicators
The annual values of the annual-scale indicator values were post-processed as follows. Indeed, the study focuses on ski resort operations and how they have been modified in the past decades, accounting for the role of snowmaking. What matters most for ski resorts is not the mean multi-annual snow reliability or other metrics characterizing their multi-annual average conditions but rather how frequently challenging operating conditions are encountered. Indeed, snow-scarce winters, primarily related to the interannual variability in meteorological and snow conditions and influenced by long-term climate trends, are a key concern for the ski tourism industry Morin et al., 2021). We therefore focus on the characteristics of challenging winter seasons and whether their frequency has changed in time. We take as a reference the value of the 20th quantiles of snow reliability indicators for the 30-year long 1960/61 to 1989/90 (i.e. 1961-1990) reference time period (referred to as Q20 values), which corresponds to the typical frequency of challenging snow conditions every 1 in 5 years and is relevant for ski resort operators. This value depends on the ski resort (because of their various elevations and elevation spans), on the snowmaking fractional coverage and on the time period of the year. We then analyse, for the 29-year time period 1990/91 to 2018/19 (i.e. 1991-2019), the frequency of winter seasons for which the snow reliability values are below the reference Q20. If the value corresponds to more than 20 % of the years (i.e. 6 years for a full 30-year time period), this means that the snow conditions have worsened compared to the reference period and vice versa. This approach enables fair comparisons across time, for a given ski resort, without comparing, in absolute terms, the mean and low-quantile values characterizing snow conditions in different ski resorts or with different snow management options. The absolute values for mean or Q20 snow reliability values for different ski resorts can be very different and need to be interpreted separately in assessments at the scale of ski resorts, together with the business model of the tourism destination encompassing each ski resort, within which the reliability of the snow cover on ski pistes can play a very different role. For example, large high-elevation ski resorts are based on high reliability values with low interannual variability and offer few alternatives to downhill skiing in case of ski resort closure, while lower-elevation ski resorts embedded in a diversified destination are less dependent on a reliable snow cover on ski pistes, which is often less guaranteed in such cases (Luthe et al., 2012;George-Marcelpoil and François, 2016). Analyses and presentations were undertaken in open-source R software (R core team, 2020), along with the ggplot2 package (Wickham, 2016). Figure 1 presents a schematic description of the method and the data processing.

Relationship between snowmaking investments and fractional coverage
The linear modelling for reconstructing snowmaking coverage, described in Sect. 2.3, was applied to 16 ski resorts in Savoie over the observed period (2010-2018). Figure 2 shows the estimated evolution of the snowmaking coverage rate over 1997-2018 for these 16 ski resorts. Estimates of snowmaking coverage are derived from the individual investments of each ski resort. The figure highlights the heterogeneity of the snowmaking deployment trajectories depending on ski resorts. Investment choices directly affect the increase in the snowmaking coverage for each ski resort: the snowmaking coverage rates range from 0 % to 37 % in 1997 and from 14 % to almost 70 % in 2018. The average of the evolution of the snowmaking coverage based on the 16 ski resorts (black line) is higher than the national average (dashed line) throughout the period considered. The figure shows that a larger difference is reached for the year 2006 with a difference of almost 10 %, while in 2018 and 1997 the difference is 2 % and 5 %, respectively. The linear modelling estimated a mean alpha coefficient of 8.48 × 10 −3 , which corresponds to EUR 118 000 per hectare with a 95 % confidence interval (CI 95 % ) of [87,180]. A bootstrap method (n = 1000 re-samplings of coefficient values) was applied to provide the 95 % confidence interval for the coefficient estimate.

Time evolution of combined (Christmas and
February) snow reliability and relationship to snowmaking coverage Figure 3 displays the time evolution of the combined snow reliability indicator and snowmaking fractional coverage for two contrasting very large ski resorts in the Savoie from 1961 to 2019, with a snowmaking coverage below (reaching ca. 24 % in 2018) and above (ca. 67 % in 2018) the national average (estimated to be on the order of 38 % in 2018), respectively. While some interannual variability is observed and both ski resorts display contrasting geographical characteristics regarding the share of their ski area below 2000 m, in general the snow reliability values have remained rather high throughout the entire time period, except for the peculiar winter of 1963/64. Compared to a situation without snowmaking, the various snowmaking dynamics lead to very similar snow reliability values, with either no increase in reliability (when the reliability with groomed snow only is already very high and close to 100 %) or an increase on the order of 5 % to 10 %. In this example, ski resort no. 6 with the larger snowmaking coverage fraction shows lower interannual variability, when snowmaking is taken into account, than ski resort no. 12 with its lower snowmaking fractional coverage. Figure 4 shows the relationship between the annual snow cover reliability difference between configurations with snowmaking and the snow cover reliability with grooming (no snowmaking, x axis), for two contrasting ski resorts (lower-than-average snowmaking coverage with 83 % of the ski area below 2000 m on the left and higher-than-average snowmaking coverage on the right and with 30 % of the ski area below 2000 m). It highlights a large range of variability for the natural (groomed) snow reliability, from 15 % to 100 % for the ski resort with a greater elevation range below 2000 m elevation and from 80 % to 100 % for the ski resort with a larger fraction above 2000 m. The gain in snow cover reliability due to snowmaking, using either the time-varying snowmaking fractional coverage or its value in 2018 for all considered seasons, is generally larger than the difference between these two configurations and generally correlates with the natural snow cover reliability -i.e. the higher the natu-   indicate the corresponding snowmaking fractional coverage. Panels (c) and (f) represent the increase in snow cover reliability compared to the grooming-only situation. Ski resort no. 6 shows a higher-than-average snowmaking equipment value, while ski resort no. 12 shows a lower-than-average snowmaking equipment value. ral (groomed) snow cover reliability, the lower the gain from snowmaking -although the snow cover reliability values including snowmaking are increasingly lower than 100 % when the natural (groomed) snow cover reliability is lower. This shows that, in these cases, the natural snow cover deficit can almost be fully compensated for (but not entirely), for this indicator focusing on the Christmas and February time periods with high snowmaking coverage, except in the cases where the deficit is too large (in particular, winter 1963/64). Figures 5 and 6 show the time evolution of the monthly snow reliability indicators for the two contrasting ski resorts introduced above, showing the large differences between reliability indicator values across the different months, in particular the lower values at the beginning and end of the winter season, and peak values, on average, in January and February. Figure 7 shows the relationships between the monthly indicator values in December to April, for the same contrasting ski resorts in terms of their snowmaking coverage. It clearly outlines the differences in snow reliability between months and between ski resorts and how snowmaking can (or cannot) compensate for the natural snow cover deficit, depending not only on the meteorological conditions but also on the snowmaking fractional coverage of individual ski resorts, which depend on individual strategic choices (and the means to implement them). The gap between snow reliability gains for the two contrasting ski resorts is larger, regardless of the snowmaking configuration, at the beginning (December and January) rather than the peak (February and March) of the winter season. For the December index, there is a maximum gain below 25 % for ski resort no. 4 (with lower-thanaverage snowmaking coverage and 83 % of its ski area below 2000 m), while ski resort no. 14 (higher-than-average snowmaking coverage and 30 % of its ski area below 2000 m) can benefit from up to 50 % snow reliability gains. With a snow coverage rate corresponding to 54 % in 2018 (brown dots), ski resort no. 14 would have increased its snow reliability with snowmaking almost always up to 25 %, whereas its grooming-only snow reliability went down below 50 %.

Change in the frequency of challenging snow conditions for the 16 ski resorts in Savoie
In this section we analyse the change in the frequency of occurrence of challenging snow conditions from 1961 to 2019 by splitting the entire time period into two periods of approximately 30 years (1961-1990 and 1991-2019) and analysing how frequently, in the second time period, the snow conditions were better or worse than the Q20 threshold of the reference period . We provide the results for the 16 ski resorts in Savoie for which a dedicated analysis of the snowmaking coverage evolution was carried out. Table 2 shows a synthesis, for the 16 ski resorts of our sample, of the monthly reliability values (reference Q20, which characterizes the upper threshold of challenging snow conditions) and the evolution of the frequency of exceedance below this threshold in a situation without snowmaking. The largest snow reliability values for the reference period are found in February and March, with Q20 values ranging from 86 % to 100 % in February and from 90 % to 100 % in March. Snow reliability values in November are low, on the order of 0 (grey cells) to a 10 % maximum, increasing through December and January to February and March and with lower values in April. The table shows that the months most affected by a change from 1961-1990 to 1991-2019 are in November (starting from low values) and February to April, with different rates of increases in the frequency of challenging conditions progressing with the month (February through April) as a function of ski resort elevations (the largest and highest-elevation ski resorts show changes mostly in November and April, with changes in a growing number of months as the elevation range starts from lower elevations). The rate of change for the Christmas time period and January shows instead a decrease in the frequency of snow conditions below the reference value. This result is fully consistent with the analysis of natural snow cover trends in the European Alps (Matiu et al., 2021), showing contrasting patterns depending on the elevation (represented here through the fraction of the ski resort below 2000 m) and on the month of the year and marked decreases in snow indicators, especially at low to intermediate elevations throughout the winter season and for ski resorts positioned at higher elevations mostly during springtime. Table 3 shows the same results as Table 2 but accounting for snowmaking with a fixed snowmaking coverage for each ski resort, corresponding to the value reached in 2018. Compared to the snow reliability values of Table 2 with grooming only, the reference Q20 value, with snowmaking, is generally higher. However, with the difference in the November situation, the general trend in the frequency of exceedance below the Q20 reference time period shows a very similar pattern to that in the natural snow situations, with a shift in snow reliability values but the same temporal pattern. This shows that the trend towards increased snow scarcity is rather similar even with snowmaking but starting from a higher snow reliability level. Changes are most pronounced for the ski resorts with a larger fraction below 2000 m and for higher-elevation ski resorts only at the beginning and end of the winter season (November and March-April), with fewer negative changes (and even increases) in the core of the winter. Table 4 shows the results in the same form as Tables 2 and 3 but takes into account the fact that each ski resort has had different snowmaking coverage dynamics over the past decades. However, due to the fact that snowmaking became generalized in the 1990s, the Q20 reference values for the period 1961-1990 are roughly similar to the figures in Table 2. In general, the figure shows that snowmaking had a positive   impact on snow reliability in the ski resorts of this sample, especially in December to January and for the entire winter season for the highest-elevation ski resorts. For smaller and lower-elevation ski resorts, snow conditions with snowmaking during the period 1991-2019 are often worse than without snowmaking between 1961 and 1990, which means that in such cases snowmaking did not counterbalance the decreasing trend in snow cover reliability.

Discussion
In this study, we have quantified the change in snow reliability in several ski resorts of the French Alps from 1961 to 2019, during a time when natural snow cover evolution has shown very strong changes especially at low to intermediate elevations, mostly attributed to the effect of atmospheric warming due to climate change (Matiu et al., 2021). Our study based on a long time period makes it possible to analyse the changes in terms of the frequency of occurrence of challenging snow conditions, defined during the reference period 1961-1990 as the value separating the 20 % worst seasons (lower snow cover indicator values) from the full record.
Based on snowmaking investment figures for 16 ski resorts in Savoie (northern French Alps), we reconstructed their snowmaking coverage trajectory since 1997, enabling us to quantify the impact of snowmaking on the changes in snow reliability over the entire time period. This study fills several gaps in the scientific literature but comes with limitations which deserve some discussion.

Limitations to the methods employed
The SAFRAN meteorological reanalysis (Vernay et al., 2019(Vernay et al., , 2022 used for this research is affected by changes in observation network density from 1961 to 2019 and changes to the numerical weather prediction models used as a guess for producing the reanalysis. This can potentially impair its ability to be used for trend analysis, as noted in Spandre et al. (2015) and Ménégoz et al. (2020). In particular, Vernay et al. (2022) indicate that the magnitude of the wintertime temperature trend is probably underestimated in the SAFRAN dataset, owing to temporal heterogeneities in the input data to SAFRAN and stronger temperature deviations from reference observations at the beginning of the time period. The magnitude of the trend on snow cover (strongly influenced Figure 7. Relationship between the monthly (December, January, February, March and April, from top to bottom) snow cover reliability difference (y axis) between configurations with snowmaking and the snow cover reliability with grooming only (no snowmaking, x axis), for two contrasting very large ski resorts and for the time period from 1961 to 2019, on the left with below-average snowmaking coverage and on the right with above-average snowmaking coverage. by air temperature trends especially at lower elevation) is therefore potentially underestimated. However, we note that our analysis here is not based on direct trend analysis but rather on how frequently snow conditions, computed at the monthly scale, fall above or below a given threshold, which depends on the meteorological conditions throughout all the winter seasons. We also note that several key features (snowscarce winter 1963/64; low-snow conditions in the late 1980s and early 1990s; challenging snow conditions in 2001/02, 2006/07 and further recent years; etc.) are aptly reproduced by the model chain, hence providing confidence in its ability to characterize the main features of the meteorological conditions of the winters of the past. While further analysis of the time heterogeneity in the SAFRAN reanalysis is needed and is currently being investigated, we consider that, due to the fact that SAFRAN uses all available observations in a robust assimilation framework, it is a usable data source for this work.
Another limitation of this work, in the snow cover modelling part, is the fact that the Crocus snow cover model not only uses the same snow management configuration for all ski resorts but also uses the same one throughout the entire time period. Yet, snowmaking technology and management strategies and tactics have evolved over the past decades (Morrison and Pickering, 2013;Wilson et al., 2018). However, taking into account other sources of uncertainty and the Table 2. Frequency of occurrence of monthly snow reliability challenging conditions, between 1991 and 2019, below the reference Q20 value for the time period 1961-1990, without snowmaking throughout the entire time period. Each row corresponds to a ski resort (ranked as a function of the increasing fraction of the domain below 2000 m). Each column corresponds to a monthly reliability indicator. The values in the matrix correspond to the Q20 reference snow reliability indicator. Red and orange cells correspond to larger frequencies during 1991-2019 than during the reference time period (more often challenging conditions) and vice versa in blue. Table 3. Same as Table 2 but for the snowmaking coverage of 2018 throughout the entire period from 1961 to 2019 (including for the reference time period . Table 4. Same as Tables 2 and 3 but with a snowmaking coverage variable in time for each ski resort. fact that the time evolution of the snowmaking coverage in ski resorts is itself largely unknown, accounting for changes in snowmaking practices and for the technological characteristics of snowmaking units appears out of reach and may be approached in future studies, requiring that further data be made available (or reconstructed from technical historical sources from ski resort internal documentation, which does not always exist). We however consider that our study, which accounts for grooming and snowmaking and brings together a large amount of original information into a consistent analysis framework, is a relevant addition to the sparse literature on observed changes in ski resort operating conditions and provides relevant and valuable results. Another limitation related to changes in time for some of the characteristics of the systems studied is the fact that we used a fixed number of ski resorts across time (even if some opened after 1961), with constant levels of ski lift equipment, corresponding to the end of the time period. This could cause some heterogeneities due to the fact that ski slopes are related, in our work, to the location of ski lifts, and if major changes in ski lifts had occurred from the early 2000s when snowmaking developed the most, inconsistencies would occur. However, we believe that our approach is relevant because the modelling framework makes it possible to simplify this complex problem, to some extent, by setting some of its dimensions to a given value (here, the geometry of ski resorts) and because the main features of ski resorts have in fact evolved little since the early 2000s, so changes in operating conditions are more influenced by the inception and increase in snowmaking than by changes in ski resort spatial organiza-tion, e.g. related to the opening of a new subdomain owing to the installation of a new ski lift.
The modelling approach employed to estimate the evolution of the snowmaking coverage rate also has some limitations. In our study, the increase in the snow reliability is only possible through a growth in the part of the ski area covered with snowmaking. However, snowmaking investments considered in our study may correspond to other strategies of ski lift operators. Our model assumes that all the snowmaking investments always resulted in an increase in the ski slope area covered by snowmaking. Therefore we removed 13 ski resorts without any snowmaking coverage growth despite investments on purpose. This indicates that snowmaking investment choices are more complex than those taken into account in our linear modelling. For instance, our approach does not represent investments made to optimize existing systems -e.g. the replacement of obsolete snow guns -or to increase the density of snow guns on slopes already equipped. Moreover, based on discussions with snowmakers, we know the ski industry has experienced an ongoing harmonization of the calculation method of the snowmaking coverage rate for many years, including through the development of tools dedicated to snowmaking management. We can speculate that the more extended snowmaking systems are, the more they require investments for their maintenance. Considering our model, the constant relationship between investments and area equipped lets us think that the increase in snowmaking coverage may be overestimated. Our modelling also provides estimations considering the technology remains constant. As indicated above on the snow cover modelling side, we did not implement the potential effects of technological progress, particularly for recent snowmaking systems, that could lead to yield increases. Similarly, we did not consider variations in investment costs in snowmaking over time, especially the potential decrease in investment costs to cover ski slopes as snowmaking systems become increasingly popular or the influence of the snowmaking industry consolidation and its consequences on the prices.
Altogether and despite its limitations, we believe that our study provides original and relevant information on observed changes in ski resort operating conditions over climate timescales, contributing to filling a critical gap in the literature. We also note that many of these limitations can only be waived by an increasing availability of data relevant to snow management in ski resorts, making it easier to quantify the impact of climate change on their operating conditions and, as a result, better understand the resorts' sensitivity to changing climate conditions and contribute to climate change adaptation in this domain (Spandre et al., 2019;Gerbaux et al., 2020;Berard-Chenu et al., 2020;Morin et al., 2021).

Quantification of the impact of climate change on snow reliability, with and without snowmaking
This study provides an original appraisal of the time evolution of snow cover reliability, based on a detailed modelling framework enabling the computation of a resort-level snow reliability indicator, spanning the time period from 1961 to 2019 for multiple ski resorts. Our results indicate that, over the past 60 years in 16 ski resorts in Savoie, the climate conditions have become increasingly challenging in terms of the natural snow cover of ski resorts, especially for ski resorts with a significant fraction of their surface area below approx. 2000 m. This is particularly the case in the early (November, December) and late (March, April) season, with contrasting patterns depending on the fraction of the ski resort above 2000 m elevation. The results also indicate that not only natural snow cover reliability and its changes but also the benefits of snowmaking in terms of snow reliability significantly vary between ski resorts. While at the core of the winter season the snowmaking deployment trajectory has counterbalanced the decrease in natural snow cover reliability for ski resorts with a low-elevation domain and in most cases increased snow cover reliability, in the early season and late season snowmaking development has not counterbalanced the snow cover declining trend. In this respect, our results fully corroborate, in the case of the northern French Alps, the findings highlighted in the IPCC SROCC Summary for Policymakers (IPCC, 2019) that "[i]n nearly all high mountain areas, the depth, extent and duration of snow cover have declined over recent decades, especially at lower elevation (high confidence)" and that "[t]ourism and recreation, including ski . . . tourism . . . have . . . been negatively impacted in many mountain regions (medium confidence)" and "[i]n some places, artificial snowmaking has reduced negative impacts on ski tourism (medium confidence)". Past impacts of climate change on natural snow cover reliability in ski resorts depend more on the indicator chosen to quantify the impacts and on resort characteristics than on their mountain area they belong to. Indeed the added value of snow management and snowmaking to reduce these impacts is also heterogeneous across ski resorts within a given mountain area. In this sense, while aggregation of results at the scale of an entire mountain range (e.g. Spandre et al., 2019, for the French Alps) makes it possible to provide a compact and general picture of the impact of climate change, past or future, on the supply side of the ski tourism industry and raises awareness at the sectoral scale of climate change impacts and risks, it provides a lumped message, which is not applicable at the scale of individual ski resorts, although this is most often the most relevant governance level for mountain tourism development. Previous studies addressing large geographical domains have generally recognized that studies at the individual scale would be most appropriate to inform deliberation and decisions at the local scale, prompting more detailed study at the scale of individual ski resorts (e.g. Gerbaux et al., 2020;Steiger et al., 2019). The present study reinforces this position and reinforces the need to exercise extreme caution when attempting to derive general statements, at the local or regional levels, about the evolution of ski resort operating conditions. Furthermore, our results indicate that the trends in natural or managed snow cover reliability strongly depend on the period of the year (at the monthly scale), which supports the need to carefully choose relevant indicators depending on the ski resort considered : depending on their business model, size, elevation range and climate settings, key periods for ski tourism can be different depending on the ski resort under consideration. The snow reliability index and its derivatives remain one indicator among others to analyse the effect of climate conditions on the operation of the ski tourism industry. For instance, Mayer et al. (2018) illustrated that a thermal comfort index could have a more significant effect on ski demand than snow depth. Several microclimatic characteristics may influence both the operation and the visitation of the ski resorts. Therefore, the search for a single snow reliability indicator universally applicable to ski resorts of all sizes and settings is probably elusive.
4.3 Can the climate change risk reduction factor due to snowmaking implementation be calculated?
Besides our analysis of the trends in natural and managed snow cover reliability in ski resorts in the northern French Alps, which we show varies across ski resorts due to differences in geographical settings (in particular the elevation distribution) and the pace of snowmaking equipment investments over the past decades, our analysis provides some of the information required to assess whether snowmaking has indeed decreased the overall climate risk faced by ski resorts, due to climate variability leading to the occasional occurrence of particularly unfavourable meteorological and snow conditions at the scale of 1 year. In a context where both the climate and the extent to which snowmaking is implemented have changed, it is virtually impossible to assess the partitioning of the climate change risk into its different factors. Still, Domaines Skiables de France, the national organization bringing together French ski resorts, and the National Association of Mayors of Mountain municipalities (ANMSM, Association Nationale des Maires des Stations de Montagne) have stated that, thanks to snow management (grooming and snowmaking) and slope preparation, the risk induced by a given level of climate hazard has been divided by a factor of 3 over the past 25 years (Domaines Skiables de France and ANMSM, 2015). The modelling approach taken in the present study allowed us to assess how the current level of snowmaking coverage would have behaved against the situation which prevailed, for example, during the notoriously challenging winter seasons in the late 1980s and early 1990s or even the winter of 1963/64, which is, for our study domain, the most challenging winter season, in terms of snow conditions, of the entire record from 1961 to 2019. Figures 3, 5 and 6 show, for example, absolute values and differences to a no-snowmaking situation, in terms of annual and monthly snow cover reliability, including these challenging winters. Figures 4 and 7 show the relationship between reliability indicators accounting for a fixed, 2018-level snowmaking coverage and the snow cover reliability values obtained without snowmaking. These figures show that challenging winters without snowmaking remain challenging when accounting for snowmaking at current snowmaking equipment levels, although in many cases the snow conditions are improved but never, especially in the early (November and December) and late (March and April) season, to the extent that snowmaking fully compensates for challenging natural snow conditions. In addition, the evidence provided in this study shows not only that the snow reliability benefits of snowmaking for challenging winter seasons, i.e. snowmaking's capacity to compensate for challenging meteorological and natural snow conditions, are strongly heterogeneous across ski resorts due to the resorts' differences in geographical settings (in particular their elevation range) but also the key periods during the winter season which are most impactful on their economic results. Providing a single factor for the climate risk reduction due to snowmaking across the entire ski industry in a given region would require that a single indicator could be applied to all ski resorts uniformly and furthermore that this indicator would be directly related to their economic performance in order to apply a risk approach meaningfully and compute the risk reduction factor related to a given snowmaking equipment level. While our modelling framework makes it possible to enable further quantifications in the future along these lines, our results suggest that the climate change risk reduction factor for the past period from 1961 to 2019 will undoubtedly be different for each ski resort, so seeking a single climate change risk reduction factor for the entire ski tourism industry is most probably elusive and would mask out the strong heterogeneity in the climate change impact drivers (hazards) and exposure and vulnerability of the individual elements (ski resorts) forming the backbone of the industry as a whole.

Past performance is not a guarantee of future results
Our results are consistent with Gonseth (2013), who highlighted snow reliability gains with snowmaking. However this study stressed diminishing returns to snowmaking investments, similarly to Falk and Vanat (2016). Gonseth (2013) also pointed out snow reliability gains expected with snowmaking remained based on a dual assumption of economic and technical feasibility under future meteorological conditions. The past increases in snow reliability due to snowmaking identified for some ski resorts in our study are also consistent with the faith in snowmaking technology emphasized by some ski industry stakeholder risk perception studies. Trawöger (2014) mentioned a "technosalvation" belief among Austrian CEOs of cable car companies, and in a similar way Bicknell and McManus (2006) showed an "overwhelming cornucopian belief" regarding snowmaking technology among Australian CEOs of ski lift companies. Steiger et al. (2019) stated that this belief is a widely held consensus in the ski industry whereby "with constantly improving snowmaking techniques and continued investment, [ski lift operators] are well prepared for climate challenges that lie ahead". According to Abegg et al. (2017) the faith placed in snowmaking by the ski tourism industry is one of the reasons why there remains a perception gap between the scientific literature and some ski tourism industry stakeholders regarding upcoming climate change challenges. The partial success of snowmaking in counteracting challenging natural snow conditions may skew ski lift operators' risk perception and lead them to consider past performance to plan future returns on snowmaking in reducing the impacts due to climate change. Yet, future success of snowmaking as an adaptive measure is not a guarantee in increasingly challenging climate conditions, and ski lift operators should not assume that snowmaking investments will continue to do well in the future simply because they have performed satisfyingly over the past decades. This is especially the case considering that the snow production demand will increase in a warming climate (Hock et al., 2019).

Conclusions and perspectives
Our contribution introduces several innovative elements in the assessment of past changes in natural and managed snow reliability in the ski industry. Our approach takes into account the changing snowmaking coverage at the resort level. In order to achieve this, we inferred the individual snowmaking dynamic of each ski resort based on its respective snowmaking investments. We then quantified snow reliability changes at the monthly timescale and focused on challenging winter seasons. We showed that the frequency of challenging conditions for ski resort operation over the 1991-2019 period increased in November and February to April compared to the reference period . In general, snowmaking had a positive impact on snow reliability, especially in December to January. If for the highest ski resorts, snowmaking improved snow reliability for the entire winter season, it did not counterbalance the decreasing trend in snow cover reliability for lower-elevation ski resorts. Contrasting outcomes of snowmaking highlight the various degrees of vulnerability across ski resorts due to their geographical settings and their business models, which affected their snowmaking deployment trajectory. Such a high heterogeneity in the ski industry advocates for climate change risk assessments to be carried out, whenever possible, at the scale of individual ski resorts and calls for caution in developing generalizing discourse about climate change impacts on the ski tourism industry at the regional level. Despite limitations on both snow cover modelling and the modelling of the snowmaking dynamics, our study is an additional element in providing more precise assessments about past changes in ski resort operating conditions. Further progress in taking into account operating condition pathways (for a larger number of ski resorts or with refinements in accounting for changes in technology), including snowmaking equipment, will depend on the availability of reliable resort-level information in various ski tourism markets. Nevertheless, our results pave the way for further studies addressing the detection and attribution of past changes and impacts of climate change on socio-ecological systems, which require analysing not only changes in climatic impact drivers but also actual impacts, taking into account all the components of climate change risk (climatic impact drivers, exposure and vulnerability). Beyond these methodological perspectives, our work opens the way to broader investigations into the consequences of the development of snowmaking facilities in mountain regions. Steiger et al. (2019) mentioned the challenge of better understanding ski tourism path dependency. We speculate that initial gains provided by snowmaking can foster the pursuit of these investments and set ski tourism stakeholders on a course to path dependency. Beyond a general trend of snowmaking investments in Savoie, we pointed out that individual situations greatly differ in terms of snowmaking equipment dynamics and snow reliability outcomes, depending on ski resort particularities. While this study explicitly considers geographical characteristics, it still underestimates the potential influence of business models, which plays a key role in climate change vulnerability, with implications for the climate change adaptation strategy at the individual ski resort level, and warrants further investigations. Code availability. The code used for the aggregation of snow cover simulations at the scale of individual ski resorts  is available at this address: https://doi.org/10.15454/COS00O. Data availability. The S2M atmospheric and natural snow cover simulations dataset is available at this DOI: https://doi.org/10. 25326/37 (Vernay et al., 2019).
Snow cover simulations accounting for grooming and snowmaking are not available publicly because of the risk of misuse and due to the high sensitivity of the dataset -note that individual ski resorts are anonymized in this study for the same reason. The corresponding datasets can be made available by the authors upon reasonable request.