Copula-Based Abrupt Variations Detection in the Relationship of Seasonal Vegetation-Climate in the Jing River Basin, China

Understanding the changing relationships between vegetation coverage and precipitation/temperature (P/T) and then exploring their potential drivers are highly necessary for ecosystem management under the backdrop of a changing environment. The Jing River Basin (JRB), a typical eco-environmentally vulnerable region of the Loess Plateau, was chosen to identify abrupt variations of the relationships between seasonal Normalized Difference Vegetation Index (NDVI) and P/T through a copula-based method. By considering the climatic/large-scale atmospheric circulation patterns and human activities, the potential causes of the non-stationarity of the relationship between NDVI and P/T were revealed. Results indicated that (1) the copula-based framework introduced in this study is more reasonable and reliable than the traditional double-mass curves method in detecting change points of vegetation and climate relationships; (2) generally, no significant change points were identified during 1982–2010 at the 95% confidence level, implying the overall stationary relationship still exists, while the relationships between spring NDVI and P/T, autumn NDVI and P have slightly changed; (3) teleconnection factors (including Arctic Oscillation (AO), Pacific Decadal Oscillation (PDO), Niño 3.4, and sunspots) have a more significant influence on the relationship between seasonal NDVI and P/T than local climatic factors (including potential evapotranspiration and soil moisture); (4) negative human activities (expansion of farmland and urban areas) and positive human activities (“Grain For Green” program) were also potential factors affecting the relationship between NDVI and P/T. This study provides a new and reliable insight into detecting the non-stationarity of the relationship between NDVI and P/T, which will be beneficial for further revealing the connection between the atmosphere and ecosystems.


Introduction
Terrestrial vegetation is the most important component of an ecosystem due to its critical roles in energy budget, water, and biogeochemical cycle of the earth ecosystem through photosynthesis and transpiration, and it can change the climate system and alter the land surface [1][2][3][4].From another perspective, the variation of vegetation coverage is a reflection of corresponding climatic conditions such as precipitation and temperature [5,6].It is well known that climate changes have been characterized by an increase of 0.85 • C global surface temperature during 1880-2010, which was documented by the Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report [7].Accordingly, this warming trend provides abundant thermal energy to regulate the inner biogeochemical processes of plants, thus affecting vegetation coverage, vegetation health, and vegetation productivity [8][9][10][11].
There is increasing evidence that vegetation growth has a relationship with many physiographic and climatic attributes such as soil moisture [12,13], leaf area index [14], and regional evapotranspiration [15,16].Among the climatic factors, precipitation (P) and temperature (T) have been proven to be the main factors influencing the distribution and composition of vegetation.Several studies have investigated the relationship between the Normalized Difference Vegetation Index (NDVI) and P/T [17][18][19].For instance, Wen et al. [20] reported that vegetation was closely related to asymmetric warming, and it responded to the asymmetric warming with nearly one year delays at a global scale.Ji and Peters [21] investigated the relationship between vegetation and rainfall in the northern and central Great Plains of the USA and found the time lag and accumulative effects of rainfall on vegetation.All the above-mentioned studies have revealed the relationships between NDVI and P/T, and have confirmed the time lag effects between NDVI and P/T.However, under the background of a changing environment, the stationarity of the relationship between NDVI and P/T might be altered, and this could strongly affect the accuracy of the prediction models for vegetation based on its correlation with climatic factors.Nonetheless, fewer comparative studies have identified the change points in the relationship between NDVI and P/T.Hence, exploring the non-stationarity of the relationship between NDVI and climatic factors is the main objective in the present study, which helps to further understand the mechanisms of vegetation dynamic.
Nevertheless, some statistical methods have been adopted to investigate the change points of the binary relationship.Among multiple statistical methods, the correlation coefficient method was commonly employed in previous studies [22][23][24].However, the correlation coefficient method which directly reflects the general trend and the strength of the bivariate relationship assumes that the relationships of the two time series are monotonous and stationary [25].Due to the complexity of ecosystems, the relationship between NDVI and climatic factors is affected by multiple factors and characterized by nonlinearity which will cause the correlation coefficient method to assume that stationarity is not applicable here, and fail to capture the change points in the relationship between variables.Besides, conceptual models were utilized in some previous studies to identify change points.Lee et al. [26] proposed 12 conceptual model structures for the regionalization of the precipitation-runoff relationship on 28 UK catchments.In the Birkinshaw and Bathurst [27] study, a spatially distributed model was used to investigate the variation of relationship between sediment yield and river basin area.However, those model-based methods require a large amount of input data, and it is difficult to determine their parameters.Copula functions, which can represent the dependency structure between variables by constructing the joint distributions of related variables and reflect the linear and nonlinear relationship between variables, overcome the shortcomings of previous methods.As copula functions can directly deal with non-stationarity, and are characterized by flexibility and less limitation on the types of marginal distributions to be connected, they have been widely used in the financial field and are now increasingly adopted in hydro-meteorology fields to detect the abrupt change of the relationships between two variables.Gu et al. [28] developed a new regional contagion detection method based on the vine copula to conduct a study on financial contagion.Huang et al. [29] detected the possibly changing relationship between precipitation and temperature in the Wei River Basin, China, by copula functions.However, the copula function has not been used to identify the change point of the relationship between vegetation and meteorological variables.Therefore, this study seeks to utilize copula functions to present the relationship between NDVI and P/T and simultaneously detect the change points of their relationships [30].
Furthermore, various studies have been conducted to investigate the possible connection between atmospheric teleconnection factors and vegetation growth and climatic factors [31][32][33].For instance, Li et al. [34] documented that the spatial patterns of NDVI have closer connections to El Niño over Eurasia.Cavazos [35] presented that changes in the circulation related to Arctic Oscillation (AO) strongly controlled the winter extreme precipitation in Balkans.Additionally, Wu and Wang [36] evaluated the impact of AO and Siberian High on the East Asian winter monsoon, and they concluded that higher winter air temperature occurs over East Asia during the positive AO phase.Nevertheless, considerable studies previously have focused on exploring the impacts of climatic and teleconnection factors on the individual components such as NDVI, P, and T, but the underlying causes for the variations in the relationship between NDVI and P/T by the climatic and atmospheric teleconnection factors (e.g., soil moisture (SM), potential evapotranspiration (PET), AO, Pacific Decadal Oscillation (PDO), El Niño-Southern Oscillation (ENSO), and sunspots) across seasons still lacks detection.Moreover, the dominant drivers leading to variations in the relationship between vegetation and climatic factors have not been revealed.Thereby, a greater understanding of how each teleconnection factor influences the relationship between NDVI and P/T will be fully explored in this present study.
The Jing River Basin (JRB), which belongs to the Loess Plateau, China, is a typical eco-environmentally vulnerable region in the world.The Loess Plateau is characterized by a highly erosive loess layer, which leads to serious degradation of the ecological environment in the Loess Plateau after centuries of unsustainable agricultural production and population explosion.Therefore, the Chinese government has carried out a series of ecological restoration projects since the 1970s to solve the serious environmental and ecological problems.Among them, the "Grain for Green" project implemented since 1999 has received remarkable feedback.The vegetation condition has been greatly improved, which might influence the relationship between NDVI and P/T.Consequently, it is necessary to investigate the relationship between NDVI and P/T in the JRB.Based on the long-term NDVI data and hydrometeorological observations, we aim (1) to introduce a copula-based framework for the change points detection and verify the reliability and superiority of this method; (2) to identify the change points of the relationship between seasonal NDVI and P/T; (3) to examine the impact of climate change on the variations of the relationships between NDVI and P/T and to determine the dominant driver of the seasonal variations.Generally, this study provides a new and reliable insight into detecting the non-stationarity of the relationship between NDVI and P/T, which will be beneficial for further revealing the connection between the atmosphere and ecosystems.

Study Area
The JRB (106.2 • E-109.1 • E, 34.8 • N-37.4 • N) is a typical arid and semi-arid region which is located in the southwest of the Loess Plateau (Figure 1).As the second level tributary of the Yellow River Basin (YRB), the JRB covers a total area of approximately 4.54 × 10 4 km 2 .Located in the transitional zone between the temperate semi-humid and temperate semi-arid regions, the JRB experiences a typical temperate continental monsoon climate and the mean annual precipitation is 545 mm.The mean temperature in the coldest month ranges from −3 to −1 • C, whilst that in the hottest month varies from 23 to 26 • C [37].Overall, the JRB is characterized by abundant precipitation and high temperature in summer and by rare precipitation and low air temperature in winter.Loessial soil and dark loessial soil are two typical types of soil which are highly erodible and widely distributed in the study areas [38,39].Nearly 60% of annual precipitation is concentrated in the flood season (from June to September), which makes water loss and soil erosion occur frequently in the JRB.Thus, the JRB is widely known as a sediment-laden basin with 2.6 × 10 8 t sediment transporting into the YRB on average [40].Accordingly, the ecological environment is extremely fragile and the vegetation coverage is sparse in the JRB.

Datasets
The NDVI is a vegetation indicator which reflects the vegetation coverage by separating vegetation from water and soil.In this study, we focused on the vegetation dynamics during 1982-2010.The NDVI remote sensing data was obtained from the US National Oceanic and Atmospheric Administration's (NOAA) Advanced Very High Resolution Radiometer (AVHRR) (https://nex.nasa.gov/nex/projects/1349/).Due to the poor vegetation coverage in winter on the Loess Plateau [41], we only focused on spring, summer, and autumn vegetation in the present study.Moreover, in order to avoid data changes caused by cloud cover and other uncertain factors, we used the maximum value of each season to represent the seasonal vegetation condition.In the current study, the gridded NDVI data within the JRB was averaged for use.
Additionally, the point daily meteorological data covering 1982-2010, including daily maximum and minimum temperature, wind speed, air pressure, sunshine duration, relative humidity, and daily mean precipitation, were recorded by nine meteorological stations and the dataset were obtained from the National Climate Center (NCC) of the China Meteorological Administration (CMA).Recommended by the World's Food and Agriculture Organization (FAO) in 1998, the Penman-Menteith method was used to calculate the point PET.Based on the point meteorological data, the regional P and PET in the JRB were calculated by the Thiessen polygon method which is a widely used technique in deriving regional P and PET [42].
Besides, in order to examine the correlations between climatic/large-scale atmospheric circulations and relationships between seasonal NDVI and P/T, the PET, SM, AO, PDO, Niño 3.4, and sunspots were also employed in this study.The gridded monthly SM data was estimated by the Variable Infiltration Capacity (VIC) model.In the current study, the gridded SM data within the JRB area was averaged for use.AO data were collected from the NOAA National Climatic Data Center (http: //www.cpc.ncep.noaa.gov/products/precip/CWlink/daily_ao_index/ao_index.html), and the monthly PDO index was downloaded from the Tokyo Climate Center (http://ds.data.jma.go.jp/tcc/tcc/products/ elnino/decadal/pdo.html).The monthly Niño 3.4 index was used to represent ENSO activities which were acquired from the NOAA Earth System Research Laboratory (http://www.esrl.noaa.gov/psd/data/correlation/nina34.data).Monthly sunspot indexes were obtained from the National Geophysical Data Center (NGDC) of the NOAA (https://www.esrl.noaa.gov/psd/gcos_wgsp/Timeseries/SUNSPOT/).Furthermore, in order to explore the correlations between human activities and the relationship between NDVI and P/T, the effective irrigated area data of the Shaanxi Province, Gansu Province, and Ningxia District were collected from the Ministry of Agriculture and Rural Affairs of the People's Republic of China (http://zzys.agri.gov.cn/nongqing.aspx).Based on the province data and the area weight of the JRB in each province, the effective irrigated area data in the JRB was obtained.Moreover, land use maps (1:100,000) for six periods (1980, 1990, 1995, 2000, 2005, and 2010) were used in the current study.They were generated from Landsat Thematic Mapper (TM) images using visual interpretation and obtained from the Resource and Environment Data Cloud Platform of Chinese academy of sciences (http://www.resdc.cn/data.aspx?dataid=176) which provides accurate, abundant, reliable, dynamic geographic information sources by collecting dynamic observations of the earth from the 'Spot' satellites.According to the national standard of classification of land use status (GB/t21010-2007), the land use types were classified into farmland, forestland, grassland, water bodies, construction land, and unused land.

Trend Analysis
The original Mann-Kendall (MK) test is a popular non-parametric method for accessing the trends of hydrometeorological variables [43].However, the MK method is based on the uncorrelated data and the results are easy to be affected by the persistence of time series.Therefore, Hamed and Rao [44] proposed a modified Mann-Kendall (MMK) method by containing the lag-i autocorrelation to remove the persistence.In this study, we applied the MMK test to explore the variations trends of vegetation cover and precipitation/temperature, and the significance level was set at 0.05.The detailed computational processes can be referred to our previous research [45,46].

A Copula-Based Framework for Identifying the Change Points of the Relationship between NDVI and P/T
To identify the change points of the relationship between seasonal NDVI and P/T, a copula-based framework is implemented in the current study.This framework contains of such steps as fitting the most appropriate marginal distributions, selecting the joint distributions, and identifying the change points.

Marginal Distribution
A marginal distribution is the projection of the joint distribution of a set of random variables onto a subspace which is defined by a subset of components.Since the seasonal NDVI, P, and T series utilized in the current study are continuous, some parametric distributions (such as Gamma distribution, Generalized Extreme Value (GEV) distribution, and lognormal distribution), which are frequently applied to fit the distribution of hydrometeorological time series, were applicable in this study [47].Besides, the parameters of these distributions were evaluated by the Maximum Likelihood Estimation (MLE) method [48].Since the errors in marginal distributions can be amplified in the determination of joint distribution, efforts should be made to fit marginal distributions as accurately as possible.Based on the research of Wilks [49], the goodness-of-fit of each individual distribution was evaluated by the Kolmogorov-Smirnov (K-S) test, thus the most suitable distribution for each individual NDVI and P/T was obtained.

Joint Distribution
The copula function can construct the joint probability distribution of the multivariate hydrologic series through integrating their corresponding marginal distributions.The formula of the copula is expressed as follows: where ϕ denotes the convex function; and u and v represent the two variables.Generally, the Archimedean copulas are widely used for the frequency analysis of multivariate hydrologic series due to their easily constructed characteristics and the numerous available copula families [50].In order to avoid over-parameterization and easily estimate the parameter, Gumbel, Frank, and Clayton copulas which were widely known as three simple Archimedean copulas, were chosen to assess the joint probability distribution of the NDVI and climate factors in the JRB.Acquiring the generating function from multivariate observation data is the crucial step to determine a copula which was shown in detail in the study of Genest and Rivest [51].In the present study, the minimum criterion of the Akaike Information Criterion (AIC) and the Root-Mean-Square Error (RMSE) were adopted to select the most suitable copula from the three Archimedean copulas [52].The Archimedean copula with the lower AIC and RMSE values indicated better fit to the joint distributions of seasonal NDVI and P/T and superior performance.

Identify Change Points
The Copula-based Likelihood-ratio method (CLR) was utilized to identify the change points of the relationship between seasonal NDVI and P/T.For a time series of (x 1, y 1 ), • • • , (x n, y n ), if only one change point of the relationship between seasonal NDVI and P/T series exists, the null hypothesis (H 0 ) and alternative hypothesis (H 1 ) are expressed as follows: If the null hypothesis is rejected, k * is the corresponding time of the change point.When k * = k is known, the logarithmic likelihood ratio statistics of copula based on the maximum likelihood estimation method can be formed as: where F(x) and G(y) are the probability distribution function of seasonal NDVI and P/T series, respectively.C 12 represents the joint distribution function of seasonal NDVI and P/T series.λ k , λ k * , and λ n denote the parameter λ of the joint distribution of the seasonal NDVI and P/T series which is estimated by the maximum likelihood method.Considering the fitting effects of the relevant marginal and joint distributions, the range of k is confined as 6 and n−5.
If k * is unknown, then When the statistic Z n , which follows the chi square distribution, is large enough to reject the null hypothesis, a change point of the relationship between seasonal NDVI and P/T on the basis of Archimedean Copula was determined.Introduced in da Costa Dias [53], the threshold of the statistic Z n is approximately 8.8.The time of the change point can be assessed as follows:

Correlation Analysis
Pearson correlation coefficients were computed to access the relationship between climate change and the relationships between seasonal NDVI and P/T (denoted by Z n series).The correlation coefficients can be written as: where R xy is the correlation coefficient, x and y represent the seasonal values of two variables, and x and y are the average of the two variables.n is the length of the time period.The overall methodology adopted in the current study is presented in Figure 2. First, the temporal trend is assessed by the MMK method.Then, the change point was identified by the copula-based framework and the double-mass curve method.Additionally, the Pearson correlation coefficient analysis was used for the attribution analysis.

Temporal Change of Seasonal NDVI, Precipitation, and Temperature
The seasonal NDVI series covering 1982-2010 in the JRB is plotted in Figure 3.It can be roughly observed that NDVI in three seasons exhibited a noticeably increasing trend with fluctuation, and the maximum increase rate was found in autumn (0.0023/year).The MMK statistics of NDVI in spring, summer, and autumn are 0.8, 1.4, and 4.1, respectively (Table 1), which means the autumn NDVI had a statistically significant upward trend at the 99% confident level, whereas spring NDVI and summer NDVI had a non-significant increasing trend.It has been documented that vegetation growth exhibited time-lag responses to climate variations [54,55].Thereby, the correlation coefficients between the NDVI and P/T were computed by considering the time-lag effects.The results revealed that spring NDVI had the greatest correlation with the previous 1-2 months accumulative precipitation and average temperature, and the coefficients were 0.49 and 0.38, respectively.Summer NDVI exhibited a one month lag with precipitation and no lag with temperature.In autumn, the time-lag of NDVI to precipitation and temperature was previous 0-1 month and 0-3 months, respectively (0 represented the current month).Consequently, the long-term precipitation and temperature series mentioned in this study were the corresponding series with time-lag effects taken into account.As shown in Table 1, precipitation in the three seasons showed a non-significant trend, while temperature had a remarkable upward trend at the 95% confident level in autumn and at the 99% confident level in spring and summer.Generally, seasonal vegetation coverage and temperature in the JRB have obviously increasing trends, while precipitation has a non-significant trend.Note: "*" and "**" represent significant at 95% and 99% confidence level, respectively.

The Selection of the Appropriate Marginal Distribution
In this study, Gamma distribution, GEV distribution, and lognormal distribution were chosen to fit the seasonal NDVI and P/T series, and then the most appropriate marginal distributions were obtained.The goodness-of-fit values of the three probability distributions were calculated by the K-S method and are shown in Table 2. Apparently, the H values of all the selected distributions are 0, indicating that all three distributions passed the K-S test in fitting seasonal NDVI and P/T series.In general, the results demonstrated that the p-values of the GEV distribution were the largest among the three distributions in fitting the seasonal NDVI and P/T, except in fitting spring temperature and autumn NDVI.Therefore, the lognormal distribution which was shown to have the maximum p-value in fitting spring temperature and autumn NDVI was selected as the most appropriate marginal distribution for spring temperature and autumn NDVI, while the GEV distribution was chosen as the most appropriate marginal distribution for the other time series.

The Selection of the Appropriate Copula Function
When the most appropriate marginal distribution was determined, the RMSE and AIC method were utilized to select the most appropriate copula from the Archimedean copulas (including Clayton, Frank, and Gumbel copulas) which were chosen to fit the joint distributions of seasonal NDVI and P/T in the current study.The values of the RMSE and AIC are presented in Table 3.According to the minimum criterion of AIC and RMSE, the Clayton copula is the most appropriate copula in fitting the joint distribution of spring NDVI and temperature series, and the joint distribution of autumn NDVI and precipitation series.Whereas, the RMSE and AIC values of the Frank copula in fitting other NDVI and climatic series are the lowest, indicating that Frank copula is the most appropriate copula for other joint distribution of NDVI and temperature/precipitation in the JRB (Table 3).Note: The bold letters represent the most appropriate copula which was chosen to fit the joint distributions of seasonal NDVI and P/T in the current study.

The Identification of Change Points in the Relationship between NDVI and P/T
The copula-based method was adopted to detect the change points of the relationship between seasonal NDVI and P/T in the JRB, and the statistics of the copula-based method are shown in Figure 4.It was found that the statistic (Z values) of the NDVI-P in the three seasons in the 1982-2010 range from 0.1 to 2.1, 0 to 2.8 and 0.5 to 5.2, respectively, and the largest Z values in spring, summer, and autumn were 2.1 (2002), 2.8 (1998), and 5.2 (1988), respectively.Evidently, the largest statistics were all less than the threshold of rejecting null hypothesis (approximately 8.8), thus, no change point was detected between seasonal NDVI and P in the JRB during 1982-2010.For the relationship between NDVI and T in the JRB, the largest Z values in spring, summer and autumn were 5.2 (1992), 1.9 (1994), and 3. 3 (1990), respectively.Similarly, no statistic was larger than the threshold of rejecting the null hypothesis.Therefore, the relationships between seasonal NDVI and P/T are still stationary under the backdrop of the changing environment.Although there were no change points between autumn NDVI and P and spring NDVI and T, their relationships might have experienced slight variations due to the relatively large Zmax value which was 5.2 for both autumn NDVI-P and spring NDVI-T.Generally, the relationships between seasonal NDVI and P/T in the JRB fluctuated with no remarkable changing point during 1982-2010.

Methodology
Climate change (P and T) has proven to be the primary driving force inducing the vegetation variation.For a long-term period, if there is no disturbance, the relationship between climate change and vegetation variation will remain in a steady-state.However, it is evidenced that temperature and precipitation are fluctuating unsteadily under a changing environment [7][8][9][10][11].Therefore, the steady-state of the relationship between climate change and vegetation tends to be unstable.Most previous studies have focused on the statistical methods to identify the relationships without considering their non-linearity and non-stationary [22][23][24].In the current study, a copula-based method was proposed to determine whether the relationship between climate change and vegetation has changed.Moreover, previous studies have adopted this copula-based method to the identification of the change points of the relationship between annual runoff and sediment series, annual rainfall and runoff [56].
In order to further verify the reliability of the results obtained by the copula-based method, the double mass curve, which is known as a relatively simple and practical method to investigate the consistency and tendency of the hydrometeorological series [56], was employed in the present study to further check the variations of the relationships between seasonal NDVI and P/T (Figure 5).It was obvious that the double-mass curve for the spring NDVI and P experienced a larger slope during 2003-2010 than during 1982-2002, which indicated that the relationship between spring NDVI and P in 1982-2002 was different from that in 1988-2010.This finding further verified that the time corresponding to the maximum Z value was the year when the relationship between spring NDVI and P began to change.Similarly, it can be observed in Figure 5 that the slope of the double mass curve for spring NDVI and T in 1992-2010 was smaller than that in 1982-1991.Thus, the relationship between spring NDVI and T began to vary in 1992, which was highly consistent with the year when Z value peaked (5.2).For the relationship between autumn NDVI and P, it can be also observed that the relationship for 1988-2010 was different from that for 1982-1987 as the slopes of the double mass curve in the two periods were different, and 1988 happened to be the year when the Z value peaked (5.2).Conversely, the slopes of the double mass curves of summer NDVI and P/T, autumn NDVI and T did not exhibit any variations during 1982-2010, indicating that the relationships of them were stationarity from 1982 to 2010.Meanwhile, the Z values of the relationships between summer NDVI and P/T and autumn NDVI and T were found to be relatively small compared to other relationships analyzed above.Therefore, the results obtained by the double-mass curve were consistent with that obtained by the copula-based method.Consequently, according to the double mass curve, the reliability and accuracy of the copula-based method applied in this study was further verified.

The Climatic Drivers for the Variations of the Relationship between NDVI and P/T
It has been proven that the large-scale atmospheric circulation anomaly and sunspots exhibited strong associations with P and T at regional and global scales [57,58].Besides, SM and PET play important roles in vegetation dynamics [59,60].Thereby, both the climatic and teleconnection factors might disturb the relationship between NDVI and P/T.To further identify the potential drivers at seasonal scale, two climatic variables (including SM and PET) and four teleconnection factors (including AO, PDO, Niño3.4,and sunspots) were selected and the Pearson correlation analysis was performed to investigate the detailed linkages between climatic/teleconnection factors and the relationship between NDVI and P/T (represented by the Z values which were obtained in Section 3.4).Specifically, the interaction between teleconnection factors and Z value series varied over time scales.Therefore, teleconnection factors at different time scales had different impacts on the Z value series due to their specifically strongly periodic fluctuations.Based on the above, the continuous wavelet transform analysis was utilized to detect the periodicity of long-term AO, PDO, Niño3.4,and sunspots series [61].Results indicated that AO, PDO, Niño3.4,and sunspots series had 14a, 7a, 18a, and 11a periodicity, respectively.In order to eliminate the possible influence of the periodicity of the teleconnection factors on the results, we thereby took 14a, 7a, 18a, and 11a as the corresponding length of the moving windows to generate new stationary series of AO, PDO, Niño3.4,and sunspots which synchronized with the Z value series coving 1987-2005.Accordingly, the Pearson correlation coefficients between the climatic/teleconnection factors and Z value series were calculated to trace the sensitivity of the Z value to different climate variations, and the results are presented in Table 4.
Obviously, the spring Z NDVI-P series had statistically significant correlations (p < 0.01) with AO, PDO, Niño3.4,as well as sunspots, while the spring Z NDVI-T series was significantly correlated with AO and PDO indices (p < 0.01).Similarly, there were certain correlations between the autumn Z NDVI-P series and all these teleconnection factors (p < 0.01), and the correlations were characterized by negative values for PDO and positive values for AO, Niño3.4,and sunspots.Moreover, for the autumn Z NDVI-T series, AO, PDO, and Niño3.4 showed remarkable correlations with it (p < 0.01).These findings indicated that relationships between NDVI and P/T were closely associated with teleconnection factors in the JRB.Nevertheless, summer Z NDVI-P series was not only significantly related to PDO (p < 0.01), but also significantly correlated with PET (p < 0.05).Additionally, the relationships between summer Z NDVI-P series and PDO/PET were positive.This implied that climatic factors such as PET could promote the variation of the relationship between NDVI and P in summer.It was notable that all the six selected factors had no significant correlation with the summer Z NDVI-T series.This finding directly demonstrated that variation of the relationship between summer NDVI and T had no significant correlations with climatic and teleconnection factors.That is, other factors might affect the variation of the relationship between NDVI and P/T during 1982-2010 due to their complexity and uncertainty which deserves further study in the future.According to this research, relationships between seasonal NDVI and P/T were sensitive to the teleconnection factors over the past three decades, especially in spring and autumn.This result is consistent with those of related studies, which have shown clearly the essential roles that teleconnection factors played on vegetation variations in spring/autumn vegetation and climate variations.Anyamba and Eastman [62] found that strong spatial and temporal teleconnection patterns exist between ENSO related climate anomalies and patterns of variation in NDVI.Cho et al. [58] revealed that 17% of the spring vegetation variance was explained by the previous winter's AO variations over the total northern high latitudes, emphasizing the importantly predictive role that winter AO on the vegetation greenness or dynamics in subsequent spring.Moreover, global climate anomalies (such as displacements in the rainfall patterns) were proven to be teleconnected with the anomalous warming of oceanic waters and changes in the Walker circulation [62].Therefore, the large-scale atmospheric circulation anomaly which coupled the ocean-atmosphere system has a significant influence on climate change.Gong and Shi [63] have addressed that large-scale climate systems such as ENSO and AO can significantly influence the inter-annual variations of regional temperature, especially the temperature in spring.Previous studies have demonstrated that temperature change could result in noticeably biological consequences, including lengthened growing season, promoted vegetation activity, and greater biomass productivity [64,65].Therefore, by changing the regional temperature, the large-scale climate fluctuations affect the vegetation coverage.Consequently, the direct impacts of teleconnection factors on climate change and indirect impacts of teleconnection factors on vegetation coverage jointly affect the relationship between NDVI and P/T.
Observations during 1982-2010 have shown a rising trend in seasonal temperature in the JRB, especially in summer.The increase in summer temperature could accelerate evapotranspiration and decrease soil moisture, resulting in drought and inhibition of vegetation growth.Additionally, the JRB is widely known as a water-limited region and the vegetation in water-limited region is sensitive to precipitation [7,66].However, the precipitation in summer covering 1982-2010 in the JRB has a non-significant trend.Therefore, in the context of poor precipitation and high temperature in summer, the relationship between NDVI and P/T in summer is more sensitive to high intensity evapotranspiration and soil moisture loss than in other seasons.
Briefly, considering the important role climatic and teleconnection factors play on vegetation and P/T variations, the deep exploration of their relationships helps better reveal and understand their evolutionary characteristics, thus contribute to more robust predictions of future climate and ecosystem.4.3.The Anthropogenic Drivers for the Variations of the Relationship between NDVI and P/T Although climatic factors have remarkable impacts on the variations of the relationship between NDVI and P/T, the important role of human activity on vegetation dynamics cannot be overlooked which tend to influence the relationship between NDVI and P/T.
Based on the land use data and map illustrated in Table 5 and Figure 6, it was easily found that farmland, forestland, and grassland were three main land types in the JRB which together accounted for more than 95% of the total area during 1980-2010.However, their respective areas have experienced fluctuated variations between 1980 and 2010.It has been reported that the population in the JRB increased rapidly during 1990-1995 [67], resulting in the inevitably increased demand for farmland and the accelerated urban expansion [68].Therefore, the grassland and water area were converted to farmland to meet the food demand in the context of population explosion.It can be obviously seen from Table 5 that grassland and water area decreased during 1990-1995, while farmland and construction land increased significantly.On the other hand, the rapid growth of population and urban expansion could also lead to overgrazing and excessive wood cutting, and thus reducing the forestland area during 1990-1995.Consequently, human activities mainly played a destructive role in vegetation, resulting in vegetation degradation in the JRB during 1990-1995.It is well known that the JRB has experienced intensive soil erosion, with an average of nearly 2.6 × 10 8 t of sediment transported into the JRB per year, accounting for 14% of the annual sediment load in the YRB [37].Hence, a series of water and soil conservation measures have been implemented on the Loess Plateau which contains the JRB to reduce the soil erosion rate, especially the "Grain for Green" program (GGP) launched by the Chinese government in 1999 [69].The It is well known that the JRB has experienced intensive soil erosion, with an average of nearly 2.6 × 10 8 t of sediment transported into the JRB per year, accounting for 14% of the annual sediment load in the YRB [37].Hence, a series of water and soil conservation measures have been implemented on the Loess Plateau which contains the JRB to reduce the soil erosion rate, especially the "Grain for Green" program (GGP) launched by the Chinese government in 1999 [69].The improvement of vegetation coverage and the alteration of land cover patterns jointly affected the land surface hydrological processes through intercepting precipitation and promoting the infiltration rate [70].The current study confirmed that the farmland area exhibited a decreasing trend during 1995-2005, while the forestland and grassland area increased during 1995-2005 (Table 5).This is in accordance with the previous research on the Loess Plateau [71,72].Xiao [71] found a significant GGP-induced increase in the leaf area index and tree coverage on the Loess Plateau.Tang et al. [72] noted a significant decrease in the farmland areas and a significant increase in the forest areas on the Loess Plateau.Therefore, human activities such as the GGP could have obviously promoted the vegetation coverage during 1995-2005.
In conclusion, both the negative and positive human disturbances affect vegetation growth, and their impacts on vegetation were more direct and efficient than climatic factors.Besides, human activities indirectly regulate the climate of a basin by affecting vegetation growth.Therefore, the variations of relationship between NDVI and P/T during 1982-2010 were closely associated with human activities in the JRB.Furthermore, based on the effective irrigated area data (EIA), the correlation coefficients between seasonal Z value series and EIA were calculated and the results are presented in Table 4. Obviously, the EIA was significantly positively correlated with the summer Z NDVI-P series and autumn Z NDVI-T series (p < 0.01), indicating that human activities such as EIA could promote the relationship changes between NDVI and P/T in summer and autumn.Consequently, considering anthropogenic disturbances will be helpful for deeper investigating the possible non-stationarity of the relationship between vegetation growth and climate change under a changing environment, thus revealing the detailed attributions.

Conclusions
Climate disturbance played a crucial role in vegetation dynamics under the background of a changing environment.Therefore, understanding the variation of their relationships and identifying the change points of the relationships was an important basis to evaluate the potential influence of climate change on ecosystems which will support regional land planning and management.
In this study, temporal variations of NDVI, P, and T were investigated in the JRB which is widely known as a typical arid and semi-arid region in China.According to the MMK test, the seasonal T and autumn NDVI exhibited a significant increasing trend, whilst P showed an insignificant trend during 1982-2010 in the JRB.Besides, a copulas-based method was introduced in present study to detect the possible change points between NDVI and P/T.Simultaneously, the double-mass curves method was adopted to verify the reliability of the results acquired by the copula-based method.The two methods consistently showed that the relationships between spring NDVI and P/T and autumn NDVI and P basin have slightly changed, while the relationships of summer NDVI-P/T and autumn NDVI-P have not strikingly changed in the JRB.Generally, no significant change point was identified from the relationship between seasonal NDVI and P/T.We further revealed the dominant drivers of the fluctuations of the relationships between NDVI and P/T at seasonal scales by considering the large-scale atmospheric circulation systems.On average, teleconnection factors have significant impacts on the relationship between seasonal NDVI and P/T.Additionally, the human-induced vegetation degradation and vegetation greenness showed a remarkable influence on the variations of the relationship between NDVI and P/T.
In conclusion, a new framework was proposed in the present study to explore the non-stationarity of the relationship between seasonal NDVI and P/T, which could help effectively and reliably investigate the sensitivity of vegetation dynamics to climate change, thus conducting ecological restoration on vegetation cover at regional scales.

Figure 1 .
Figure 1.Location of the Jing River Basin.(a) The basin map of the Jing River Basin; (b) location of the Jing River Basin in the Yellow River Basin, China.

Figure 2 .
Figure 2. Flowchart of the framework for the change point identification.

Figure 4 .
Figure 4.The statistics of copula-based method in detecting the change points of the relationship between seasonal NDVI and P/T during 1982-2010 in the JRB.

Figure 5 .
Figure 5.The double mass curves of seasonal NDVI and P/T covering 1982-2010 in the JRB.(a) and (b) represent the double-mass curve of spring NDVI-P/T; (c,d) represent the double-mass curve of summer NDVI-P/T; (e,f) represent the double-mass curve of autumn NDVI-P/T.

Figure 6 .
Figure 6.Land use map in different periods of the Jing River Basin.

Figure 6 .
Figure 6.Land use map in different periods of the Jing River Basin.

Table 1 .
The results of the modified Mann-Kendall (MMK) trend test.

Table 2 .
The goodness-of-fit values of the three probability distributions in fitting the distribution of individual seasonal NDVI and P/T series.

Table 2 .
Cont.H represents the hypothesis test result, returned as a logical value, if H = 0, indicating a failure to reject the null hypothesis at 95% confident level, if H = 1, indicating the rejecting of the null hypothesis at 95% confident level; P denotes the p-value of the test, returned as a scalar value in the range (0, 1), small values of P cast doubt on the validity of the null hypothesis.The bold letters represent the most appropriate marginal distributions. Note:

Table 3 .
The Root-Mean-Square Error (RMSE) and Akaike Information Criterion (AIC) value of Clayton, Frank, and Gumbel copulas in fitting the joint distribution of seasonal NDVI and P/T.

Table 4 .
The correlations between climatic/teleconnection factors and Z value series in the JRB.Note: "*" and "**" represent significant at 95% and 99% confidence level, respectively.The bold numbers represent the significant correlations between climatic/teleconnection factors and Z value series.

Table 5 .
Characteristics of inter-annual variations of land-cover in the JRB.