Increasing Winter Baseflow in Response to Permafrost Thaw and Precipitation Regime Shifts in Northeastern China

Rapid permafrost thaw and precipitation regime shifts are altering surface and subsurface hydrological processes in arctic and subarctic watersheds. Long-term data (40 years) from two large permafrost watersheds in northeastern China, the Tahe River and Duobukuer River watersheds, indicate that winter baseflows are characterized by significant positive trends of 1.7% and 2.5%·year−1, respectively. Winter baseflows exhibited statistically significant positive correlations with mean annual air temperature and the thawing index, an indicator of permafrost degradation, for both watersheds, as well as the increasing annual rainfall fraction of precipitation for the Duobukuer River watershed. Winter baseflows were characterized by a breakpoint in 1989, which lagged behind the mean annual air temperature breakpoint by only two years. The statistical analyses suggest that the increases in winter baseflow are likely related to enhanced groundwater storage and winter groundwater discharge caused by permafrost thaw and are potentially also due to an increase in the wet season rainfall. These hydrological trends are first apparent in marginal areas of permafrost distribution and are expected to shift northward towards formerly continuous permafrost regions in the context of future climate warming.


Introduction
Climate warming has been driving hydrologic changes across the globe, but especially in high latitude and high altitude regions [1].Obvious links between climate warming and hydrology, such as increasing evaporation [2], changing precipitation distribution [3], and melting glaciers [4], have been widely investigated and identified as important drivers affecting hydrological processes [5].However, more recently, other indirect impacts of climate warming on hydrological processes have been considered.For instance, permafrost thaw due to climate warming [6,7] alters how water is routed and stored in watersheds, and thus impacts both surface and subsurface hydrology [8,9].Also, a temperature-induced shift of precipitation from snow towards rain and seasonal shifts of precipitation can influence hydrologic processes and regimes [10,11].These topics have recently received increased attention in both research and management communities, as discussed in more detail below.
A warming climate results in less precipitation falling as snow and an earlier onset of snowmelt, and these changes in turn influence the timing of the intra-annual streamflow distribution [10].Advances in the timing of peak spring season flows due to earlier snowmelt have been observed in the western United States [12,13] and Siberia [14].Also, seasonal shifts in precipitation and runoff impact groundwater recharge rates and, consequently, baseflow that is predominantly sourced from the release of groundwater storage [15].However, the response of baseflows to precipitation regime changes varies regionally due to the influence of climate conditions.Godsey et al. [16] found that summer baseflows decreased due to the reduction of peak snow water equivalent in California's Sierra Nevada Mountains, where baseflows are primarily sustained by groundwater recharged during snowmelt.Conversely, Jyrkama and Sykes [17] found that warmer winter air temperatures reduced the extent of ground frost and shifted the spring snowmelt event earlier in the year, allowing more water to infiltrate into the ground and thus more groundwater recharge and presumably higher baseflow [18,19].For example, Yang et al. [14] observed significant increases (25%-90%) in cold season (October to April) baseflows sourced from groundwater discharge in Siberia.Like Siberia, the watersheds in the present study are characterized by more precipitation in the summer (rain) than in the winter (snow); however, the study watersheds are located considerably further south in a region that is sensitive to climate warming.The response of winter baseflows to climate warming-induced precipitation regime shifts in this region in northeastern China and at similar latitudes (~50 • N) is still unknown.
Permafrost thaw is thought to alter surface and subsurface hydrology by enhancing infiltration, deepening groundwater flow paths, and increasing groundwater discharge [20].These changes are caused by the difference in the hydraulic properties of frozen and unfrozen ground [21,22].Several recent studies have related observed increases in winter baseflow to permafrost thaw.Walvoord and Striegl [23] analyzed long term (>30 years) streamflow records of the Yukon River Basin, and attributed increases in groundwater contributions (0.7% to 0.9%•year −1 ) to streamflow to climate warming and concomitant permafrost thaw that enhanced infiltration and opened groundwater flowpaths.In eastern Siberia, Brutsaert and Hiyama [24] related changes in baseflow to the rate of the active layer thickening resulting from permafrost thaw, and inferred that the active layer thickness was increasing at average rates of 0.3 to 1 cm•year −1 from 1950 to 2008 in areas with discontinuous permafrost.Although the physical tenability of enhanced groundwater flow rates and increased baseflow due to the thawing of continuous and discontinuous permafrost is supported by numerical modeling results (e.g., [25][26][27]), as well as baseflow recession analyses [28,29], more observational evidence in different regions is required to upscale these results to the pan-Arctic basin and to inform regional water management decisions.In particular, increased attention should be placed in watersheds within the discontinuous permafrost zone, as this region is hydrologically and ecologically sensitive to atmospheric warming [6] because the permafrost temperature is close to 0 • C.
The southeastern margin of the Eurasian Cryolithozone in northeastern China has experienced rapid climate warming over the past five decades, with an accelerated trend after the late 1980s [30].Such climate warming can result in significant permafrost degradation, and changes in the precipitation regimes [31], and permafrost warming or degradation has already been observed in this region [32][33][34].For instance, the active layer thickness increased from 1.2 m in the 1990s to more than 1.8 m in 2004 in the Yituli'he Permafrost Observatory (121 • 20 E, 50 • 56 N) [33].The results of direct observations from boreholes (CK38, CK3, 123 • 11 E, 52 • 50 N) showed that the permafrost thickness decreased from 7 m in 1979 to 2 m in 1991, resulting in an average thinning rate of 0.4 m/year during this period [33].Also, inactive ice wedges that were observed in the Wuma River basin (121 • 45 E, 52 • 58 N) in 1987 [34] had disappeared by 2007 due to recent climate warming [32].Although permafrost degradation is well-documented in this region, the hydrological response to permafrost thaw and precipitation regime shifts in this region has received little attention compared to similar permafrost regions in North America and Siberia.In this study, hydrologic and meteorological data from two large permafrost watersheds in northeastern China are used to (1) assess any significant trends in winter baseflows over the past four decades and (2) determine any correlations and temporal delays between the long-term winter baseflow and permafrost thaw metrics (air temperature and thawing index) and characteristics of the precipitation regime (annual rainfall fraction of precipitation and annual summer precipitation).These analyses will help identify how regional hydrological processes in northeastern China are being influenced by climate warming and permafrost thaw and determine if hydrologic trends observed in northern watersheds in other countries are also manifested in northeastern Asia.Such studies are lacking in China and are required to understand the sensitivity of water resources to future climate conditions and to develop effective water management strategies in other locations that are more densely populated.

Study Area
Two watersheds in the permafrost environment of northeastern China (Figure 1a), the Duobukuer River watershed (DBKE, 3094 km 2 ) and Tahe River watershed (TH, 6581 km 2 ) which both feed the Heilongjiang River, were selected for this investigation.They are located in the southernmost distribution of permafrost that extends from the Arctic region of Eurasia and span discontinuous (defined as 50%-90% permafrost coverage), sporadic (10%-50%), and isolated (0%-10%) permafrost zones [35].The fractions of permafrost distribution in the study watersheds gradually increase from south to north with more discontinuous permafrost in TH and more isolated permafrost in DBKE (Figure 1b).The elevations of the two watersheds range from 155 to 1401 m above sea level, with an average of 722 m.However, most (79%) of the area for both watersheds is concentrated in the elevation range of 500 to 900 m.There are no zones with permanent snow cover in either watershed.The study watersheds are dominated by gentle hills covered by forest with fens in a relatively flat valley bottom [36].According to the long-term forest data observed using the Forestry Standards "Observation Methodology for Long-term Forest Ecosystem Research" of China (LY/T1952-2011) and supplied by the local Forestry Bureau, the two study watersheds are characterized by high forest coverage (~75% over entire study period).The soil in this region is dominated by coniferous forest soils with thickness ranging between 20 and 30 cm with an overlying 0.05 to 0.2 m thick organic layer [37].Below the soil layer is a 1.0 to 3.0 m thick mixture of sandy clay with gravel [38].The thickness of permafrost ranges from 0 to 30 m with a supra-permafrost thawed zone (active layer plus potential taliks) ranging from 1.0 to 3.0 m [39].Although the forest vegetation and organic soil can thermally insulate the ground [38], the permafrost in the study regions has been experiencing rapid thaw in the past several decades [33].
Water 2017, 9, 25 3 of 15 countries are also manifested in northeastern Asia.Such studies are lacking in China and are required to understand the sensitivity of water resources to future climate conditions and to develop effective water management strategies in other locations that are more densely populated.

Study Area
Two watersheds in the permafrost environment of northeastern China (Figure 1a), the Duobukuer River watershed (DBKE, 3094 km 2 ) and Tahe River watershed (TH, 6581 km 2 ) which both feed the Heilongjiang River, were selected for this investigation.They are located in the southernmost distribution of permafrost that extends from the Arctic region of Eurasia and span discontinuous (defined as 50%-90% permafrost coverage), sporadic (10%-50%), and isolated (0%-10%) permafrost zones [35].The fractions of permafrost distribution in the study watersheds gradually increase from south to north with more discontinuous permafrost in TH and more isolated permafrost in DBKE (Figure 1b).The elevations of the two watersheds range from 155 to 1401 m above sea level, with an average of 722 m.However, most (79%) of the area for both watersheds is concentrated in the elevation range of 500 to 900 m.There are no zones with permanent snow cover in either watershed.The study watersheds are dominated by gentle hills covered by forest with fens in a relatively flat valley bottom [36].According to the long-term forest data observed using the Forestry Standards "Observation Methodology for Long-term Forest Ecosystem Research" of China (LY/T1952-2011) and supplied by the local Forestry Bureau, the two study watersheds are characterized by high forest coverage (~75% over entire study period).The soil in this region is dominated by coniferous forest soils with thickness ranging between 20 and 30 cm with an overlying 0.05 to 0.2 m thick organic layer [37].Below the soil layer is a 1.0 to 3.0 m thick mixture of sandy clay with gravel [38].The thickness of permafrost ranges from 0 to 30 m with a supra-permafrost thawed zone (active layer plus potential taliks) ranging from 1.0 to 3.0 m [39].Although the forest vegetation and organic soil can thermally insulate the ground [38], the permafrost in the study regions has been experiencing rapid thaw in the past several decades [33].Both watersheds are characterized by a typical continental monsoon climate with a wet, hot summer and a cold, dry winter.Mean annual air temperature (MAAT) is −1.5 • C. The average maximum air temperature is 22 • C (July), while the average minimum is −29 • C (January).The summer is heavily influenced by the Southeast Monsoon from the Pacific Ocean, which results in a wet season from June to September.The average precipitation in the wet season is 440 mm, accounting for 85% of the annual total.During the winter from November to March, precipitation occurs in the form of snow, which usually contributes to runoff during the snowmelt season (late April to May).Native vegetation in the study area consists of forest communities dominated by Larch (Larix gmelinii), along with some broadleaf species, such as birch (Betula platyphylla) and Mongolian oak (Quercus mongolica).

Data
Daily river discharge data from 1973 to 2012 for DBKE and TH were inferred from stage measurements and rating curves at the Sgongling and Tahe hydrological stations, respectively (Figure 1c).These stations also contained rain and snow gauges.No dams or hydropower stations regulate flow in the study watersheds.All daily discharge data (volumetric flow rate, m 3 /s) were transformed into runoff depth (mm) by dividing by the watershed area, and the runoff depth was then aggregated to obtain winter baseflow.Winter baseflow is herein defined as the cumulative runoff (discharge/area) between November and March.During this period, precipitation is temporarily stored in the snowpack, the rivers are covered by ice, and streamflow is predominantly sourced from groundwater discharge.The changes in winter baseflow cannot be affected by a shift in the snowmelt timing as this event occurs later in the year.
The MAAT and air temperature thawing index (integral of air temperature over time during thawing period, • C•day) were calculated from the daily mean air temperature data from Jiagedaqi and Xinlin climate stations for DBKE and TH watersheds, respectively (Figure 1c).The thawing index is the standard parameter used to drive the Stefan equation which has been used to predict seasonal ground freeze-thaw as well as permafrost thaw in both hydrological and large earth system models [40][41][42].This index has a physical basis, and the Stefan equation is the most common approach used in large-scale permafrost modeling [43].The thawing index ( • C•days) for each year was calculated by summing up the mean daily air temperatures for all of the days that it exceeded 0 • C. The annual precipitation, summer precipitation (June to August), and annual rainfall fraction (depth of annual rain divided by depth of annual precipitation) data for both watersheds were calculated from the daily precipitation data at Songling and Xinlin from rain and snow gauges (Figure 1c).According to the National Meteorological Information Center, which provided these data, these precipitation data are corrected for snow under-catch.Data from these meteorological stations are assumed to be representative of each basin; however, we acknowledge that the elevation difference across both basins is significant.Although the MAAT and precipitation would not be uniform across the site, we expect that the trends in MAAT and precipitation obtained at both stations would be reasonably representative of the trends of the watershed-averaged values for both sites.The similarities of the climate trends observed between the two watersheds (Section 4.1) support this inference.

Trend Analysis
The non-parametric trend test, Mann-Kendall tau [44,45], was employed to detect whether significant trends exist in the long-term winter baseflow, MAAT, thawing index, annual precipitation, summer precipitation, and annual rainfall fraction.The magnitude of the trend, i.e., the slope (change per unit time), was estimated using the Kendall-Theil Robust Line method [46,47].The average change rate (%•year −1 ) of flows was calculated as the slope (mm•year −1 ) divided by the average flows (mm) for the entire period.One approach to visualize long-term trends in hydrologic data is to average out the time series data over constant intervals; herein we use the common interval length of 11 years [48].

Time Series Correlation Analysis
Cross-correlation analysis is an effective approach to detect statistical relationships between time series of hydrological and climatic variables and has the advantage of being able to determine lagged effects [49,50].Before cross-correlation analysis, all hydrological data series and climatic data series were pre-whitened to eliminate autocorrelations existing in the data series by using the best fitting Autoregressive Integrated Moving Average (ARIMA) models [51][52][53].Model residuals were then applied for cross-correlation analysis [54].

Breakpoints Analysis
The Pettitt's test [55] was applied to detect if there were any breakpoints in winter baseflow, MAAT, thawing index, or annual rainfall fraction.The nonparametric change-point test method developed by Pettitt [55] is widely used for detecting breakpoints in long-term hydrometeorological data series [56,57], and details can be found in Pettitt [55] and Zhang et al. [58].The difference in the timing of the breakpoints in the winter baseflow and permafrost thaw metrics yields the phase shift between these variables.

Trends in Long-Term Winter Baseflow, Permafrost Thaw Metrics, and Precipitation
The winter baseflow for each watershed exhibited significant (p < 0.01) positive trends over the full study period , with average increasing rates of 2.5%•year −1 and 1.7%•year −1 and total increases of 98.3% and 66.9% for DBKE and TH, respectively (Table 1).The MAAT and thawing index for both watersheds exhibited similar and significant positive trends over the past 40 years, with MAAT increases of 0.043 • C•year −1 for both watersheds and thawing index increases of 7.1 and 7.2 • C•days•year −1 for DBKE and TH, respectively.Although annual precipitation and summer rain were relatively stable with no significant (p > 0.05) trend during the entire study period, the annual rainfall fraction of DBKE showed a significant (p < 0.05) positive trend of 0.16%•year −1 .There was a strong temporal correlation among winter baseflow, MAAT, thawing index, and annual rainfall fraction (symbols, Figure 2).The 11-year moving windows of winter baseflow, MAAT, thawing index, and annual rainfall fraction (lines, Figure 2) also indicate generally consistent trends among these variables over the entire period, during which there was a positive trend from the beginning of records to the 1990s followed by relatively stable regimes for winter baseflow and MAAT.Trends in monthly precipitation were also detected as shown in Figure 3. Athough no significant (Mann-Kendall tau test, p > 0.05) trends were found in any monthly precipitation series, May, June, July, and October exhibited increasing patterns, while April, August, and September exhibited decreasing patterns.
Water 2017, 9, 25 5 of 15 data is to average out the time series data over constant intervals; herein we use the common interval length of 11 years [48].

Time Series Correlation Analysis
Cross-correlation analysis is an effective approach to detect statistical relationships between time series of hydrological and climatic variables and has the advantage of being able to determine lagged effects [49,50].Before cross-correlation analysis, all hydrological data series and climatic data series were pre-whitened to eliminate autocorrelations existing in the data series by using the best fitting Autoregressive Integrated Moving Average (ARIMA) models [51][52][53].Model residuals were then applied for cross-correlation analysis [54].

Breakpoints Analysis
The Pettitt's test [55] was applied to detect if there were any breakpoints in winter baseflow, MAAT, thawing index, or annual rainfall fraction.The nonparametric change-point test method developed by Pettitt [55] is widely used for detecting breakpoints in long-term hydrometeorological data series [56,57], and details can be found in Pettitt [55] and Zhang et al. [58].The difference in the timing of the breakpoints in the winter baseflow and permafrost thaw metrics yields the phase shift between these variables.

Trends in Long-Term Winter Baseflow, Permafrost Thaw Metrics, and Precipitation
The winter baseflow for each watershed exhibited significant (p < 0.01) positive trends over the full study period , with average increasing rates of 2.5%•year −1 and 1.7%•year −1 and total increases of 98.3% and 66.9% for DBKE and TH, respectively (Table 1).The MAAT and thawing index for both watersheds exhibited similar and significant positive trends over the past 40 years, with MAAT increases of 0.043 °C•year −1 for both watersheds and thawing index increases of 7.1 and 7.2 °C•days•year −1 for DBKE and TH, respectively.Although annual precipitation and summer rain were relatively stable with no significant (p > 0.05) trend during the entire study period, the annual rainfall fraction of DBKE showed a significant (p < 0.05) positive trend of 0.16%•year −1 .There was a strong temporal correlation among winter baseflow, MAAT, thawing index, and annual rainfall fraction (symbols, Figure 2).The 11-year moving windows of winter baseflow, MAAT, thawing index, and annual rainfall fraction (lines, Figure 2) also indicate generally consistent trends among these variables over the entire period, during which there was a positive trend from the beginning of records to the 1990s followed by relatively stable regimes for winter baseflow and MAAT.Trends in monthly precipitation were also detected as shown in Figure 3. Athough no significant (Mann-Kendall tau test, p > 0.05) trends were found in any monthly precipitation series, May, June, July, and October exhibited increasing patterns, while April, August, and September exhibited decreasing patterns.

Relationships between Winter Baseflow, Permafrost Thaw Metrics, and Precipitation
According to the results of the linear regression analyses (Figure 4), the winter baseflows exhibited stronger linear relationships with the permafrost thaw metrics than the precipitation regime characteristics.Significant (p < 0.05) linear relationships were detected between winter baseflows and MAAT or thawing index for both watersheds (Figure 4a-d).In the past 40 years, the winter baseflows in the two watersheds increased by 0.83 (TH) and 1.80 (DBKE) mm per • C of MAAT, and 0.72 and 1.46 mm per 100 • C•days of thawing index.For the precipitation characteristics, a significant (p < 0.05) relationship was found between winter baseflow and annual rainfall fraction for DBKE (Figure 4i), but no significant linear relationships were found among any of the other precipitation metrics and winter baseflow (Figure 4e-h,j).

Relationships between Winter Baseflow, Permafrost Thaw Metrics, and Precipitation
According to the results of the linear regression analyses (Figure 4), the winter baseflows exhibited stronger linear relationships with the permafrost thaw metrics than the precipitation regime characteristics.Significant (p < 0.05) linear relationships were detected between winter baseflows and MAAT or thawing index for both watersheds (Figure 4a-d).In the past 40 years, the winter baseflows in the two watersheds increased by 0.83 (TH) and 1.80 (DBKE) mm per °C of MAAT, and 0.72 and 1.46 mm per 100 °C•days of thawing index.For the precipitation characteristics, a significant (p < 0.05) relationship was found between winter baseflow and annual rainfall fraction for DBKE (Figure 4i), but no significant linear relationships were found among any of the other precipitation metrics and winter baseflow (Figure 4e-h,j).The time series cross-correlation analyses revealed that all the winter baseflows for both of the study watersheds exhibited significant and positive correlations (p < 0.05) with MAAT and thawing index over the entire study period (Table 2).Both watersheds exhibited the same lag of two years between winter baseflow and MAAT, and the same lag of four years between winter baseflow and the thawing index.The winter baseflow of DBKE also showed a significantly positive correlation (p < 0.05) with annual rainfall fraction, but no significant correlation (p > 0.05) was found between winter baseflow and annual rainfall fraction for TH (Table 2).Cross-correlation analyses were also conducted for the winter baseflows with the annul precipitation and summer rainfall.However, no significant correlations were found between them, and thus the results are not shown in Table 2.

Phase Shifts in Winter Baseflow
There were obvious phase shifts in the winter baseflows compared to the MAAT, thawing index, and annual rainfall fraction over the entire period for both watersheds.Table 3 demonstrates that the winter baseflows, MAAT, and thawing index series exhibited significant (p < 0.01) breakpoints in 1991, 1989, and 1988, respectively (Table 3).The annual rainfall fraction only showed a significant (p < 0.05) breakpoint in 1986 for DBKE.According to these results, there was a lag of two years between winter baseflows and MAAT for both watersheds, which is consistent with the results of the cross-correlation analyses for the winter baseflow and MAAT.There was a lag of three years between thawing index and winter baseflows.The statistically significant differences (p < 0.01) in winter baseflow, MAAT, and thawing index for both watersheds and annual rainfall fraction for DBKE between the periods before and after breakpoints were further confirmed by the Kolmogorov-Smirnov Z-test (Table 3).The time series cross-correlation analyses revealed that all the winter baseflows for both of the study watersheds exhibited significant and positive correlations (p < 0.05) with MAAT and thawing index over the entire study period (Table 2).Both watersheds exhibited the same lag of two years between winter baseflow and MAAT, and the same lag of four years between winter baseflow and the thawing index.The winter baseflow of DBKE also showed a significantly positive correlation (p < 0.05) with annual rainfall fraction, but no significant correlation (p > 0.05) was found between winter baseflow and annual rainfall fraction for TH (Table 2).Cross-correlation analyses were also conducted for the winter baseflows with the annul precipitation and summer rainfall.However, no significant correlations were found between them, and thus the results are not shown in Table 2. Notes: 1 Bold italics values indicate the statistical significance of p < 0.05.The first, second, third, and fourth model columns are the Autoregressive Integrated Moving Average (ARIMA) models to pre-whiten time series data of flows [44], MAAT, TI, and rain fraction for each watershed, respectively.

Phase Shifts in Winter Baseflow
There were obvious phase shifts in the winter baseflows compared to the MAAT, thawing index, and annual rainfall fraction over the entire period for both watersheds.Table 3 demonstrates that the winter baseflows, MAAT, and thawing index series exhibited significant (p < 0.01) breakpoints in 1991, 1989, and 1988, respectively (Table 3).The annual rainfall fraction only showed a significant (p < 0.05) breakpoint in 1986 for DBKE.According to these results, there was a lag of two years between winter baseflows and MAAT for both watersheds, which is consistent with the results of the cross-correlation analyses for the winter baseflow and MAAT.There was a lag of three years between thawing index and winter baseflows.The statistically significant differences (p < 0.01) in winter baseflow, MAAT, and thawing index for both watersheds and annual rainfall fraction for DBKE between the periods before and after breakpoints were further confirmed by the Kolmogorov-Smirnov Z-test (Table 3).Based on the analyses above, the trend in winter baseflow was closely related to MAAT and the thawing index.To quantify changes between time periods, the entire study period was divided into two sub-periods based on the breakpoints: 1973 to 1989 (warming period) and 1990 to 2012 (stable period).Similar increases in average MAAT (0.95 and 0.79 • C) and the thawing index (164 and 172 • C days) were found between the warming period and the stable period for both DBKE and TH (Figure 5b,c), but there was a higher increase in the average annual rainfall fraction in DBKE (4.2%) than in TH (2.3%) between these two periods (Figure 5d).The average winter baseflow for DBKE also exhibited a higher increase (101%) than that for TH (65%) (Figure 5a), and these findings are consistent with the results of the trend analyses which determined the total changes over the period of record to be 67% (TH) and 98% (DBKE) (Table 1).Based on the analyses above, the trend in winter baseflow was closely related to MAAT and the thawing index.To quantify changes between time periods, the entire study period was divided into two sub-periods based on the breakpoints: 1973 to 1989 (warming period) and 1990 to 2012 (stable period).Similar increases in average MAAT (0.95 and 0.79 °C) and the thawing index (164 and 172 °C days) were found between the warming period and the stable period for both DBKE and TH (Figure 5b,c), but there was a higher increase in the average annual rainfall fraction in DBKE (4.2%) than in TH (2.3%) between these two periods (Figure 5d).The average winter baseflow for DBKE also exhibited a higher increase (101%) than that for TH (65%) (Figure 5a), and these findings are consistent with the results of the trend analyses which determined the total changes over the period of record to be 67% (TH) and 98% (DBKE) (Table 1).

The Impacts of Climate Warming-Induced Permafrost Thaw on Winter Baseflow
Northeastern China has experienced a pronounced warming trend over the past half a century [59], which has resulted in rapid warming and thawing of permafrost as shown by active layer thickening, thinning permafrost, increasing ground temperatures, expanding taliks, and disappearing permafrost islands in northeastern China [33].Such changes have specifically been documented in the region of the study watersheds.For instance, the active layer thickness in the Amu'er area (90 km northwest from TH) increased by about 50 cm from the 1970s to the 1990s [60,61].Also, ground temperature at 13 m depth in an undisturbed meadow on the first terrace of the northern bank of the Yituli'he River (180 km southwest of the study watersheds) increased by about

The Impacts of Climate Warming-Induced Permafrost Thaw on Winter Baseflow
Northeastern China has experienced a pronounced warming trend over the past half a century [59], which has resulted in rapid warming and thawing of permafrost as shown by active layer thickening, thinning permafrost, increasing ground temperatures, expanding taliks, and disappearing permafrost islands in northeastern China [33].Such changes have specifically been documented in the region of the study watersheds.For instance, the active layer thickness in the Amu'er area (90 km northwest from TH) increased by about 50 cm from the 1970s to the 1990s [60,61].Also, ground temperature at 13 m depth in an undisturbed meadow on the first terrace of the northern bank of the Yituli'he River (180 km southwest of the study watersheds) increased by about 0.28 • C from 1984 to 1997 [62].Jin et al. [33] found that the observed permafrost degradation in this region can be mainly attributed to the persistent warming that occurred from 1970 to 2000.Other researchers have demonstrated a trend between hydrologic trends and ground thermal regimes in Eurasia.Smith et al. [63] noted increasing patterns of minimum daily flows throughout Russia, and these were partly attributed to the reduced influence of seasonal ground freezing.Increased winter flow in the Aldan subbasin of the Lena Basin in Siberia was thought to be related to permafrost thaw [64].However, the hydrologic response to permafrost thaw in the southeastern margin of the Eurasian Cryolithozone (i.e., the study region) has not been previously considered.
Our results indicate that there was a significant increase in MAAT of 0.043 • C•year −1 and the thawing index of approximately 7 • C•days•year −1 in the study watersheds over the past 40 years.Based on the observed thawing and warming of permafrost in the study region and the well-known linkage between that and winter baseflow, we postulate that the significant positive trends in winter baseflows are at least partly related to permafrost thaw and increased groundwater discharge.These results are also supported by the strong significant positive correlations (p < 0.05) between winter baseflow and both the MAAT and the thawing index in the study watersheds.Winter baseflow is a standard indicator of groundwater contributions to streamflow in watersheds where winter precipitation is stored in the snowpack [65].Thus, the positive trends in winter baseflow indicate increased winter groundwater discharge, which likely was caused by enhanced groundwater recharge and groundwater storage due to ground thaw [23].In general, permafrost thaw is expected to increase surface and subsurface hydrologic connectivity of watersheds and thereby lead to enhanced interactions between surface water bodies and aquifers [66].It should be noted that the analyses presented herein do not disentangle the hydrologic changes due to permafrost thaw from those caused by seasonal frost changes, and the latter warrants further investigation.
Our results also indicate that there are higher increasing rates of winter baseflow in the study watersheds (1.7%•year −1 and 2.5%•year −1 ) than those rates (ranging from 0.4%•year −1 to 0.9%•year −1 ) observed in permafrost watersheds with similar average winter baseflow (>300 m 3 •s −1 ) from the Canadian Northwest Territories [20] and the Yukon River Basin in Alaska [23].These differences likely arise because (1) the watersheds in the present study span discontinuous, sporadic, and isolated permafrost systems that are more sensitive to climate warming than continuous permafrost environments [6,66] and (2) the climate warming rate in northeastern China (0.32 • C•decade −1 from 1954 to 2005) [67] is higher than that in Canadian Northwest Territories (with the greatest temperature change occurring during the winter with a positive trend of 0.23 • C•decade −1 from 1950 to 2010) [68] and in Alaska (with the highest increases in winter in the interior region of 0.22 • C•decade −1 from 1949 to 1998) [69].The spatial distribution of our results also points to the importance of permafrost thaw as a driver of hydrologic change.The southern DBKE watershed, which contains the most climate-sensitive sporadic and isolated permafrost zones (Figure 1b), was characterized by more pronounced changes in winter baseflow (0.16 mm•year −1 or 2.5%•year −1 ) over the period of record, while the more northern watershed (TH) only exhibited increasing rates of 0.09 mm•year −1 (1.7%•year −1 ).

Impact of Climate Warming-Induced Precipitation Regime Changes on Winter Baseflow
Temperature-induced shifts of precipitation from snow towards rain and earlier snowmelt are two of the most profound and widely anticipated changes in the hydrological cycle in response to climate warming [70,71].In the snow-dominated catchments of California's Sierra Nevada mountains, Godsey et al. [16] found that for every 10% decrease in peak snow water equivalent, annual dry season minimum flows decreased from 9% to 22% due to the reduction in groundwater recharged from snowmelt.However, the study watersheds are rain-dominated (the average rain fraction of precipitation is 87.2% over the entire period of record).Thus, the surface and subsurface water are mainly sourced by the rainfall in the wet season (summer) [72], during which the groundwater is recharged from both rainfall and surface water.Although the amount of snowpack decreased due to the increased rain fraction and earlier snowmelt, the duration and extent of infiltration increased, allowing for more groundwater recharge [17], which can then sustain winter baseflows [18,19].
Our results showed that, although no significant (p > 0.05) increase in summer precipitation was found, there were positive trends in monthly precipitation in May, June, July, and October for DBKE (Figure 3).Also, there was a significant (p < 0.05) positive trend (0.16%•year −1 ) for rain fraction over the period of record and a total increase of 4.2% in rain fraction between the warming period (1973)(1974)(1975)(1976)(1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989) and the stable period .Several studies have demonstrated that there is a positive correlation between annual rainfall and groundwater recharge (e.g., [73,74]).Thus, as demonstrated by others [18,19], the increased rainfall during the spring and summer seasons in our study watersheds likely resulted in enhanced infiltration and annual groundwater recharge that in turn sustained winter baseflow.Thus, this precipitation regime shift is another possible source for the increased winter baseflows in the study watersheds.This hypothesis is supported by the fact that DBKE exhibited a higher total increase (98%) in winter baseflow than TH (67%), and DBKE also exhibited a significant positive trend (p < 0.05) in annual rainfall fraction and a high increasing rate of wet season precipitation, while TH did not (Table 1 and Figure 3).

Lagged Effects of Winter Baseflow in Response to Climate Warming
The surface hydrological response is naturally phase shifted from climate warming because permafrost is not in equilibrium with the climate due to the thermal inertia associated with pore water phase change [9].Some modeling studies have also demonstrated a phase shift between base flow response and permafrost thaw, but the results vary depending on the hydrogeologic system in question.For instance, Bense et al. [75] found that the process of uptake of water into aquifer storage following an increase in sub-permafrost hydraulic heads delays base flow increases by decades to centuries.However, the modeling results of Ge et al. [26] indicated changes to the maximum shallow supra-permafrost groundwater discharge lagged the maximum subsurface temperature by less than one year.
In this study, two statistical methods were applied to identify any breakpoints in the observed long-term winter baseflow, MAAT, and thawing index data.The breakpoint analysis indicated that the breakpoint for winter baseflow occurred in 1991, which was two years after the MAAT breakpoint.These results are consistent with those of the cross-correlation analysis, which indicated that the increasing winter baseflow was lagged relative to the MAAT and thawing index by two and four years, respectively.The observed lags between MAAT or thawing index and winter baseflow tentatively support the inference that permafrost thaw is a driver of hydrologic change in these study watersheds.To our knowledge, this is the first field study examining the phase shifts between river discharge conditions in permafrost regions and MAAT or thawing index.The relatively rapid response of winter baseflow to the increasing MAAT likely indicates the occurrence of shallow permafrost thaw and demonstrates that groundwater and surface water exchanges in the southeastern marginal area of the Eurasian Cryolithozone are sensitive to increased air temperature.A recent study in the Lena River basin of Northern Eurasia showed that the daily minimum flows also exhibited a significant positive trend with an abrupt change in the 1990s and early 2000s [76].The timing of these changes are lagged relative to those calculated in this study, which we attribute to the fact that the present study is focused on more southern watersheds with predominantly discontinuous permafrost.Given the expected northward shifting of the southern permafrost limit in the context of future climate warming, more pronounced hydrologic regime shifts may occur in environments that are presently characterized as continuous permafrost.Such changes have implications for water resources distribution and availability in northeastern China and across the pan-Arctic basin, and suggest that flexible infrastructure and management policies should be implemented to accommodate future hydrologic change.

Conclusions
Climate warming-induced permafrost thaw and precipitation regime shifts are altering hydrological processes at high latitudes and elevations.Our study presents historical long-term (1973 to 2012) data and demonstrates that winter baseflows exhibited a significant positive trend of 2.5% and 1.7%•year −1 over the past 40 years in two watersheds in northeastern China.Based on the statistical analyses between the climate data, permafrost thaw metrics, and hydrologic data, the increases in winter baseflow are postulated to be related to permafrost thaw and likely also partly due to the increased rainfall in the wet season during the study period.Furthermore, there was a breakpoint in the winter baseflow in 1989, which lagged the mean annual air temperature breakpoint by only two years, suggesting that the permafrost and hydrologic regimes in these environments respond quickly to increasing air temperature.These findings emphasize the dynamic nature of future hydrologic processes and water resources management in northeastern China and other similar cold regions around the world.

Figure 2 .Figure 3 .
Figure 2. Annual and 11-year moving windows of (a) winter baseflow; (b) mean annual air temperature (MAAT); (c) thawing index (TI); (d) annual precipitation; (e) annual rainfall fraction; and (f) the summer rainfall in DBKE (blue) and TH (red).For the 11-year moving windows, the year indicated is the median year (i.e., 1990 represents the period from 1985 to 2005).

Figure 3 .
Figure 3. Trends of monthly precipitation in DBKE and TH.

Figure 4 .
Figure 4. Winter baseflow vs. the MAAT (panels a and b), thawing index (panels c and d), annual precipitation (panels e and f), annual summer rainfall (panels g and h), and annual rainfall fraction (panels i and j) for the DBKE and TH watersheds.S represents the slope (the magnitude of the trend).

Figure 5 .
Figure 5. Box plots for the (a) winter baseflow; (b) MAAT; (c) thawing index; and (d) annual rainfall fraction in the warming period and stable period for both study watersheds.

Figure 5 .
Figure 5. Box plots for the (a) winter baseflow; (b) MAAT; (c) thawing index; and (d) annual rainfall fraction in the warming period and stable period for both study watersheds.

Table 1 .
Trend analyses resultsfor winter baseflow, mean annual air temperature (MAAT), thawing index, and precipitation metrics.Bold italics and bold italics underlined values indicate trends with statistical significance of p < 0.05 and 0.01, respectively; 2 TC and AC represent total change over the period of record and average change per year, respectively.

Table 1 .
Trend analyses resultsfor winter baseflow, mean annual air temperature (MAAT), thawing index, and precipitation metrics.Bold italics and bold italics underlined values indicate trends with statistical significance of p < 0.05 and 0.01, respectively; 2 TC and AC represent total change over the period of record and average change per year, respectively.
Notes:1Bold italics and bold italics underlined values indicate trends with statistical significance of p < 0.05 and 0.01, respectively; 2 TC and AC represent total change over the period of record and average change per year, respectively.

Table 2 .
Cross-correlations between winter baseflow and climatic variables in the different study watersheds and lags between the climatic variable MAAT or thawing index and the winter baseflow.
[44]s:1Bold italics values indicate the statistical significance of p < 0.05.The first, second, third, and fourth model columns are the Autoregressive Integrated Moving Average (ARIMA) models to pre-whiten time series data of flows[44], MAAT, TI, and rain fraction for each watershed, respectively.

Table 2 .
Cross-correlations between winter baseflow and climatic variables in the different study watersheds and lags between the climatic variable MAAT or thawing index and the winter baseflow.

Table 3 .
Results of breakpoint analysis from Pettitt's test and the Kolmogorov-Smirnov Z-test.

Table 3 .
Results of breakpoint analysis from Pettitt's test and the Kolmogorov-Smirnov Z-test.