Analysis of the Variability and Future Evolution of Snowfall Trends in the Huaihe River Basin Under Climate Change

Global warming changes the characteristics of regional climate and crop growth environments and affects human production and life in the long term. In this study, historical snowfall observation data from 199 conventional meteorological stations in the Huaihe River Basin were analyzed with the climate trend method, the nonparametric Mann-Kendall (M-K) test and the sliding t test to determine the variability and abrupt changes in climate trends in the Huaihe River Basin during 1951–2018. Then, the Coupled Model Intercomparison Project phase 6 datasets were used to analyze the evolutionary trend of snowfall under four climate change scenarios in 2015–2065. The results show the following: 1) Due to climate warming, both the historical data and the simulated data under different climate change scenarios show declining snowfall in the Huaihe River Basin. 2) Based on the analysis of historical data from meteorological stations, the snow events had a clear latitudinal distribution—the higher the latitude, the lower the occurrence of snow events. In contrast, as shown by the analysis results obtained under different climate change scenarios, the snow mass flux had a clear longitudinal distribution—the higher the longitude, the higher the snow mass flux. 3) In Scenario SSP370, 2050 is important because the snowfall changes from decreasing to increasing, and in 2065, the snowfall in most areas of the basin is greater than that in 2050. Furthermore, the snowfall along the southeastern coast of the basin has the greatest variability.


INTRODUCTION
Global warming and the frequent occurrences of extreme meteorological disasters and heavy snow events in the Northern Hemisphere during winter have caused large-scale traffic disruption, damage to power facilities and severe crop damage (Khattak and Knapp, 2001;Mccabe et al., 2007;Gao, 2009;Feng and Chen, 2016). Additonally, climate change alters the characteristics of the regional climate and environment and directly affects the productivity and safety of humans (Austnes et al., 2008;Johansson et al., 2011). The Huaihe River, with a total length of more than 1,000 km, is one of the seven major rivers in China and flows through the three provinces of Henan, Anhui, and Jiangsu; the occurrence and evolution trend of snowfall in the Huaihe River Basin can provide an important reference for the prevention and mitigation of local disasters.
The study of snowfall study is very important for safety and property protection, and the prediction, and modeling of the evolution of snow events are hot research fields. Previous studies on snowfall have mainly focused on the following areas: 1) the use of remote sensing images to retrieve the seasonal variability in snow cover (Pu et al., 2007;Dietz et al., 2013) and to simultaneously assess the variability in snow cover in grasslands, farmland and forests; 2) the analysis of the weather conditions during snow events and the causes and effects of snow events by investigating snowfall over large areas (Wang et al., 2009;Shi et al., 2010)-for example,  and Gao et al. (2008) analyzed the climate characteristics, influences, and causes of heavy snowfall over large areas of South China in January 2008; 3) the investigation of the historical occurrences of snow events in an administrative region or in a mountain range-for example, Sun et al. (2010) analyzed the characteristics of snowfall variability in Northeast China during 1960China during -20154) the exploration of the occurrence and evolution of snowfall and snow events by analyzing the records of relevant snow events based on in meteorological data (Li, 1993); and 5) the analysis of the characteristics and formation mechanism of snow events by analyzing the circulation patterns associated with snowfall (Birsan and Dumitrescu, 2014). Snowfall occurs over large areas in Central China every year and has resulted in severe traffic disruption and crop damage. However, few studies have been performed on winter snowfall in Central China, and studies on historical snow events and comprehensive analyses of snowfall evolution under climate change scenarios are lacking.
In this study, the Huaihe River Basin in Central China was selected as the study area to analyze the evolution of snowfall trends obtained from historical data and under climate change scenarios. This study finds that as more effort is made to mitigate climate change, the snowfall variability caused by climate warming gradually weakens, starting in 2050.

Data
In this study, the Huaihe River, one of the seven largest rivers in China, was selected as the study area. The Huaihe River Basin is located between the Yangtze River Basin and the Yellow River Basin. It has an area of 270,000 km 2 . The Huaihe River is in the climate transition zone between northern and southern China and is located in the warm semihumid monsoon climate zone. The Huaihe River Basin is an important grain production region in China. The snowfall variability in the basin has a great impact on the growth of winter wheat and other crops, large-scale transportation, and power facilities. This study used snowfall data and weather phenomenon data from 199 conventional meteorological stations in the Huaihe River Basin, basin elevation data and meteorological station distribution data. Figure 1 shows the topography and Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 594704 2 meteorological station distribution in the Huaihe River Basin. In addition, the Coupled Model Intercomparison Project phase 6 (CMIP6) datasets and the WATCH Forcing Data methodology applied to ERA5 (WFDE5) reanalysis data (Cucchi et al., 2020) were used to verify the applicability of multiple models in the CMIP6 datasets. The CMIP6 data of 2015-2065 were used to analyze the snowfall evolution in different scenarios.
Before using the multimodel data in CMIP6 for snowfall analysis, it is necessary to evaluate the applicability of multiple models and choose suitable models to analyze snowfall evolution. In this study, the WFDE5 reanalysis data were used to evaluate the multiple models in the CMIP6 data for the years 2015-2018. WFDE5 is based on the ERA5 reanalysis data published by the European Center for Medium-Range Weather Forecasts, which provides high-accuracy global snowfall data (Cucchi et al., 2020). The CMIP data are obtained from climate change reports of the Intergovernmental Panel on Climate Change and provide important references. They have been widely used in the fields of meteorology (Priestley et al., 2020;Zhu et al., 2020), hydrology (Almazroui et al., 2020), ocean research (Grise and Davis, 2020) and others. The CMIP6 datasets included 11 models obtained under four climate change scenarios in 2015-2065 ( Table 1). The ensemble member r1i1p1 was run for all models. In SSP126, the lowest radiative forcing scenario, the radiative forcing is stabilized at 2.6 W/m 2 , and it supports studies with an objective to control the temperature rise within 2°; in SSP245, the lowto-moderate radiative forcing scenario, the radiative forcing is stabilized at 4.5 W/m 2 ; in SSP370, the moderate-to-high radiative forcing scenario, the radiative forcing is stabilized at 7 W/m 2 ; and SSP585, the high radiative forcing scenario, is the only scenario that can intentionally achieve a radiative forcing of 8.5 W/m 2 by the end of the present century . Since different spatial data have different spatial resolutions, prior to the assessment of model quality, the spatial resolutions of all models must be bilinearly interpolated to be consistent with the spatial resolution of WFDE5, i.e., a 0.5°× 0.5°grid.

Introduction and Verification of Model Data
Climate simulations can provide powerful information on climate change trends. However, the quality of the information needs to be assessed before application. In this study, the mean correlation coefficient, root mean square error (RMSE), and standard deviation (STD) were used to assess the climate simulation data quality. The average correlation coefficient reflects the linear correlation between the predicted value and the reference value of the model. The higher the absolute value of the correlation coefficient is, the higher the simulation quality of the model is. The optimal value of the correlation coefficient is 1. The formula is as follows (Taylor, 2001):

R
Cov prsn x,t,AOGCM , prsn x,t,WFDE5 Var prsn x,t,AOGCM Var prsn x,t,WFDE5 (1) where x represents a grid point in the study area, x ∈ (1, N), and t represents a time step in the evaluation period, t ∈ (1, n), where n represents the total number of time steps studied. Additionally, prsn x,t,AOGCM represents the predicted model ensemble results of grid point x at time step t, and prsn x,t,WFDE5 represents the WFDE5 value of grid point x at time step t, which is used as the reference value. Cov is the covariance of the predicted value and the reference value within the same grid point. The RMSE shows the magnitude of the error between the simulated and observed values at the station. The lower the RMSE is, the closer the response pattern is to the reference value, and the optimal value of RMSE is 0. The formula is as follows (Zhang et al., 2018): The STD response model simulates its own stability. The lower the STD is, the better the model stability is. The optimal value of the STD is 1 (relative to the reference value). The formula is as follows (Rosero et al., 2009):

Trend Analysis Method
The climate trend can well express the long-term trends of change in meteorological elements. It is the simplest and most effective method for analyzing the relationship between snowfall and time.
A positive climate trend indicates that snowfall increases over time, while a negative climate trend indicates that snowfall decreases over time. The calculation formula is as follows : where k is determined by the least square method, k × 10 mm/10a is defined as the climate trend, x is the time series, and y is the snowfall.

Mann-Kendall Test for Abrupt Changes
The nonparametric Mann-Kendall (M-K) test is mainly used to test the change trends and abrupt changes of a certain factor. In the time series test, the nonparametric M-K test was first proposed by Mann and Kendall. The advantage of the M-K test is its low requirement of samples the samples do not need to follow a specific distribution pattern and can have a small number of outliers. Therefore, the M-K test is widely used in meteorology, hydrology, and other disciplines. The calculation steps are as follows (Gocic and Trajkovic, 2013): For a time series with n sample sizes, x 1 , x 2 , x 3 , /, x n , let us construct an ordered series: where R i represents the cumulative number of times when X i is greater than X j (1 ≤ j ≤ i). Under the assumption that the time series is random and independent, UF k is defined as follows: In the formula, when k 1, UF 1 0. E(S k ) and Var(S k ) are the mean and variance in cumulative number S k , respectively. When X 1 , X 2 , /, X n are mutually independent and have the same continuous distribution, they can be calculated with the following formulas: Var(S k ) n(n − 1)(2n + 5) 72 (8) Time series X is reversed to obtain series X n , X n−1 , /, X 1 . The above process is repeated while making UB k −UK F (k n, n − 1, /, 1). If the UB k and UF k curves intersect and the intersection point is between the critical lines, then the time corresponding to the intersection point is the time when an abrupt change occurs.

Sliding T-test
The sliding t test is applied to test whether there is a significant difference between the means of two subsequences in the climate factor sequence based on whether there is a significant difference in the two overall means. If the test results exceed a certain level, a sudden change is considered to have occurred.
For time series x with a sample size of n, a moment is specified as the reference point, and the samples of the two subsequences x 1 and x 2 before and after the reference point are n 1 and n 2 respectively. The average values of the two sequences are x 1 and x 2 , and the variance is S 2 1 and S 2 2 . The statistic t is defined as follows (Du et al., 2019): where S is: S n 1 s 2 1 + n 2 s 2 2 n 1 − n 2−2 (10)

Spatial Distribution Characteristics of Snowfall
A spatial distribution map of snowfall in the Huaihe River Basin can be obtained using the inverse distance weighted spatial interpolation of the annual mean snowfall over 68 years from 1951 to 2018. As shown in Figure 2, the overall annual mean snowfall in the Huaihe River Basin differs significantly from south to north and gradually increases from northeast to southwest. The snowfall is lowest in the eastern Shandong Peninsula and highest in southeastern Anhui Province. The snowfall at the provincial boundaries between Shandong Province and Henan, Anhui, and Jiangsu Provinces is approximately 3.7 mm.
Distribution Characteristics of the Climate-related Snowfall Trend Using the annual mean snowfall data for the Huaihe River Basin from 1961 to 2018, we calculated the climate trend at each station. The climate distribution map was obtained using inverse distance weighted spatial interpolation ( Figure 3). As shown in Figure 3, the positive changes in snowfall in the Huaihe River Basin are mainly distributed in western Shandong Province, southwestern and northwestern Jiangsu Province, central Anhui Province, and central and western Henan Province. The highest trend is 11.06 mm · 10 a −1 . Negative changes in snowfall are mainly distributed in northern Hubei Province, southern Henan Province, and southwestern Anhui Province. The lowest trend is −7.22 mm · 10 a −1 . The mean climate trend in the basin is −0.85 mm · 10 a −1 , indicating that the snowfall in the basin showed an overall decrease.

Analysis of the Snowfall Trend
The changes in the number of snow events and changes in snowfall in the Huaihe River Basin between 1952 and 2018 were statistically analyzed. Figure 4A shows that high numbers of snow events occurred in 1956, 1969, 1985, and 2012, and the corresponding snowfall reached 147.8, 242.9, 175.9, and 129.6 mm, respectively. The number of snow events is positively correlated with the annual snowfall, and the number of snow events in a year has gradually decreased with time. Figure 4B shows that high snowfall values occurred in 1993 and 1969, and the annual snowfall gradually decreased with time. Figure 4A,B show that the number of snow events and the amound of snowfall both decrease. According to the snow classification criteria, snow events can be divided into light snow, moderate snow, heavy snow, snowstorms, heavy snowstorms, and severe snowstorms. The time dependence of the number of days with snow of each grade in the Huaihe River Basin from 1951 to 2018 was plotted, as shown in Figure 5. Figure 5A shows that large numbers of snow days occurred in 1971, 1985, and 2012, and the number of snow days in a year gently declined overall. Figure 5B shows that large numbers of days with moderate snow occurred in 1984 and 2010, and the number of days with moderate snow in a year decreased overall. Figure 5C shows that large numbers of days with heavy snow occurred Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 594704 in 1974 and 2000, and the decrease in the annual number of days with heavy snow was more significant than of the decrease in the number of snow days and the number of days with moderate snow. Figure 5D shows that large numbers of days with snowstorms occurred in 1954, 1964, 1985, 1987, 1989, and 1998. Figures 5E,F show that the interannual changes in the number of days with heavy snowstorms and the number of days with severe snowstorms were not significant. Overall, the number of days with light snow, moderate snow, heavy snow, snowstorms, heavy snowstorms, and severe snowstorms decreased over time. The declines in the number of days with heavy snow and the number of days with snowstorms were the most significant. Figure 6 shows the changes in the annual number of snow events, the annual average of the lowest temperature and the annual lowest of the mean temperature in the Huaihe River Basin from 1951 to 2018. The annual average temperature and the annual minimum temperature show an upward trend over time, while the number of annual snowfalls shows the opposite trend, showing a downward trend. It can be seen that there is a negative correlation between temperature and snowfall times.
The nonparametric M-K test and the sliding t test were used to analyze the abrupt changes in snowfall and snowfall days at weather stations in the Huaihe River Basin from 1951 to 2018, respectively, as shown in Figures 7 and 8. In the analysis of the sudden change in the number of snowfall days, Figure 7A shows a sudden change in the number of snowfall days in 1998, while the sudden change in Figure 8A occurred in the five years of 1962, 1979, 1990, 1991, and 2012. Therefore, we take the intersection of the nonparametric M-K test and the sliding t test as our final result. From 1951 to 2018, there was no obvious sudden change in the number of snowfall days. In the snowfall analysis in Figure 7B, changes occurred in 1972, 1986, and 2012, and in Figure 8B, changes occurred in 1986, 1988, and 2010. Therefore, by comparing the two analysis methods, it is found that a sudden change in snowfall occurred in 1986 and that this change in 1986 corresponded to a sudden switch from decreasing to increasing snowfall.

Quality Assessment and Ensemble of Climate Change Scenarios
The quality of the CMIP6 multiscenario data was assessed using the WFDE5 climate reanalysis dataset of 2015-2018. A Taylor diagram was used to compare the quality of the data under different scenarios. It can well describe the closeness of a model (or a group of models) to the observed values (Taylor, 2001). Figure 9 shows the quality assessment results of snowfall predicted by the models under different emission scenarios. In  Figure 9D), and GFDL-ESM4 ( Figure 9D) to the STD of the observed values are less than 1. This result indicates that the changes in the snowfall predicted by these models are smaller than the reference value and that the predicted results are stable. The ratio of the STD of the values simulated by the remaining models to the STD of the  Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 594704 6 observed values is greater than 1. The results indicate that the changes in snowfall predicted by these models are larger than the reference value. The ratio between the RMSE of the results simulated by these models and the RMSE of the reference value shows that models EC-Earth3-Veg, FGOALS-f3-L, and NorESM2-MM have the lowest quality in the corresponding scenarios.
Based on the above quality assessment results, in each scenario, the model with the poorest quality performance (model GFDL-ESM4 in scenario SSP126, model GFDL-ESM4 in scenario SSP245, model GFDL-ESM4 from scenario SSP370, and model EC-Earth3 in scenario SSP585) was excluded. The remaining models in each scenario are high-quality. After the model with the worst performance is excluded from each scenario, the screened models are obtained. The selected model data are grouped according to the characteristics of the corresponding scenario. In each time step and each grid point, the median snowfall is selected and called the ensemble median, which is used for subsequent simulations of the trend of snow mass flux.

Analysis of the Evolution of Snowfall Under Climate Change Scenarios
In this study, the ensemble median method was adopted to obtain the multimodel snowfall (the snowfall under the climate change scenarios is the snow mass flux), and the interannual evolution of the annual mean snowfall of 2015-2065 was plotted, as shown in Figure 10. There are certain differences between the evolution of snowfall in different climate scenarios. In the SSP126 scenario with the lowest emissions, snowfall is expected to reach 3/4 of its current value by 2060. In the SSP545 scenario with the highest emissions, the snowfall will decrease to half of its current value by 2060. In other scenarios with intermediate emissions, the  predicted snowfall is similar and somewhat lower than the current value. Through the analysis of the four climate change scenarios, the overall trends are relatively consistent. Relevant studies mostly use SSP370 as a representative model to analyze climate change (Shao et al., 2019;Butchart et al., 2020). This study also used the SSP370 model to analyze the snowfall trend of the Huaihe River Basin. As shown in Figure 10, the annual mean snowfall will change abruptly in 2050 from decreasing to rapidly increasing and will reach approximately 125% of the value in 2015 by the end of 2065.
In the study of long-term climate change, the interdecadal variability in the climate in China is correlated with the interdecadal variability in large-scale circulation patterns. Studying the variability in meteorological elements can help reveal the changes in circulation at a large scale (Nidheesh et al., 2017). Through mean filtering for interdecadal variability, the influence of low-frequency fluctuations can be eliminated. As shown in Figure 11, the changes in the future trend of snowfall are not significant in the SSP370 scenario and reach the minimum value in 2050, and the linear trend rate is −0.2e −7 kg · m −2 · s −1 · a −1 . The snowfall increases rapidly after 2050, and the linear trend rate reaches 0.5e −7 kg · m −2 · s −1 · a −1 . In 2065, snowfall returns to the snowfall value observed in 2018.

Spatial Variability Characteristics of Snowfall in the Climate Scenarios
Through the accumulation of monthly snowfall, we can obtain the snowfall distribution map of 2015-2065 in future climate scenarios, as shown in Figure 12. The snowfall gradually increases from the coastal areas to the inland mountains and differs significantly from east to west. The cumulative snowfall is highest (reaching 6 e −4 kg · m −2 · s −1 ) in the northwestern part of the basin of Tongbai Mountain, where the Huaihe River originates, and lowest (less than 1 e −4 kg · m −2 · s −1 ) in the eastern part of the basin.
By calculating the climate trend of annual snowfall in each grid over time, the spatial distribution of the climate trend is obtained, as shown in Figure 13. In most areas, the snowfall tendency rate is negative, and the snowfall decreases over time. In particular, snowfall declines slowly in the eastern coastal regions, whereas the snowfall trend is especially low in the western hilly areas, where snowfall declines most significantly. The rate of annual  Based on previous studies, the percentage change in snowfall in 2050 and 2065 relative to 2015 was calculated. Figure 14 shows the spatial distribution of the percentage change in snowfall. Figure 14A shows that by 2050, the snowfall in most areas decreases, with percentage changes of −60% to −20%. Additonally, the snowfall in the southern regions, southeastern coast regions, and northern regions decreases by more than 80%, and the snowfall in other regions decreases by approximately 40%. In contrast, snowfall increases in a small portion of the central regions. Figure 14B shows that by 2065, relative to 2015, snowfall has increased greatly (approximately 175% of that in 2015) in the central and southern regions, and the areas where the snowfall has doubled are mainly concentrated in the southeastern regions. Based on Figure 14A,B, the snowfall variability in the southeastern region is more intense than that in other regions; the snowfall in the northern region decreases significantly, reaching −50%; and the snowfall in most regions in 2065 is higher than that in 2050.

DISCUSSION
The uncertainty in this work is mainly concentrated in the following aspects. First, the spatial resolution of WFDE5 is 0.5°. However, the spatial resolutions of different models used in the different scenarios in CMIP5 are different. The bilinear interpolation method was used to unify the spatial resolution of the above data, to match that of WFDE5, and there is uncertainty in the spatial distribution accuracy. Second, independent of whether CMIP6 or CMIP5 is selected, many assessments and comparisons of the quality of CMIP5 and CMIP6 are underway (Cucchi et al., 2020;Freund et al., 2020;Kim et al., 2020;Nie et al., 2020). Roach et al. found that compared with the previous CMIP5 models, most CMIP6 models are superior in quality in terms of characterizing the East Asian summer monsoon (EASM)-related interannual abnormal precipitation model in East China (Roach et al., 2020). The overall conclusion is that CMIP6 is better than CMIP5. Last, based on the snowfall data obtained at meteorological stations, the snowfall gradually increased from northeast to southwest, and the conclusion derived from the model data obtained under different scenarios is that the snowfall gradually increases from east to west (from coastal to inland regions). Although the two conclusions are slightly different, overall, there is less snowfall in coastal regions than in inland regions. The study found that high emissions will lead to an extreme reduction in snowfall in the basin, which is related to the reduction in snow events caused by climate warming. The underlying mechanism needs to be further investigated. Sun et al. (2010) analyzed the characteristics of snowfall changes in Northeast China from 1960 to 2015, and found Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 594704 differences in the amount of snowfall in time, years in which major changes occurred, and space. In the research of this article, we not only study the historical evolution but also add the trend of snowfall in the coming decades based on climate scenario data. In snowfall research based on satellite imagery, there are numerous studies on snow cover and the changes in snow cover in middleand high-latitude areas and high-elevation areas using CMIP5 and CMIP6 data to predict snow cover in future climate scenarios (Birsan and Dumitrescu, 2014). This research method is similar to the research in this study, and the Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 594704 11 research conclusions are of great significance to regional climate change research, disaster prevention and disaster mitigation.

CONCLUSIONS
From the analysis of the climate trends based on meteorological observation data and the analysis of abrupt changes in snowfall, the total number of snow days, and the snowfall evolution in the Huaihe River Basin under climate change scenarios, the following conclusions are drawn.
(1) As shown by the analysis of the historical snowfall data from the meteorological stations, the snowfall amount and the number of snow events in the basin decreased over time in the period 1951-2018. Based on the analyze with the nonparametric M-K test and sliding t test, 1986 was the year in which a sudden change occurred in the average annual snowfall amount, but there was no obvious sudden change in the number of snow events.
(2) The analysis of snowfall evolution in the SSP370 scenario shows that 2050 is an important year because the snowfall amount will change from decreasing to increasing and that in 2065, the snowfall amount in most areas of the basin will be significantly larger than that in 2050. Additionally, the variability in snowfall is highest in the southeastern coastal region of the basin, indicating high instability in this region.
Under the background of climate change, snowfall will generally decrease over time. The snowfall shows a decline from 2015 to 2050 and reaches its lowest value in 2050. With continuous attention paid to climate change issues, the annual snowfall will be in the recovery stage in 2050-2065 and will rebound significantly.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
HD, QL, and XZ: research conceptualization, methodology, formal analysis, and writing first draft preparation. XZ, ZZ, JS, and YH: supervision. All authors contributed to manuscript revision, read, and approved the submitted version.