Nonlinear Response of Streamflow to Climate Change in High-Latitude Regions : A Case Study in Headwaters of Nenjiang River Basin in China ’ s Far Northeast

Assessment of the response of streamflow to future climate change in headwater areas is of a particular importance for sustainable water resources management in a large river basin. In this study, we investigated multiscale variation in hydroclimatic variables including streamflow, temperature, precipitation, and evapotranspiration in the Headwater Areas of the Nenjiang River Basin (HANR) in China’s far northeast, which are sensitive to climate change. We analyzed 50-year-long (1961–2010) records of the hydroclimatic variables using the ensemble empirical mode decomposition (EEMD) method to identify their inherent changing patterns and trends at the inter-annual and inter-decadal scales. We found that all these hydroclimatic variables showed a clear nonlinear process. At the inter-annual and inter-decadal scales, streamflow had a similar periodic changing pattern and transition years to that of precipitation; however, within a period, streamflow showed a close association with temperature and evapotranspiration. The findings indicate that the response of streamflow in headwater regions to climate change is a nonlinear dynamic process dictated by precipitation at the decadal scale and modified by temperature and evapotranspiration within a decade.


Introduction
Global change is expected to alter hydrological processes, redistributing water resources in time and space as the frequency and intensity of climatic extremes increase in response to a warmer world [1,2].The change in regional water resources, particularly streamflow, will affect local agriculture, environmental conditions, and socio-economics.Headwater streams are the smallest parts of river and stream networks, but they have the longest cumulative length within a river basin [3,4] and therefore govern hydrologic connectivity and supplying and preserving water sources [5].For instance, the headwaters of the Colorado River Basin produced 75% of the annual streamflow from 25% of the area, while less than 10% of the annual streamflow was contributed by the lower basin [6].These characteristics of a drainage network make headwater regions more vulnerable to global climate change in a river basin [7,8].Despite headwaters spatial dominance and importance, they are often not well considered in management policies [9], especially in high-latitude regions.Hence, understanding of streamflow response to future climate change in headwaters is urgently required.
There have been an increasing number of studies on headwater streamflow to climate change [3,4,9,10].These studies employed methods designed under the hypothesis that the time series of hydroclimatic variables are stationary and therefore can be considered as a linear system (such as moving average or polynomial, linear regression, empirical or spline function fitting) [11].However, studies also found that hydroclimatic processes are complex and hydroclimatic variables, including runoff, precipitation, temperature and evapotranspiration, often behave nonlinearly and nonstationarily [12,13].Therefore, traditional (i.e., linear interpretation) methods cannot satisfy for capturing nonlinear, nonstationary response patterns of streamflow to climatic variability.Wu and Huang [14] developed a time-series signal processing method, ensemble empirical mode decomposition (EEMD), that has a stronger self-adaptability and local variation characteristics based on the signal [13,15,16].The method can gradually separate the oscillations at different scales in hydroclimatic change research (Intrinsic Mode Function-IMF) and the trend components from the original signal.Although the method has some drawbacks (for instance, the residual noise and the mode mixing problem), it has been found to be a robust approach suitable for non-stationary and nonlinear signal detection [17][18][19].In recent years, there has been a growing application of the EEMD method in the field of hydroclimatic variables change, and some meaningful results have been achieved.In this study, we used the method to investigate nonlinear responses of streamflow in a headwater region of China's far northeast to climate change in the past half century.
As mentioned above, headwater regions play a crucial role in quantity and variability of streamflow for downstream areas [20] because headwater streams dominate the cumulative drainage length (60-80%) of a typical river basin [3,4].This is especially true for the headwater region of the Nenjiang River Basin in China's far northeast that is dominated by snow hydrology.The lower Nenjiang River Basin is part of China's largest commercial grain production base as well as the country's largest inland wetland region (including three Wetlands of International Importance and homed many endangered wildlife) [21].It has been reported that, during the past 50 years, the annual air temperature increased by 0.38 • C per decade across Northeast China, while the annual precipitation showed a decreasing trend, especially in the summer and autumn [22].Concurrently, streamflow in the Nenjiang River Basin showed a remarkable decline, especially during the spring and summer months [23].Several studies have looked at the temporal variability of streamflow in the river basin under warming and drying climate conditions [21,[24][25][26][27].These studies have found that climate change was the dominant factor for the reduction in the runoff in the river basin, especially in the upper basin region.However, all these studies use a linear approach to interpret the role of climate conditions on streamflow, giving little consideration to nonlinear and nonstationary characteristics of hydroclimatic variables.
With this mind, we employed a non-linear approach in this study to understand the complexity of streamflow responses to temperature, precipitation and evapotranspiration changes at multiple time scales.The study utilized observational data from seven hydrological stations and 15 meteorological stations across the Nenjiang River Basin to pursue the following objectives: (1) investigate nonlinear nature of streamflow and climate variability in the basin's headwater region by analyzing inter-annual and inter-decadal trends using the ensemble empirical mode decomposition approach; and (2) determine nonlinear relationships between streamflow and climate.The ultimate goal of this study was to understand the response of headwater streamflow in a high-latitude region to variable changes in temperature, precipitation and the corresponding evapotranspiration.

Study Area
The Nenjiang River is the largest tributary of the Songhua River in Northeast China, and drains a total land area of 297,000 km 2 , geographically covering the area of latitude 43 • 36 -51 • 11 N and longitude 118 • 25 -128 • 81 E (Figure 1).The river originates from the Greater Khingan Mountains, traveling southward through a semiarid region, the western Songnen Plain, and entering the Songhua River.The river flow has been the most important water source for agricultural irrigation and the large area of wetlands in the region that support diverse ecosystems and many endangered wildlife species [21].The headwater region of the Nenjiang River, composing of the upper basin of several tributaries, has a drainage area of 55,702 km 2 , accounting for 26.03% of the entire river basin.The average annual streamflow (or water yield) is 78.26 × 10 8 m 3 (computed from 1951 to 2010) that it contributes to 36.53% of total flow at the last gauge station before the Nenjiang River's confluence with the Songhua River, the Dalai Station (G7 in Figure 1).In recent years, multipurpose water resource management is gradually being carried out to make the Songhua River Basin become China's future major grain production base and wetland reserve in the downstream region of Nenjiang River Basin.The streamflow in the basin plays an important role in maintaining ecology security of wetlands and feeding commodity grain production bases in downstream basins.
Water 2018, 10, x FOR PEER REVIEW 3 of 17 traveling southward through a semiarid region, the western Songnen Plain, and entering the Songhua River.The river flow has been the most important water source for agricultural irrigation and the large area of wetlands in the region that support diverse ecosystems and many endangered wildlife species [21].The headwater region of the Nenjiang River, composing of the upper basin of several tributaries, has a drainage area of 55,702 km 2 , accounting for 26.03% of the entire river basin.
The average annual streamflow (or water yield) is 78.26 × 10 8 m 3 (computed from 1951 to 2010) that it contributes to 36.53% of total flow at the last gauge station before the Nenjiang River's confluence with the Songhua River, the Dalai Station (G7 in Figure 1).In recent years, multipurpose water resource management is gradually being carried out to make the Songhua River Basin become China's future major grain production base and wetland reserve in the downstream region of Nenjiang River Basin.The streamflow in the basin plays an important role in maintaining ecology security of wetlands and feeding commodity grain production bases in downstream basins.The headwater region has a cold climate with long and chilly winter as well as widespread snow and seasonal frozen soil and permafrost soil.With rolling hills and deep river valleys, its topography changes greatly and the vertical drop of the headwater is high.Average monthly temperatures range from −21.46 °C in January to 20.38 °C in July with an annual mean of 1.52 °C.Main land covers are grassland and forest.Moreover, human activity is limited ahgnd so it has little impact on local streamflow.The eco-environment of HANR remains in a semi-natural or pristine state, so it is an ideal place to study the response of streamflow to climate change.

Hydro-Meteorological Data
Monthly streamflow records from 1961 to 2010 were collected from six river gauge stations: Shihuiyao (G1), Beian (G2), Xiaoergou (G3), Zhanlantun (G4), Suollun (G5) and Tuliemaodu (G6).They are located at the ecotone of the Khingan Mountains and Songnen Plain, representing a natural state of streamflow (i.e., with minimal anthropogenic activities) over the HANR (Figure 1 and Table 1).Daily meteorological data (including precipitation, air temperature, relative humid, sunshine duration and wind speed) covering fifteen weather stations and spanning 1961-2010 were obtained The headwater region has a cold climate with long and chilly winter as well as widespread snow and seasonal frozen soil and permafrost soil.With rolling hills and deep river valleys, its topography changes greatly and the vertical drop of the headwater is high.Average monthly temperatures range from −21.46 • C in January to 20.38 • C in July with an annual mean of 1.52 • C. Main land covers are grassland and forest.Moreover, human activity is limited ahgnd so it has little impact on local streamflow.The eco-environment of HANR remains in a semi-natural or pristine state, so it is an ideal place to study the response of streamflow to climate change.

Hydro-Meteorological Data
Monthly streamflow records from 1961 to 2010 were collected from six river gauge stations: Shihuiyao (G1), Beian (G2), Xiaoergou (G3), Zhanlantun (G4), Suollun (G5) and Tuliemaodu (G6).They are located at the ecotone of the Khingan Mountains and Songnen Plain, representing a natural state of streamflow (i.e., with minimal anthropogenic activities) over the HANR (Figure 1 and Table 1).Daily meteorological data (including precipitation, air temperature, relative humid, sunshine duration and wind speed) covering fifteen weather stations and spanning 1961-2010 were obtained from the Climate Data Center, China Meteorological Administration.The data quality for these forcing has been controlled by China Meteorological Administration and has been proven to be unmistakable by means of the extreme value examination and time synchronization detection.Few stations have missing data which were filled using linear regression method.Spatial distribution of the gauge/meteorological stations and their subordination to which sub-basins are shown in Figure 1 and Table 1, respectively.Temperature and precipitation are the main driving forces to the variant streamflow in Northeast China [23,27], especially facing a warmer trend with less precipitation since 1980s [21,28].Evapotranspiration (ET 0 ) also plays an important role in controlling the water budget over those water-deficit regions [29].The variability of climatic factors is a nonlinear and complex process.Therefore, nonlinear response of streamflow to temperature, precipitation and potential evapotranspiration were analyzed.The Penman-Monteith (PM) equation recommend by the United Nations Food and Agriculture Organization [30] was used to calculate daily potential evapotranspiration (ET 0 ).

Ensemble Empirical Mode Decomposition
We applied the ensemble empirical mode decomposition (EEMD) method [14,31] to decompose the streamflow, temperature, precipitation and evapotranspiration during the period of 1961-2010 into different time-scales and its nonlinear trend.The EEMD method was evolved from the EMD method to alleviate or reduces the scale-mixing problem.In the EEMD method, a time series x(t) is decomposed in terms of adaptively obtained, amplitude-frequency modulated oscillatory components C i (i = 1, 2, ..., n) and a residual R n , a curve that is either monotonic or contains only one extremum from which no additional oscillatory components can be extracted from: where the ith IMF C i can be written in polar coordinates as , where r(i) is the ith time-dependent amplitude, C i is the ith time-dependent frequency, and R is the residual [31].An IMF is different from Fourier modes for which both r i and C j are time independent, and is defined by the following two properties: (1) each IMF C j has exactly one zero crossing between two consecutive local extremes (i.e., a sequence of maxima and minima); and (2) the local mean of each IMF C j is zero.R n (RES) follows a priori shape and changes with the time after the intrinsic variability at various time scales is removed.This trend also has low sensitivity to the extension (addition) of new data.The latter property guarantees that the physical interpretation at specified time intervals does not change with the addition of new data, which is to say, to be consistent with a physical constraint that the subsequent evolution of a physical system cannot alter the reality that has already happened [13,14,31].Here, the ensemble number is set as 100 and the added noise has an amplitude that is 0.2 times the standard deviation of the corresponding data.

Application of EEMD to Hydro-Meteorological Data
First, we analyzed annual anomaly series to show nonlinear characteristics of streamflow in six gauge stations (Figure 2).Then, we decomposed annual streamflow anomaly of six gauge stations into four Intrinsic Mode Function components (IMF1, IMF2, IMF3 and IMF4) and one trend component (RES), which can adequately reflect multiscale nonlinear variation characteristics of original streamflow series (Figure 3).Each IMF components reflect the fluctuation of streamflow from high frequency (shorter than 10-year period) to low frequency (longer than 10-year period) at different time scales, and whereas the RES component indicates the general change tendency associated with streamflow throughout the whole time period.We obtained the mean period of each IMF by counting the number of peaks of time series curves [14].Next, we investigated multiscale variations of streamflow by reconstructing inter-annual and inter-decadal series (Figure 4 and Table 2).The inter-annual runoff was obtained by IMF1 and IMF2, which represent the inter-annual runoff variation with trend component, while the inter-decadal runoff was obtained by IMF3 and IMF4, which are representatives of the inter-decadal runoff variation plus trend component.Finally, from the view of fluctuation, in-phase or anti-phase and variety of scales or periodic oscillations, we comparatively analyzed reconstructed inter-annual, inter-decadal series and trend components of streamflow, temperature, precipitation and potential evapotranspiration to identify the multiscale response of streamflow to climate change at different timescales in six headwater areas of HANR (Figures 5-7 and Table 3).The Augmented Dickey-Fuller (ADF) unit root test method [32] was used to analyze the stationarity of annual streamflow series for six gauge stations.
Water 2018, 10, x FOR PEER REVIEW 5 of 17

Application of EEMD to Hydro-Meteorological Data
First, we analyzed annual anomaly series to show nonlinear characteristics of streamflow in six gauge stations (Figure 2).Then, we decomposed annual streamflow anomaly of six gauge stations into four Intrinsic Mode Function components (IMF1, IMF2, IMF3 and IMF4) and one trend component (RES), which can adequately reflect multiscale nonlinear variation characteristics of original streamflow series (Figure 3).Each IMF components reflect the fluctuation of streamflow from high frequency (shorter than 10-year period) to low frequency (longer than 10-year period) at different time scales, and whereas the RES component indicates the general change tendency associated with streamflow throughout the whole time period.We obtained the mean period of each IMF by counting the number of peaks of time series curves [14].Next, we investigated multiscale variations of streamflow by reconstructing inter-annual and inter-decadal series (Figure 4 and Table 2).The inter-annual runoff was obtained by IMF1 and IMF2, which represent the inter-annual runoff variation with trend component, while the inter-decadal runoff was obtained by IMF3 and IMF4, which are representatives of the inter-decadal runoff variation plus trend component.Finally, from the view of fluctuation, in-phase or anti-phase and variety of scales or periodic oscillations, we comparatively analyzed reconstructed inter-annual, inter-decadal series and trend components of streamflow, temperature, precipitation and potential evapotranspiration to identify the multiscale response of streamflow to climate change at different timescales in six headwater areas of HANR (Figures 5-7 and Table 3).The Augmented Dickey-Fuller (ADF) unit root test method [32] was used to analyze the stationarity of annual streamflow series for six gauge stations.

Nonlinear Variations of the Streamflow in HANR
All six headwater areas showed a decreasing trend in annual streamflow from 1961 to 2010 (Figure 2), with varying magnitude, timing and fluctuation.The Beian (Figure 2a) and Shihuiyao stations (Figure 2b) had relatively low streamflow in two periods (1970-1985 and 1995-2010): Beian, starting from the late 1960s to mid-1980s; and Shihuiyao, from the mid-1990s to 2000s.Relatively high water periods occurred in the early or mid-1960s and mid-1980s until early 1990.However, the streamflow at the Xiaoergou (Figure 2c), Zhalantun (Figure 2d) and Suolun (Figure 2e) stations were relatively high during 1983-1994.The streamflow at Tuliemaodu (Figure 2d) was comparably high during 1969-1970 and 1986-1994, but was low during 1970-1984 and 1995-2010.In addition, 1998 witnessed an unusually high flow, causing flood stage at the Xiaoergou, Zhalantun, Suolun and Tuliemaodu stations.The five-year moving average and transformation between and low water periods indicate that the streamflow anomaly varied with time nonlinearly and unstably.The Augmented Dickey-Fuller (ADF) unit root tests also suggest that annual streamflow of six headwater areas is non-stationary for all six stations (not shown here).
By applying the EEMD decomposition method, the annual streamflow anomaly at six gauge stations were decomposed into four IMF components and one trend component (denoted by "IMF1-4" and "RES", respectively) in HANR for 1961-2010 (Figure 3).As is shown in Figure 3, each IMF component reflects frequency and amplitudes on streamflow variability: IMF1 has the highest frequency and largest amplitudes than the IMF components, which allows it to present inter-annual amplitudes (at two-or three-year time intervals) of the streamflow in HANR.When moving from IMF1 to IMF4, the frequency and amplitude degrade whereas the periodic increase obviously.RES plots deliver approximate nonlinear trends among six gauge stations during 1961-2010.RES showed a decreasing trend at the earlier stage while an increasing trend at the later stage at the Beian.However, an opposite trend was found at other stations (Figure 3b-f).
Figure 4 demonstrated the reconstructed IMF1 + IMF2 (inter-annual variations) and IMF3 + IMF4 (inter-decadal variations) sequences plus RES (overall trend) component decomposed by EEMD method from original streamflow anomalies at six gauge stations.The reconstructed inter-annual sequences, basically coinciding with the variation trend of original streamflow anomaly series, can portray the fluctuations of the original runoff anomaly series in the study period, which illustrate that the inter-annual oscillations play an important role in the overall streamflow variation of six sub-basins in HANR.Beian and Shihuiyao, to some extent, shared similar inter-decadal patterns during the high-flow periods (1960-1971 and 1983-1995) as well as during the low-flow periods (1975-1980 and 2000-2010).However, the amplitude of the streamflow at Shihuiyao on decadal scales were much larger than Beian's.For the four stations, the amplitudes of the reconstructed inter-decadal signals were flatter before of extreme hydrological events happened in later period (Figure 4c-f).The streamflow in HANR's six sub-basins has a 3-4-year period.
The RES component reflects nonlinear trend at six gauge stations between 1961 and 2010: strong decreasing-slight rising in Beian, strong rising-slight decreasing in Shihuiyao, slight rising-strong decreasing in Xiaoergou and Zhalantun, slight rising-strong decreasing in Suolun and slight rising-strong decreasing in Tuliemaodu (Table 3).Table 2 and RES plot in Figure 3 also reveal a transition year for each sub-basin, i.e., 1990 at Beian, 1991 in Shihuiyao, 1982 at Xiaoergou, 1970 at Zhalantun, 1992 at Suolun and 1984 at Tuliemaodu.A signal transition between an increment and a decrement actually reflects a transition in streamflow between high level and low over a longer time scales.A larger decrement (from 3.29 to −1.31) before 1990 and smaller amplification (from −1.31 to 0.74) after that reveal the overall decreasing trend of streamflow in Beian (Figure 3a).Notably, the streamflow variation of Beian may change in a rising trend stage in future.However, the streamflow of other stations will continuously experience a decreasing trend.Moreover, the inter-decadal variation of streamflow has significant periodic oscillation.As shown in Table 2, the streamflow has the longest period in Beian and Shihuiyao (24.50 years), followed by Tuliemaodu (21.00 years), and by Zhalantun, Suolun and Xiaoergou (16.33 years).

Multiscale Response of Streamflow to Temperature Variations
As depicted in Figure 5, six gauge stations all showed an overall increase (see "Trend" panel) in temperature and approximate quasi-periodicity (see 1st and 2nd panels).Moreover, the temperature signals share similar patterns with streamflow at inter-annual and inter-decadal scales at six gauge stations.Compared with the variations in temperature at six gauge stations, their streamflow has complicated in-phase and anti-phase alternatively at inter-annual time steps.The streamflow has the periods ranging from 3.06 to 3.82 years, which roughly coincides with the periods temperature's ranging between 3.27 and 3.76 years (Table 2).For inter-decadal variation, the streamflow shows an obvious in-phase response before 1980 and anti-phase after 1980 to temperature at Beian, Shihuiyao, Xiaoergou and Zhalantun.With respect to the periodic oscillations, the streamflow at Xiaoergou, Zhalantun and Suolun showed a similar value (i.e., 16.3 years) against the temperature's (16.3 years), but not at Beian and Shihuiyao and Tuliemaodu (Table 2).According to Figure 5 and Table 3, the streamflow and temperature can be divided into three patterns: (1) whole anti-phase at Xiaoergou, Zhalantun and Suolun; (2) anti-phase to in-phase at Beian; and (3) in-phase to anti-phase at Shihuiyao and Tuliemaodu.However, the transition year of trend variation of streamflow and temperature are different.

Multiscale Response of Streamflow to Precipitation
Overall, the precipitation shows an increasing to decreasing trend at all gauge stations during 1961-2010, against an opposite trend at Beian (Table 3 and Figure 6).The precipitation at six gauge stations has an inter-annual periods ranging from 2.88 to 3.27 years, which is similar to the values of 3.06 to 3.82 in streamflow (Table 2).When comparing to inter-decadal periods of precipitation, Suolun has 16.33 years in streamflow and longer periods at remaining gauge stations.Moreover, the streamflow reveals some similar inter-annual patterns (in-phase) with the precipitation at most gauge stations except for Beian (Figure 6a) in 1960s, mid-1980s and 2000s.The trend patterns of precipitation show a consistence with non-liner trend in streamflow at six gauge stations.Note that, although both streamflow and precipitation presented similar nonlinear characteristics in trend, the streamflow have earlier transition years than precipitation's (except Suolun).From the analysis above, precipitation likely plays a major role in variation of streamflow in HANR.

Multiscale Response of Streamflow to ET 0 Variations
Over the past 50 years, ET 0 in the six headwater areas of HANR showed an overall increasing trend with inter-annual periods ranging from 2.88 to 3.06 (Figure 7 and Table 2).The streamflow varied with ET 0 from in-phase to anti-phase at inter-annual scales.At inter-decadal scales, the streamflow shows an anti-phase with ET 0 , especially at Shihuiyao, Xiaoergou, Zhalantun and Suolun with a period of 16.33 years.The streamflow has a longer change period than ET 0 at Beian and Tuliemaodu, which could be expected because ET 0 has less contribution on the change in streamflow.Overall, the trend of ET 0 at six gauge stations all experienced a transition from a relative larger to smaller stage in mid-1980s.Moreover, streamflow response to ET 0 in following patterns: from anti-phase to in-phase at Beian, in-phase to anti-phase at Shihuiyao, Xiaoergou, Zhalantun, Suolun and Tuliemaodu.However, the transition year of ET 0 at Suolun is not obvious, with a strong increase to slight decrease in 2005.

Discussion
This paper investigated nonlinear response of streamflow to temperature, precipitation and potential evapotranspiration in the HANR.Instead of performing a traditional method, we focus on nonlinear and nonstationary characteristics of hydro-meteorological variables based on ensemble empirical mode decomposition method.We found that the response of streamflow to climate change is a nonlinear dynamic process as well as a complicated inconsistence multiscale process.The temperature, precipitation, potential evapotranspiration as well as streamflow all exhibited obviously nonlinear characteristics.This is inconsistent with previous studies conducted by Feng et al. [26], Dong et al. [25] and Li et al. [21].They only found that the streamflow showed a decreasing trend over the past 50 years.In addition, the trend component of streamflow reveals that the change trend of Shihuiyao (Figure 4b) has a transition year in 1991 (strong rising to slight decreasing), which is inconsistent with previous studies conducted by Dong et al. [25] who found two change points (1975 and 1981) based on the non-parametric Pettitt test and Li et al. [21] who detected an abrupt change in 1990 based on Mann-Kendall test.Probably, this is because the trend component focuses on changing process, while the non-parametric Pettitt test and Mann-Kendall test mainly consider rapid change from a statistic characteristic to another.In addition, the test only can assess the statistical significance of trend, but cannot determine its specific shape, no matter the trend is linear or nonlinear [33].Hence, nonlinear analysis of hydroclimatic variables can capture more information about change characteristics than traditional statistical methods (e.g., non-parametric Pettitt test and Mann-Kendall test).
It should be noted that the temperature and ET 0 showed strong nonlinear rising trend, while the streamflow and precipitation experienced nonlinear decreasing trend, which indicated a warm-drying trend in HANR.Yu et al. [34] showed that drought events has become ever more serious both in severity and in extent during the crop-growing season in recent decades over Northeast China.Moreover, wetlands have great sensitivity and vulnerability to the changes in water resources [26,35,36].The warm-drying trend has contributed to considerable wetland degradation and agricultural hazards in Northeast China [37].The continue decreasing in streamflow may possess a large pressure on surface and subsurface water resources [38][39][40].This calls for more authority basin-wide regulations to manage water resources to meet the balance among agricultural sector, ecosystem and domestic water utility.
There are some uncertainties in this study because many other factors can directly or indirectly impact on the response of streamflow to climate change.Here, we only analyzed temperature, precipitation and potential evapotranspiration without consideration of spatial heterogeneity at the watershed and sub-watershed scales [41].Changes in land use/cover is a crucial factor triggering spatiotemporal variability in hydrologic regulation services of ecosystem [41,42].Hydrologic regulation is a significant component of watershed services, and exerts a dynamic force on hydrological circulation and processes that form an ecosystem, leading to temporal and spatial changes in precipitation-streamflow relationship and water distribution [43].For example, the catchment of Shihuiyao Gauge Station (G1 in Figure 1) are primarily covered by forests and wetlands.Wetlands and forests provide important hydrological services including storage and slow release of water, flood mitigation and low flow support, which can alter the response of streamflow regimes to climate change [44][45][46].
Extensive and widespread seasonally frozen and permafrost soils exist in the HANR, which is essential to preserve the water resources and streamflow generation.However, frozen soil experienced a significant degradation trend coupled with a significant climate warming [47,48].Due to warmer temperature, the increase in active layer thickness of permafrost has two effects on hydrological processes as well as the response of streamflow to climate change: (1) more water for evaporation and hence fewer water yield for streamflow; and (2) more water storage in the soil layers and less surface streamflow generation.Although it is accepted that climate change is one of the major drivers for hydrological changes, the effects of permafrost degradation on runoff processes are poorly understood and remain controversial.Thus, more studies are need to address on the question what extent permafrost degradation might impact streamflow changes in the HANR.

Conclusions
This study investigated the nonlinearity of streamflow and its responses to changes in temperature, precipitation, and the corresponding evapotranspiration at the multiple time scales.An Intrinsic Mode Function approach with EEMD was applied to identify the inherent changing patterns and trends at multiple time scales to examine the nonlinear response of streamflow to hydroclimatic variables.We found that, from 1961 to 2010, streamflow, temperature, precipitation and potential evapotranspiration showed apparent multiscale change characteristics, i.e., inter-annual scale (around three-year periodic), inter-decadal scale (around 12-, 16-and 24-year periodic) and nonlinear change trend.The comparatively analysis revealed that streamflow showed corresponding consistent variation characteristics with precipitation, while streamflow experienced a slight consistent different variation change characteristic with temperature and evapotranspiration initially and intensified inconformity variation change characteristic afterwards.The findings suggest that the response of streamflow in headwater regions to climate change is a nonlinear dynamic process dictated by precipitation at the decadal scale and modified by temperature and evapotranspiration within a decade.The results provide insights into how headwater hydrometeorology variables vary at different time scales and how dynamic response processes of streamflow to climate change.

Figure 1 .
Figure 1.Geographical location of the Neijiang River Basin in Northeast China and the spatial distribution of gauge stations (black dots) and weather stations (blue-red symbol) in the basin's headwater region used in this study.

Figure 1 .
Figure 1.Geographical location of the Neijiang River Basin in Northeast China and the spatial distribution of gauge stations (black dots) and weather stations (blue-red symbol) in the basin's headwater region used in this study.

Table 2 .
Streamflow (S), temperature (T), precipitation (P) and potential evapotranspiration (ET 0 ) at six gauge stations on different time scales during 1961-2010 in the headwater region of the Nenjiang River Basin, Northeast China.

Table 1 .
Summary of river gauge stations, tributaries, representing drainage area (DA), annual streamflow (AS), and the identification numbers of meteorological stations (MS No.) in the headwater regions of the Nenjiang River Basin, Northeast China.

Table 2 .
Streamflow (S), temperature (T), precipitation (P) and potential evapotranspiration (ET0) at six gauge stations on different time scales during 1961-2010 in the headwater region of the Nenjiang River Basin, Northeast China.

Table 3 .
Trends (Tre)and transition years (TraY) for temperature (T), precipitation (P), evapotranspiration (ET0), and streamflow (S) variations during 1961-2010 for six gauge stations in the headwater region of the Nenjiang River Basin, Northeast China.

Table 3 .
Trends (Tre) and transition years (TraY) for temperature (T), precipitation (P), evapotranspiration (ET 0 ), and streamflow (S) variations during 1961-2010 for six gauge stations in the headwater region of the Nenjiang River Basin, Northeast China.