Trends in seasonal snowpack and their relation to climate variables in mountain catchments in Czechia

ABSTRACT This study investigated trends in snowpack for the period 1965–2014 in 40 catchments located in five mountain regions in Czechia. We analysed daily series of air temperature, precipitation, and snow water equivalent (SWE) that were simulated with a conceptual model. The Mann-Kendall test showed strong increasing trends in air temperature at all elevations, mostly at the end of the cold season. This increase caused a decrease in snowfall fraction and SWE. Maximum SWE decreased mainly in western parts of Czechia (by up to −45 mm/decade). The length of the snow-covered period decreased by up to −6.8 days/decade, mainly due to earlier melt-out. Snowpack was more sensitive to changes in air temperature at elevations below 900 m a.s.l., while precipitation had a larger impact on snowpack at elevations above 1200 m a.s.l. The relative importance of air temperature for snow variability increased at all elevations in the last few decades.


Introduction
Seasonal snowpack is an important component of the hydrological cycle. During the cold season, a large amount of water is temporarily stored in snowpack and thus cannot contribute to runoff. During spring snowmelt, accumulated water is released and streamflow increases. Snowmelt water is a key source for groundwater recharge (Tague and Grant 2009) and groundwater storage affects streamflow during summer. Therefore, snowpack variability has an impact on catchment runoff during summer (Jenicek and Ledvinka 2020a). Recent decreases in snow water equivalent (SWE) and the shift in the date of the snow cover melt-out results in lower spring and summer streamflow (Godsey et al. 2014, Jenicek et al. 2016, Langhammer and Bernsteinová 2020. The amount of accumulated snow is primarily affected by air temperature and precipitation amounts during the cold season. Recent studies have found that snowpack has been decreasing worldwide due to the increase in air temperature (Lemke et al. 2007). The air temperature increase causes a decrease in the snowfall fraction (S f ), a rate of snowfall water equivalent to the total precipitation (Barnett et al. 2005, Knowles et al. 2006. In the last few decades, the decrease in S f ranged from 0.8%/decade in high-elevation catchments in Switzerland to 5.5%/decade for catchments in Czechia and Slovakia (Blahušiaková et al. 2020). According to (Hynčica and Huth 2019), for the region of Czechia, the phase of precipitation in the cold season shifted from solid to mixed precipitation over the last 35 years, resulting in a decrease in S f by on average 6%/decade. The decrease in S f was the strongest in February and January (Hynčica and Huth 2019).
The increase in air temperature also results in a decreasing number of days with snow cover (Marty 2008, Brown andMote 2009). The snow season has shortened over the last few decades, mainly due to a shift in the date of melt-out to an earlier calendar date (Klein et al. 2016, Zhang andMa 2018). In catchments in the Swiss Alps, the snow season decreased on average by 8.9 days/decade in the period 1970-2015. The major cause of this decrease was an earlier date of melt-out, which shifted by 5.8 days/decade (Klein et al. 2016). The shortening of the snow-covered period was reported also from Czech mountain catchments over the last few decades (Langhammer et al. 2015, Blahušiaková et al. 2020. Shortening of the snow season in Czechia was reported especially for climate stations located at low elevations (Kliment et al. 2011). The shift of the snowmelt season to an earlier date causes less energy available for snowmelt which may result in lower snowmelt rates (Musselman et al. 2017, Wu et al. 2018. The sensitivity of snow cover to changes in air temperature depends on elevation. The most sensitive are areas with mean air temperature around 0°C in the cold season (Adam et al. 2009, Brown andMote 2009). In areas where mean seasonal air temperature is well below 0°C, changes in snow cover are controlled mainly by changes in winter precipitation. In these areas, SWE increases despite the increase in air temperature (Mote et al. 2005, López-Moreno et al. 2009. A linear relationship between air temperature, precipitation and snow cover was described by Morán-Tejeda et al. (2013) for the region of the Swiss Alps. These authors proposed 1400 m a.s. l. as the threshold elevation below which air temperature is a better predictor of SWE. Above this elevation threshold, the snow cover variability is controlled mainly by precipitation amount. The location of the elevation threshold changes for different regions. For example, an elevation of 850 m a.s.l. was proposed for the region of Norway (Skaugen et al. 2012). However, such threshold elevations are often increasing with recent increases in air temperature (Morán-Tejeda et al. 2013).
A decrease in SWE during recent decades was reported for the United States (Mote et al. 2018, Zeng et al. 2018, Asian high mountains (Smith and Bookhagen 2018), and German low mountain ranges (Dong and Menzel 2020). Zhang and Ma (2018) described a decrease in mean and maximum SWE over large parts of Eurasia during the period 1979-2004. Recent changes in SWE in the Alps were investigated by Marty et al. (2017). These authors found a decreasing trend in SWE, which was independent of latitude or longitude. The decrease in SWE was stronger in spring than in winter, and the loss in SWE was greater at lower elevations (up to 80% for 60 years) than at higher elevations (about 10%). Negative trends in snow depth in mountain catchments in central Europe were reported by Blahušiaková et al. (2020). According to this study, snow depth decreased by 1% at higher elevations (~2300 m a.s.l.) and by 63% at lower elevations (~800 m a.s.l.). A decrease in snow depth in Czechia was also described by Kliment et al. (2011). With an ongoing increase in air temperature, it is expected that the SWE will continue to decrease, and the shortening of the snow season will continue in the future (Bavay et al. 2013, Fyfe et al. 2017, Jenicek et al. 2018. In past decades, snow cover has decreased worldwide mainly due to the increase in air temperature. Although these changes have a large impact on the runoff regime, snowpack trends have rarely been systematically investigated in central Europe outside the Alpine region. This is especially important since continued warming is predicted for central Europe in the 21st century while future changes in precipitation remains uncertain (Svoboda et al. 2016). Therefore, it is essential to answer the question of how snow cover responds to the mutual interaction of changes in air temperature and precipitation in different regions and at different elevations.
The above points guided the two main objectives of our study, specifically (1) to analyse long-term changes in selected climate and snow cover characteristics in mountain catchments in Czechia and (2) to investigate the snowpack sensitivity to climate characteristics at different elevations. We investigated the changes in selected snow indices at a multicatchment level for 40 catchments located at different elevations. The long-term trends were analysed using simulated time series resulting from a bucket-type catchment model.

Data
The study area consists of 40 catchments located in five mountain ranges in Czechia ( Fig. 1 and Table 1). Nine catchments are located in the Bohemian Forest (BF), two are in the Ore Mountains (OM), 12 are in the Western Sudetes (WS), 12 are in the Eastern Sudetes (ES) and five are in the Western Carpathians (WC). These catchments were selected for their hydrological regime, which is affected by seasonal snow cover. The runoff regime in the study catchments is not affected by any major human activity (dams, irrigation etc.). Each study catchment covers an area from 1.8 to 383 km 2 , and the mean catchment elevation ranges from 493 to 1297 m a.s.l. The mean annual maximum SWE for individual catchments ranges from 40 mm for the lowest catchments to 689 mm for the highest catchments (Table 1).
Time series of mean daily air temperature and daily sum of precipitation (period 1965-2014), mean daily discharge and weekly SWE (period 1980SWE (period -2014 were available from climatological and hydrological stations operated by the Czech Hydrometeorological Institute. Station data were further used to calibrate and validate a hydrological model and to simulate the individual components of the water cycle.

Model calibration
The conceptual semi-distributed hydrological model HBVlight (Seibert 2000) was used to simulate individual components of the rainfall-runoff process. HBV-light is a bucket type model, and it consists of four main routines: (1) the snow routine which simulates snow accumulation and snowmelt based on the degree-day method; (2) the soil routine, which simulates soil moisture, actual evapotranspiration and groundwater recharge; (3) the response routine, which simulates runoff from two groundwater boxes; and (4) the routing routine, which simulates the runoff propagation to the catchment outlet (Seibert and Vis 2012). Each catchment was divided into sub-zones according to elevation at intervals of 100 m. Precipitation, air temperature, SWE, amount of snowmelt water, and soil moisture were simulated separately for these elevation zones.
The HBV model was calibrated for each catchment using a genetic algorithm procedure (Seibert 2000). The model was calibrated against both observed runoff and SWE. Since the observed SWE was available only for the period 1980-2014, the model was calibrated for the period 1980-1997 and validated for the period 1998-2014. Results of the calibration were evaluated using a combination of three objective functions: (1) logarithmic Nash-Sutcliffe efficiency for runoff (Nash and Sutcliffe 1970; 60% weight), (2) volume error for runoff (20% weight), and (3) Nash-Sutcliffe efficiency for SWE (20% weight). This multicriteria calibration of the model against both runoff and SWE enabled us to better control the simulation of the SWE in individual catchments. One hundred calibration runs were performed, to reduce the parameter uncertainty resulting from the fact that the same simulation results might be reached by several combination of model parameter values (equifinality). The results of the model calibration are described in Jenicek and Ledvinka (2020a), and resulting simulations are publicly available (Jenicek and Ledvinka 2020b). Objective functions for all catchments are shown in the Appendix (Fig. A1). A comparison between observed and simulated SWE for selected catchments is shown in Fig. A2.
Based on 100 optimized parameter sets, we created 100 simulations for the entire period 1965-2014, based on observed data of precipitation and air temperature. A median simulation resulting from these 100 parameter sets was used for further analysis. Figure A3 in the Appendix shows the variability of trends based on the 100 parameter sets and representativeness of the median values used in analyses.

Snow and climate characteristics and data analysis
Several snow and climate characteristics were selected for the analysis (Table 2). To investigate regional variability, we used mean catchment values as well as values representing individual elevation zones defined by the model (100 m step). Changes in snow and climate characteristics were analysed at seasonal (October-May) and monthly levels.
The seasonal mean air temperature (T cold ) and precipitation amounts (P cold ) were analysed either in terms of the catchment mean or separately for individual elevation zones. Snowfall fraction (S f ), a ratio of snowfall water equivalent to total precipitation, was calculated using a single threshold temperature approach. The threshold temperature was calibrated by the HBV-light model separately for each catchment. Monthly mean values of these characteristics were also investigated.
To describe the overall snow storage, maximum seasonal SWE (SWE max ) and mean SWE were calculated. Mean SWE was calculated as the mean of daily SWE values for the specific month or for the whole season (SWE mean ). SWE max was calculated as the absolute daily maximum of SWE for the whole cold season. Sums of daily melt water from snowpack were calculated for individual months and for the whole cold season (S melt ).
In addition to the changes in amount of snow, we also analysed changes in length of snow cover duration. First, we defined the day of the year of snow cover onset (DOY onset ) as the first day in a season after which SWE is equal to or higher than 5 mm for at least 7 days. Second, we defined the day of the year when snow cover melts out (DOY melt ) as the first day after the day with SWE max when SWE drops below 5 mm. The snow cover duration (S dur ) was then calculated as the number of days between DOY onset and DOY melt . Since S dur does not account for temporary snow melt-out during the season, we additionally defined the number of days with snow (N dws ), as a sum of days in which SWE exceeded 5 mm. Since the values of the above characteristics may be affected by thresholds used for their calculation, we tested a variety of thresholds and found no major impact on results or their interpretation.
The Mann-Kendall test (Mann 1945, Kendall 1975) was used to determine the statistical significance of potential trends in time series of individual snow and climate characteristics.
The Theil-Sen's slope estimator (Sen 1968) was used to estimate the slope of the monotonic trend. Analyses were done using the Python package pyMannKendall version 1.4.1 (Hussain and Mahmud 2019). Since the Shapiro-Wilk test showed that not all study characteristics are normally distributed, the Spearman's correlation coefficient r s was calculated to analyse the strength of relationships between selected variables. Three p values were used as thresholds of statistical significance: .01, .05 and .1. All analyses were done separately for individual catchments and elevation zones.

Snow cover changes between two time periods
The study period 1965-2014 was divided into two sub-periods, 1965-1989 and 1990-2014, to detect potential changes in mean values of the snow characteristics in the study catchments between the two periods. In the period 1990-2014, the mean daily SWE in all catchments was lower compared to the period 1965-1989 ( Fig. 2(c)). An exception was a short period at the beginning of the snow-covered season when the differences were negligible. In general, mean annual SWE max was lower by 28 mm (21%) in the second sub-period  compared to the first sub-period . The mean date of SWE max occurred 9 days earlier in the second sub-period and the mean DOY melt occurred 12 days earlier. The mean DOY onset did not change considerably between the two sub-periods (Fig. 2(c)). Since the amount of snow decreased considerably in the second sub-period, maximum values of S melt also decreased, and the period with the highest S melt values occurred earlier as a response to higher air temperatures and thus earlier snowmelt period ( Fig. 2(a) and (d)). The mean daily precipitation amount was similar in the two sub-periods ( Fig. 2(b)).
The changes in climate and snow characteristics between the two sub-periods showed different patterns for individual elevation zones (Fig. 3). The difference in cold season air temperature (T cold ) between the two sub-periods was around 0.6°C for all elevation zones. Precipitation was moderately higher (3-6%) in the second sub-period at all elevations except the highest elevation zone, which showed a small decrease in cold season precipitation (P cold ), by 1% ( Fig. 3(b)). The absolute decrease in S f was consistent across all elevations, since it followed the decrease in T cold which did not differ across individual elevation zones ( Fig. 3(a) and (c)). The relative decrease in SWE max for the second sub-period compared to the first sub-period was the greatest at lower elevations, where air temperature was close to the freezing point in the first sub-period but increased for the second sub-period ( Fig. 3(a) and (d)). In the first subperiod, SWE max was higher than in the second sub-period by 18-26%, at elevations between 700 and 1000 m a.s.l. The difference decreased with elevation and it was always negative, i.e. the decrease in SWE max in the second sub-period was smaller for higher elevations. The variability in SWE max was higher at higher elevations for both sub-periods, due to overall higher snow accumulations. The variability in SWE max at all elevations increased in the second sub-period as a result of a decrease in lower quartile values, while the upper quartile was more or less the same in both sub-periods.   In the second sub-period, 1990-2014, S dur was shorter in most of elevation zones (Fig. 4). The decrease in S dur was mainly caused by the shift in DOY melt to an earlier calendar date. The largest decrease in S dur was found at elevations below 1000 m a.s.l. (10-24 days), caused mainly by earlier DOY melt . Interestingly, changes in DOY onset were small or even negligible at these low elevations. In contrast, the DOY onset occurred earlier at elevations between 1100 and 1400 m a.s.l. in the second sub-period compared to the first sub-period. In general, changes in DOY onset were smaller compared to changes in DOY melt in most elevation zones, except elevations from 1100 to 1300 m a.s.l., where both DOY onset and DOY melt occurred earlier, and thus S dur did not change.

Seasonal trends in snow and climate characteristics
The most significant trends were found for T cold , S dur , DOY melt , and S f . Increasing trends in T cold were found at all elevations (p values mostly lower than .05; Fig. 5). The change in T cold reached about +0.2°C/decade at elevation zones between 600 and 1300 m a.s.l. Above 1300 m a.s.l., the temperature increase was slightly smaller (+0.15°C/decade). Regionally, increasing trends were found in most catchments (p < .01), except for some catchments located in the Western Carpathians and most of the catchments in the Eastern Sudetes (Fig. 6).
No significant trends in P cold were found across individual elevation zones (Fig. 5). Nevertheless, significant trends were found for 15 catchments located mostly in western parts of the Eastern Sudetes (up to +45 mm/decade) and in the Ore Mountains (up to +33 mm/decade). In contrast, decreasing trends were found in catchments in western parts of the Western Sudetes (Fig. 6).
At elevations from 700 to 1100 m a.s.l., S f decreased significantly. The strongest trends were found at elevations from 1100 to 1300 m a.s.l. (from −0.015 to −0.018/decade; Fig. 6). Regionally, significant decreases in S f occurred in catchments in the Western Sudetes, Bohemian Forest, and Ore Mountains (Fig. 6).
At elevations from 700 to 1000 m a.s.l., SWE mean decreased significantly. SWE mean also significantly decreased above 1400 m a.s.l., where the decrease was the most apparent and reached −14 mm/decade. In contrast, at elevations from 1000 to 1400 m a.s.l.,the decrease in SWE mean was not significant. Decreasing trends were found in 16 catchments in the Western Sudetes and Bohemian Forest, and sporadically in the Eastern Sudetes. The strongest trends were found in the Western Sudetes (up to −27 mm/decade). No significant increasing trends in SWE mean were found in the western part of the Eastern Sudetes, which showed an increasing trend in P cold (around +40 mm/decade, p < .01), and which could partly compensate for an increase in T cold . No significant trends were found in the Western Carpathians (the most eastern part of Czechia).
Significant decreasing trends in SWE max were found mostly in the same catchments that showed trends in SWE mean , except for catchments located in the Bohemian Forest, where only one significant trend was identified (Fig. 6). Increasing trends in SWE max were found in two catchments with a strong increase in P cold and no significant trend in T cold . The decrease in SWE mean and SWE max was followed by a decrease in S melt .
In all elevation zones, except the lowest one, DOY melt shifted significantly to an earlier date. The most evident shift occurred at elevations around 800 m a.s.l.(−5.6 days/decade). The trends in DOY melt were significant at elevations above 800 m a.sl., with p values lower than .01, giving the highest confidence of derived trends among all the study snow characteristics. Decreasing trends in DOY melt were detected in most of the catchments located in the Western Sudetes (up to −9 days/decade) and in about half of the catchments in the Bohemian Forest, Ore Mountains, and Western Carpathians.
At elevations from 700 to 1400 m a.s.l., S dur decreased significantly (up to −6.8 days/decade). The decrease in S dur was strongest below 1000 m a.s.l. A similar pattern was detected in DOY melt . At the catchment level, trends in S dur corresponded to trends in DOY melt . Catchments without significant trends in S dur are located mainly in the Eastern Sudetes. No significant trends in DOY onset were found at either catchment or elevation level.  1965-1989 (red) and 1990-2014 (blue). Grey, green and orange bars show differences between medians of S dur , DOY melt and DOY onset , respectively.

Trends in snow and climate characteristics for individual months in the cold season
The strongest increasing trends in mean air temperature were detected in April (from +0.47 to +0.58°C/decade), followed by May and November (about +0.3°C/decade). Trends were significant at the .01 level at all elevations in April, and they were stronger at lower elevations ( Fig. 7(a)). In contrast, the analyses found almost no trends in monthly precipitation. Weak trends in precipitation were found only for April in a few elevation zones. The strongest (and only one with p < .05) decrease is −7.5 mm/decade ( Fig. 7(b)).
The strongest significant trends in monthly S f were found mainly in January, March, and April. The strongest decrease was detected at elevations above 800 m a.s.l. In April, S f decreased at all elevations (Fig. 7(c)). The most pronounced decreases were found at the highest elevations in April (up to 0.065/decade). The largest decrease in mean SWE occurred at the end of the cold season and was caused by air temperature increase together with situations when the temperature was close to 0°C. Such combinations had a large impact on snow storage. It resulted in a significant decrease in mean SWE in April at lower elevations and in May at higher elevations. The largest decrease (−29 mm/decade) was found in April at elevations above 1400 m a.s.l. (Fig. 7(d)).
The main period with the main spring snowmelt was shifted from May to April (Fig. 7(e)). A decreasing trend in the snowmelt rate was found in the majority of elevation zones for May with values up to −47 mm/decade for elevations above 1400 m a.s.l. In contrast, an increasing trend (p value < .01) occurred in the two highest elevation zones in April (+31 mm/ decade). Strong decreasing trends in N dws were found at the end of the cold season ( Fig. 7(f)). This is also the case for mean SWE, with N dws decreasing in April at elevations below 1000 m a.s.l., and in May at elevations above 1000 m a.s.l. In both months, the decrease in N dws was about 2-3 days/decade. Significant trends were also found for December, but the slope of this decrease was rather small.

Correlation between climate and snow characteristics
To analyse mutual interactions between snow and climate characteristics at different elevations, the Spearman's correlation coefficient r s was calculated for all possible pairs of the study characteristics. The individual correlations were calculated separately for the whole period and the two sub-periods for three elevation zones representing low, middle and high elevations: 600-900, 900-1200, and 1200-1500 m a.s.l. (Fig. 8).
In general, the lower T cold is for a specific year, the higher is S f (Fig. 8(a)). Due to the increase in T cold , the correlation between T cold and S f was considerably stronger for low elevations in the second sub-period than in the first sub-period (−0.64 and −0.45, respectively). While no significant correlation between T cold and S f was found at mid and high elevations in the first sub-period, there was strong correlation in mid elevations in the second sub-period (r s = −0.49; p < .05). A weaker correlation (r s = −0.34; p < .1) was found for high elevations in the second sub-period.
Correlations between SWE mean , SWE max , and T cold have similar patterns to S f , although the correlation coefficients are lower. With the increasing air temperature, the relationship of S f to SWE mean and SWE max became stronger. In the second sub-period, the values of the correlation coefficients between S f and related snow characteristics rapidly increased. For example, the correlation coefficients between S f and SWE max at mid elevations are twice as high in the second sub-period compared to the first sub-period. Low or no impact of T cold on S f , SWE mean and SWE max at higher elevations was probably caused by lower air temperature, which was well below the freezing point, even in warmer episodes. Therefore, air temperature had a smaller effect on precipitation phase and consequently on snow storage.
The variability in P cold has a larger impact on SWE mean at higher elevations than does the variability in T cold . The positive correlation between P cold and SWE mean was found at all elevations for the whole study period . The impact of P cold on SWE mean was weakest at low elevations (r s = 0.26), and it increased with elevation to 0.68 ( Fig. 8(a)). The correlation between SWE mean and P cold was stronger for the first subperiod compared to the second sub-period (Fig. 8).
At higher elevations, the relationship between S dur and DOY onset was strong, with a Spearman's correlation coefficient of −0.91. This correlation decreased with decreasing elevation. In contrast, DOY melt was strongly positively correlated with S dur at lower elevations (r s = 0.79), and values of correlation coefficients were lower at mid and high elevations, where the values were almost the same (0.66 and 0.69, respectively). This indicates that S dur at lower elevations was more affected by changes in DOY melt than DOY onset for the whole study period. The correlation between S dur and DOY onset was similar for both sub-periods, while the correlation between S dur and DOY melt-out was stronger in the second sub-period. The latter indicates an acceleration of the processes, causing an earlier onset of the snowmelt period and thus earlier snow melt-out.
At elevations from 600 to 900 m a.s.l., values of SWE max were mainly affected by T cold (Fig. 9(a)). SWE max was higher in years with lower T cold even for years with low P cold , as indicated by black circles located mostly in both left quadrants of Fig. 9(a). In contrast, at elevations above 1200 m a.s.l., P cold was an essential factor that affected SWE max . At these higher elevations, higher SWE max occurred more often in years with high P cold regardless of T cold (black circles in both top quadrants; Fig. 9(c)). The above pattern is found in both sub-periods. In general, SWE max was lower and the T cold was higher in the second sub-period than in the first sub-period ( Fig. 9  (a-c)).

Long-term changes in climate variables, snowfall fraction, and snow storage
Our results show that air temperature increased the most in April, which corresponds to the end of the snow-covered period. During this time, air temperature is near the freezing point and thus the snowfall fraction is sensitive to changes in air temperature. The changes in air temperature and phase of precipitation led to a decrease in SWE and shift of snow cover melt-out to an earlier date, which is the main reason for the shortening of the snow-covered period.
The decrease in S f and consequent changes in snow storage shown in our study are consistent with previous studies from regions with similar climate, although with a larger elevation range (Knowles et al. 2006, Jenicek et al. 2016). Significant trends in S f were detected at mid elevations between 700 and 1300 m a.s.l. Comparing different regions, the strongest decreasing trends in S f were found for catchments in the Bohemian Forest (up to −3.5%/decade). The magnitude of the slopes of the trends agrees with results by Blahušiaková et al. (2020), who focused on a similar study period and area in central Europe. S f decreased mostly in April at higher elevations (up to −6.3%/decade) and in March at middle elevations (up to −6.3%/decade), which corresponds to the findings of Blahušiaková et al. (2020). In contrast, Hynčica and Huth (2019) reported that S f decreased the most in February in Czech climate stations (−10%/decade). A possible explanation for the differences is that rather high elevation areas were included in our study, and thus the changes in S f were not detected for February due to overall low air temperatures at these high elevations. Additionally, Hynčica and Huth (2019) investigated the period 1983-2018, which is known to be a period with overall acceleration of trends in air temperature and S f compared to earlier periods (Langhammer and Bernsteinová 2020).
The results showed that SWE was lower at all elevations in the second sub-period. The difference was largest at low elevations and it decreased with increasing elevation (Fig. 3). Significant long-term trends were found mainly for lower elevations. These findings corroborate studies from Switzerland (Marty et al. 2017), Norway (Skaugen et al. 2012), and the United States (Mote et al. 2018). In contrast with the results of Skaugen et al. (2012) or Mote et al. (2005), we did not find significant increasing trends in SWE mean even in catchments with significant increasing trends in P cold . Increasing trends in SWE max were found only for two catchments in the Eastern Sudetes, which were probably caused by the increase in precipitation. Otherwise, trends in SWE max were similar to trends in SWE mean .
At low and middle elevations, mean SWE decreased mostly in April, and at high elevations in May, which was caused by the fact that increasing trends in monthly air temperature were the strongest in these spring months. The absolute decrease in mean SWE was strongest at higher elevations. A larger decrease in SWE in April compared to February was reported by Marty et al. (2017) from regions of the Alps with investigated trends between 500 and 3000 m a.s.l. This study also showed that the highest absolute SWE decrease is related to higher elevations. In the elevation zones below 800 m a.s.l., significant decreasing trends in mean SWE were found from January to March, as also reported from German mountain areas below 820 m a.s.l. (Dong and Menzel 2020).
Since the precipitation data were adjusted by the HBV model, they may differ from measured data used for model calibration. Specifically, a snowfall correction factor (S FCF ) was used in the model to correct solid precipitation for the undercatch (e.g. due to wind). Similarly, precipitation and air temperature were adjusted for elevation using lapse rate, a model parameter calibrated separately for each catchment. Consequently, the trends detected with these adjusted data may differ from trends in the observed station data. Therefore, we compared trends in observed and simulated seasonal precipitation and air temperature and did not find any differences; the trends are nearly identical (Appendix and Fig. A4).
Similarly, we compared seasonal and monthly trends in observed and simulated SWE despite the fact that only weekly observed SWE data were available from climate stations, making any evaluation rather complicated. Nevertheless, the direction of the trends for both observed and simulated SWE was the same in most of the catchments, although the magnitude and significance of the trends partly differed (Fig. A4). However, we found no statistically significant trend in both observed and simulated SWE with opposite slope directions. The differences in trend magnitude mentioned above were probably caused by the comparison of trends based on weekly observed SWE at climate stations on the one hand, and trends based on simulated SWE from 100 m wide elevation zones on the other hand. Additionally, the compared time series were relatively short for making a convincing trend analysis .

Long-term trends in the length of the snow-covered period
The results showed that S dur decreased significantly in most elevation zones after 1965 due to the shift in DOY melt to an earlier calendar date. Trends in both variables were stronger at lower elevations. The shortening of the period with snow cover and the shift of melt-out day was also reported by Klein et al. (2016). These authors highlighted the shift in melt-out day as the main factor causing the shortening of the snow-covered period. In contrast to our results, Klein et al. (2016) found significant trends in the shift in snow onset day to a later calendar date, which was not found in our study catchments. Therefore, an earlier start of the snowmelt period and thus earlier melt-out is the main cause of the shortening of the snow-covered period in our study region.
A shortening of the snow-covered period caused mainly by an earlier melt-out was also described by Zhang and Ma (2018) at a sub-continental scale in Eurasia and by Zeng et al. (2018) for the western United States. However, the latter study additionally reported later snow cover onset which also contributed to the shortening of the snow-covered period. The snow-covered period has been shrinking, with average rates of 8.9 days/decade in the Swiss Alps and 10 days/decade in the western United States, over the last few decades (Klein et al. 2016, Zeng et al. 2018. In our study catchments, the average rate of shortening of the snowcovered season reached 4.8 days/decade with a maximum of 10 days/decade. Our results show that the strongest decreases in S dur are mostly situated in the western part of Czechia (Bohemian Forest, Western Sudetes), while only a few catchments showed a decrease in S dur in the eastern regions of Czechia. This might be caused by an increasing continentality from west to east and by different synoptic situations in the west and east causing snowfall and partial snowmelt during the winter season (Steirou et al. 2017).
Strong decreasing trends in N dws were detected at the end of the cold season in our study catchments. Similar to SWE, N dws decreased in April at lower elevations and in May at higher elevations. The decrease in N dws in spring was strongly connected to the increase in air temperature. Langhammer et al. (2015) identified a step-like drop of N dws in the 1980s. After this drop, N dws slightly increased. A step-like drop in N dws in the 1980s due to an increase in air temperature was also reported in Switzerland (Marty 2008).
In this study we calculated only simple monotonic trends for the whole period 1965-2014, in a similar way to other studies, such as Hynčica and Huth (2019) and Blahušiaková et al. (2020). Therefore, we might have overlooked step-like changes in time series described in the previous paragraph, and it might thus be beneficial to look for step-like changes in the time series in future studies.

Mutual influence of air temperature and precipitation on snow storage
Snow variables were negatively correlated with air temperature and positively correlated with precipitation. The strength of the relationship between snow characteristics and temperature (precipitation) decreased (increased) with elevation. This indicates that snow storage at higher elevations is controlled mostly by precipitation amounts, while air temperature is the main controlling factor of snow storage at low elevations. The results of our study are in good agreement with findings from the Alps (Morán-Tejeda et al. 2013, Marty et al. 2017) and the western United States (Mote et al. 2005, Sospedra-Alfonso et al. 2015. In our study, the relationship between SWE and precipitation was stronger than the relationship between SWE and temperature at elevations above 1200 m a.s.l. From 900 to 1200 m a.s.l., the impact of air temperature and precipitation on SWE seems to be of equal importance. Below 900 m a.s.l., the SWE was controlled mainly by air temperature. The threshold below (above) which the snow cover variability is controlled mainly by air temperature or precipitation was found to be around 1560 m a. s.l. for the central Rocky Mountains, and around 1400 m a.s.l. for the Swiss Alps (Morán-Tejeda et al. 2013, Sospedra-Alfonso et al. 2015, Marty et al. 2017. According to Mote (2003), spring SWE in the western United States decreased due to the increase in air temperature at elevations below approximately 1800 m a.s.l. in the second half of the 20th century. Similar findings, although with a lower elevation threshold of 850 m a.s.l., were reported in Norway (Skaugen et al. 2012).
The threshold elevation separating the different effects of climate variables to control snow storage may also increase with an ongoing increase in air temperature (Morán-Tejeda et al. 2013). This elevation increase was indicated also by our results, showing an increasing correlation of temperature and snow variables at higher elevations in the sub-period 1990-2014 compared to the sub-period 1965-1989 (which showed significant correlations only at low elevations). In addition, correlations increased at lower elevations and they became significant at lower levels (p < .01) in the second sub-period. However, the results from elevations above 1400 m a.s.l. may be affected by the uneven distribution of these highest elevations among individual mountain ranges. Elevations above 1400 m a.s.l. are present in only two of the five studied mountain ranges.

Conclusions
We analysed changes in selected snow cover characteristics in 40 catchments in five mountain ranges in Czechia for the period 1965-2014. We were mostly focused on the sensitivity of snow cover at different elevations to changes in air temperature and precipitation.
We compared climate and snow characteristics for two subperiods, 1965-1989 and 1990-2014. In the second sub-period, SWE decreased at all elevations. The largest decrease in SWE and the shortening of the snow-covered period were found at elevations around 800 m a.s.l., where mean seasonal air temperature was close to the freezing point. At these elevations, the length of the snow-covered period was shorter by up to 24 days in the second sub-period compared to the first subperiod and the SWE max decreased by up to 28 mm.
The long-term trend analysis showed an increase in seasonal air temperature at all elevations in our study catchments. The increase in air temperature was strongest at the end of the cold season (April, May). Consequently, the April and May snowfall fraction and mean SWE decreased. Additionally, the period with spring snowmelt shifted from May to April at higher elevations and from April to March at lower elevations. It resulted in an earlier day of year of melt-out, which was shifted by up to 5.6 days/decade at elevations from 700 to 1500 m a.s.l.
Earlier melt-out caused a shortening of the snow-covered period despite the fact that no significant trends were found for snow cover onset. The length of the snow-covered period significantly decreased at elevations between 700 and 1400 m a.s.l. The decrease ranged from 4.2 days/decade at low elevations to 6.8 days/decade at higher elevations.
In contrast to lower mean SWE in April and May, trends in mean December to March SWE were not detected at most elevations except elevations below 800 m a.s.l. The mean seasonal SWE and annual SWE max decreased significantly in the Western Sudetes, partly in the Bohemian Forest, and in a few catchments in the Eastern Sudetes (up to −27 mm/decade for SWE mean and −45 mm/decade for SWE max ). In the case of two catchments in the Eastern Sudetes, increasing trends in SWE max were detected owing to increased winter precipitation. In general, the results indicate that decreasing trends in winter snow storage and SWE max occurred more often in western parts of Czechia than in eastern parts of Czechia.
The results for our study catchments show that for changes in snow characteristics, trends in both seasonal air temperature and seasonal precipitation are important, although their importance is highly dependent on elevation. The impact of air temperature on snow variability was stronger at elevations below 900 m a.s.l., while seasonal precipitation was more important at elevations above 1200 m a.s.l. At elevations from 900 to 1200 m a.s.l., the impact of the two characteristics on SWE was of equal importance. With the increase in air temperature in the last few decades, the importance of air temperature in controlling snow storage and variability increased at all elevations. Figure A2. Comparison of observed and simulated snow water equivalent (SWE) for five selected catchments (one catchment per mountain range; see Table 1). The catchments were selected to represent a typical elevation of the respective mountain range. Simulated values represent the nearest elevation zone to stational observed values. Note that observed SWE might not be fully comparable to simulated SWE due to point and weekly SWE observations measured in open area at the nearest climate station. Figure A3. Variability of Theil-Sen's slopes for 100 simulations based on 100 optimized parameter sets for five selected catchments (one catchment per mountain range; see Table 1). Red dots represent values resulting from analyses based on median simulations. The catchments were selected to represent a typical elevation of the respective mountain range. Figure A4. Observed and simulated long-term trends in seasonal and monthly air temperature (T), precipitation (P), SWE and runoff (Q). The values represent Theil-Sen's slopes. Significant Mann-Kendall trends are highlighted in black (p < .1), in black bold (p < .05) or with asterisks (p < .01). Increasing trends in shades of red; decreasing trends in shades of blue. Catchments are sorted according to mountain ranges from west to east.