Analysis on the Variation of Hydro-Meteorological Variables in the Yongding River Mountain Area Driven by Multiple Factors

Abstract

In recent decades, strong human activities have not only brought about climate change including both global warming and shifts in the weather patterns but have also caused anomalous variations of hydrological elements in different basins all around the world. Studying the mechanisms and causes of these hydrological variations scientifically is the basis for the management of water resources and the implementation of ecological protection. Therefore, taking the Yongding River mountain area as a representative watershed in China, the changes of different observed and simulated hydro-meteorological variables and their possible causes are analyzed on an inter-annual scale based on ground based observations and remotely sensed data of hydrology, meteorology and underlying surface characteristics from 1956 to 2016. The results show that the annual natural runoff of Guanting hydrological station in the main stream of the Yongding River, Cetian hydrological station and Xiangshuibao hydrological station in the tributary of the Yongding River all have a significant decreasing trend and abrupt changes, and all the abrupt change points of the annual natural runoff series of the three hydrological stations appear in the early 1980s. On the inter-annual scale, the water balance model with double parameters is unable to effectively simulate the natural surface runoff after the abrupt change points. The annual average precipitation after the abrupt change points decreases by no more than 10%, compared with that before the abrupt change points. However, the precipitation from July to August, which is the main runoff-production period, decreases by more than 25%, besides the intra-annual temporal distribution of precipitation becoming uniform and a significant decrease in effective rainfall, which is the source of the runoff. Meanwhile, the NDVI in the basin show an increasing trend, while the groundwater level and land water storage decrease significantly. These factors do not lead only to the continuous reduction of the annual natural runoff in the Yongding River mountain area from 1956 to 2016, but also result in significant changes of the hydro-meteorological relationship in the basin.

1. Introduction

Sustainable water supply is the basis of the ecosystem and human survival, and the climate system and characteristics of the underlying surface are the dominant factors for water resources in a basin [1,2]. In the past few decades, strong human activities have been rapidly changing the natural environment in different regions, having a profound impact on the water resources [3,4,5]. Due to the strong human activities, hydrological elements in many regions show an increasing or decreasing trend. This leads to the sharp decrease in river runoff in some arid areas, which will bring challenges to the management of water resources and ecological protection in the basin at present or in the future [6,7,8,9] Therefore, it is of great significance to analyze the variability and causes for the changes in different hydro-meteorological variables in order to understand the climate change induced changes of the hydrological cycle, as well as to improve the hydrological prediction, water-resource evaluation methods, water-resource planning and management, etc.

Climate change has a large influence on runoff and other hydro-meteorological variables by changing precipitation, potential evapotranspiration and snowmelt through changes in air temperature [10,11], while human activities affect runoff through land use/cover change, urbanization process, reservoir operation and direct exploitation of surface water and groundwater [12,13,14]. At present, with the increasing scarcity of water resources, especially in arid and semi-arid areas, the impact of climate change and land use/cover change on hydrological elements is more and more obvious, which has attracted the attention of hydrological researchers, decision makers and scholars [15]. The International Association of Hydrological Sciences (IAHS) has taken the evolution and attribution of river runoff under natural influences and human activities as one of its key research topics in recent years [16]. A large number of literatures are devoted to the study of hydrological variation and attribution identification of the runoff variation in different basins or regions. Pankaj summarized the main methods of attribution analysis for the runoff variation at present stage from four aspects of hydrological modeling approaches, conceptual approaches, analytical approaches and methods based on field observations, and then sorted out the conclusions of the runoff change contribution rate in some regions of the world in relevant literature [17]. Kuk conducted an attribution analysis on runoff changes through four methods and found that human activities are the main reasons for the runoff reduction in each state, which are mainly due to urbanization, agricultural development, rainwater utilization measures and water conservancy project construction in Indiana, New York, Arizona and Georgia in USA [18]. Wang quantitatively analyzed the causes of runoff change of main rivers in China and concluded that human activities are the main causes of the runoff reduction in northern China, while the runoff change in southern China is mainly affected by climate change, based on the RCCC-WBN model [19]. Yang analyzed that the significant decrease in water resources in the Yellow River Basin of China in the past half century is mainly due to the decrease in precipitation, the increase in evapotranspiration and the increase in irrigation water consumption [20]. Gao analyzed the attribution of runoff change in the Loess Plateau by the Budyko theory and concluded that the sensitivity of hydrological response to land use is higher than that to climate change, and the contribution rate of land use to runoff reduction in most areas of the Loess Plateau is about 60% [21]. Bao quantitatively identified the runoff change in Haihe River Basin based on the four element driving theory, and the results showed that the contribution rates of climate change and human activities to runoff reduction are 32–38% and 62–68%, respectively [22], which is similar to Xu’s research conclusion [23]. Although human induced changes observed for different elements of the hydrological cycle and their attribution have been widely analyzed, it is still not fully understood about the variation and its causes of some hydro-meteorological variables. At present, the researchers mainly focus on the diagnosis and attribution analysis of the observed runoff, but natural runoff which is obtained by adding the water amount utilized by human beings to the observed runoff is less discussed in literatures. As the key factor to characterize the precipitation redistribution and water production capacity, natural runoff is very important for analyzing the change of water resources from the source. Moreover, even though many papers have combined analyzing the relationship of runoff and precipitation, evapotranspiration, land use and other elements, few discussed the in-depth reasons of runoff change sufficiently with a variety of observed data and remotely sensed data, which restricts the systematic understanding and correct attribution of the driving mechanism for the runoff change.

The Yongding River is one of the seven major water systems of the Haihe River and the mother river of Beijing. Maintaining the health of the Yongding River and creating a flowing, clean and green Yongding River are the important foundation for ensuring the water safety in the capital and promoting the coordinated development of Beijing, Tianjin and Hebei. However, since the 1980s, the situation of the Yongding River has changed, and the problems of river cut-off and dry-up have become increasingly prominent, and the ecological corridor function of the river has been lost [24,25]. The development and the changes observed for different hydro-meteorological variables for the watershed of the Yongding River is representative for Northern China. Considering the declining trend of the surface runoff in the Yongding River mountain area (YRMA) since 1980s, the Haihe River Water Conservancy Commission proposed a “recover” treatment method for the natural runoff seriously affected by underlying surface [26]. Wang used the geomorphology-based hydrological model (GBHM) to analyze the impact of the land use/cover change on the runoff production in the Yongding River Basin (YRB). The results showed that the increase in woodland and grassland would reduce the runoff production in the basin [27]. Zhang used the SWAT model to simulate the monthly runoff process in the upper reaches of the Yongding River from 1957 to 2000 and pointed out that climate change is the dominant factor for the runoff change in the basin from 1980 to 2000 [28]. Based on the Budyko hypothesis, Hou used the climate elasticity coefficient method to analyze the driving factors for the runoff change in the YRB and pointed out that the human activities have been the main factors for the significant decline of the observed runoff in the YRB since the 1980s [29]. It can be seen that these studies mainly focus on the observed runoff changes of the Yongding River, but the conclusions are different or even contradictory. Considering that there still exist many factors not yet properly understood, as well as the important geographical location, resource ecology and historical and cultural significance of the Yongding River, the aim of this paper is to analyze the changes in natural runoff and other hydro-meteorological variables for the YRMA in recent decades, as well as to understand the causes of the observed changes using a large variety of ground-based and remotely sensed observations. First, the inter-annual variation of surface runoff, precipitation and evaporation capacity in the YRMA are explored. Then, the changes in the relationship between precipitation, vegetation cover and runoff are analyzed, also by utilizing a water balance model. Finally, the causes of the observed changes are investigated by analyzing the temporal distribution of the precipitation, change of vegetation cover and groundwater exploitation.

2. Materials

2.1. Diagnosis Methods for the Variability of Hydrological Elements

Time series of different hydrological variables are often non-stationary due to trend and abrupt changes. At present, there are several statistical tools to analyze the variability of these time series, however, most of them cannot detect the trend and abrupt changes at the same time, and the results obtained by different methods are not the same. In order to improve the robustness and reliability of the diagnosis methods for the variability of hydrological elements, a process called “Hydrological Alteration Diagnosis System” based on multiple methods is proposed by Xie [30], including three steps, that is preliminary diagnosis, detailed diagnosis and comprehensive diagnosis.

Preliminary diagnosis: preliminary diagnosis checks the variability of a sequence according to the Hurst method [31] to determine whether there is an abrupt change in the sequence.

Detailed diagnosis: if the preliminary diagnosis confirms that there is a change in the sequence, the detailed diagnosis will be used for further analysis, that is, multiple statistical tests are used to further diagnose the trend and abrupt changes of the sequence. Among them, slide F-test, slide T-test, L-H test, ordered clustering test, R/S test and optimal segmentation are used for analyzing abrupt changes in the time series, while Pearson, Spearman and Kendall correlation coefficients are used for trend analysis [30]. If a test method judges that the trend or abrupt change of the sequence is significant, the test result is recorded as +1; if it is found to be non-significant, then the test result is recorded as –1; if it cannot be judged, then the test result is recorded as 0.

Comprehensive diagnosis: comprehensive diagnosis means comprehensively analyzing the detailed diagnosis results. When the sum of all the results obtained by the trend test or abrupt change test in detailed diagnosis is greater than 1, the change of the sequence is considered significant, otherwise it is not significant. If the abrupt change points as a results of the detailed diagnosis are different, the point of year with the most occurrences can be regarded as the year in which the abrupt change of natural runoff occurred.

After finding the abrupt change points in the series of hydrological elements, the time series before the abrupt change points can be regarded as the reference period. While the series after the abrupt change points can be regarded as the variation period.

2.2. Water Balance Model for River Basins

Basin water balance model is an empirical or semi-empirical hydrological model with the advantages of a simple structure and few parameters. At present, there are a series of water balance models, which are widely used in the simulation of the hydrological process and the attribution analysis of the runoff variability from daily to multi-year time-scales in basins with different topography and climates. In this study, the water balance model with double parameters (WBM-DP) [32] is used to simulate the annual runoff of the YRB. The WBM-DP is composed of three equations. The first is the water balance equation, which is used to describe the relationship between hydro-meteorological variables at adjacent times. The second is the coupled water–energy balance equation, used to describe the actual evapotranspiration under the control of the evaporation capacity and the available water supply such as precipitation. The last is the storage–discharge equation, used to describe the relationship between the runoff at the outlet section and the available water supply such as precipitation. The principle and methods of the model are described as follows.

Water balance equation: taking the water storage of the basin as the state variable, the water balance relationship at the beginning and end of the period is expressed as follows.

$$S_t=S_{t-1}+P_t+I_t-E_t-Q_t$$

(1)

where $S_{t-1}$ and $S_t$ are the water storages of the basin at the beginning and end of a certain period, respectively; $P_t$, $E_t$ and $Q_t$ represent the precipitation, actual evapotranspiration and runoff in this period, respectively; $I_t$ is the inflowing water in this period.

The coupled water–energy balance equation is shown as follows. The actual evapotranspiration in the basin is calculated by the Choudhury–Yang formula [33].

$$E_t=\frac{E{P}_t{H}_t}{{(H_t^n+E{P}_t^n)}^{1/n}}$$

(2)

where EP is the atmospheric evapotranspiration capacity in a certain period; H is the amount of available water for evapotranspiration in this period; n is a parameter reflecting the influence of underlying surface characteristics on evapotranspiration; the larger n is, the more sensitive EP is to the change of P, especially in arid area. For a basin, H is the sum of the water storage (S) in the basin at the beginning of this period, the precipitation (P) during this period and the inflow water (I), i.e.,

$$H_t=P_t+I_t+S_{t-1}$$

(3)

The storage–discharge equation adopts the hyperbolic tangent function to express the relationship between the runoff and the amount of water deducting evapotranspiration.

$$Q_t=(H_t-E_t)\tanh (\frac{H_t-E_t}{SC})$$

(4)

where SC is the storage–discharge coefficient of the basin. The value of SC is small in humid area and large in arid area.

By combining Equations (2)–(4), the times series of annual runoff, evapotranspiration and water storage can be obtained on the annual time scale. The WBM-DP can be used to simulate the runoff and evapotranspiration processes on an annual or monthly scale. In the paper, the only two parameters of the model were calibrated through the shuffled complex evolution method developed by the University of Arizona [34]. The annual natural runoff in reference period and in variation period was simulated separately, with the model parameters calibrated for each period separately. Then the water balance error and Nash–Sutcliffe efficiency can be obtained:

$$Re=1-{\displaystyle \sum _{i=1}^{i=N}{Q}_t^s}/{\displaystyle \sum _{i=1}^{i=N}{Q}_t^o}$$

(5)

$$NSE=1-{\displaystyle \sum _{i=1}^{i=N}{(Q_t^s-Q_t^o)}^2/}{\displaystyle \sum _{i=1}^{i=N}{(Q_t^o-{\overline{Q}}_{}^o)}^2}$$

(6)

where Re is water balance error; NSE is the Nash–Sutcliffe efficiency; $Q_t^o$ and $Q_t^s$ are the observed and simulated runoff in the t^th month, respectively; ${\overline{Q}}_{}^o$ is the average value of $Q_t^o$ in the period; N is the total number of months in the period.

2.3. Traditional Runoff Restoration Method

For watersheds influenced by human activities caused obvious changes in the runoff, the runoff restore should be performed. The runoff restore is calculated by using the traditional runoff restoration method [35]. The natural runoff of the hydrological stations of the YRMA can be obtained by adding the water amount utilized by human beings to the observed runoff. The formula of natural runoff calculation is as follows.

$$W={\displaystyle \sum _{i=1}^7{W}_i}$$

(7)

where W is the natural runoff after restoration; W₁ is the observed runoff; W₂ is the agricultural irrigation water consumption; W₃ is industrial water consumption; W₄ is domestic water consumption; W₅ is the variation amount of the storage capacity of water storage project; W₆ is the runoff influence by soil and water conservation measures; W₇ is the runoff inter basin water transfer. Section 3.2. gives a detailed overview on the data necessary for the calculations.

2.4. The Transfer Matrix of Land Use Types

The land-use transfer matrix is an effective tool to represent the transfer area of different land use types in a certain period and can present the overall transfer situation of land use [36]. On the ArcGIS platform, the land use data in 1980 and 2015 are statistically analyzed to obtain the spatio-temporal pattern variations of land use from 1980 to 2015. The transfer-out area of one land use type mainly reflects the area of a certain land use type transforms into another land use type from 1980 to 2015. While the meaning of transfer-in area of one land use type reflects the area of a certain land use type transformed into by another land use type from 1980 to 2015. The proportion of different land use types refers to the proportion of different land use types in the total land area in a certain year. The calculation formulas is as follows.

$$P_i=\frac{T_{ij}}{{\displaystyle \sum _{j=1}^n{T}_{ij}}}\times 100\%$$

(8)

In the formulas, P_i is the proportion of different land use types in the total land area; T_ij is the land area of one land use type; i is the certain year of 1980 or 2015; j is the type of land use, including cultivated land, woodland, grassland, water area, construction land and unused land.

2.5. Slope Analysis Method

In this study, an ordinary least squares (OLS) method based on grid size is used to analyze the climate and vegetation trend on interannual time scales. The calculation equation for trend analysis is:

$${\theta }_{slope}=\frac{n\times {\displaystyle \sum _{i=1}^{n}i\times Y_i-{\displaystyle \sum _{i=1}^{n}i{\displaystyle \sum _{i=1}^n{Y}_i}}}}{n\times {\displaystyle \sum _{i=1}^n{i}^2-}{({\displaystyle \sum _{i=1}^{n}i})}^2}$$

(9)

where ${\theta }_{slope}$ is the slope of the linear regression, indicating the changing trend of NDVI; n represents the total number of the studied years, and it is 35 in this study; and Y_i is the value of NDVI in the ith year.

3. Study Area and Data

3.1. Study Area

The YRB is located in the Middle East of China and the northwest of the Haihe River Basin (111°58′ E–116°22′ E, 38°50′ N–41°16′ N), flowing through five provincial administrative regions of Inner Mongolia, Shanxi, Hebei, Beijing and Tianjin. Its catchment area is about 47,000 km², accounting for 14.7% of the Haihe River system. In 2016, the total population in the basin was about 8.13 million, and the urbanization rate was 58.3%. The primary tributaries of the YRB upper reaches are the Sanggan River and the Yanghe River, which are called Yongding River after the confluence. The YRB has a semi-humid and semi-arid monsoon climate, with an annual average precipitation of 360–650 mm. The precipitation from June to August accounts for about 80% of the annual average precipitation.

Guanting hydrological station is located at the downstream of Guanting Reservoir. The basin in the upper reaches of the Guanting Reservoir is a mountainous area, Yongding River mountain Area (YRMA), with an area of 45,100 km², accounting for 95.8% of the whole basin, which is the main runoff-producing area of the basin. The Cetian hydrological station and Xiangshuibao hydrological station are located at the downstream of Cetian Reservoir of Sanggan River and Xiangshuibao Reservoir of Yanghe River, with the catchment area of 17,100 and 14,500 km², respectively (Figure 1).

3.2. Data

In this paper we utilize a large dataset of various ground-based and remotely sensed observations, such as meteorology and hydrology, exploitation and utilization of water resources, land use, vegetation cover and land water storage. The following data are used in the study.

(1)

Meteorological data: the daily surface observation data (monthly precipitation, maximum/minimum temperature, relative humidity, average wind speed, sunshine) at 14 meteorological stations (Beijing, Huailai, Zhangjiakou, Zhangbei, Jining, Yuxian, Yangquan, Datong, Youyu, Wutaishan, Youyu, Wuzhai, Hequ, and Yuanping) in the study area from 1956 to 2016, as shown in Figure 1 (http://cdc.cma.gov.cn/home.do/, accessed on 5 August 2020).

(2)

Areal precipitation data: the areal precipitation data in the YRMA from 1956 to 2016 based on the precipitation data of the meteorological stations shown in Figure 1, obtained by the Thiessen polygon method.

(3)

Potential evapotranspiration data: the potential evapotranspiration at each station is calculated by the Penman–Monteith [37] formula based on the meteorological data of meteorological stations, and the averaged surface potential evapotranspiration in the study area is obtained by the Thiessen polygons method.

(4)

Observed runoff data: according to the Hydrological Data of Haihe River Basin of Annual Hydrological Report of China [38], the observed runoff of Cetian, Xiangshuibao, Guanting hydrological stations from 1956 to 2016 are selected.

(5)

Water resources exploitation and utilization data: according to the water resource communiques of Shanxi, Beijing, Hebei and evaluation of water resources (http://slt.shanxi.gov.cn/zncs/szyc/szygb/; http://slt.hebei.gov.cn/dynamic/search.jsp; http://swj.beijing.gov.cn/zwgk/szygb/, accessed on 20 October 2020) and Evaluation of water resources in Haihe River Basin [26], the exploitation and utilization data of water resources in YRMA are selected, including the total amount of water supply and consumption, the consumption of different water users and the water supply from surface and underground.

(6)

Land use data: the remotely sensed data of land use in the YRB in 1980 and 2015 are obtained, mainly including crop land, forest land, grassland, water area, land for industrial, construction land, as well as unused land, with the data resolution of 1 km. The land use data are based on Landsat TM/ETM remotely sensed image in 1980 and 2015, interpreted by ArcGIS platform. (https://www.resdc.cn/data.aspx?DATAID=184/, accessed on 8 November 2020).

(7)

Vegetation data: the monthly normalized vegetation index (NDVI) data from 1982 to 2016 are obtained. The monthly NDVI data from 1982 to 2015 is from the Global Inventor Modeling and Mapping Studies (GIMMS) NDVI data set, with a resolution of 8 km and a temporal resolution of 15 days (https://ecocast.arc.nasa.gov/data/pub/gimms/3g.v1/, accessed on 14 November 2020). The monthly NDVI values were obtained by using the maximum value composite (MVC) method [39]. The monthly NDVI data in 2016 is from the website of resources and environment science and data center, which is based on the SPOT/VEGETATION, MODIS and other satellite remotely sensed images, with a resolution of 1 km and a temporal resolution of one month (https://www.resdc.cn/DOI/doi.aspx?DOIid=50, accessed on 14 November 2020)). It could be used after unifying the resolution of GIMMS NDVI and MODIS NDVI data [40].

(8)

GRACE data: The GRACE data from April of 2002 to December of 2016 in the YRMA are used to study the long-term spatiotemporal changes of terrestrial water storage anomaly (TWSA) [41]. The TWSA data is released by the University of Texas with a spatial resolution of 0.25° (http://www2.csr.utexas.edu/grace/RL06_mascons.html/, accessed on 20 December 2020). Among them, 33 months of data are missing, which are interpolated by the cubic spline interpolation method [42]. According to the topographic map of the basin, the area weighted average is carried out, and the monthly average land water storage change of the basin is obtained.

4. Results and Discussion

4.1. Changes of the Annual Runoff and other Hydrological Elements

The inter-annual variations of natural and observed runoff in the YRMA from 1956 to 2016 are shown in Figure 2. It can be seen that there is an obvious decreasing trend in the annual observed runoff of the three hydrological stations downstream of the Cetian, Xiangshuibao, Guanting Reservoirs in the YRMA. The variabilities of the annual natural runoff series of the Guanting hydrological station, the Cetian hydrological station and the Xiangshuibao hydrological station are tested at the 95% confidence level by using the diagnosis system for hydrological variability. The results in

Table 1. Diagnostic results of annual natural runoff variation in the mountain area of Yongding River from 1956 to 2016.

Hydrological Stations			Cetian	Xiangshuibao	Guanting
Diagnostic Process
Preliminary diagnosis		Hurst index	0.982	1.000	0.975
		Overall variation	Strong	Strong	Strong
Detailed diagnosis	abrupt	Slide F-test	2010(+)	2014(+)	2000(+)
		Slide T-test	1982(+)	1983(+)	1983(+)
		L-H test	1982(0)	1983(0)	1983(0)
		Ordered clustering test	1982(0)	1983(0)	1983(0)
		R/S test	1982(0)	1986(0)	1983(0)
		Optimal segmentation	1982(+)	1983(+)	1998(+)
	trend	Pearson correlation test	(+)	(+)	(+)
		Spearman correlation test	(+)	(+)	(+)
		Kendall correlation test	(+)	(+)	(+)
Compre -hensive diagnosis	abrupt	Abrupt change point	1982	1983	1983
		Comprehensive significance	2(+)	2(+)	1(+)
	trend	Trend Comprehensive significance	decline 3(+)	decline 3(+)	decline 3(+)
Conclusion			The natural runoff of all hydrological stations showed significant decreasing trend and abrupt changes

show that the annual natural runoff of the three hydrological stations all has a significant decreasing trend and abrupt changes, and the abrupt change points are in 1983, 1982 and 1983, respectively. Considering that the abrupt change points of natural runoff in the three hydrological stations are almost the same, 1983 is regarded as the abrupt change point for the hydrological stations, and the periods of 1956–1982 and 1983–2016 are recorded as the reference period and the variation period, respectively. The mean values of annual natural and observed runoff of the Guanting hydrological station, the Cetian hydrological station and the Xiangshuibao hydrological station before and after the abrupt change point are shown in

Table 2. Annual mean values of meteorological and hydrological elements in the basins of Yongding River mountain area in different periods.

Basins	Period	P/mm	EP/mm	EP/P	R1/mm	R2/mm
Cetian Reservoir	1956–1982	422.9	1012.8	2.39	37.8	25.1
	1983–2016	384.6	988.3	2.57	22.5	3.7
Xiangshuibao Reservoir	1956–1982	427.0	1033.1	2.42	51.0	34.9
	1983–2016	389.0	971.5	2.50	21.9	10.2
Guanting Reservoir	1956–1982	432.2	1004.4	2.32	43.4	28.5
	1983–2016	395.5	962.9	2.43	21.2	6.0

Remarks: R₁ and R₂ in the table represent annual natural runoff and observed runoff respectively.

, as well as the mean values of annual precipitation and evaporation capacity at reference period and variation period. Since the 1980s, the surface natural runoff of the above three hydrological stations in the YRMA has declined significantly, indicating that the water production capacity of the basin has greatly decreased. The annual average natural runoff of the Guanting hydrological station is 43.4 mm in the reference period and 21.2 mm in the variation period, decreased by 51.1%. For the Cetian hydrological station and the Xiangshuibao hydrological station, they are 37.8 mm, 22.5, 40.5% and 50.1 mm, 21.9 mm, 57.1%, respectively. Compared with the natural runoff, the reduction of the observed runoff in each hydrological station is more severe. From the reference period to the variation period, the reduction ratios of annual average natural runoff of Guanting hydrological station, Cetian hydrological station and Xiangshuibao hydrological station are 78.9%, 85.3% and 70.8%, respectively.

The annual average precipitation in the Guanting Reservoir basin is 432.2 mm in the reference period and 395.5 mm in the variation period, with a decrease in about 8.5%. For the basins of the Cetian Reservoir and the Xiangshuibao Reservoir, the reduction rates are 9.1% and 8.9%, respectively. At the same time, the potential evapotranspiration in the basin of each hydrological station in the variation period also decreases in different degrees compared with the reference period. The annual potential evapotranspiration in the Guanting Reservoir basin is 1004.4 mm in the reference period and 962.9 mm in the variation period, with a decrease in about 4.1%. The reduction rates in the basins of Cetian and Xiangshuibao reservoirs are 2.4% and 6.0%, respectively. As the reduction of annual precipitation is higher than that of potential evapotranspiration, the average drought index (the ratio of annual potential evapotranspiration to precipitation) of the YRMA in the variation period increases compared with that in the reference period, and the annual average drought indexes in the basins of the Guanting Reservoir, the Cetian Reservoir and the Xiangshuibao Reservoir increase by 0.11, 0.18 and 0.08, respectively.

4.2. Variations of Hydro-Meteorological Relationships in the Basin

The variation of the hydrological elements in the YRMA from 1956 to 2016 does not show only the sharp decrease in the natural runoff, but also the significant change of the hydro-meteorological relationship and the coupled water–energy balance on an annual scale.

Figure 3 shows the scatter plot and linear regression relationship between the annual precipitation and natural runoff in the main and tributary stream sections in the YRMA in the reference period and variation period. It can be seen that the annual precipitation–runoff scattered points in the variation period are basically below those in the reference period, and the slope of linear regression line in the variation period is significantly lower than that in the reference period. Thus the runoff production capacity of the YRMA decreases greatly under the same annual precipitation conditions. The average annual runoff coefficients at the Guanting Reservoir station before and after the abrupt change point are 0.10 and 0.05, respectively, and that at the Xiangshuibao Reservoir station are 0.12 and 0.06, respectively. The relative reductions of the natural runoff of the Guanting and Xiangshuibao hydrological stations are both close to 50%. For the Cetian hydrological station, they are 0.09, 0.06 and 33.3%, respectively. This means that before the 1980s, about 10% of the annual precipitation in the YRMA could be converted into natural water resources. Since the 1980s, this proportion has reduced to about 5%. At the same time, the scatter of annual precipitation–runoff of the three hydrological stations in the variation period is more dispersed than that in the reference period. The determination coefficient of corresponding linear regression equation of annual precipitation–runoff decreases from 0.57, 0.63 and 0.59 in the reference period to 0.10, 0.17 and 0.07 in the variation period respectively, which means that the contribution of the annual precipitation to the natural runoff is greatly reduced.

The natural water production capacity of the YRMA is greatly reduced, and the reduction of natural runoff is significantly higher than the annual precipitation. This reflects that the natural runoff converted from precipitation supply is decreasing, while the water consumption is increasing. It shows that the relationship among precipitation, runoff and evapotranspiration has changed, which means that the precipitation redistribution and the relationship of coupled water–energy balance have changed. At the same time, it should be recognized that the change of the coupled water–energy balance in the YRMA does not show only the increase in the water-consumption proportion in the water supply, but also may involve the significant change of the water–energy response process and mechanism. In order to illustrate this problem, the WBM-DP is used to simulate the annual natural runoff of the three hydrological stations. According to the abrupt change points of annual natural runoff in Section 4.1, the models for the reference period and the variation period are separately calibrated. Figure 4 shows the simulation results of the natural runoff at each hydrological station in 1956–1982 and 1983–2016. In the reference period, the accuracy of WBM-DP is acceptable, and the simulation result is good. The water balance error at each hydrological station is close to zero, and the Nash–Sutcliffe efficiency (NSE) coefficients are 0.78, 0.82 and 0.89, respectively. The model parameters are also within the normal limits. This indicates that the relationship between the actual evapotranspiration and natural runoff with the annual precipitation and atmospheric evaporation capacity in the YRMA basically conforms to the Budyko hypothesis and the storage–discharge equation. However, the simulation result of WBM-DP in the variation period is poor, and the NSE values of the three hydrological stations are only 0.46, 0.58 and 0.55, respectively. Especially since the mid-1990s, the difference between the simulated value and the actual value of the natural runoff has been particularly significant. Taking Guanting hydrological station as an example, the WBM-DP fails to accurately simulate the natural annual runoff process after 2000, and the simulated value is obviously higher than the actual value. This indicates that the influencing factors and mechanisms of the natural runoff and evapotranspiration in the basin are more complex during the variation period, so the relationship with the annual precipitation and atmospheric evaporation capacity can no longer be described by the Budyko hypothesis and the storage–discharge equation. Only by changing the parameters of the WBM-DP, it is impossible to adapt to the significant changes of the coupled water–energy balance.

4.3. Driving Factors for the Change of Hydro-Meteorological Relationship

Understanding the dominant driver factors of hydrological change is essential for water resources management. Similar to natural runoff, non-stationary trends also exists in hydrological drivers. Currently, the hydrological drivers of variation in streamflow have been widely discussed all around the world [43,44]. In undisturbed basins, changes in runoff are likely primarily due to climatic factors [45]. While in human-impacted basins, development and land cover change can be the most important drivers for runoff variation [46,47]. The relative importance of hydrologic drivers varies spatially across different types of basins. As a representative arid watershed in Northern China, the variation and its causes of hydro-meteorological variables of the YRB can be analyzed from the aspects of precipitation, land use, water utilization, et al. The detailed analysis and discussion are shown below.

4.3.1. Intra-Annual Temporal Distribution of Precipitation

Precipitation is the basic driving factor for the hydrological process in a river basin, and its spatio-temporal distribution characteristics have a profound impact on the runoff. Previous studies have pointed out that under the premise of a slight decrease in the atmospheric evaporation capacity, the annual average precipitation in the variation period of the YRMA is reduced by less than 10% relative to the reference period, but the annual average natural runoff in the main stream is reduced by more than 40% and that in the tributaries is reduced by more than 50%. Although the elastic response coefficient of annual runoff to precipitation is significantly greater than 1 in many cases [25], it is still intriguing. Since the YRMA is located in North China influenced by the monsoon climate, the precipitation in a year is highly concentrated in the flood season (from June to September), and thus the river runoff is also mainly concentrated in the flood season. The precipitation in the flood season of the YRMA determines the scale of the annual runoff to a considerable extent, so the change of the intra-annual temporal distribution of precipitation and its impact on water production should be paid attention to. Figure 5 shows the intra-annual precipitation distribution in the basins of the Guanting, Cetian and Xiangshuibao reservoir control basin in the reference and change periods. The precipitation difference between the variation period and the reference period is significant mainly from July to August, while there is little difference in other months. From 1983 to 2016, the precipitation in the basins of the Guanting Reservoir, Cetian Reservoir and Xiangshuibao Reservoir from July to August decreases by 26.3%, 27.8% and 26.5%, respectively, compared with that from 1956 to 1982. As a result, the nonuniformity coefficient of monthly precipitation in the three control basins decreases from 1.12, 1.16 and 1.14 to 0.95, 0.96 and 0.96, respectively, and the precipitation concentricity in a year decreases.

It can also be seen from Figure 5 that, except for individual cases, the monthly average natural runoff in the three control basins of the Guanting Reservoir, the Cetian Reservoir and the Xiangshuibao Reservoir during the variation period decreases with varying degrees, compared with that in the reference period. From July to September, the natural runoff of the three hydrological stations decreases by 10.67, 9.68 and 12.79 mm, respectively (the relative reductions are 56.7%, 51.0% and 61.2%, respectively), accounting for 47.9%, 59.2% and 43.9% of the annual reduction, respectively. The nonuniformity coefficient of the natural runoff in the three hydrological stations decrease from 0.67, 0.56 and 0.60 during 1956–1982 to 0.57, 0.55 and 0.50 during 1983–2016, respectively. It shows that the intra-annual distribution of the natural runoff in each hydrological station tends to be even. Obviously, the decline of the natural runoff from July to September in the YRMA has a clear correspondence with the substantial decrease in precipitation from July to August.

Furthermore, the daily precipitation from July to August in the reference period and in the variation period in the YRMA are analyzed, respectively, and then the frequency histogram of days with different-level daily precipitation in the control basin of each hydrological station from July to August is plotted (Figure 6). Compared with 1956–1982, the number of days with daily precipitation less than 0.1 mm, that is, the number of days without rainfall, increases significantly in 1983–2016; the proportion of days with daily precipitation between 0.1 and 2 mm does not change significantly; the ones with daily precipitation in 2–10 mm and 10–50 mm both decrease significantly, while the one with daily precipitation over 50 mm is basically zero. This shows that compared with the reference period, the number of effective precipitation days for the runoff production in the YRB has been significantly reduced during the variation period. The daily average precipitation in the Guanting Reservoir basin decreases from 4.3 mm during 1956–1982 to 2.9 mm during 1983–2016, the variation coefficient decreases from 2.3 to 1.9, and the kurtosis coefficient decreases from 73.3 to 18.3 (

Table 3. Daily precipitation statistics of different meteorological stations from July to August.

Hydrological Station	1956–1982			1983–2016
	Average	Variation Coefficient	Kurtosis Coefficient	Average	Variation Coefficient	Kurtosis Coefficient
Cetian	4.4	2.5	53.4	2.8	2.5	24.6
Xiangshuibao	4.2	2.5	75.3	2.8	2.0	21.4
Guanting	4.3	2.3	73.3	2.9	1.9	18.3

). The other two control basins also show similar trends. It shows that compared with 1956–1982, the daily average precipitation in the YRMA in 1983–2016 decreases significantly, the concentration decreases and the distribution becomes uniform. The changes of the precipitation from July to August in the YRB obviously have an important influence on the runoff production. The decreases of precipitation intensity and concentration will be beneficial for increasing the vegetation interception and reducing the excess runoff by throughfall and surface infiltration. As a result, the rainfall flowing into the vadose zone will move more slowly, but the proportion of returning to the atmosphere through vegetation transpiration and soil evapotranspiration will also increase, so the groundwater infiltration recharge by rainfall may not increase or even decrease. Therefore, the total result is an increase in evapotranspiration and a decrease in runoff production. As a result of the decrease in total precipitation and its uniform distribution, as well as the increase in evaporation in the basins, the natural runoff of the basins could decrease.

4.3.2. Changes of Vegetation Cover

Land use and vegetation cover are important para meters to reflect the characteristics of the underlying surface, and their changes have an important impact on the water cycle and the associated processes. The land use types in the YRMA mainly include crop land, forest land, grassland, water area, construction land and unused land, in which the crop land, forest land and grassland cover most of the land. Figure 7 shows the distributions of land use types in the study area in 1980 and 2015, which are used to show the land use conditions in the reference period and the variation period, respectively. Compared with 1980, the forestland area increases by 2364.4 km² in 2015, with a relative increase in 34.6%. The areas of crop land and grassland respectively decrease by 1610.1 and 1056.8 km², with a relative decrease in 7.7% and 7.5%, respectively. Urban and rural construction land increase by 926.9 km², with a relative increase in 116.2%. According to the transfer matrix of land use area in the YRMA shown in

Table 4. Transfer matrix of land use types in the mountain area of Yongding River from 1980 to 2015 unit: km².

2010	Cultivated Land	Woodland	Grassland	Water Area	Construction Land	Unused Land	Total	P2010
1980
cultivated land	13139.1	2093.1	4160.8	373.1	1111.5	142.5	21,020.1	47.3%
woodland	822.7	4560.1	1347.2	39.5	55.3	16.5	6841.2	15.4%
grassland	4078.7	2378.6	7184.8	137.6	328.1	47.9	14,155.6	31.9%
water area	612.3	100.1	224.3	117.8	71.7	21.5	1147.6	2.6%
construction land	453.7	49.1	137.0	12.3	140.9	4.7	797.6	1.8%
unused land	303.6	24.6	44.7	28.7	17.0	37.0	455.7	1.0%
total	19,410.0	9205.6	13,098.9	708.9	1724.5	270.0	44,417.9	100.0%
P1980	43.7%	20.7%	29.5%	1.6%	3.9%	0.6%	100%	-

, from 1980 to 2015, there is mainly the mutual transfer of crop land, forest land and grassland, as well as the conversion from parts of crop land and grassland into construction land. As a result, the proportion of the area of cultivated land to the whole basin area decreases from 47.3% in 1980 to 43.7% in 2015, but the proportion of the total area of forest land and grassland increased from 47.3% to 50.2%, and that of the construction land increases from 1.8% to 3.9%. The proportion of conversion area of different land use types to the whole basin area is small and varied between 0.4% to 5.3%. Therefore, in general, the land use pattern of the YRMA changed little compared 1980 with 2015.

In order to further study the change law of the vegetation in the YRMA, the NDVI data during 1982–2016 are analyzed. Figure 8 shows the variation of the annual average and maximum NDVI in the YRMA during 1982–2016. The annual average NDVI in the YRMA during 1982–2016 showed a significant increasing trend (at the confidence level of 99%), and it increased by 0.0013 per year, with the minimum value of 0.272 (1985) and the maximum value of 0.341 (2016). The fluctuation range of the spatial change rate of annual average NDVI obtained by the slope analysis method is −0.08% a⁻¹–0.33% a⁻¹. The areas where the annual average NDVI show increasing and decreasing trends account for 95.36% and 4.64% of the total area, respectively. Among them, the area with a faster growth (slope ≥ 2 × 10⁻³ a⁻¹) accounts for about 12.32%. From 1982 to 2016, the annual maximum NDVI in the YRMA showed a significant increasing trend (at the confidence level of 99%), with an average increase in 0.0023 a⁻¹; the minimum value is 0.459 (1983), the maximum value is 0.607 (2016); the fluctuation range of spatial change rate is −0.18% a⁻¹–0.6% a⁻¹. In addition, in only 8.12% of this area the annual maximum NDVI shows a decreasing trend.

Comparing the data of land use in 1980 and 2015 as well as the NDVI data in 1982 and 2015, the changes of NDVI in crop land, woodland and grass are analyzed. The annual average NDVI of crop land, forest land and grassland increases from 0.282, 0.367 and 0.254 in 1982 to 0.327, 0.431 and 0.296 in 2016, with the growth rates of 13.2%, 17.4% and 16.3%, respectively. The annual maximum NDVI of crop land, forest land and grassland increase from 0.488, 0.627 and 0.404 in 1982 to 0.587, 0.801 and 0.506 in 2016, with growth rates of 20.2%, 27.8% and 25.2%, respectively. It can be seen that the NDVI all increases significantly under different land use types, and the most obvious increase in NDVI is in woodland, followed by grassland, and finally crop land. To a certain extent, due to meteorological factors such as precipitation and temperature and human activities such as the vegetation conservation and the soil and water conservation in the basin, as well as the promotion of irrigation measures and the management improvement of agriculture, forestry and grassland, the vegetation cover under various land use types increased, so that the annual and seasonal average NDVIs in the growth period are all increasing.

Therefore, from 1982 to 2016, the vegetation in the YRMA increased significantly, showing a clear trend of greening, which is consistent with the land cover change after the policy implementation of “returning farmland to forest and grassland” in the North China Plain and the economic and social development, and it is also consistent with the development trend of land use [48,49,50].

The vegetation change in the YRMA is simultaneously affected by strong human activities. The changes of precipitation and temperature will inevitably affect the vegetation growth. At the same time, as one of the key areas for soil erosion control in China, the vegetation-cover changes in the study area are also affected by the large-scale soil and water conservation measures. As for the impact of strong human activities on the vegetation-cover change, that which one is dominant is not the focus of this paper. In any case, since the 1980s, the vegetation cover in the YRMA has been increasing significantly in general. The vegetation-cover change causes the change of the underlying surface conditions, which will inevitably affect the hydrological cycle in the basin and then lead to the change of the river runoff, which has been confirmed by a large number of studies [51,52,53,54,55,56,57]. However, whether the increase in vegetation cover increases or decreases the river runoff is closely related to the geographical and climatic characteristics. Some studies have suggested that because of the humid climate in the upper reaches of the Yangtze River, the increase in forest vegetation can strengthen the water conservation capacity, and will not increase the actual evapotranspiration, and then increase the runoff. However, there are more reports in northern China about that the increase in vegetation cover leads to the decrease in water production in the basin. The main physical mechanism is that the increase in vegetation cover strengthens the ability of vegetation canopy to intercept the precipitation, increases the evapotranspiration of forest and grass vegetation and strengthens the consumption of precipitation, soil water and groundwater, which is not good for the production of the surface runoff and subsurface flow. For example, Liu X. et al. pointed out that in the semi-humid or semi-arid loess area, the runoff coefficient and flood yield coefficient will be reduced with the improvement of forest and grass vegetation, and the more arid the climate is, the more the runoff or flood will decrease [54]. It is found that the large-scale vegetation restoration aggravates the water resource consumption in the Loess Plateau, resulting in a decrease in the runoff production. Zhao et al. integrated the Gravity Recovery and Climate Experiment (GRACE) land water storage data and dynamic vegetation model simulation, and it was found that the large-scale artificial vegetation restoration project consumes the land water storage capacity of Maowusu Sand Land in the northern China at the rate of 17 mm per year [58]. The YRMA is adjacent to the middle reaches of the Yellow River, and the climate characteristics in these two areas are similar. Therefore, in this study the annual average NDVI and annual maximum NDVI series are compared with the natural runoff series from 1982 to 2016 in the YRMA (Figure 9), and they show opposite inter-annual variation trends, with a significant negative correlation. This indicates that the change of the vegetation cover may be an important reason for the decrease in the natural runoff in the YRMA.

4.3.3. Groundwater Exploitation

In addition to the intra-annual temporal distribution characteristics of precipitation and the vegetation cover change, another important reason for the decrease in the natural runoff and the change of hydro-meteorological relationship in the YRMA is the high-intensity utilization of water resources. Although the water consumption caused by the economic and social water consumption has been taken into account in calculating the restored natural runoff, because the economic and social water consumption in the Yongding River basin is highly dependent on groundwater resources, the increase in economic and social water consumptions will not only reduce the observed runoff, but also affect the natural water production of the basin due to its impact on the vadose zone and groundwater aquifer.

It can be seen from Figure 10 that the exploitation of groundwater in the YRB is relatively weak before the 1970s, with an annual extraction volume of less than 200 million m³. In general, the groundwater depth is 1.5–3.0 m, with a thin vadose zone. Therefore, in the case of large-scale precipitation, it is easier for the excess infiltration of rainfall to replenish the groundwater through unsaturated aquifer, and it is also possible to saturate the vadose zone. However, since the 1970s, the economic and social water demand in the basin has been increasing. It can be seen from Figure 11 that the social water consumption in the basins of the Guanting, Cetian and Xiangshuibao reservoirs from 1983 to 2016 is significantly higher than that from 1956 to 1982. According to statistics analysis, in the control basins of the three hydrological stations, the annual average economic and social water consumption in the variation period increases by 209 million m³, 154 million m³ and 22 million m³ compared with that in the reference period. At the same time, the amount of surface water resources in the YRMA continues to decline, so the amount of groundwater exploitation is increasing. The annual exploitation amount was 710 million m³ in 1980, which reached 1.2 billion m³ in 1990, further increased to 1.36 billion m³ in 2012, and then began to decline. It was 1.13 billion m³ in 2016. According to the “Haihe River Basin Water Resources Bulletin” and other data during 2010–2016, under current conditions, the proportion of the groundwater supply in the total water supply is as high as 56.5–67.1% in the YRMA. There are varying degrees of groundwater overexploitation in the urban areas of Datong City and Zhangjiakou City and the surrounding areas of Shanyin, Huairen and other cities, and in some irrigation areas the groundwater is also exploited. As a result, the groundwater level continues to sink, and the maximum buried depth reached more than 60 m in the 2010s [59]. Although the time series of the anomaly changes of GRACE land water storage data obtained is relatively short, it can still provide important evidences for the analysis of the groundwater overexploitation in the YRMA. Figure 12 shows that the overall decline trend of land water storages in the YRMA in 2002–2016 is very obvious and sustainable. According to statistical analysis, the land water storage in the YRMA decreased by about 0.71 mm per month and 8.4 mm per year on average. From January of 2002 to December of 2016, the land water storage of the YRMA decreased by 172.5 mm. The decrease in the land water storage in the YRMA is not only closely related to the groundwater exploitation, but also may be affected by the increase in the vegetation cover.

In the YRMA, the long-term overexploitation of groundwater has many influences on the hydro-meteorological relationship. First of all, the drastic decline in the groundwater level means that the thickness of the current vadose zone may be several times or even tens of times thicker than that before 1980s, which greatly improves the water holding capacity of the vadose zone. It is difficult for the vadose zone to be filled with water even under the condition of large-scale rainfall. In addition, it is difficult to replenish the groundwater through infiltrations when the regional precipitation intensity is decreasing, and therefore the precipitation mostly stays in the vadose zone and is consumed through the vegetation evapotranspiration or soil evapotranspiration. Secondly, the continuous decline of the groundwater level will change the recharge and discharge relationship between groundwater and water in river channels, and the lateral recharge of groundwater to the river runoff will weaken, disappear, or even be recharged by the river runoff. To this end, in general, the long-term overexploitation of groundwater will increase the proportion of water consumption, resulting in the reduction of natural runoff production. Zhang believes that the groundwater exploitation is the main factor for the runoff reduction in the YRB [59].

5. Conclusions

This study analyses the changes of natural runoff in the main and tributaries of the Yongding River from 1956 to 2016, reveals the change of the hydro-meteorological relationship in the basin around the 1980s, and comprehensively analyzes the important driving mechanism leading to the significant decline of the natural runoff and the change of the hydro-meteorological relationship from three aspects, including the intra-annual temporal distribution of precipitation, the change of vegetation cover and the groundwater exploitation. The main conclusions are as follows.

(1)

From 1956 to 2016, the annual natural runoff of the Guanting hydrological station located in the main stream, Cetian hydrological station and Xiangshuibao hydrological station located in the tributary of Yongding River mountain area all had significant decreasing trends and abrupt changes. The abrupt change points were in 1983, 1983 and 1982, respectively. Compared with the reference period, the annual average precipitation in the basins of the three hydrological stations during the variation period decreased by less than 10%, the annual average potential evapotranspiration decreased slightly, but the annual average natural runoff decreased sharply, and the relative reduction was more than 40%.

(2)

The variation of the hydrological elements in the Yongding River mountain area from 1956 to 2016 presented a sharp decrease in the natural runoff and change of the relationship between the hydro-meteorological variables since 1980s. With the precipitation–runoff relationship becoming weaker, the WBP-DP based on the Budyko hypothesis and storage–discharge equation cannot adapt to the significant changes of the coupled water–energy balance in the basin and cannot effectively describe the complex relationship of hydro-meteorological variables on an annual scale.

(3)

Compared with 1956–1982, the intra-annual temporal distribution of precipitation in the Yongding River mountain area underwent important changes during 1983–2016. In the main runoff generation period from July to August, the precipitation amount dropped significantly, the concentration of precipitation distribution decreased significantly, and the number of effective precipitation days decreased significantly. This is one of the important reasons for the decrease in the natural runoff and the change of hydrological responses.

(4)

From 1980 to 2015, the overall land use pattern in the Yongding River mountain area did not change much, but the NDVI had a significant increasing trend, and there was a significant negative correlation between the NDVI and the annual natural runoff. The evapotranspiration increased and the consumption of underground water storage caused by the increase in the vegetation cover was likely to be another important reason for the substantial reduction of the natural runoff since the 1980s.

(5)

Due to the long-term and large-scale exploitation of the groundwater in the Yongding River mountain area, the thickness of vadose zone may increase by several times or even tens of times compared with that before the 1980s. The groundwater level in the basin has drastically dropped and the land water storage has been continuously consumed. It has changed the mechanism of runoff production and the relationship between groundwater recharge and discharge, and also has a significant impact on the runoff decreasing and the adjustment of the hydro-meteorological relationship in the YRB.

Under the combined influences of climate change and human activities, the sharp natural-runoff decrease and the change of the hydro-meteorological relationship in the Yongding mountain area are remarkable. The continuous and substantial decline of natural runoff indicates that the hydrological cycle in the basin has been significantly changed, which intensifies the contradiction and competition between the economic-social water use and the ecological water use, making the sustainable utilization of water resources and ecological stability recovery in the basin face severe challenges. In order to realize the coordinated benign evolution of the water resources, economy, society and ecological environment system, we should strive to build a benign water cycle system in the comprehensive management of the river basin, including strengthening the total water consumption control, improving the level of efficient and intensive utilization of water resources, compressing the groundwater exploitation and optimizing the structure of water supply sources. In addition, we should systematically evaluate the hydrological and ecological effects of large-scale vegetation restoration, reasonably allocate vegetation restoration measures and strengthen the natural vegetation restoration.