Climate Change Impacts on Cold Season Runo ﬀ in the Headwaters of the Yellow River Considering Frozen Ground Degradation

: Climate change has e ﬀ ects on hydrological change in multiple aspects, particularly in the headwaters of the Yellow River (HWYR), which is widely covered by climate-sensitive frozen ground. In this study, the annual runo ﬀ was partitioned into four runo ﬀ compositions: winter baseﬂow, snowmelt runo ﬀ , rainy season runo ﬀ , and recession ﬂow. In addition, the e ﬀ ects of global warming, precipitation change, and frozen ground degradation were considered in long-term variation analyses of the runo ﬀ compositions. The moving t -test was employed to detect change points of the hydrometeorological data series from 1961 to 2013, and ﬂow duration curves were used to analyze daily runo ﬀ regime change in di ﬀ erent periods. It was found that the abrupt change points of cold season runo ﬀ , such as recession ﬂow, winter baseﬂow, and snowmelt runo ﬀ , are di ﬀ erent from that of the rainy season runo ﬀ . The increase in winter baseﬂow and decrease in snowmelt runo ﬀ at the end of 1990s was closely related to global warming. In the 21st century, winter baseﬂow presented a larger relative increase compared to rainy season runo ﬀ . The correlation analyses indicate that winter baseﬂow and snowmelt runo ﬀ are mainly controlled by water-resource-related factors, such as rainy season runo ﬀ and the accumulated precipitation in cold season. To analyze the global warming impacts, two runo ﬀ coe ﬃ cients—winter baseﬂow discharge rate ( R w ) and direct snowmelt runo ﬀ coe ﬃ cients ( R s )—were proposed, and their correlation with freezing–thawing indices were analyzed. The increase of R w is related to the increase in the air temperature thawing index (DDT), but R s is mainly controlled by the air temperature freezing index (DDF). Meanwhile, the direct snowmelt runo ﬀ coe ﬃ cient ( R s ) is signiﬁcantly and positively correlated to DDF and has decreased at a rate of 0.0011 / year since 1980. Under global warming, the direct snowmelt runo ﬀ (runo ﬀ increment between March to May) of the HWYR could decrease continuously in the future due to the decrease of accumulative snow in cold season and frozen ground degradation. This study provides a better understanding of the long-term runo ﬀ characteristic changes in the HWYR.


Introduction
Hydrological impacts of climate change have drawn worldwide concern, especially in cold regions due to the existence of the climate-sensitive cryospheric hydrogeology environment, such as frozen ground, seasonal snowpack, etc. [1][2][3][4][5]. Climate change impacts the hydrological change of cold regions, both directly through changes in climate forcing (precipitation and evaporation) and indirectly This study analyzed the long-term runoff changes in the HWYR, with the aims to (1) characterize the long-term runoff variations of different seasons in different change periods, (2) analyze the potential causes of runoff change in considering global warming, precipitation change, and frozen degradation, (3) analyze frozen ground degradation impacts on cold season runoff; and (4) provide a better understanding of runoff response to climate change in the HWYR.

Study Area
The Yellow River originates from the east of the Tibetan Plateau. The catchment above Jimai station, possessing a large permafrost area and a cold climate, was selected to analyze hydrological change of the headwaters of the Yellow River (HWYR). The study area, with a drainage area of about 45,000 km 2 , is embraced by the Bayan Har mountains, the Buqing mountains, and Animaqing mountain from the south, west, to northeast, with elevation ranges from 3951 to 5345 m (Figure 1a). Two major large lakes, Gyaring and Ngoring, are located in the Northwest of the catchment. Two meteorological stations, Maduo and Dari, are located in this area to observe climate change. The Jimai hydrological station is located at the outlet of the catchment and gauges the runoff change. As no water-supply-related glacier existed in the study area, the cryosphere mentioned in this study focused on frozen ground and seasonal snowpack. The study area is fully underlain by frozen ground, including seasonally frozen soil and permafrost. The distribution of the frozen ground in the study area in the 1980s (Figure 1b) and the 1990s (Figure 1c) were obtained from the Map of Snow, Ice and Frozen Ground in China (1988) and the Map of Permafrost on the Qinghai-Tibetan Plateau [35], respectively. The two maps were provided by the Environmental and Ecological Science Data Center for West China, National Natural Science Foundation of China (http://westdc.westgis.ac.cn). According to these two maps, frozen soil of the study area was briefly classified as seasonally frozen soil and permafrost. As shown in Figure 1b,c, permafrost distribution in the study area has shrunken significantly since the 1980s.

Data Processing
Two kinds of datasets-hydrometeorological data and remote sensing data-were used in this study. Note that we defined a freezing-thawing year and a hydrological year in the following data processing and analyses, since the time cycles of the freeze-thaw process and hydrological process are different from each other. When calculating the accumulative snow in winter, using the freezing index (DDF) and thawing index (DDT), we regarded the period from September to the following August as a freezing-thawing year (FTY). However, we took March to the following February as a hydrological year (HY) when analyzing the runoff changes in different seasons.

Hydrometeorological Data
The daily precipitation and temperature data from HY 1961 to 2013 at Dari and Maduo meteorological stations were provided by the China Meteorological Administration (CMA). The daily runoff from 1961 to 2013 (lack of 1990 daily runoff data) and monthly runoff from 1961 to 2013 at the Jimai outlet hydrological station were provided by the Yellow River Conservancy Commission (YYCC). The annual precipitation (P, mm), annual mean air temperature (T, • C), and precipitation intensity (I, mm/day) were obtained from the daily meteorological data and used in the following analyses. The precipitation intensity (I) is defined as the ratio of the total annual precipitation amount to the total number of days of rain in a hydrological year.

Snow and Frozen Ground
A simple threshold method was used to divide the precipitation into rain, snow, and sleet in this study (described in Section 3.1). The rain, snow, and sleet were calculated separately with the respective daily temperature and precipitation of the two meteorological stations by the threshold method. The accumulative snow from September to May was used to indicate the seasonal snowpack. The Moderate Resolution Imaging Spectroradiometer (MODIS) snow product-the MOD10A2 dataset (the MODIS/Terra Snow Cover 8-day L3 Global 500 m Grid) is the most widely used snow product. The MOD10A2 dataset is produced and distributed by the National Aeronautics and Space Administration (NASA) Distributed Active Archive Center (DAAC) and is available at the website of the National Snow and Ice Data Center (NSIDC) (http://nsidc.org/). The MOD10A2 dataset was used to validate the obtained snow by the proposed method in this study. The MOD10A2 dataset from 2001 to 2013, with 8-day temporal resolution and 500 m spatial resolution, was collected and used to detect the snow cover seasonal variations in the HWYR. As MOD10A2 is an 8-day composite dataset, the cloud cover effects were effectively minimized. Many studies have proved the validity and efficiency of the MOD10A2 in mapping snow cover [36,37]. The percentage snow cover area of each image was defined as the ratio of snow cover area to the whole basin area. The mean value of the snow cover area of every same date image from 2001 to 2013 was used to represent the snow cover distribution of the date.
The air temperature freezing index (DDF) and air temperature thawing index (DDT), as commonly used and easily obtained indices, were used to indicate the frozen depth and thawed depth variations of frozen ground in this study. The DDF and DDT were defined as the absolute accumulative negative air temperature and the accumulative positive air temperature in a freeze-thaw year (FTY) (from September to August), respectively. The physical basis of using DDF and DDT to indicate the frozen ground is derived from the commonly used active layer frozen-thawed depth simulation model Stefan equation [38]. The DDF and DDT were calculated separately with the daily air temperature obtained from the respective station.
The mean value of DDF (or DDT) of the two meteorological stations was used to represent the frozen (or thawed) ground depth of the study area. All climate factors, such as P, T, I, rain, sleet, and snow, used in the following analyses were the mean value of the two stations.

Methods
In this study, we proposed a simple but efficient way to distinguish snow, sleet, and rain by using only daily precipitation and daily air temperature (Section 3.1). We used the moving t-test (Section 3.2) to detect abrupt change points and used flow duration curves (Section 3.3) to detect runoff characteristics changes. Then correlation analyses were applied to reveal relationships and interactions between runoff, climatic factors, and freezing-thawing indices.

Identification of Precipitation Types
The snow, sleet, and rain records are available only from 1956 to 1979 in the two meteorological stations. We analyzed the snow, sleet, and rain distribution from the data, along with air temperature, and the results are shown in Figure 2. Snow rate increased as temperature decreased, and rain rates increased as temperature increased. Sleet rate showed a normal distribution in relation to temperature ( Figure 2). Using 50% (the horizontal red dashed line in Figure 2) as a critical value, we were able to find the critical temperature for snow, sleet, and rain. According to the obtained results, we defined a simple method to distinguish snow, sleet, and rain at the two stations, as indicated in Table 1. Table 1. Critical temperature of snow, sleet, and rain.

Snow Sleet Rain
Jimai The obtained critical temperature (1 • C) of snow is a commonly used threshold value in distinguishing snow and rain in previous studies in the HWYR [39,40]. The obtained critical temperature of snow, sleet, and rain was used to distinguish precipitation types in the two meteorological stations, respectively, and the mean values of snow, sleet, and rain of the two stations were used in the following analyses. The correlation coefficient between obtained snow and sleet, and observed snow and sleet is 0.75 and 0.56 at an annual scale.

The Moving t-Test
The original t-test is an efficient and valid way to assesses whether the mean values of two subsamples are significantly different from each other. For the original t-test, there are two major statistic values, t and s, which are defined as follows [29]: As shown in the above Equations (1) and (2), x i , s, and n i are the mean value, standard deviation, and the size of the two independent samples (i = 1, 2). The moving t-test applies the original t-test to a subseries with a given length n before and after the potential change points (n = n 1 = n 2 ). Therefore, the length of the moving window is set as 2n to detect if the middle point is an abrupt change point. If the calculated t-value is larger than the critical value t α , then the middle point of the moving window is accepted as a change point at the significant level α. In this study, n = 10, and the critical value was t α = 2.1 for α = 0.05. The moving t-test was employed to detect the abrupt change points of the long-term runoff and climate factors series.

Flow Duration Curves
FDC has been generally used since approximately 1915, which provides a simple and efficient way to show the runoff distribution in different percentage of exceedance (P e ) ranges [41,42]. To plot the FDC, firstly, the runoff in the corresponding period should be sorted from large values to small values. Then, the percentage of exceedance (P e ) of the corresponding runoff is calculated by the following equation: where j is the sequence number of the sorted runoff series in the corresponding period and m is the total number of the runoff series. In this study, the daily runoff from 1961 to 2013 was divided into three periods according to the long-term runoff variations pattern (described in Section 4.1). The daily runoff in the corresponding period was integrated into one-year mean daily runoff to represent the daily runoff change regime in the corresponding period. In this study, to maintain the consistence of the length of daily runoff series in each year, the last daily runoff of a leap year was removed and the m in Equation (3) was set as 365.

Runoff Compositions and Their Correlations
The moving t-test was employed to detect the change points of the long-term annual mean runoff (AMR) variations. As shown in Figure 3a, there are two abrupt change points, 1989 and 2002, dividing the long-term runoff change into three periods: I. pre-change period ; II. 1990s (1990II. 1990s ( -2002; and III. 2000s (2003III. 2000s ( -2013. The abrupt change points obtained in this study are similar to previous studies in the Yellow River sub-basins using other approaches [43,44]. The intra-annual variations of daily runoff in the three periods are displayed in Figure 3b. Compared to the pre-change period, period II is a low-flow period, and period III is a recovery period. Although runoff volume changed significantly in the three periods, the seasonal patterns showed similarity. Daily runoff intra-annual variations can be generally divided into four stages: winter baseflow (WB) (January-February), snowmelt runoff (SR) (March-May), rainy season runoff (RSR) (June-September), and recession flow (RF) (October-December) account for 4%, 15%, 61%, and 20% of annual runoff, respectively. As runoff variations from October to the following February can be integrated into a full recession process, we defined March to the following February as a hydrological year (HY) in following runoff change analyses ( Figure 3b). The correlation analyses were applied to reveal the four runoff composition interactions ( Table 2). Rainy season runoff is a major runoff composition before cold season runoff and has continuous effects on recession flow (RF), winter baseflow (WB), and snowmelt runoff (SR). The correlation coefficients between rainy season runoff (RSR) and RF, WB, and SR are 0.85, 0.80, and 0.54, respectively, reducing along with the hydrological processes. The correlation coefficient between RSR and SR is relatively lower due to depletion of recession processes and disturbance of snowmelt recharge. However, SR still showed a high correlation with RSR (0.54), RF (0.68), and WB (0.60), as the continuous groundwater discharge from recession processes affects SR. To remove groundwater discharge impacts on SR, we defined direct snowmelt runoff (DSR) as the difference of SR and WB. As shown in Table 2, the correlations of DSR and RSR, and RF and WB are obviously lower compared to the correlation of SR and RSR, and RF and WB. The impacts of RSR, RF, and WB are effectively removed. The correlation coefficient between DSR and WB is negative; however, the value is low (−0.08), thus the correlation is negligible.

The Changes of Runoff Compositions
Flow duration curves (FDC) of the three periods are shown in Figure 4a, and the relative runoff variations of period II (1990s) and period III (2000s) to period I (the pre-change period) are displayed in Figure 4b. Points in red, yellow, blue, and purple indicate spring (SR), rainy season runoff (RSR), recession flow (RF), and winter baseflow (WB), respectively. Dots in the blue line, square points in the black line, and triangle points in the red line represent FDC of period I, period II, and period III, respectively. Percentage of exceedance (P e ) of RSR and WB range from 0% to 30% and 80% to 100%, respectively. P e values of SR and RF range from 30% to 80%. The runoff of period II is lowest among the three periods in almost all ranges of P e . Runoff of period III returned as a same level as period I, but WB showed a relatively larger increase (Figure 4a). Using period I as a benchmark, the relative changes of period II and period III are presented in Figure 4b. For period II, both high flow (P e < 30%) and low flow (P e > 70%) decreased by 30-40%, medium flow (30% < P e < 70%) decreased by 10-20%. The runoff of period III showed a recovery; the high flow and a part of the medium flow increased by 0-20%. However, the low flow, which was dominated by winter baseflow (WB), increased by more than 40%. Meanwhile, a part of the medium flow (40% < P e <60%), which was dominated by SR showed a 10% reduction. The abrupt change points of snowmelt runoff (SR), rainy season runoff (RSR), recession flow (RF), winter baseflow (WB), annual mean runoff (AMR), and direct snowmelt runoff (DSR) are detected by the moving t-test ( Figure 5). As mentioned above, variations of AMR from 1961 to 2013 can be divided into three periods by two abrupt change points: 1989 and 2002. RSR, as the major runoff composition, has two change points (1989,2002) which are exactly the same as AMR. WB has two abrupt change points: 1989 and 1999. RF has two abrupt change points: 1989 and 2001. SR and DSR have only one abrupt change point: 1997. These abrupt change points are circled in Figure 5.
Runoff long-term inter-annual variations can be partitioned into different periods by these obtained abrupt change points. The mean values, relative changes to reference periods, and change rates in different periods of different runoffs are shown in Table 3. In the 1990s, AMR and RSR decreased by 28.15% and 29.57% relative to the reference period , respectively. In the 21st century, RSR and AMR increased by 3-5% relative to the reference period. The mean value of WB decreased by 20% in the 1990s, but presented a significant increase trend (at 5% level) with a rate of 1.10 m 3 /s/year. In the 21st century, WB continuously increased at a rate of 3.66 m 3 /s/year and presented a larger increase (42.03%) compared to the recovery of RSR and AMR. The RF showed a reduction in both the 1990s (6.9%) and 21st century (6.4%), however, increased significantly (p < 0.01) in the 21st century at a rate of 6.04 m 3 /s/year. SR deceased by 8.64% in the period 1998-2013, however, DSR showed a more considerable decrease (31.21%). Both SR and DSR showed no significant change in trend during the period 1998-2013.

Precipitation Changes
In this section, we investigate precipitation changes, including the changes of precipitation amount, precipitation intensity, and precipitation types. For precipitation types, we used a threshold method mentioned in Section 3.1 to distinguish snow, sleet, and rain. The intra-annual snow variations and inter-annual snow variations in the HWYR are shown in Figure 6a,b, respectively. The monthly snow of Maduo and Jimai stations is derived from daily precipitation from HY 1961 to 2013. The 8-day snow cover variations in one year was obtained from 2001 to 2013 from the MOD10A2 dataset. The time span of snow obtained from the observation station is different from the remote sensing dataset. As indicated in Figure 6a, however, the intra-annual variations of snow and snow cover area are consistent with each other, mutually proving the efficiency of the data. As indicated by Figure 6a, snow starts accumulating in September and begins melting in March, and is completely melted during June to July, and no permanent snowpack exists in this area. The percentage snow cover area ranges from 30% to 50% in winter. The accumulative snow from September to May is defined as accumulative snow in the cold season (ASC) and is taken as seasonal snowpack in this study. The long-term variations of ASC are displayed in Figure 6b; the ASC showed a slight increase from 1980 to 1997 and then decreased during the period 1998-2013.
We focused our study on analyzing snow, sleet, and rain in the cold season, which are more related to cold season runoff; moreover, 98% of snow happened in the cold season (September-May). To investigate the snowmelt runoff and winter baseflow changes, we summarized sleet and rain from September to February, and March to May, separately. The moving t-test results of precipitation amount, precipitation intensity, and precipitation types are shown in Figure 7, and the respective mean values, relative changes, and change rates are summarized in Table 4.  The long-term variations of precipitation and rain are divided into two periods. Precipitation and rain increased by 12.12% and 16.09% in the 21st century. Although a decrease in precipitation in the 1990s in the HWYR was reported by previous studies [31,32], the decrease is not significant, as there is only one change point (2002) detected by the moving t-test (Figure 7). However, precipitation intensity decreased by 8.36% in period 1988-1995 and increased by 4.49% in period 1996-2013. To consider the impact of changes in precipitation types on cold-season runoff, we identified cold-season precipitation into three groups: (Snow)  , (Sleet + Rain)  , and (Sleet + Rain) (3)(4)(5) , indicating snow happened in the cold season, sleet and rain happened between September and the following February, and sleet and rain happened between March and the following May.

The Changes of Temperature and Freezing-Thawing Indices
In this study, we used annual mean air temperature to indicate global warming and freezing-thawing indices (DDF, DDT) to indicate frozen ground degradation following global warming.
As shown in Figure 8, air temperature (T) and thawing index (DDT) have two abrupt change points, and freezing index (DDF) has only one abrupt change point. As seen in Table 5, both T and DDT have two increase periods. T showed a slight increase (0.52 • C) in the period 1986-1996 but a substantial increase (1.27 • C) in the period 1997-2013, with a significant increase rate 0.07 • C/year (p < 0.05). Similar to T, DDT increased by 3.26% in the period 1986-1991 and increased by 14.29% with a significant increase rate of 11.42 • C·day/year in the period 1992-2013 (p < 0.01). The abrupt change point (1997) of DDF is relatively late compared to DDT and air temperature. DDF showed a 14.99% reduction in the period 1998-2013 compared to the period 1961-1997.

Discussion
In previous studies, the decrease of annual mean runoff (AMR) and rainy season runoff (RSR) in the low-flow period (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002) is generally attributed to decreases in precipitation, and increases in evaporation [29,31]. Meanwhile, the significantly decreased precipitation intensity is also a crucial factor that resulted in an evaporation increase in the 1990s. The increase in precipitation after 2002 offset the effects of increased evaporation, and is the main reason leading to a recovery of the runoff in the recovery period (2003-2013) [31,32].
This study focused on the cold season runoff changes. The cold season runoff, recession flow (RF), and winter baseflow (WB) shared one abrupt change point (1989) with rainy season runoff (RSR). It is a rational result, as RF and WB are results from the recession of RSR. In addition, RSR, RF, and WB are highly correlated with each other (Table 2). However, the second abrupt change points of RF and WB have a three-year (1999) deviation and a one-year (2001) deviation compared to RSR (2002), indicating the runoff supplement and runoff processes changed during the 1990s ( Figure 5). Both cold season precipitation and freezing-thawing indices showed abrupt changes in the 1990s (Tables 4  and 5). The potential reasons for the inconsistent variations of RSR, WB, and RF are the changes of cold season precipitation and frozen ground. Meanwhile, previous studies reported that the existence of frozen ground impedes infiltration, resulting in higher direct flow ratio [23][24][25]. However, the ways in which the snowmelt runoff changes along with cold season precipitation changes and frozen ground degradation are still unknown. In this section, the precipitation change impacts on cold season runoff and the potential hydrological effects of frozen ground degradation are discussed.

Precipitation Change Impacts on Cold Season Runoff
The correlations between cold season runoff and cold season precipitation are shown in Table 6. By comparing Tables 2 and 6, we found that the winter baseflow (WB) is mainly controlled by recession flow (RF) and rainy season runoff (RSR) but is also affected by (Sleet + Rain)  . Snowmelt runoff (SR) is mainly controlled by (Snow  + (Sleet + Rain)  ), RF, and WB. After removing the WB impacts on SR, the correlation coefficient between direct snowmelt runoff (DSR) and Snow  increased to 0.54. The correlation coefficient between direct snowmelt runoff (DSR) and Snow  + (Sleet + Rain)  is 0.55. The correlation coefficient between direct snowmelt runoff (DSR) and Snow  + (Sleet + Rain) (3)(4)(5) is 0.29. Therefore, DSR is mainly fed by the snowmelt from the accumulative snow in cold season (ASC), therefore, the changes of precipitation types between September and May have little impacts on direct snowmelt runoff change. The recession flow (RF) has a relatively weak correlation with cold season precipitation, which is mainly controlled by rainy season runoff (Table 2).  (3)(4)(5) 0.07 −0.15 −0.24 0.03 (Sleet + Rain)  0.36 0.49 0.32 0.13 Snow (9-5) + (Sleet + Rain) (3)(4)(5) −0.06 0.19 0.29 0.07 Snow  + (Sleet + Rain)  0.25 0.60 0.55 0.14 The increased (Sleet + Rain)  , caused by the increase of air temperature in the 1990s, is an extra supplement for RF and WB compared to the pre-change period . RSR and RF showed a decrease trend in the 1990s, however, WB increased significantly at a rate of 1.10 m 3 /s/year (Table 3). In addition, (Sleet + Rain)  showed no significant change until 1999, indicating that the change of supplement could not fully explain the gradually increase phenomenon of WB in the 1990s. Therefore, another important factor, frozen ground degradation, should be involved in the explanation of WB changes. In the 2000s, both RF and WB showed a significant increase trend, partly because the extra supplement, (Sleet + Rain)  , increased by 6.35% from 2000 to 2013. The continuous effects of frozen ground degradation on WB increase in the 2000s are further discussed in the following section.

Frozen Ground Degradation Impacts on Winter Baseflow
As indicated by the obtained results in Table 2, winter baseflow is significantly correlated to rainy season runoff (RSR) and recession flow (RF). Recession flow is a depletion process of rainy season runoff, in which surface runoff gradually decreased and was replaced by groundwater discharge. Winter baseflow is mainly fed by groundwater discharge, which is related to the water stored in the active layer and aquifer (WRSQ) in the rainy season [45]. Due to the lack of groundwater observation data, RSR was used as an index to indicate volume of WRSQ in this study. Additionally, assume that RSR is linearly related with WRSQ when the aquifer characteristic does not evidently change.
To investigate the variations of winter baseflow (WB), the related correlation analyses are plotted in Figure 9a,b. The rainy season runoff (RSR) and winter baseflow (WB) showed similar variations, as shown in Figure 9a. As shown in Figure 9b, WB is significantly positively correlated with RSR (R 2 = 0.65, p < 0.01) which supports the theory that winter baseflow is fed by water resources stored in aquifers in the rainy season. If we assume that aquifer characteristics, such as hydrological conductivity and storage capacity, have not undergone dramatic change, WB and RSR could show the same variation. As analyzed in Section 4.1, however, the moving t-test results indicate the abrupt change point of WB appeared three years earlier than RSR, indicating the change of the relationship between WB and RSR at the end of 1990s. To analyze the effects of frozen ground degradation to winter baseflow (RSR) change, we defined winter groundwater discharge rate as the ratio of WB to RSR. The long-term variations of freezing-thawing indices and winter groundwater discharge rate (R w ) are shown in Figure 9c. The labelled years indicate the abrupt change points of freezing-thawing indices and R w detected by the moving t-test. DDT showed a continuous significant increase (R 2 = 0.70, p < 0.01) since 1986, DDF remained stable between 1985 and 1997 but showed a significantly decrease (R 2 = 0.22, p < 0.01) in the period 1997-2013. The results indicate permafrost thawed depth significantly increased since 1985, and frozen depth significantly decreased since 1997. The enlarged thawed depth provided more room for water storage in the active layer in the rainy season. Meanwhile, the decrease of frozen depth leads to taliks formation which provided extra groundwater discharge channels. The combination effects of thawed depth and frozen depth changes lead more grounder water discharge in winter. The R w remains about 0.1-0.2 but gradually significantly increased (R 2 = 0.34, p < 0.01) since 1985 and reached about 0.2-0.3 in the 21st century.
The correlation of R w and freezing-thawing indices after temperature significantly increased (1985) are plotted in Figure 9d. The R w is significantly and negatively (R 2 = 0.14, p < 0.01) correlated with DDF, and is significantly and positively (R 2 = 0.33, p < 0.01) correlated with DDT. This finding indicates the increase in thawed depth and decrease in frozen depth will increase winter groundwater discharge rate, which is consistent with the understanding that under the warming climate the deeper permafrost table and the reduction frozen depth of active layer will lead more groundwater discharge in the winter [12,46]. The higher correlation coefficient between DDT and R w (R 2 = 0.33, p < 0.01) indicated that the deeper thawed depth played a more important role in influencing groundwater discharge rate change. However, the correlation coefficient between thawing-freezing indices and R w (R 2 = 0.33, R 2 = 0.14) is lower than the correlation coefficient of RSR and WB (R 2 = 0.65, p < 0.01), which indicates water resources stored in the rainy season is the main controlling factor of WB, and global warming induced frozen ground degradation altered the relationship between WB and WRSQ. Meanwhile, R w in the period 2000-2013 presented a more distinct fluctuation compared to period the 1985-2000 (Figure 9c). This may be related to the increased sleet and rain from September to February since 1999, which provided an extra supplement to the winter baseflow.

Frozen Ground Degradation Impacts on Direct Snowmelt Runoff
The correlation analyses matrix between cold season runoff and cold season precipitation indicates that direct snowmelt runoff is mainly fed by accumulative snow in the cold season (ASC) ( Table 6). The direct snowmelt runoff (DSR) and respective runoff coefficient changes are plotted in Figure 10a,b, respectively. The direct snowmelt runoff coefficient (R s ) was defined as the ratio of direct snowmelt runoff depth (mm) to accumulative snow in the cold season (ASC) (mm). As shown in Figure 10a, the DSR slightly increased in the period 1961-1997 and decreased in the period 1998-2013, presenting a similar variations pattern with ASC plotted in Figure 8b. The change pattern can be generally divided into two periods; a pre-change period ) and a decrease period (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013). The DSR remains at the same level (mean value 52.76 m 3 /s) during the period 1986-1997 but sharply decreased to 36.30 m 3 /s during the period 1998-2013. The decrease in the DSR (31.21%) is three times larger than the decrease in SR (8.64%) after removing the increasing winter baseflow (WB) impacts (Table 3). Due to the impacts of gradually increased WB in the 1990s, SR presented a slight increase trend during the 1990s, although snow decreased.
Notably, the R s of DSR shown in Figure 10b displayed a different variation pattern; the R s was approximately 0.16 from 1961 to 1980 but decreased with a rate of 0.0011/year from 1981 to 2013 (significant at 10% level (p < 0.1)). The variation in the direct snowmelt runoff coefficient (R s ) is different from direct snowmelt runoff, indicating there is another factor controlling snowmelt runoff generation besides accumulative snow in cold season (ASC). Since frozen ground is characterized as having low permeability, one possible reason for the R s decrease is the decreased frozen depth and shortened frozen duration, leading to more snowmelt infiltration in the ground in spring rather than directly contributing to the runoff. As indicated by previous studies, the hillslopes produce more runoff than those without permafrost, moreover, snowmelt runoff in cold regions has a higher direct runoff ratio than the storm runoff in watersheds without frozen ground [22,47,48].
To analyze frozen ground degradation impacts on direct snowmelt runoff generation, DDF and DDT were used to represent the frozen depth and thawed depth of the frozen ground. The correlations between freezing-thawing indices and direct snowmelt runoff coefficient (R s ) are shown in Figure 11. As indicated by the comparison of Figure 11a,b, the direct snowmelt runoff coefficient (R s ) is more related to DDF (0.42, p < 0.01) rather than to DDT (0.16, p < 0.01). This indicates that ground frozen depth in winter controlled direct snowmelt runoff generation, and with the deeper frozen depth and the longer frozen duration in the winter, more direct snowmelt runoff will be produced in the spring.
The results indicate frozen ground degradation will reduce direct snowmelt runoff generation.

Frozen Ground Degradation Impacts on Water Balance
In this study, we ignored frozen ground degradation impacts on water balance in the HWYR and focused on its impacts on winter baseflow and snowmelt runoff. Nevertheless, the frozen ground degradation impacts have been considered in long-term water balance analyses lately [7,32,34]. Duan [7] used a sensitive-based method to analyze permafrost effects on runoff change of watersheds in north-eastern China and concluded that permafrost thaw increased runoff, but it is difficult to separate the permafrost degradation impacts from complex driving factors. Based on a modified decomposition method, Wu [32] also found that the permafrost degradation presented a positive effect, which is contrary to temperature impacts on runoff change in the headwaters of the Yellow River (HWYR). Sun [49] used a water balance simulation model to simulate and quantify the runoff change in the HWYR. The study found that total permafrost degradation will increase discharge at an annual scale. However, Wang [34] reported a different conclusion based on the Budyko framework, indicating that frozen ground degradation decreases the runoff of the HWYR by about 34%. In addition, Zheng [50] and Zhang [51] used a physical based model with different freeze-thaw process scenarios to simulate runoff in cold regions, and the results showed that the simulated runoff generally reduced when the soil freeze-thaw mechanism is not considered.
As described above, the conclusions of frozen ground degradation's contribution to water balance are still debatable. Meanwhile, as indicated by the moving t-test and the statistical results, the change in precipitation, precipitation intensities, and temperature can explain the annual mean runoff changes. The results agree with most previous studies that the long-term annual runoff change in the HWYR is mainly controlled by climate changes (decrease in precipitation and increase in evaporation) [29][30][31][32]. Furthermore, it is difficult to separate the frozen ground degradation impacts on water balance, because the increase in temperature affects actual evaporation and the frozen ground at the same time [32]. Considering frozen ground degradation impacts on water balance by a statistical method is out the scope of this study. However, frozen ground impacts on annual water balance is still an unsolved question which deserves further analyses in future research.
In period II, the decrease in high flow (P e < 30%) and low flow (P e > 70%) presented the same levels (about 30%) as period I. In period III, however, low flow increased by 40% but high flow increased by less than 20% relative to period I. (2) The intra-annual variations of runoff in the HWYR can be divided into winter baseflow (January to February), snowmelt runoff (March to May), rainy season runoff (June to September), and recession flow (October to December), which account for 4%, 15%, 61%, and 20% of annual runoff, respectively. (3) The winter baseflow is mainly controlled by rainy season runoff, but the frozen ground degradation altered the relationship between rainy season runoff and winter baseflow. Frozen ground degradation increased ground water discharge rate in winter. (4) Direct snowmelt runoff remained at the same level during the period 1961-1997 but decreased 31.21% in the period 1997-2013 relative to the period 1961-1997. The direct snowmelt runoff coefficient (R s ) has been linearly decreasing at a rate of 0.0011/year since 1980, which is significantly controlled by DDF.
According to the obtained conclusions summarized above, we can deduce that with the decreasing accumulative snow in the cold season and the degrading frozen ground, the direct snowmelt runoff of the HWYR may decrease continuously under a warming climate in the future.