Understanding the Main Causes of Runoff Change by Hydrological Modeling : A Case Study in Luanhe River Basin , North China

In the traditional point of view, if there is a significant decreasing trend for a runoff time series, while no significant trend for a precipitation series is present, then an unreliable conclusion will be made that the land surface change is the main contributor to the runoff change. To test it, we selected four sub-watersheds in the Luanhe river basin as the study areas where land use has changed severely. We first detected the long-term rainfall and runoff trend by the Mann–Kendall test, Sen’s slope, and the moving average method, and found that the runoff had a decreasing trend at the 0.05 significance level, while the rainfall had no significant trend in all sub-watersheds. Then an orderly cluster analysis and moving T test method were used to detect the change point of the runoff series. We quantified the contributions of the land surface change and climate variability based on Soil and Water Assessment Tool (SWAT), and the contribution of climate variability accounted for more than 50%, which implies that climate change is the main factor of runoff decrease in the study areas. To further test this, a trend analysis of a reconstructed annual runoff time series under undisturbed conditions has been done. The results showed that in some sub-watersheds, although rainfall series had no significant decreasing trend, the runoff series had significant downward trend. This can be explained by the nonlinear relationship between rainfall and runoff. This study came to a different conclusion from the common view, which observes that runoff decrease is mainly caused by land surface change if rainfall series lacks a significantly decreasing trend.


Introduction
Hydro-climatic variables are subjected to random trends and periodic components, and the hydrologic variability might be attributed to the effect of global climatic change and man-induced disturbances [1].Prominent research of recent years has investigated the existence of non-stationarity in time series, and identified the most probable reason for the non-stationarity around the world.For example, Wagesho et al. (2012) used parametric and non-parametric tests to evaluate the temporal variation of streamflow, lake level, and rainfall at Rift Valley lakes basin of Ethiopia, and further found that the observed trend in these variables is the result of the combined effect of climate change and altered catchment condition [2].Murphy and Ellis (2014) reported that under the potential impacts of climate change, statistically significant temperature increases were detected in eight watersheds of the Colorado River Basin, but their precipitation and runoff remain stationary processes [3].
In the past 50 years, the measured runoff of the six major rivers in China (Yangtze River, Yellow River, Haihe River, Huaihe River, the Pearl River and Songhua River) has decreased [4,5].Especially under the influence of climate change and anthropogenic activities, the temporal and spatial distribution characteristics of water resources have changed significantly.According to the results of the second water resources evaluation of China, it was found that the precipitation of the four river basins (Haihe River, Huaihe River, Yellow River, and the Liaohe River) in North China had decreased 6% on average while the surface water resources had reduced by 17% from 1980 to 2000.Among them, the decrease in Haihe river basin is the most severe [6].
As one of the primary partitions of the Haihe River basin, the Luanhe River not only bears the water demand of the social development of the basin, but also supports the development of the Tianjin and Tangshan as the source of Luanhe-Tianjin Water Diversion Project [7].Therefore, much effort has been made to understand the characteristics of the runoff change in this basin.Among them, Li et al. (2014) pointed out that the characteristics of runoff yield had changed markedly around the 1980s, based on the analysis of the runoff yield characteristics and its change trend [8].Wang et al. (2013) investigated the hydrological effects of the water resources development in Luanhe River basin, in which the results showed that the annual runoff had obviously declined [9].Moreover, Wang et al. (2014) evaluated the spatial and temporal variations in hydro-climatic variables and runoff in response to climate change [10].Feng et al. (2008) studied the change trend of the inflow water resources of Panjiakou Reservoir, which showed an obvious decreasing trend.Their results indicated that the runoff had a decreasing trend in the past half century [11].
Runoff change could be the product of the combined effects of climate change and human activities [12,13].Climate change usually involves both the decrease of precipitation and increases in temperature in the Luanhe river basin, which could lead to a reduction of water resources.Human activities, which affect the runoff, mainly include the land use change, construction of conservancy projects, and water use condition, etc. [14].As a consequence of dramatic changes in the global scale in recent years, it has become a difficult problem to quantify the influence on runoff impacted by the two main factors, surface change and climate variability, in both environmental and hydrological research.
Precipitation is the primary factor that impacts runoff.Increasing precipitation is closely related to the increasing of flow in summer and the increasing of flow in summer impacts the increasing of annual flow.The ascending of air temperature, especially in autumn, causes greater snow melting at high mountain elevations, and therefore causes flow increase during the non -flood period, which accelerates river basin evaporation, especially in summer, and also makes the depletion of ground water resource increasing, and has an elimination function for the annual flow increase caused by the precipitation.Water conservation of different types of land use is significantly varied, thus having an indirect effect on the runoff.What is more, land use changing to land weakens the effect of the rainfall intensity and the antecedent soil moisture conditions on the rainfall-runoff relationship.
Liu et al. [15] built a Soil and Water Assessment Tool (SWAT) model in the Luanhe River basin to estimate the runoff at disturbed and undisturbed stages with different climatic data and human activity effects, thus analyzed the quantitative impact analysis of climate change and human activities on runoff.From the common perspective, land surface change (or human activities) plays a prominent role in runoff decrease if there is no significant decrease trend in the rainfall series.The study as mentioned above draws a similar conclusion that human activities, with a contribution rate of 60-80%, are the most important factors to the runoff decrease [9,16,17].However, taking the primary effect on runoff into consideration, different conclusions have been drawn by some other researchers.For example, Liu et al. [18] found that the influence of climatic factors and human activities on runoff amount respectively accounted for 55% and 45% of the total reduced runoff; i.e., the climatic factors play the primary role in the runoff change.
To further understand the main causes of runoff change in Luanhe river basin, this paper carried out a systematical analysis based on the long-term rainfall and runoff series.Firstly, we analyzed the variation trend of rainfall and runoff by the Mann-Kendall (M-K) test, and used Sen's slope method and the moving T method to detect the change point of the runoff series.Then, the runoff series under undisturbed (before the change point) conditions for the disturbed period (after the change point) was reconstructed using a SWAT model, and the contributions of land surface change and climate variability were quantified based on the modeling results.

Study Area
The Luanhe river basin, located in the northeast part of China, has a drainage area of 44,750 km 2 , and lies between 115 • 40 -119 • 20 E longitude and 39 • 10 -42 • 35 N latitude with the elevation ranging from 2 to 2205 m (the average elevation is 766 m).The topography significantly descends from northwest to southeast and about 98% of the drainage area is mountainous region while about 2% is plain.The watershed receives an average precipitation of 560 mm, mostly in the summer (70-80%), with an average runoff of 46.94 × 10 8 m 3 per year.The average temperature ranges from −0.3 • C to 11 • C for the whole area and presents a diminishing trend during the periods concerned.With an average potential evapotranspiration of 950-1150 mm per year, the Luanhe river basin sometimes could reach a high value of 1801 mm annually.The average absolute humidity is 9.6 g/m 3 , and the average relative humidity is 60-70%.
The Luanhe River basin is well known for water supply function for the Tianjin city, an important metropolis of China.It was planned to introduce water of a billion cubic meters to Tianjin per year from the basin.However, the actual amount of water transferred to Tianjin was less than the plan in several years, especially in the last decade.The deficient water supply to Tianjin is mainly due to the decrease of annual water storage in the Panjiakou reservoir in the 21st century.Due to rainfall reduction, land use change and construction of many small check dams for soil and water conservation, the average annual runoff has decreased by about 30% after 1980 [7].The annual runoff was below the mean of the time series in the last few years The long-lasting water shortage aggravated the water crisis in Tianjin city.In this study, four sub-watersheds were selected, which were the Luanhe (Sandaohezi station), Yixunhe (Hanjiaying station), Wuliehe (Chengde station), and Liuhe (Liying station) with the area of 17,100, 6761, 2460 and 626 km 2 , respectively.
The total area of the four sub-watershed accounts for more than 60% of the entire Luanhe River basin.The locations of the four hydrological stations seem to be close to each other, but their upstream control watersheds distribute quite uniformly in the entire basin.Specifically, the watershed controlled by Sandaohezi station covers themost area of the upstream Luanhe River basin.The watersheds of Hanjiaying and Chengde stations are located in the middle part of the basin, and the watershed of Liying station lies in the downstream part.Therefore, the four typical sub-watersheds can represent the overall condition of the Luanhe River basin.

Data
Provided by Hydrology and Water Resource Survey Bureau of Hebei Province, daily rainfall and monthly runoff data recorded at 4 hydrological stations were adopted for the selected sub-watersheds.The locations of the four stations are depicted in Figure 1.The data were available from 1956-2010 in Luanhe, Wuliehe, and Liuhe sub-watersheds and 1956-2003 in Yixunhe sub-watershed.
On the basis of various data, such as station information, precipitation, temperature, solar radiation, daily average wind speed and relative humidity, a meteorological database in a SWAT model has been built to simulate the effects on runoff.Daily precipitation data ranging from 1957 to 2010 is acquired from 10 gauge stations (Sandaohezi, Chengde, Kuancheng, Liying, Xiabancheng, Xiahenan, Xuanjiangyingzi etc.).We appreciate the China Meteorological Data Sharing Service System for providing highly valuable documents of the meteorological data, which contained temperature, wind speed, sunshine duration, and relative humidity.Daily meteorological feature values measured by the Duolun, Weichang, Chengde, and Qinglong meteorological stations from 1957 to 2010 were adopted.Monthly flow data, recorded by Sandaohezi, Chengde, Liying, Pingquan stations with a time series of 1957 to 2010 were selected, and a relative shorter time series measured in Hanjiaying station during 1957-2003 was also applied in the SWAT model.Locations of the gauge stations are shown in Figure 1.
stations with a time series of 1957 to 2010 were selected, and a relative shorter time series measured in Hanjiaying station during 1957-2003 was also applied in the SWAT model.Locations of the gauge stations are shown in Figure 1.

Methods for Trend Identification in Hydrological Time Series
Methods for trend identification in hydrological time series include the linear regression method, cumulative departure method, moving average method, and the Mann-Kendall rank correlation test method, etc.
The Mann-Kendall (M-K) test is one of the most widely used non-parametric tests for trend detection in hydrologic time series, and Sen's slope estimator is usually employed for identifying the true slope of trend analysis [19][20][21].The Sen's procedure is applied following the M-K test to measure the magnitude of any significant trend found in the M-K test.
Thus, this study adopted the Mann-Kendall rank correlation test and Sen's Slope method to diagnose the trend of the runoff and rainfall change of four sub-basins in Luanhe River basin, furthermore moving average method was used to verify the results.

Mann-Kendall Rank Correlation Test
The Mann-Kendall rank correlation test is widely used in the trend analysis of continuous time series [22,23].It does not require the test sequence to obey a certain distribution.The test method is as follows: For time series X = {x1, x2, …, xn}, Xi (i = 1, 2, …, n) are variables for analysis, where n is the length of the time series, and its standard normal variable T can be calculated.
The M-K test statistic T is generally used to measure the degree to which a trend is consistently increasing or decreasing.If T is positive, the sequence has an increasing (rising) trend; if T is negative, it shows a decreasing (descending) trend.For samples that are trend-free, the T value is close to zero and the p-value is close to 0.5 [24].Compared with the standard normal variate (Tα/2) at the desired

Methods for Trend Identification in Hydrological Time Series
Methods for trend identification in hydrological time series include the linear regression method, cumulative departure method, moving average method, and the Mann-Kendall rank correlation test method, etc.
The Mann-Kendall (M-K) test is one of the most widely used non-parametric tests for trend detection in hydrologic time series, and Sen's slope estimator is usually employed for identifying the true slope of trend analysis [19][20][21].The Sen's procedure is applied following the M-K test to measure the magnitude of any significant trend found in the M-K test.
Thus, this study adopted the Mann-Kendall rank correlation test and Sen's Slope method to diagnose the trend of the runoff and rainfall change of four sub-basins in Luanhe River basin, furthermore moving average method was used to verify the results.

Mann-Kendall Rank Correlation Test
The Mann-Kendall rank correlation test is widely used in the trend analysis of continuous time series [22,23].It does not require the test sequence to obey a certain distribution.The test method is as follows: For time series X = {x1, x2, . . ., xn}, Xi (i = 1, 2, . . ., n) are variables for analysis, where n is the length of the time series, and its standard normal variable T can be calculated.
The M-K test statistic T is generally used to measure the degree to which a trend is consistently increasing or decreasing.If T is positive, the sequence has an increasing (rising) trend; if T is negative, it shows a decreasing (descending) trend.For samples that are trend-free, the T value is close to zero and the p-value is close to 0.5 [24].Compared with the standard normal variate (T α/2 ) at the desired significance level α, the significance of the trend is evident if |T| > |T α/2 |, and the trend is not statistically significant otherwise.

Sen's Slope Method
In order to further estimate the degree of change in time series, this study calculates the Sen's Slope [25] of each meteorological and hydrological variable.Sen's Slope is able to avoid the defect in time series and eliminate the outliers' interference to the time series.Positive(negative) represents a rising (descending) trend.

Moving Average Method
The moving average method can avoid the impact of the sequence fluctuations to some extent and make the tendency and periodicity of hydrological time series more intuitive and evident [26].

Method for the Change Point Identification of Time Series
The aim of the detection of change point is to test the most possible change time and amplitude of variation of the studied statistical characteristic values.There are many methods to determine the change point, such as moving F test method [27], moving T test method [28], orderly cluster analysis method [29], R/S analysis method [30], moving rank sum test method [31], Lee-Heghinian test method, Brown-Forysthe test method, Bayesian test method [32] and optimal information segmentation method [33], etc.Since the hypothesis and applicable conditions of these methods for determining the change point in hydrological time series are different, the results of these various test methods are quite different.Therefore, the hypothesis and applicable conditions should be paid attention on when one selects the test method for determining the change point.The drawback of the moving F test is that the selection of the two terminal lengths before and after the change point is not clear.The R/S analysis method is not robust for short-term memory and heteroscedasticity in time series, so it cannot accurately distinguish between short-term memory and long-term memory in time series.At the same time, the orderly cluster analysis and the moving T test method have high efficiency for the mutation sequence test, while the R/S test, the optimal information dichotomy method and the moving F test are not sensitive to the mutation sequence, and the test efficiency is low and the effect is poor.Based on this, this study adopted an orderly cluster analysis and moving T test methods to test the change point of the time series and carry out a comparative analysis.

Ordered Cluster Analysis (OCA) Method
The ordered cluster analysis method is used to estimate the possible significant interference points of the sequence after ordered clustering.In this method, the order of the time series cannot be disturbed and the ordered classification is used for deriving the possible interference point.Its essence is to ascertain the optimal segmentation point, i.e.; the sum of the squared deviations of the same class is required to be minimal, while that of the different class should be large.The principle is as follows [34]: where, x τ and x n−τ are average value of the sequences before and after the disturbance point τ.Sum of the squared deviation is calculated by: which should satisfy S * n (τ) = min|S n (τ)|, 1 ≤ τ ≤ n and τ is the optimal segmentation point, namely the disturbance point (change point).This method applies in sequences of various types of distribution, but it is difficult to deal with the case that the interference point is located at the end of the sequence.

Moving T Test (MTT) Method
The moving T test method was proposed, aiming to solve the problem of the traditional T test, which is able to test but not search for the change point.Assume the distribution function before and after the change point τ as F 1 (x) and F 2 (x), and extract two samples with the capacity of n 1 and n 2 from them respectively, then construct a T statistic as follows: where x i , S i are the means and variances of the samples.T obeys the distribution of t(n 1 + n 2 − 2).At a confidence level of α, if T > t α/2 , the sequence has significant difference.Select the point that T reaches the maximum among all the points τ which satisfy T > t α/2 as the most possible change point τ 0 .
The first three most significant change points should be excluded as required by these test methods.Considering this condition and the causes of the change, a most significant change point can be determined.

Method for Defining the Impacts of Climate Change and Human Activities on Runoff Using SWAT
Because the hydrological modelling method can take the complexity of hydraulic problems into consideration [35], as well as describe the hydrological variables and their evolution laws accurately [16], it is popular in the quantitative identification of climate-and human-induced impacts on runoff.
SWAT is a physically based distributed hydrological model running on a daily time step [36].In SWAT, a study watershed is first divided into multiple sub-basins and then subdivided into hydrologic response units (HRUs).These HRUs are defined as homogeneous spatial units of unique soil, land use, and slope characteristics, capable of describing spatial heterogeneity in land cover and soil types within the watershed.For each HRU, the hydrologic cycle is simulated based on the water balance equation, and the relevant hydrologic components (such as precipitation, interception, surface runoff, evapotranspiration, percolation, and lateral flow from the soil profile, and return flow from shallow aquifers) are estimated within the model.Flows summed across all HRUs in a sub-basin are routed through the stream system to the watershed outlet.Several options are provided by the SWAT for the simulation of hydrologic processes.Details of the model description and the related calculations are given by Neitsch et al. (2011) and Arnold et al. (2012) [37,38].
Based on this, SWAT distributed hydrological model has been selected in this study.The following steps are involved in order to construct a spatial database in SWAT: (1) generate the water system and divide the basin to several sub-basins; (2) define the land use type, soil type, slope data, etc. in the study area; (3) create hydrological response units (HRUs) in the watershed.Specifically, divide the study area to 33 sub-basins using the SWAT model based on the DEM data in the upstream Luanhe River basin firstly, then check the National Soil Survey using Chinese Soil Database and transfer the soil texture, map the soil type distribution and establish the soil type index table, thus build up the soil database.In addition, divide the slope to five grades according to the "Technical specification for investigation of land use status" issued by National Agricultural Zoning Committee in 1984, and set the critical values of the land use type, soil type and slope as 20%, 10%, 20% respectively, and finally use Multiple HRUs to divide the 33 sub-basins into 857 HRUs.
The SWAT model considers not only the change of precipitation and temperature, but also other factors such as the wind speed, humidity, and sunshine duration.Using SWAT to identify the impact of climate change and human activities on runoff includes two procedures.Firstly, divide the hydrological time series to two stages (undisturbed stage or natural stage, and disturbed stage).During the natural stage, measured runoff and natural runoff are almost the same, and the runoff change is mainly because of climate change.Secondly, simulate the natural runoff during the disturbed stage according to the measured runoff during the undisturbed stage.
The different parts of impact can be calculated by the following equations: where, ∆Q is the total amount of runoff change; ∆Q C and ∆Q H are the runoff change amount due to climate change and human activities; Q HR is the observed runoff during the disturbed stage; Q HN is the simulated natural runoff during the disturbed stage and Q HB represents the runoff during the undisturbed stage.
According to the Equation ( 7), the part of impact caused by climate change can be calculated by the natural runoff during the undisturbed stage and simulated natural runoff during the disturbed stage, and based on Equation ( 8), the other part of impact due to human activities is the difference between the simulated natural runoff and the observed runoff during the disturbed stage.Then the contribution rates of the two factors can be computed by the following formulae: in which, Cs and Hs represent the contribution rates of climate change and human activities to the runoff change respectively.

Analysis of Rainfall Variation in Luanhe River Basin
The rainfall data are mainly from 10 rainfall gauge stations in the Luanhe River basin that can well reflect the rainfall variation in the basin and are distributed relatively uniformly.The 5-year moving average curve can be depicted according to the annual average rainfall amount data from 1957 to 2010 measured by the 10 stations (see Figure 2), and the characteristic values of the annual average rainfall amount of each decade in the Luanhe River basin are listed in Table 1.
Water 2018, 10, x FOR PEER REVIEW 7 of 17 where, ∆Q is the total amount of runoff change; ∆QC and ∆QH are the runoff change amount due to climate change and human activities; QHR is the observed runoff during the disturbed stage; QHN is the simulated natural runoff during the disturbed stage and QHB represents the runoff during the undisturbed stage.
According to the Equation ( 7), the part of impact caused by climate change can be calculated by the natural runoff during the undisturbed stage and simulated natural runoff during the disturbed stage, and based on Equation ( 8), the other part of impact due to human activities is the difference between the simulated natural runoff and the observed runoff during the disturbed stage.Then the contribution rates of the two factors can be computed by the following formulae: (7) in which, Cs and Hs represent the contribution rates of climate change and human activities to the runoff change respectively.

Analysis of Rainfall Variation in Luanhe River Basin
The rainfall data are mainly from 10 rainfall gauge stations in the Luanhe River basin that can well reflect the rainfall variation in the basin and are distributed relatively uniformly.The 5-year moving average curve can be depicted according to the annual average rainfall amount data from 1957 to 2010 measured by the 10 stations (see Figure 2), and the characteristic values of the annual average rainfall amount of each decade in the Luanhe River basin are listed in Table 1.In Figure 2, it can be seen that the 5-year moving average curve shows a downward trend, indicating that the average annual rainfall amount from 1957-2010 in Luanhe River basin has a decreasing trend.In Table 1, the decades of the 1960s, 1970s and 1990s are wet as their annual rainfall amounts are higher than the multi-year average level (521 mm), while the decades of the 1980s and 2000s are dry with large anomalies (−4.4% and −7.1%, respectively).
In addition, when comparing the annual rainfall amount of each decade, it is found that since the 1960s the annual rainfall firstly decreased, then increased and again decreased in the central part of the basin.Especially in the decade of the 1990s, the rainfall was relatively abundant.
In order to test if the rainfall amount of the basin has changed temporally with an obvious trend, the Mann-Kendall (M-K) correlation test method is adopted to test the change trend of the annual rainfall amount of each sub-basin (see Table 2).From the test results in Table 2 it can be found that: (1) the T test values of the annual rainfall amount of all the sub-basins are lower than the critical value at the significance level of α = 0.05 and negative, along with the negative Sen's Slope estimates, which demonstrates a decreasing but not significant trend of the rainfall in the Luanhe River basin; (2) it can be seen that the change trend of rainfall calculated by Sen's Slope method in the Hanjiaying gauge station is relatively evident with a decreasing rate of around 1.26 mm/a, while it is least significant in Chengde station with a decreasing trend of 0.65 mm/a.

Analysis of Runoff Variation in Luanhe River Basin
Hereon, runoff change trend and change point of the four selected stations are analyzed.Among them, the observed data of Sandaohezi station, Chengde station, and Liying station are from 1957 to 2010, while the data length of Hanjiaying station is from 1957 to 2003.

Trend Analysis of the Runoff Time Series
The 5-year moving average method and M-K trend test method are used to analyze the trend of the annual average runoff of the four selected hydrological stations in Luanhe River basin.The results of the two methods are shown in Figure 3 and Table 3  In Figure 2, it can be seen that the 5-year moving average curve shows a downward trend, indicating that the average annual rainfall amount from 1957-2010 in Luanhe River basin has a decreasing trend.In Table 1, the decades of the 1960s, 1970s and 1990s are wet as their annual rainfall amounts are higher than the multi-year average level (521 mm), while the decades of the 1980s and 2000s are dry with large anomalies (−4.4% and −7.1%, respectively).
In addition, when comparing the annual rainfall amount of each decade, it is found that since the 1960s the annual rainfall firstly decreased, then increased and again decreased in the central part of the basin.Especially in the decade of the 1990s, the rainfall was relatively abundant.
In order to test if the rainfall amount of the basin has changed temporally with an obvious trend, the Mann-Kendall (M-K) correlation test method is adopted to test the change trend of the annual rainfall amount of each sub-basin (see Table 2).From the test results in Table 2 it can be found that: (1) the T test values of the annual rainfall amount of all the sub-basins are lower than the critical value at the significance level of α = 0.05 and negative, along with the negative Sen's Slope estimates, which demonstrates a decreasing but not significant trend of the rainfall in the Luanhe River basin; (2) it can be seen that the change trend of rainfall calculated by Sen's Slope method in the Hanjiaying gauge station is relatively evident with a decreasing rate of around 1.26 mm/a, while it is least significant in Chengde station with a decreasing trend of 0.65 mm/a.

Analysis of Runoff Variation in Luanhe River Basin
Hereon, runoff change trend and change point of the four selected stations are analyzed.Among them, the observed data of Sandaohezi station, Chengde station, and Liying station are from 1957 to 2010, while the data length of Hanjiaying station is from 1957 to 2003.

Trend Analysis of the Runoff Time Series
The 5-year moving average method and M-K trend test method are used to analyze the trend of the annual average runoff of the four selected hydrological stations in Luanhe River basin.The results of the two methods are shown in Figure 3 and Table 3 respectively.From Figure 3 it can be observed that the change of annual average runoff of all the four typical hydrological stations is similar and show a decreasing trend.In Table 3, the T test values of annual runoff of all these stations are negative and larger than the corresponding critical value at a confidence level of 95%, indicating that the annual runoff of these stations have a significant decreasing trend.The Sen's Slope of runoff is maximum in Chengde station, (−1.21 m 3 /(s•year)) and minimum in Liying station (−0.05 m 3 /(s•year)).From Figure 3 it can be observed that the change of annual average runoff of all the four typical hydrological stations is similar and show a decreasing trend.In Table 3, the T test values of annual runoff of all these stations are negative and larger than the corresponding critical value at a confidence level of 95%, indicating that the annual runoff of these stations have a significant decreasing trend.The Sen's Slope of runoff is maximum in Chengde station, (−1.21 m 3 /(s•year)) and minimum in Liying station (−0.05 m 3 /(s•year)).

Change Point Analysis of the Time Series
The change point of the runoff time series is taken as the transition point of the undisturbed stage and disturbed stage.In order to guarantee the reliability of the analysis results, the ordered cluster analysis method (for simplicity, log(Sn) denote as S) and moving T test method have been adopted to analyze the change point of the annual runoff time series (see Figures 4 and 5).The detailed results of the change point test of all the selected stations are listed in Table 4.

Change Point Analysis of the Time Series
The change point of the runoff time series is taken as the transition point of the undisturbed stage and disturbed stage.In order to guarantee the reliability of the analysis results, the ordered cluster analysis method (for simplicity, log(Sn) denote as S) and moving T test method have been adopted to analyze the change point of the annual runoff time series (see Figures 4 and 5).The detailed results of the change point test of all the selected stations are listed in Table 4.    1959, 1964, 19791959, 1964, 19791959, 1964, 1979Yixunhe Hanjiaying station. 1959, 1979, 19641959, 1979, 19621959, 1979, 1962, 1964Wuliehe Chengde station. 1962, 1979, 19981962, 1979, 19981962, 1979, 1998Liuhe Liying station. 1962, 1979, 19961962, 1979, 19961962, 1979, 1996 From Figures 4 and 5, and Table 4, it can be observed that the possible change points of the annual runoff time series of the sub-basins detected by ordered cluster analysis method and moving T test method are 1959, 1962, 1964, 1965, 1979, 1991, 1996 and 1998.In order to determine the most possible change point, one needs to take the causes of the change into account.Among them, in 1962Among them, in , 1963Among them, in , 1966Among them, in , 1991 and 1996 extreme floods occurred during the flood season, which were caused by

Change Point Analysis of the Time Series
The change point of the runoff time series is taken as the transition point of the undisturbed stage and disturbed stage.In order to guarantee the reliability of the analysis results, the ordered cluster analysis method (for simplicity, log(Sn) denote as S) and moving T test method have been adopted to analyze the change point of the annual runoff time series (see Figures 4 and 5).The detailed results of the change point test of all the selected stations are listed in Table 4.    1959, 1964, 19791959, 1964, 19791959, 1964, 1979Yixunhe Hanjiaying station. 1959, 1979, 19641959, 1979, 19621959, 1979, 1962, 1964Wuliehe Chengde station. 1962, 1979, 19981962, 1979, 19981962, 1979, 1998Liuhe Liying station. 1962, 1979, 19961962, 1979, 19961962, 1979, 1996 From Figures 4 and 5, and Table 4, it can be observed that the possible change points of the annual runoff time series of the sub-basins detected by ordered cluster analysis method and moving T test method are 1959, 1962, 1964, 1965, 1979, 1991, 1996 and 1998.In order to determine the most possible change point, one needs to take the causes of the change into account.Among them, in 1962,1963,1966,1991 and 1996 extreme floods occurred during the flood season, which were caused by   1959, 1964, 19791959, 1964, 19791959, 1964, 1979Yixunhe Hanjiaying station. 1959, 1979, 19641959, 1979, 19621959, 1979, 1962, 1964Wuliehe Chengde station. 1962, 1979, 19981962, 1979, 19981962, 1979, 1998Liuhe Liying station. 1962, 1979, 19961962, 1979, 19961962, 1979, 1996 From Figures 4 and 5, and Table 4, it can be observed that the possible change points of the annual runoff time series of the sub-basins detected by ordered cluster analysis method and moving T test method are 1959, 1962, 1964, 1965, 1979, 1991, 1996 and 1998.In order to determine the most possible change point, one needs to take the causes of the change into account.Among them, in 1962Among them, in , 1963Among them, in , 1966Among them, in , 1991Among them, in and 1996 extreme floods occurred during the flood season, which were caused by extreme storms in the Luanhe River basin [39], therefore, it is difficult to assert that the change of runoff in these years is caused by the land use change.In addition, the runoff data of the Luanhe, Wuliehe, Liuhe River basins are from 1957-2010, and from the perspective of statistics, the points close to the head (1959) and end points (1998) are not reliable to be taken as change points.
Remotely sensed land use data were provided by the Chinese Academy of Sciences.The land use in the Luanhe River Basin was classified into six types: forest (FRST), grassland (PAST), arable land (AGRL), water body (WATR), urban area (URHD), and bare land (URLD).Remotely sensed land use data of 1970 and 1980 were used to quantify the effects of land use change on annual runoff in the Luanhe River basin.The area of each land use type was obtained by ArcGIS software and the results are shown in Figure 6 and Table 5.
From 1970 to 1980, the great transformation of these land use types resulted in the variation of their generation mechanisms.In addition, soil and water conservation mainly constructed after 1980 in the area held more water in land surface, so less runoff flow was yielded in the basin.
Water 2018, 10, x FOR PEER REVIEW 11 of 17 extreme storms in the Luanhe River basin [39], therefore, it is difficult to assert that the change of runoff in these years is caused by the land use change.In addition, the runoff data of the Luanhe, Wuliehe, Liuhe River basins are from 1957-2010, and from the perspective of statistics, the points close to the head (1959) and end points (1998) are not reliable to be taken as change points.
Remotely sensed land use data were provided by the Chinese Academy of Sciences.The land use in the Luanhe River Basin was classified into six types: forest (FRST), grassland (PAST), arable land (AGRL), water body (WATR), urban area (URHD), and bare land (URLD).Remotely sensed land use data of 1970 and 1980 were used to quantify the effects of land use change on annual runoff in the Luanhe River basin.The area of each land use type was obtained by ArcGIS software and the results are shown in Figure 6 and Table 5.
From 1970 to 1980, the great transformation of these land use types resulted in the variation of their generation mechanisms.In addition, soil and water conservation mainly constructed after 1980 in the area held more water in land surface, so less runoff flow was yielded in the basin.Thus, according to the actual underlying surface change condition of the Luanhe River basin, 1979 is determined as the change point of all the sub-basins.

Calibration and Verification of SWAT Model
As stated above, the change point of the runoff series in the Luanhe river basin occurred in 1979, thus we consider 1956-1979 to be the reference period (undisturbed period) and 1980-2010 to be the disturbed period.In the application of SWAT model, the observed runoff data of the undisturbed stage, from 1965-1975, are used for calibration while data from 1976-1979 are adopted for validation.
In order to avoid the influence of initial conditions on the simulation results, the warm-up period of the model was selected from 1962 to 1964.The time step is set as 1 day.Monthly runoff data from 1965 to 1975 of Sandaohezi, Hanjiaying, Chengde, and the Li Ying hydrological stations are used to calibrate model parameters.The results can be seen in Table 6.Thus, according to the actual underlying surface change condition of the Luanhe River basin, 1979 is determined as the change point of all the sub-basins.In order to avoid the influence of initial conditions on the simulation results, the warm-up period of the model was selected from 1962 to 1964.The time step is set as 1 day.Monthly runoff data from 1965 to 1975 of Sandaohezi, Hanjiaying, Chengde, and the Li Ying hydrological stations are used to calibrate model parameters.The results can be seen in Table 6.We selected four sub-watersheds with hydrological stations to calibrate and verify the SWAT model in the reference period.The Nash-Sutcliffe coefficient of efficiency (E ns ) [40] and coefficient of determination (R 2 ) [41] in all the stations are greater than 0.55, illustrating a good performance of SWAT in these sub-watersheds [39].Especially in Liying St., the E ns and R 2 are both 0.76 in the calibration period.During the verification period, the E ns and R 2 reach as high as 0.82 and 0.83 respectively.This reflects the capability of SWAT for modeling the runoff processes in these regions.
The hydrological module in SWAT is calibrated by four sub-basins, and the monthly runoff data from 1976 to 1978 are selected as the validation period.The calibration and verification results of each sub basin are shown in Table 7.In order to distinguish the effects of different factors, we simulated monthly discharges in the disturbed period under 1970 land use condition.We obtained the simulated annual runoff time series by summing up the simulated monthly discharge in the disturbed period .The mean annual runoff during the undisturbed period was calculated and defined as the benchmark value (42.92 mm).The value was considered as the representative of runoff magnitude under the natural situation without impact of human activities.Based on this, taking the Luanhe sub-watershed for example, the contributions of climate change and human activities to runoff can be calculated using the Formulas ( 6)-( 9), and the results are shown in Table 5 [18] where the HBV model has been used for the runoff time series of 1960-2007.In the same way as stated above, the impact of climate change and human activities on runoff in all the selected sub-basins can be calculated, see Table 9.As can be seen in Table 9, climate change is the main cause of the runoff reduction in all the sub-basins.In the Liuhe sub-basin, the contributions of climate change and human activities to the runoff reduction accounted for 84.11% and 15.89%, respectively.The relative impact of climate change in the Wuliehe sub-basin is less than the other sub-basins, namely 60.93%.
As mentioned above, there were many different results of the impacts of climate and land use changes in the Luanhe river basin, and different conclusions were drawn in these studies.Some studies concluded that human activities are the most important factor in the runoff decrease [15,42], which is consistent with the common view that "runoff decrease was mainly caused by land surface change (or human activities) if rainfall series had no significant decrease trend".However, there are also some other researches that have drawn a similar conclusion as this study [18,43].They found that climate change was the dominant factor of runoff decrease.The different conclusions of these studies may be due to the selection of different research periods, hydrological stations, and hydrological models.Furthermore, as stated in [44], one could not separate the contributions of land use change and precipitation variability accurately before understanding the mechanism of the concurrent effects of them.However, one thing is for sure: with the enhancement of human activities, the impact of climate change on runoff is gradually decreased, while that of human activities is gradually increased.

Analysis of Runoff Variation Without the Effect of Human Activities
The simulated annual runoff the disturbed period (1980-2010) under 1970 land use conditions could be combined with the annual runoff in the undisturbed period  to form a long annual runoff time series (1956-2010, we can call it "reconstructed annual runoff"), in which the effect of human activities has been excluded.To further understand the main cause of the runoff change, a trend analysis of the reconstructed annual runoff time series has been analyzed using the M-K trend test method.Since there was a catastrophic flood in 1959, which may influence the trend analysis, we tested the time series under two conditions: with 1959 and without 1959.
From Table 10 it can be seen that the test values of the annual runoff in Sandaohezi St. is higher than the critical value with the significance level of α = 0.05 and negative, which demonstrates a significant decreasing trend, although the annual rainfall of this station showed a "not significant" decreasing trend.Therefore, we can conclude that even without the effect of human activities, a "not significant decrease" of annual rainfall may result in a significant decrease of runoff, so the common view that "runoff decrease was mainly caused by land surface change if rainfall series had no significant decrease trend" is arbitrary and unreliable.By comparing Tables 2 and 10, it can be found that the T values of the reconstructed annual runoff ("With 1959" column) are much larger than those of the annual rainfall in all the stations except the Hanjiaying Station.This implies that a reduction of rainfall usually leads to a greater degree of reduction in runoff, which can be explained by the nonlinear relationship between rainfall and runoff.The nonlinear relationship has been found in many previous works [45,46].The study of the impact of climate change and human activities on runoff is based on the assumption that the two factors are independent.In fact, CO 2 emissions from human activities will also have an impact on climate change, thus climate change and human activities are interacting and restricting each other.In addition, this study has analyzed uncertainties that may exist in quantitative identification.Therefore, a more accurate quantitative identification of climate change and the impact of human activity on runoff is needed to further study the impact of many factors.
In our study, the effect of human activities was evaluated by the difference between simulated and observed runoff during the disturbed period.However, because of the SWAT model error, a part of the difference might be caused by the model error.Although the model error was ignored in this study, which might lead to overestimated results, the SWAT model performed well and satisfied the accuracy requirement.In future study, we will pay more attention on the SWAT model error and parameters uncertainty to improve the results.

Conclusions
In order to understand the causal mechanism of hydrological droughts in the Luanhe river basin, this paper adopted appropriate statistical methods to analyze the variational characteristics of rainfall and runoff in the Luanhe river basin based on long term hydrological observation data.A SWAT model was used to quantitatively identify the impact of climate change and human activities on runoff.The main conclusions obtained in this work are as follows: 1.
Through M-K trend test, the precipitation showed a decrease trend but not significant, while the runoff had a significant reduction trend.Using the moving T test method and the orderly cluster analysis method to identify the change point, and considering the actual change of the underlying surface in the Luanhe river basin, the change point was determined as 1979.

2.
With the assumption that the impacts of climate change and human activities on runoff are independent, a SWAT model was adopted to simulate the "reconstructed runoff" after the change point (disturbed period, 1980-2010) using the 1970 land use condition, thus distinguishing the impact degree of climate change and human activities on runoff change.The results showed that climate change was the main cause of runoff decrease in the selected sub-basins.

3.
The simulated annual runoff sequence during the disturbed period under the 1970 land use condition was combined with the annual runoff time series in the undisturbed period to form a "reconstructed runoff time series", which excluded the impact of land use change on runoff change.Through the M-K trend test, without the effect of land use change, the rainfall decrease with a "not significant trend" may lead to a significant runoff decrease, and rainfall reduction to a certain degree usually results in a greater reduction in runoff.This is due to the nonlinear relationship between rainfall and runoff.
From the common point of view, if the rainfall sequence has no significant trend, the runoff decrease is mainly caused by land surface change.This paper draws a different conclusion from that perspective.The nonlinear relationship, which has been found in many previous studies, may make the "not significant rainfall reduction" to be the main cause of the "significant runoff reduction".

Figure 1 .
Figure 1.Location of the Luanhe River basin and the selected stations.

Figure 1 .
Figure 1.Location of the Luanhe River basin and the selected stations.

Figure 2 .
Figure 2. 5-year moving average annual rainfall amount in the Luanhe River basin.

Figure 3 .
Figure 3. 5-year moving average annual runoff hydrograph of each station in Luanhe River basin.

Figure 3 .
Figure 3. 5-year moving average annual runoff hydrograph of each station in Luanhe River basin.

Figure 4 .
Figure 4.The ordered cluster test results of the annual runoff of each hydrological station.

Figure 5 .
Figure 5.The moving T test results of the annual runoff of each selected hydrological station.

Table 4 .
Detailed change point test results of the runoff time series of each selected hydrological station in the upstream of Luanhe River basin (significance level = 0.05).

Figure 4 .
Figure 4.The ordered cluster test results of the annual runoff of each hydrological station.

Figure 4 .
Figure 4.The ordered cluster test results of the annual runoff of each hydrological station.

Figure 5 .
Figure 5.The moving T test results of the annual runoff of each selected hydrological station.

Table 4 .
Detailed change point test results of the runoff time series of each selected hydrological station in the upstream of Luanhe River basin (significance level = 0.05).

Figure 5 .
Figure 5.The moving T test results of the annual runoff of each selected hydrological station.

Table 4 .
Detailed change point test results of the runoff time series of each selected hydrological station in the upstream of Luanhe River basin (significance level = 0.05).

4. 3 .
Quantitative Identification of the Impact of Climate Change and Human Activities on Runoff Using the SWAT Model 4.3.1.Calibration and Verification of SWAT Model As stated above, the change point of the runoff series in the Luanhe river basin occurred in 1979, thus we consider 1956-1979 to be the reference period (undisturbed period) and 1980-2010 to be the disturbed period.In the application of SWAT model, the observed runoff data of the undisturbed stage, from 1965-1975, are used for calibration while data from 1976-1979 are adopted for validation.

Table 1 .
Characteristic values of the annual average rainfall amount of each decade in the Luanhe River basin.

Table 2 .
Annual rainfall variation trend test results of each sub-basin. respectively.

Table 2 .
Annual rainfall variation trend test results of each sub-basin.

Table 3 .
Test results of annual average runoff change trend of each hydrological station.

Table 3 .
results of annual average runoff change trend of each hydrological station.

Table 5 .
Comparison of land use types in the study area in 1970 and 1980.

Table 5 .
Comparison of land use types in the study area in 1970 and 1980.

Table 6 .
Results of the Soil and Water Assessment Tool (SWAT) model calibration.

Table 7 .
Results of calibration and validation.
. It can be found that: (1) Since 1980s, in comparison with the benchmark value (42.92 mm) the observed runoff of all the decades have decreased to different extents.Especially in the decade of 2000-2010, the observed runoff is 15.08 mm, which has a difference of 27.84 mm compared with the value.The contributions of climate change and human activities are 56.86% and 43.14%, respectively.(2) In comparison with the benchmark value, both the modeling and observed runoff have a large reduction in the 1980s (a reduction of 16.64 mm for the modeling runoff and 21.57mm for the observed runoff).During this period, climate change is the main factor of the runoff reduction, with a contribution of 77.15%.(3) Compared with the undisturbed stage (1956-1979), the runoff reduction of the disturbed stage (1980-2010) caused by climate change accounts for 64.74%, which indicates that the climate change is the main cause.Human activities contribute 35.26% for the runoff reduction.With the improvement of humanity's ability to change the environment, the impact of human activities on basin runoff is gradually enhanced while the effect of climate change is gradually weakened.This can be found from Table 8; that from the 1980s to 2000s the contribution of human activities increases from 22.86% to 43.14%, while that of climate change decreases from 77.15% to 56.86%.The results of the impact of climate change and human activities on runoff shown in Table 8 are similar as that of Liu et al.

Table 8 .
Impact of climate change and human activities on the runoff in Luanhe sub-basin.

Table 9 .
Impact of climate change and human activities on runoff of each sub-basins.

Table 10 .
Test results of the reconstructed annual runoff change trend of each hydrological stations.