Attribution Analysis of Climate and Anthropic Factors on Runoff and Vegetation Changes in the Source Area of the Yangtze River from 1982 to 2016

Abstract

Analyzing the temporal variation of runoff and vegetation and quantifying the impact of anthropic factors and climate change on vegetation and runoff variation in the source area of the Yangtze River (SAYR), is of great significance for the scientific response to the ecological protection of the region. Therefore, the Budyko hypothesis method and multiple linear regression method were used to quantitatively calculate the contribution rates of climate change and anthropic factors to runoff and vegetation change in the SAYR. It was found that: (1) The runoff, NDVI, precipitation, and potential evaporation in the SAYR from 1982 to 2016 all showed an increasing trend. (2) The mutation year of runoff data from 1982 to 2016 in the SAYR is 2004, and the mutation year of NDVI data from 1982 to 2016 in the SAYR is 1998. (3) The contribution rates of precipitation, potential evaporation and anthropic factors to runoff change of the SAYR are 75.98%, −9.35%, and 33.37%, respectively. (4) The contribution rates of climatic factors and anthropic factors to vegetation change of the SAYR are 38.56% and 61.44%, respectively.

1. Introduction

The source area of the Yangtze River (SAYR) is located in the hinterland of the Qinghai Tibet Plateau, which is an important ecological barrier in China [1]. It has the ecological functions of maintaining the ecological security of the Yangtze River Basin [2]. At the same time, it is also one of the most sensitive and fragile regions in the global ecological environment [3,4]. Due to the joint impact of climate change and anthropic factors (indiscriminating felling of trees, disafforestation, overgrazing, mining, and digging), the ecological environment of the SAYR is deteriorating, and the function of water conservation is weakening [5,6]. Additionally, the problems of grassland degradation, wetland shrinkage, glacier retreat, soil erosion, and biodiversity decline are becoming increasingly prominent. These issues highlight the threats to the water resources and ecological security of the whole basin.

The impact of climate change has been embodied in the components of the hydrological cycle, including precipitation [7], evapotranspiration [8], and runoff [9]. Most of the literature has focused on precipitation and evapotranspiration worldwide [10,11], and runoff has taken less attention. Nevertheless, runoff remains one of the most important water resources and has a great influence on the formation of geomorphology, the development of soil, and the growth of plants [12]. It is an important condition for the regional industrial and agricultural water supply and plays an important role in the development of social economy.

In recent years, under the joint influence of climate change and anthropic factors, the problems of climate, ecology, and hydrology in the SAYR have been a cause for concern. Therefore, numerous programs have been carried out to improve the eco-hydrological environment, among which Grain for Green proposed by the Chinese government in 1998 is the most well-known. In this context, monitoring the dynamic changes of runoff and vegetation and assessing the impacts of anthropic factors and climate change on runoff and vegetation change are of great practical value for formulating rational adaptive management countermeasures [13,14,15,16,17,18,19].

Previous studies have focused on the relationship between runoff and vegetation change and its influencing factors in the SAYR. Chen et al. analyzed the temporal and spatial change of vegetation cover in the SAYR from 1982 to 2003 and analyzed the influence of topography and human factors on the vegetation change [20]. Yao et al. analyzed the spatiotemporal variation characteristics of vegetation net primary productivity in the SAYR from 1959 to 2008 and analyzed the relationship between vegetation change and climate factors (temperature, relative humidity, precipitation, wind speed, and sunshine hours) [21]. Liu et al. analyzed the spatiotemporal variation characteristics of vegetation coverage in the SAYR from 1997 to 2012 and believed that the improvement of vegetation status in the SAYR benefited from the joint influences of climate humidification and ecological engineering [22]. Qian et al., Li et al., and Luo et al. studied the variation regularities of runoff in the source area from 1957 to 2009, 1961 to 2011, and 1961 to 2016, respectively, and pointed out that the annual surface water resources in the SAYR all showed an increasing trend [23,24,25]. Xu et al. analyzed the runoff variation characteristics and the degree of human impact in the SAYR from 1956 to 2004 [26]. Liu et al. separated the influences of climate and anthropic factors on runoff in the SAYR based on the Budyko hypothesis [27]. These research results have important scientific value for understanding the regularities of runoff and vegetation change in the SAYR, but there are few studies on the assessment and quantitative analysis of climate change and the contribution rate of anthropic factors to runoff and vegetation change in the SAYR [27,28].

The objective of this study was to quantify the impact of anthropic factors and climate change on vegetation and runoff variation in the SAYR by following three steps: (1) analyzing the temporal variation regularities of runoff and vegetation in the SAYR; (2) distinguishing the mutation year of runoff and NDVI data with the Mann–Kendall mutation analysis method; and (3) quantifying the contribution rate of climatic and anthropic factors to vegetation and runoff in the SAYR.

2. Study Area and Data

2.1. Study Area

The source region of the Yangtze River is located between 90°43′ and 96°45′ E, 32°30′ and 35°35′ N, with an area of about 13.77 × 10⁴ km² (Figure 1). Its annual average runoff is about 129.17 × 10⁸ m³, and it affects the interannual fluctuation, long-term evolution trend, and the sustainable utilization of water resources in the Yangtze River Basin. The landform is mainly high plains and hills, with an average altitude of more than 4000 m. Its climate type belongs to the semi-humid and semi-arid region of the plateau frigid zone, with abundant sunshine and large temperature difference between day and night. Temperatures are generally low throughout the year in the Yangtze River source area. The month with the lowest precipitation is December, the month with the highest precipitation is July, and the annual precipitation is concentrated in May–September, which accounts for more than 85% of the annual precipitation [29].

2.2. Data Sources

The data used in this study consist of three parts: (1) The annual runoff observation data of Zhimenda station from 1982 to 2016 were all from the Yangtze River Water Conservancy Commission (http://www.cjw.gov.cn/, accessed on 1 January 2021); (2) the climate station data in and around the source region of the Yangtze River Basin from 1982 to 2016 were obtained from the China Meteorological Administration (http://www.cma.gov.cn, accessed on 1 January 2021). First, the daily potential evaporation of 25 meteorological stations was calculated using the Penman–Monteith equation, and then the monthly precipitation and potential evaporation of 25 meteorological stations were calculated. Finally, the monthly precipitation and potential evaporation were interpolated from data of 25 meteorological stations in and around the SRYR by kriging. Annual precipitation and potential evaporation were obtained by adding monthly scale data. (3) The NDVI data of the source region of the Yangtze River Basin from 1982 to 2016 were obtained from the NOAA CDR AVHRR NDVI: Normalized Difference Vegetation Index, Version 5 (https://www.ncdc.noaa.gov/cdr/terrestrial/normalized-difference-vegetation-index, accessed on 1 January 2021). The time resolution is daily, and the spatial resolution is 0.05° × 0.05°. The monthly NDVI and annual NDVI were extracted with the maximum combination method.

3. Research Methods

3.1. Trend Analysis Method

In this study, the univariate linear regression method (Slope) was used to analyze the variation trend of runoff, precipitation, potential evaporation, and NDVI in the study period; time t was taken as the independent variable. The calculation equation is as follows [30,31,32]:

$$Slope=\frac{n{{\displaystyle \sum }}_{i=1}^{n}i{X}_i-{{\displaystyle \sum }}_{i=1}^{n}i{{\displaystyle \sum }}_{i=1}^n{X}_i}{n{{\displaystyle \sum }}_{i=1}^n{i}^2-{\left({{\displaystyle \sum }}_{i=1}^{n}i\right)}^2}$$

(1)

where Slope is the simple linear regression coefficient of runoff, precipitation, potential evaporation, and NDVI from 1982 to 2016. If Slope is positive, it means that the variable has an upward trend in the study period, whereas a negative value represents a downward trend; n is the number of years of the study period, i is 1 to n, n is 35, and X represents the annual runoff, precipitation, potential evaporation, and NDVI in the study area.

3.2. Mann–Kendall Mutation Analysis Method

Mann–Kendall (MK) mutation analysis has been widely used to test the abrupt point of hydro-meteorological elements due to its advantages of minimized interference from outliers and simple calculation [33,34].

For X with n data, a rank sequence is constructed, as shown in Equations (2) and (3):

$$Sk={{\displaystyle \sum }}_{i=1}^{k}ri (k=2, 3, \dots , i)$$

(2)

$$ri=\left\{\begin{array}{cc}+1& xi>xj\\ 0& xi\le xj\end{array}\right. (j=2, 3, \dots , i)$$

(3)

The UF_k statistics can be calculated by Equation (4):

$$U{F}_k=\frac{\left[sk-E\left(sk\right)\right]}{\sqrt{Var\left(sk\right)}}(k=1, 2, \dots , n)$$

(4)

where UF₁ = 0, E(sk) and Var(sk) represent the mean and variance of sk, respectively. The calculation equations are as follows:

$$E\left(sk\right)=\begin{array}{cc}\frac{n\left(n+1\right)}{4}& \left(2\le k\le n\right)\end{array}$$

(5)

$$Var\left(sk\right)=\begin{array}{cc}\frac{n\left(n-1\right)\left(2n+5\right)}{72}& \left(2\le k\le n\right)\end{array}$$

(6)

The inverse order of X is calculated again, and at the same time, UB_k = UF_k (k = n, n − 1, …,1), UB₁ = 0, then the UB_k curve can be calculated. When the two lines (UF curve and UB curve) intersect, and the intersection point is within the range of the 0.05 significance level, the intersection point can be considered as the mutation year.

3.3. Budyko Hypothesis

There are three assumptions in the attribution analysis of runoff change using the Budyko hypothesis: (1) Human activities and climate change do not affect each other and are independent factors. (2) For multi-year water balance, the change of water storage is usually negligible compared with the average annual precipitation depth. Therefore, the change of catchment storage water for multi-year water balance is assumed as 0; (3) The base period is only affected by climate change. Therefore, except for climate change, other factors affecting runoff change are classified as human activities.

The water balance equation of the basin is expressed as follows:

$$R=P-ET-\Delta S$$

(7)

R is the runoff depth; P is precipitation; ET is the actual evapotranspiration; and ΔS is the change of water storage. For multi-year water balance, the change of water storage is usually negligible compared with the average annual precipitation depth. Therefore, in the study of long-term hydrological data, it is generally considered that ΔS is 0.

The Budyko hypothesis, based on an assumption that the change of catchment storage water for multi-year water balance is considered as 0, is widely applied to quantitatively calculate the impact of climate change and anthropic factors on long-term annual runoff change [35,36,37,38,39]. ET can be calculated according to the Budyko hypothesis, the calculation equation of which is as follows [40,41]:

$$ET=\frac{ P\times ET0}{{\left(P^{\omega }+ET{0}^{\omega }\right)}^{1/\omega }}$$

(8)

ET₀ is the annual average potential evapotranspiration (mm); ω indicates characteristic parameters of underlying surface.

ET₀ can be calculated using the Penman–Monteith equation.

$$E{T}_0=\frac{0.408\Delta \left(R_n-G\right)+\gamma \frac{900}{T+273}{U}_2\left(e_a-e_d\right)}{\Delta +\gamma \left(1+0.34U_2\right)}$$

(9)

Equation (7) can be transformed into Equation (10):

$$R=P-\frac{ P\times ET0}{{\left(P^{\omega }+ET{0}^{\omega }\right)}^{1/\omega }}$$

(10)

where R, P, and ET₀ are known, and parameter ω can be calculated by the “fsolve” function of MATLAB.

Based on Equation (10), the elastic coefficients of precipitation (${\epsilon }_P$), potential evaporation (${\epsilon }_{ET0}$), and underlying surface parameters (${\epsilon }_{\omega }$) to runoff can be calculated by the following equations [42]. The concept of elastic coefficient is shown in Appendix A.

$${\epsilon }_P=\frac{{\left(1 + {\varphi }^{\omega }\right)}^{1/\omega +1}-{\varphi }^{\omega +1}}{\left(1 + {\varphi }^{\omega }\right)\left[{\left(1 + {\varphi }^{-\omega }\right)}^{1/\omega }-\varphi \right]}$$

(11)

$${\epsilon }_{ET0}=\frac{1}{\left(1 + {\varphi }^{\omega }\right)\left[1-{\left(1 + {\varphi }^{-\omega }\right)}^{1/\omega }\right]}$$

(12)

$${\epsilon }_{\omega }=\frac{\ln \left(1 + {\varphi }^{\omega }\right) + {\varphi }^{\omega }\ln \left(1 + {\varphi }^{\omega }\right)}{\omega \left(1 + {\varphi }^{\omega }\right)\left[1-{\left(1 + {\varphi }^{-\omega }\right)}^{1/\omega }\right]}$$

(13)

$$\varphi =\frac{ET0}{P}$$

(14)

The runoff variation amount caused by the annual average precipitation change ($\Delta RP$), annual average potential evaporation change ($\Delta RET0$), and underlying surface parameter change ($\Delta R\omega$) can be calculated by the following Equations (15)–(17):

$$\Delta R\mathrm{P}=\epsilon P\frac{R}{P}\times \Delta P$$

(15)

$$\Delta R\mathrm{ET}0=\epsilon ET0\frac{R}{ET0}\times \Delta ET0$$

(16)

$$\Delta R\omega =\epsilon \omega \frac{R}{\omega }\times \Delta \omega$$

(17)

On this basis, the contribution rate of precipitation ($\eta R_p$), potential evaporation, ($\eta R_{ET0}$) and anthropic factor ($\eta R_H$) to runoff changes can be calculated using Equations (19)–(21):

$$\Delta R=\Delta RP+\Delta RET0+\Delta R\omega$$

(18)

$$\eta R_P=\Delta RP/\Delta R\times 100\%$$

(19)

$$\eta R_{ET0}=\Delta RET0/\Delta R\times 100\%$$

(20)

$$\eta R_H=\Delta R\omega /\Delta R\times 100\%$$

(21)

3.4. Attribution Analysis of Climate and Anthropic Factors to Vegetation Change

There are two assumptions in the attribution analysis of vegetation change using this method: (1) Human activities and climate change do not affect each other and are independent factors. (2) The base period is only affected by climate change. Therefore, except for climate change, other factors affecting vegetation change are classified as human activities.

According to the mutation analysis results of NDVI time series data, the NDVI data are divided into two parts: base period (T₁) and change period (T₂). On this basis, the change of average NDVI in the two periods ($△\overline{NDVI}$) can be calculated as follows:

$$△\overline{NDVI}=\overline{NDV{I}_{T2}}-\overline{NDV{I}_{T1}}$$

(22)

where $\overline{NDV{I}_{T1}}$ and $\overline{NDV{I}_{T2}}$ are the average NDVI values in the base period (T₁) and change period (T₂), respectively. This method assumes that NDVI in the base period is only affected by climate change. Therefore, the NDVI difference between the base period and change period can be attributed to climatic factors ($△\overline{NDV{I}_C}$) and anthropogenic factors ($△\overline{NDV{I}_H}$).

The monthly NDVI values are quantitatively correlated with monthly precipitation and potential evapotranspiration [43,44]. Therefore, a multiple linear regression equation between monthly NDVI, monthly precipitation, and monthly potential evapotranspiration in the base period (T₁) can be established:

$$\overline{NDV{I}_{T1}}=a\ast PT1+b\ast ET0T1+c$$

(23)

The monthly precipitation (P_T2) and potential evapotranspiration (ET0_T2) in the T₂ period are known, and the simulated NDVI in the T₂ period ($\overline{NDV{I}_{T2,S}}$) can be calculated by Equation (24).

$$\overline{NDV{I}_{T2,S}}=a\ast PT2+b\ast ET0T2+c$$

(24)

This method assumes that NDVI in the base period ($\overline{NDV{I}_{T1}}$) is only affected by climate change. Therefore, $\overline{NDV{I}_{T2,S}}$ in the T₂ period calculated by Equation (24) is only affected by climate change. $\overline{NDV{I}_{T2}}$ is the average value of NDVI observation in the T₂ period, which is affected by the combination of climate change and human activities.

The contribution rates of anthropic factor ($\eta NDV{I}_H$) and climate change ($\eta NDV{I}_C$) to NDVI can be calculated separately by Equations (25)–(28):

$$△\overline{NDV{I}_H}=\overline{NDV{I}_{T2}}-\overline{NDV{I}_{T2,S}}$$

(25)

$$△\overline{NDV{I}_C}=\overline{NDV{I}_{T2,S}}-\overline{NDV{I}_{T1}}$$

(26)

$$\eta NDV{I}_H=△\overline{NDV{I}_H}/△\overline{NDVI}$$

(27)

$$\eta NDV{I}_C=△\overline{NDV{I}_C}/△\overline{NDVI}$$

(28)

4. Results and Analysis

4.1. Trends Analysis of Runoff, NDVI, and Climate Factors

Figure 2 displays the change trend of the annual runoff in Zhimenda hydrological station and NDVI in the SAYR. It can be found that the annual runoff of Zhimenda hydrological station shows an increasing trend from 1982 to 2018, with an average annual increase of 0.8831 × 10⁸ m³. The NDVI value in the SAYR showed a significant growth trend from 1982 to 2016, with an average annual growth of 0.0018.

Figure 3 displays the change trend of average annual precipitation and annual potential evaporation in the SAYR. It can be found from Figure 3a that the annual average precipitation in the SAYR showed an increasing trend from 1982 to 2016, with an average annual increase of 1.2394 mm. Figure 3b displays that the average annual potential evaporation in the SAYR also showed an upward trend from 1982 to 2016, with an average annual increase of 0.5965 mm.

4.2. Mutation Analysis of Runoff and NDVI

In this study, the Mann–Kendall method was used to analyze runoff and NDVI data in the SAYR, and the results are shown in Figure 4 and Figure 5.

It can be found from Figure 4 that the UF and UB curve calculated from the runoff time series data in the SAYR intersected in 2004, and this intersection point was within the range of the 0.05 significance level, indicating that 2004 was the abrupt change year of runoff in the SAYR from 1982 to 2016.

It can be found from Figure 5 that UF and UB curve calculated from the NDVI time series data in the SAYR intersected in 1998, and this intersection point was within the range of the 0.05 significance level, indicating that 1998 was abrupt change year of NDVI in the SAYR from 1982 to 2016. This may be because Grain for Green was proposed by Chinese government in 1998, and the impact of human activities on vegetation change has increased since 1999.

4.3. Assessment of the Contribution Rates of Climate and Anthropic Factors to Runoff and Vegetation Changes

Identifying the abrupt change year of runoff and NDVI was the premise for dividing time series data into different periods (i.e., base period and changing period). If there was no abrupt change year of runoff and NDVI data, we did not divide time series data into different periods.

According to mutation analysis results of runoff time series data at Zhimenda hydrological station from 1982 to 2016, the Zhimenda hydrological station 1982–2016 data were divided into two periods: the base period (1982–2004) and the change period (2005–2016), and the eigenvalues of meteorological and hydrological variables for different periods in the Zhimenda hydrological station were calculated (

Table 1. Eigenvalues of meteorological and hydrological variables in the SAYR.

Hydrological Station	Period	ET0/mm	R/mm	P/mm	$\mathit{\omega }$	R/P	ET0/P
Zhimenda	1982–2004	831.74	90.57	386.35	1.25	0.23	2.15
	2005–2016	852.93	115.06	425.29	1.18	0.27	2.01

).

In order to further evaluate the temporal evolution characteristics of the impact of climate factors and underlying surface parameters on runoff, according to Equations (11)–(13), the elastic coefficients of precipitation, potential evaporation, and underlying surface characteristic parameters in different years on runoff change of Zhimenda hydrological station were calculated (Figure 6).

It can be seen from Figure 6 that the absolute values of $\epsilon \mathrm{P}$, $\epsilon \mathrm{ET}0$, and $\epsilon \omega$ in for Zhimenda hydrological station all showed a decreasing trend, indicating that the sensitivity of runoff change in the SAYR to climate factors and underlying surface parameters is decreasing.

According to Table 1, the differences of average annual precipitation, potential evaporation, and underlying surface parameter in the base period (1982–2004) and change period (2005–2016) were calculated and are displayed in

Table 2. Attribution analysis of runoff variation in the SAYR.

Hydrological Station	εP	εET0	εω	$\Delta \mathit{R}$	ΔP	ΔET0	Δω	ΔRP	ΔRET0	ΔRL	ηRP	ηRET0	ηRH
Zhimenda	1.90	−0.90	−1.54	24.49	38.94	21.19	−0.07	18.34	−2.26	24.14	75.98%	−9.35%	33.37%

. Then, according to Equations (15)–(17), the runoff variation amount for Zhimenda hydrological station caused by climate factors (precipitation, potential evaporation) and underlying surface condition change were calculated, as displayed in Table 2. Finally, the contribution rate of climate (precipitation, evaporation) and anthropic factors to runoff change of Zhimenda hydrological station was calculated according to Equations (18)–(21), as shown in Table 2.

The contribution rates of climate and anthropic factors to runoff variation in Zhimenda hydrological station are 66.63% and 33.37%, respectively. In general, climate change plays a major role in increasing runoff, and precipitation has a more significant effect on runoff increase than the reference evapotranspiration.

This conclusion is similar to that of Liu et al. [27]. However, there are some differences in the specific contribution rate value, and this may be due to the following factors: (1) The time ranges of runoff data are different; (2) different methods are used to calculate the contribution rate of climatic (precipitation, potential evaporation) and anthropic factors to runoff variation.

According to the abrupt analysis results of NDVI data of the SAYR from 1982 to 2016, the NDVI data in the SAYR from 1982 to 2016 can be divided into two periods: the base period (1982–1998) and the change period (1999–2016). Previous studies have found that monthly NDVI is closely related to monthly precipitation and potential evaporation [43,44], and NDVI change can be used to characterize vegetation change. Referring to the method proposed by Li et al. [45], this study first established the multiple linear regression equation between monthly NDVI, monthly precipitation, and monthly potential evapotranspiration in the period of 1982–1998, and then the contribution rate of climate change and anthropic factors to vegetation change was quantitatively analyzed according to Equations (25)–(28) (

Table 3. Attribution analysis of vegetation change in the SAYR.

Fitting Equation	$△\overline{\mathit{N}\mathit{D}\mathit{V}\mathit{I}}$	$\overline{\mathit{N}\mathit{D}\mathit{V}{\mathit{I}}_{\mathit{T}1}}$	$\overline{\mathit{N}\mathit{D}\mathit{V}{\mathit{I}}_{\mathit{T}2,\mathit{S}}}$	$\overline{\mathit{N}\mathit{D}\mathit{V}{\mathit{I}}_{\mathit{T}2}}$	$\mathit{\eta }\mathit{N}\mathit{D}\mathit{V}{\mathit{I}}_{\mathit{C}} (\%)$	$\mathit{\eta }\mathit{N}\mathit{D}\mathit{V}{\mathit{I}}_{\mathit{H}} (\%)$
NDVI = 1.7912 × 10−3P + 7.1204 × 10−4ET0 + 0.1057 (R2 = 0.68)	0.0197	0.2122	0.2198	0.2319	38.56	61.44

). The results showed that compared to the period of 1982–1998, anthropic factors played a major role in increasing vegetation coverage in the period of 1999–2016, with a contribution rate of 61.44%, and the contribution rate of climate factors was 38.56%.

5. Discussion

Although the change trend and mutation characteristics of runoff and vegetation in the SAYR were analyzed, and the impact of climate change and anthropic factors on runoff and vegetation change were calculated, there are still some uncertainties. The data of precipitation and potential evaporation in the study area are interpolated from low-density meteorological station data in and around the SAYR, but there are still some deviations between the interpolation results and their actual distribution. The change of catchment storage water for multi-year water balance is assumed as 0, and we assumed that human activities and climate change do not affect each other and are independent factors. This study also assumed that the base period is only affected by climate change. All of these factors will lead to some uncertainties in the research results.

6. Conclusions

This study revealed the change trend and mutation characteristics of runoff and vegetation in the SAYR, and the contribution rates of climate change and anthropic factors to runoff change in the SAYR were quantitatively calculated using the Budyko hypothesis and multiple linear regression method, which provides theoretical support for water resource management and ecological protection in the SAYR.

The results showed that the runoff, NDVI, precipitation, and potential evaporation in the SAYR from 1982 to 2016 all showed an increasing trend. The abrupt change years of runoff and NDVI data in the SAYR from 1982 to 2016 were 2004 and 1998, respectively. Anthropic factors play a major role in annual runoff and vegetation change, with contribution rates of 75.98% and 61.44%. In the follow-up study, we will try to quantitatively analyze the contribution rate of climate and human factors to vegetation and runoff changes in different seasons.