Edinburgh Research Explorer Snow cover duration trends observed at sites and predicted by multiple models

trends observed at sites and predicted by multiple models' The Cryosphere Discussions, Copernicus Publications. Abstract. Thirty-year simulations of seasonal snow cover in 22 physically based models driven with bias-corrected meteorological reanalyses are examined at four sites with long records of snow observations. Annual snow cover durations differ widely between models but interannual variations are strongly correlated because of the common driving data. No signiﬁcant trends are observed in starting dates for seasonal snow cover, but there are signiﬁcant trends towards snow cover ending earlier at two of the sites in observations and most of the models. A simpliﬁed model with just two parameters controlling solar radia- tion and sensible heat contributions to snowmelt spans the ranges of snow cover durations and trends. This model predicts that


Introduction
The extensive seasonal snow cover of Northern Hemisphere land is sensitive to climate warming and strongly influences snow cover by climate models has improved, CMIP5 simulations underestimated significant reductions observed in spring snow cover extent (Brutel-Vuilmet et al., 2013) and had a wide spread in predictions of snow-albedo feedback strength (Qu and Hall, 2014). In preparation for the sixth IPCC assessment report, climate modelling centres have now performed CMIP6 coupled land-atmosphere-ocean simulations with their latest models. Mudryk et al. (2020) report an overall better representa-20 tion of Northern Hemisphere snow cover extent in the CMIP6 multi-model ensemble than in CMIP5, but a large spread remains in simulated trends.
In addition to coupled model experiments, snow simulations by stand-alone land surface models have been driven with prescribed meteorological variables on global grids in the Global Soil Wetness Project (Dirmeyer et al., 2006) and at individual sites in the Project for Intercomparison of Land-surface Parameterization Schemes (Slater et al., 2001)

25
Intercomparison Project (Etchevers et al., 2004;Essery et al., 2009). These studies have invariably found wide ranges in simulations and inconsistencies in model performance. The Earth System Model-Snow Model Intercomparison Project (ESM-SnowMIP; Krinner et al., 2018) includes simulations driven with both in situ meteorological measurements and bias-corrected reanalyses at ten well-instrumented snow study sites; simulations were first evaluated with between seven and twenty years of in situ driving data, but using reanalyses allows longer simulations for investigating trends. This paper examines observed 30 trends in seasonal snow cover duration and simulations driven with 1980-2010 bias-corrected reanalyses at four of the ESM-SnowMIP sites selected because they had at least 27 years of daily snow observations up to 2010. The locations of the sites are given in Table 1. Reflecting motivations for the establishment of snow study sites by national organizations, Col de Porte (France), Reynolds Mountain East (USA) and Weissfluhjoch (Switzerland) are at high elevations in mid-latitude mountains, whereas Sodankylä (Finland) is a low elevation Arctic site. All of the sites typically have between five and eight months of 35 continuous winter snow cover and can have shorter periods of ephemeral snow cover at other times of year.
Simple empirical models of snowmelt are still often used for hydrological and glaciological applications, but all of the models participating in ESM-SnowMIP are physically based and calculate coupled mass and energy balances for snow on the ground. Eighteen groups submitted simulations by 22 models and model variants driven with a common set of bias-corrected SPONSOR, SWAP and Veg3D), and snow physics models (Crocus and SNOWPACK); references for all of these models can be found in Table 1 of Krinner et al. (2018). Although snow models are much less complex than comprehensive Earth System 45 Models, they have sufficient complexity and large enough parameter spaces to make it difficult to interpret why they behave in the ways that they do. For Earth System Models, Randall et al. (2019) concluded that "we must work to create much simpler models that can semiquantitatively reproduce the key results of the comprehensive models". In that spirit, a highly simplified two-parameter energy balance model ("2PM" hereafter) is used to interpret the results of the ESM-SnowMIP models.

50
All of the meteorological variables required to drive physically based mass and energy balance snow models (air temperature, humidity and pressure, snowfall and rainfall rates, shortwave and longwave radiation fluxes, and wind speed) for 1980-2010 at the ESM-SnowMIP sites were extracted from the GSWP3 dataset and interpolated from three-hourly to hourly timesteps.
Because coupling to an atmosphere model was not required, snow models that are not part of an Earth System Model were also able to participate in this component of ESM-SnowMIP. For GSWP3, the 20th Century Reanalysis was used to nudge the 55 dynamics of the Global Spectral Model for downscaling from 2 • to 0.5 • resolution. Biases in monthly means of temperature, diurnal temperature range, precipitation and radiation fluxes relative to Climate Research Unit Time-Series (CRUTS), Global Precipitation Climatology Centre and Surface Radiation Budget datasets were then removed. Additional bias corrections had to be applied for ESM-SnowMIP site simulations because the mountain sites are at much higher elevations than the 0.5 • GSWP3 grid cells in which they lie (Table 1). Biases relative to in situ measurements for overlapping periods at each site were simply 60 removed for all driving variables, thus preserving distribution shapes, seasonal cycles and trends from the GSWP3 dataset (Menard et al., 2019). The meteorological variables extracted from GSWP3, interpolated to hours and bias-corrected to the sites are referred to as the driving data for the ESM-SnowMIP models hereafter.
The simplified model that will be used for interpreting the ESM-SnowMIP results below has two fixed dimensionless parameters: a snow albedo α and a surface-atmosphere turbulent exchange coefficient C H . A large simplification comes from neglecting heat required to warm snow to the melting point in comparison with heat required to melt snow. Snowmelt rate M (kg m −2 s −1 ) is predicted by the energy balance equation with latent heat of melting λ m (0.334×10 6 J kg −1 ), latent heat of sublimation λ s (2.835×10 6 J kg −1 ), surface temperature and 75 for air pressure p (Pa), temperature T a , specific humidity q a , heat capacity c p (1005 J K −1 kg −1 ) and density ρ (kg m −3 ); U (m s −1 ) is wind speed and q sat is the specific humidity of saturated air. Equations (1) to predict changes in snow mass S (kg m −2 ), which is limited to be greater than or equal to zero and is converted to depth using a fixed snow density of 300 kg m −3 .
3 Results Figure 1 shows monthly means and trends in air temperatures measured at the sites and in the driving data; averages and ranges 85 of observed start and end dates for continuous seasonal snow cover are also shown. The seasonal temperature cycle and trends in the driving data match observations closely at Sodankylä because the station there was included in the CRUTS database used for correcting GSWP3 temperatures. Weissfluhjoch is 60 km from the nearest CRUTS station at Säntis but only 50 m higher.  Open circles show average dates of maximum snow depth. (e-h): Monthly temperature trends calculated by the Theil-Sen method (Sen, 1968). Vertical bars for measurements and grey bands for driving data show 95% confidence intervals.
rapid December warming at Sodankylä will not directly influence simulated snow cover durations because it corresponds with a reduction in the occurrence of very low temperatures at times when snow is not melting. Other warming trends at Reynolds Mountain East and Sodankylä occur during snow-free months, but warming trends at Col de Porte and Weissfluhjoch overlap 95 the normal periods of snowmelt.
Solid precipitation is notoriously difficult to measure accurately, and quality-controlled measurements of snowfall are not available for all years back to 1980 at all of the sites. Annual snowfall amounts derived from precipitation gauge measurements are therefore only shown for comparison with the driving data in Figure 2, and snowfall trends will only be investigated in the driving data. Weissfluhjoch is the only site with a significant downward trend in snowfall at the 95% confidence level, although Weissfluhjoch has the highest snowfall, coolest summer temperatures and longest seasonal snow cover.  Start and end dates for continuous seasonal snow cover were found by searching for the first and last dates with snow depths exceeding 2 cm before and after the dates of maximum snow depth in each year. Figure 3 shows averages and trends for start 110 and end dates in observations and simulations at all of the sites (annual time series from which these were calculated are shown in supplementary Figure A1). Simulated start dates are largely determined by snowfall in the driving data and show relatively little spread between the models, except that some models will melt early snowfall at Col de Porte and others will retain it on the ground. Trends towards later start dates are observed at all sites and in most model simulations, but none of these trends are found to be significant with 95% confidence. Simulated end dates are influenced by differences in how models 115 respond to increasing air temperatures and solar radiation in spring, leading to larger spreads between models. The spread is particularly large for Weissfluhjoch; two of the models melt snow consistently earlier than the others, and three models retain year-round snow cover in some years (which has never been observed in measurements going back to 1936 at Weissfluhjoch).
Years in which a model does not melt the snow are excluded from calculations of end dates. Significant trends towards earlier snow disappearance are observed at Col de Porte and Weissfluhjoch but not at Reynolds Mountain East or Sodankylä, and 120  Lejeune et al. (2019) and at Weissfluhjoch by Marty and Meister (2012). The remaining discussion here will focus on model behaviour at those two sites.
2PM was run 10,000 times for each site with snow albedos ranging from 0.5 to 1 and turbulent exchange coefficients ranging 125 from 10 −4 to 10 −2 ; these unrealistically wide parameter ranges give results that encompass and extend beyond the ESM-SnowMIP model results in Figure 3. Increasingly negative trends in continuous snow cover end dates for 2PM simulations that melt snow later in the year at Col de Porte and Weissfluhjoch are striking features of Figure 3. The same behaviour is seen most clearly for Weissfluhjoch in the ESM-SnowMIP models.
Snow albedo and turbulent exchanges between the surface and the atmosphere vary with time in reality and in realistic 130 models, but 2PM results can be plotted as contours or a colour scale on the fixed α − C H parameter space. Figure 4 overlays contours for snow cover end dates on colour maps of end date trends. Snowmelt becomes independent of air temperature as exchange coefficients approach zero and independent of solar radiation as albedos approach 1. Lower albedos and larger exchange coefficients lead to earlier melt at Col de Porte, as might be expected. At Weissfluhjoch, however, low albedos can cause radiation-driven melt in May, when solar radiation is high but air temperatures are still often below 0 • C ( Figure 1d); 135 larger exchange coefficients then delay melt by cooling the snow, so the May and June contours curve downwards in Figure 4b.
Even in the absence of net solar radiation and sensible heat (α = 1, C H = 0), there is sufficient longwave radiation in the driving data to melt the snow at Col de Porte each year, but the 2PM parameter space includes simulations that develop permanent snow cover at Weissfluhjoch (upper left corner of Figure 4b) if the previous winter's snow has not melted by mid-August.
Average observed and ESM-SnowMIP model snow cover end dates at Col de Porte fall in April or May; 2PM can produce 140 a wide range of end date trends for snow melting in those months, seen as a bulge in Figure 3a corresponding with a region where trend and end date contours cross in Figure 4a. ESM-SnowMIP models that have average end dates close to the start of May for Col de Porte have trends at the lower end of the 2PM range in Figure 3a, consistent with small exchange coefficients characteristic of low roughness and high atmospheric stability over snow.
Trends in snow cover end date show two areas of the 2PM parameter space in Figure 4 with enhanced negative trends.

145
Strong trends for snow melting in July have already been noted in Figure 3 and will be discussed again later. Enhanced trends also occur for snow cover ending in months with warming trends, provided that exchange coefficients are large enough for simulations to be sensitive to air temperature. This is apparent in Figure 4 as protrusions of stronger trends along the end date contours for late April at Col de Porte and mid-June at Weissfluhjoch. The average end date of continuous snow cover (dotted contour) is close to the date of maximum temperature trend (dashed contour) at Col de Porte, as expected for a positive 150 feedback on warming with decreasing snow cover duration. Snow disappears at Weissfluhjoch about two weeks later than the date of maximum warming, however. It may be that warming trends at Weissfluhjoch are dominated by advection from lower surrounding areas with earlier snowmelt; Col de Porte is at the fifty-seventh elevation percentile and Weissfluhjoch is at the ninety-fourth elevation percentile for 10 km × 10 km areas centred on the sites. Warming is also expected to vary with elevation in mountain regions (Pepin et al., 2015).

155
Annual snow cover duration (SCD) depends on the timing of snowfall, how much snow has to be melted and how much energy is available to melt it. Figure 5a for Weissfluhjoch and Table 2 for all sites show that modelled interannual variations in SCD are highly correlated with annual snowfall, except at Sodankylä; low snowfall and rapid temperature increases from April to May at Sodankylä limit variations in the end date of snow cover, both between years and between models ( Figure   3c). Beyond the range of the ESM-SnowMIP models in Figure 5a, correlations in 2PM simulations inevitably decrease as 160 the model undergoes a transition from seasonal to permanent snow cover at Weissfluhjoch independent of annual snowfall.
Incoming solar radiation in the driving data for Weissfluhjoch peaks around the summer solstice in late June, whereas energy available to melt snow from longwave radiation and sensible heat peak in late July (measured solar radiation actually peaks in May at Weissfluhjoch because of seasonal variations in cloud cover and multiple reflections between high albedo snow and clouds). Snow persisting after the peak in available energy will melt more slowly, so additional snowfall increases SCD more 165 for simulations that retain seasonal snow cover later. The sensitivity obtained by linear regression of SCD against snowfall, shown for Weissfluhjoch in Figure 5b, therefore increases for late lying snow. Because SCD is highly correlated with snowfall, increased sensitivity to snowfall in simulations with late lying snow and decreasing snowfall combine to amplify trends in SCD, as seen in Figures 3 and 4 for both Col de Porte and Weissfluhjoch.   (Mudryk et al., 2017). Large-scale simulations are required for predicting large-scale trends in snow cover extent, but simulations at well-instrumented sites allow more insight into modelling of snow processes and impacts that are experienced on small scales.
Interannual variations in modelled snow cover duration are strongly correlated with annual snowfall in the driving data at 180 three of the four sites, which means that the models are also strongly correlated with each other (supplementary Figure A2) because they all shared the same driving data. This inter-model correlation will not be preserved when snow models are coupled to different atmosphere models. Coupling also allows feedbacks that are suppressed when snow models are driven with prescribed meteorology. Coupled simulations with prescribed snow conditions are proposed in ESM-SnowMIP to evaluate the effects of snow feedbacks (Krinner et al., 2018). Because water will not be conserved if snow mass is prescribed indepen-185 dently of snowfall and melt, these should be land-atmosphere simulations with prescribed sea surface temperatures to avoid perturbations of the ocean by runoff that would occur in coupled land-atmosphere-ocean simulations.
A simple two-parameter snowmelt model shows that the response of snow models to warming depends on the timing of Author contributions. RE prepared the manuscript with substantial contributions from all co-authors. HK provided the global reanalysis data, which LW extracted and interpolated for the study sites. All other co-authors either performed model simulations or provided field data.