Response of Vegetation to Changes in Temperature and Precipitation at a Semi-Arid Area of Northern China Based on Multi-Statistical Methods

: Hydrothermal and climatic conditions determine vegetation productivity and its dynamic changes. However, the legacy e ﬀ ect and the causal relationships between these climatic variables and vegetation growth are still unclear, especially in the dry regions. Based on multi-statistical methods, including bivariate correlation analysis and composite Granger causality tests, we investigated the correlation, causality, and lag length between temperature / precipitation and the vegetation growth (Normalized Di ﬀ erence Vegetation Index, NDVI) in three typical sub-watersheds in the Luanhe River Basin, China. The results show that: (1) Precipitation and temperature are the Granger causes of NDVI variation in the study catchment; (2) temperature and precipitation are not strictly positively correlated with NDVI during growing seasons along with the whole sequence, and excessive warmth and precipitation inhibits vegetative growth; (3) the lag length of vegetation growth in response to temperature / precipitation was shorter in agriculture areas (~2 months) than the forest-dominant area, which have indicated 3–4 months lag length; and (4) anthropogenic disturbance did not result in notable negative e ﬀ ects on vegetation growth at the Luanhe River Basin. Our study further suggests that use of these multi-statistical methods could be a valuable approach for comprehensively understanding the correlation between vegetation growth and climatic variations. We have also provided an avenue to bridge the gaps between stationary and non-stationary sequence, as well as to eliminate pseudo regression problems. These ﬁndings provide critical information for developing cost-e ﬃ cient policies and land use management applications for forest conservation in arid and semi-arid area.


Introduction
Precipitation and temperature are recognized as the two major climatic factors determines the vegetation biophysical processes [1]. To quantify how climatic factors affect vegetation growth at larger scale, the Normalized Difference Vegetation Index (NDVI), derived from infrared and near-infrared spectral band [2], has been widely applied as a proxy of vegetation growth [3]. The climate-vegetation-runoff correlation has been discussed based on the NDVI at the semi-arid watershed, Luanhe River Basin, China [4], where the vegetative activities are strongly affected by seasonal precipitation and temperature [5,6]. However, the relative importance and causal relationships between precipitation and temperature on vegetation growth are still unclear. Therefore, understanding the relationship between the changes of climate factors and vegetative growth is essential for further improvement of our understanding on ecosystem response to climate change and for the forest conservation management [5].
Previous studies reported that the annual vegetation growth showed a positive relationship with climate warming in most parts of China [6], but the response sensitivity differed by seasons [7] and vegetation types [8]. Similar findings have been found around the world. For example, Lamchin has pointed out that in most of Asia temperature is the main factor promotes vegetation growth [9], and the same relationship has been found for the majority of the Nigeria [10]. Admittedly, these studies have mainly focused on the correlative relationships between climatic factors and NDVI using classical mathematical methods, for example Pearson and Spearman correlation coefficients. However, these mathematical methods were introduced to reveal the general trend of the whole time sequence data, and require the sequence to be monotonous [11]. Therefore, whether the relationship between vegetation and climatic factors is positive or negative among the whole sequence or only within a specific range have not been fully discussed.
Because of the complexity of the ecosystem and its biochemical processes, the NDVI and climatic factors are often characterized as nonlinear or non-stationary sequence. Therefore, the classical mathematical method is not sufficiently applicable to reveal the nonlinear characters of ecosystem. Researchers have applied algorithms such as the Mann-Kendall test [12] that divide the non-stationary time sequence into several stationary stages to study the nonlinear trend of ecosystem. While the Copula method, which has been widely applied in econometrics to analyze the correlation between different variables, can be used to analyze the non-stationary sequence [13], the approaches and attempts for the researchers who are discerning more general principles of the holistic relationship between the two sequences is still lacking.
Additionally, correlation analysis is not identical with the causal relationship, and the legacy effect of climatic variables on the corresponding vegetation growth needs to be further explored. For example, temperature, precipitation and relative humidity have been recognized to be associated with vegetation changes of 26%, 49%, and 23%, respectively, in China, on the loess plateau [14]. However, the legacy effect of the response time has not been widely discussed. The Granger causality test, which was pioneered by Clive W. J. Granger [15], provides an approach to test the causal relationship between climate variables and its consequences. The Granger causality test is concerned with not only the correlation between different sequences but also the causal relationships of one factor to the other [16].
In this study, the Luanhe River Basin (LRB), which is located at a semi-arid climate region and characterized as an ecological conservation area, was selected as the study catchment. The Copula method and a composite Granger causality test (with Engle-Granger cointegration test and vector autoregression model) were used to fully describe how vegetation responded to the changes in temperature and precipitation. The aims of this study are to (1) address the reliability of the multi-statistical method in ecology; (2) quantify the correlation and causal relationship between vegetation and climatic variables; (3) test for difference of these correlations and causality in different seasons and vegetation types. Based on these findings, some suggestions are proposed for policymakers and stakeholders who have focused on ecosystem restoration and forest conservation.

Study Area
The Luanhe River Basin (LRB), covering a total area of 44,793.51 km 2 , as shown in Figure 1, was selected as the study catchment. The LRB is located northeast of Beijing, China (115 • 30 E-119 • 45 E, 39 • 10 N-42 • 40 N), and is characterized as an important ecological barrier to alleviate the effect of sandstorms from Mongolia and conserve water resources for the Beijing-Tianjin-Hebei region. The Luanhe River is not a tributary of the Hai River, but an independent river following the Hai River into the sea. The elevation ranges from 3 to 2241 m and decreases from the mountainous areas in the northwest to a plain in the southeast. The LRB is situated in a typical temperate continental monsoonal zone with semi-arid climate [17], where is more sensitive to climatic and anthropogenic disturbance. The study catchment has an annual temperature fluctuating from 0 • C to 9.0 • C and 478.5 mm of average annual precipitation. Heavy rainfall and relatively high temperature occur in summer, and less rainfall and low temperature occur in winter.
Forests 2020, 11, x FOR PEER REVIEW 3 of 15 mm of average annual precipitation. Heavy rainfall and relatively high temperature occur in summer, and less rainfall and low temperature occur in winter. To provide more specific ecological conservation management strategies based on the characteristics of vegetation types in the study catchment, sub-watersheds with different landscapes should be considered. Therefore, three sub-watersheds were identified based on the landcover characters: Upstream (6405 km 2 , 15.25% of the total area), midstream (18,394 km 2 , 43.57% of the total area), and downstream (19,807 km 2 , 41.18% of the total area) ( Figure 1).
In particular, agricultural land in the upstream (31.33% in 1980 and 32.38% in 2015) had the greatest contribution among the sub-watersheds. As part of the ecological conservation area, the midstream and downstream were mostly covered by forest and herbaceous vegetation and did not change during 1980-2015. The forest covers about 45% and 44% at midstream and downstream, respectively, and the herbaceous covers about 29% in both these two sub-watersheds. The urban area was mainly situated in the downstream area, with 1.64% in 1980 and 2.27% in 2015 ( Figure 2). These three sub-watersheds can be categorized by their typical land use types as agriculture area (upstream), forest conservation area (midstream), and forest-urban area (downstream). To provide more specific ecological conservation management strategies based on the characteristics of vegetation types in the study catchment, sub-watersheds with different landscapes should be considered. Therefore, three sub-watersheds were identified based on the landcover characters: Upstream (6405 km 2 , 15.25% of the total area), midstream (18,394 km 2 , 43.57% of the total area), and downstream (19,807 km 2 , 41.18% of the total area) ( Figure 1).
In particular, agricultural land in the upstream (31.33% in 1980 and 32.38% in 2015) had the greatest contribution among the sub-watersheds. As part of the ecological conservation area, the midstream and downstream were mostly covered by forest and herbaceous vegetation and did not change during 1980-2015. The forest covers about 45% and 44% at midstream and downstream, respectively, and the herbaceous covers about 29% in both these two sub-watersheds. The urban area was mainly situated in the downstream area, with 1.64% in 1980 and 2.27% in 2015 ( Figure 2). These three sub-watersheds can be categorized by their typical land use types as agriculture area (upstream), forest conservation area (midstream), and forest-urban area (downstream). In our research, the growth seasons include late spring (April to May), summer (June to August) and early autumn (September to October). The winter has been ignored, because of the low vegetation growth rate during chilly climate in the study catchment [18]. The average NDVI value during summer, early autumn and late spring was 0.611, 0.459, and 0.354, respectively. The NDVI for the three sub-watersheds from higher to lower is in the order of downstream, midstream and upstream, which is directly opposite to the elevational change. Moreover, most of the NDVI in the LRB from 1982 to 2015 has shown an increasing trend in different sub-watersheds during the growing seasons, except for the upstream during late spring ( Figure 3). Meanwhile, the upstream also showed the smallest increasing tendency compared to that of the other sub-watersheds during all three seasons ( Figure 3).  In our research, the growth seasons include late spring (April to May), summer (June to August) and early autumn (September to October). The winter has been ignored, because of the low vegetation growth rate during chilly climate in the study catchment [18]. The average NDVI value during summer, early autumn and late spring was 0.611, 0.459, and 0.354, respectively. The NDVI for the three sub-watersheds from higher to lower is in the order of downstream, midstream and upstream, which is directly opposite to the elevational change. Moreover, most of the NDVI in the LRB from 1982 to 2015 has shown an increasing trend in different sub-watersheds during the growing seasons, except for the upstream during late spring ( Figure 3). Meanwhile, the upstream also showed the smallest increasing tendency compared to that of the other sub-watersheds during all three seasons ( Figure 3). The monthly climatic data from ten meteorological stations, including precipitation and temperature, from 1982 to 2015 was provided by the National Meteorological Administration of China (http://data.cma.cn). The monthly Normalized Difference Vegetation Index (NDVI) data (0.05 degree resolution) from 1982 to 2015 was obtained from the Advanced Very High Resolution Radiometer (AVHRR) of US National Oceanic and Atmospheric Administration (NOAA) (https://nex.nasa.gov) to quantify the variation of vegetative cover changes in the LRB.
In our research, the growth seasons include late spring (April to May), summer (June to August) and early autumn (September to October). The winter has been ignored, because of the low vegetation growth rate during chilly climate in the study catchment [18]. The average NDVI value during summer, early autumn and late spring was 0.611, 0.459, and 0.354, respectively. The NDVI for the three sub-watersheds from higher to lower is in the order of downstream, midstream and upstream, which is directly opposite to the elevational change. Moreover, most of the NDVI in the LRB from 1982 to 2015 has shown an increasing trend in different sub-watersheds during the growing seasons, except for the upstream during late spring ( Figure 3). Meanwhile, the upstream also showed the smallest increasing tendency compared to that of the other sub-watersheds during all three seasons ( Figure 3).

Data Pre-Treatment
In our study, we address the climate factors (precipitation and temperature) and the vegetation growth (NDVI) at monthly scale in the three sub-watersheds. To achieve this research goal, the data pre-treatment was procedure as follows: first, for the climate factors, we interpolate the monthly precipitation and temperature data from the meteorological stations and re-summarized the climate data as the average value for each sub-watershed; second, in order to match the format of climate data, the raster NDVI was also categorized as the average value at monthly scale and calculated into different sub-watershed using the GIS (zonal statistics tool).

Methods
A multi-statistical method is proposed in our study to characterize how vegetation changes in response to temperature and precipitation. As the flowchart shows in Figure 4, the procedure is as follows: (a) The collection and pre-treatment of vegetation and climatic data, (b) bivariate correlation assessment by the tail features based on the different fitted copulas, (c) using composite Granger causality test to detect both the causal relationship and legacy effect between NDVI-T/P, (d) take the different sub-watersheds and growth seasons into consideration, we address several implications for policymakers and stakeholders. The detailed instructions for the statistical methods are described in the "Methods" section of supplementary information. In our study, we address the climate factors (precipitation and temperature) and the vegetation growth (NDVI) at monthly scale in the three sub-watersheds. To achieve this research goal, the data pre-treatment was procedure as follows: first, for the climate factors, we interpolate the monthly precipitation and temperature data from the meteorological stations and re-summarized the climate data as the average value for each sub-watershed; second, in order to match the format of climate data, the raster NDVI was also categorized as the average value at monthly scale and calculated into different sub-watershed using the GIS (zonal statistics tool).

Methods
A multi-statistical method is proposed in our study to characterize how vegetation changes in response to temperature and precipitation. As the flowchart shows in Figure 4, the procedure is as follows: (a) The collection and pre-treatment of vegetation and climatic data, (b) bivariate correlation assessment by the tail features based on the different fitted copulas, (c) using composite Granger causality test to detect both the causal relationship and legacy effect between NDVI-T/P, (d) take the different sub-watersheds and growth seasons into consideration, we address several implications for policymakers and stakeholders. The detailed instructions for the statistical methods are described in the "Methods" section of supplementary information.   1) and (2) are the two components for multi-statistical method).

Copula Based Bivariate Correlation Assessment
The most commonly used Copula functions, which were designed to assess bivariate joint dependence structures, have been divided into two families, elliptical (Gaussian and Student's t-Copula) and Archimedean (Clayton, Gumbel, and Frank). Because elliptical Copulas are derived from the well-known distribution types and widely used Pearson's correlation, the elliptical Copula has become popular in most fields. However, the elliptical is only able to reproduce an expression  1) and (2) are the two components for multi-statistical method).

Copula Based Bivariate Correlation Assessment
The most commonly used Copula functions, which were designed to assess bivariate joint dependence structures, have been divided into two families, elliptical (Gaussian and Student's t-Copula) and Archimedean (Clayton, Gumbel, and Frank). Because elliptical Copulas are derived from the well-known distribution types and widely used Pearson's correlation, the elliptical Copula has become popular in most fields. However, the elliptical is only able to reproduce an expression that is symmetrical [19]. The Archimedean Copulas (Table 1), including the Gumbel, Frank, and Clayton Copulas, have an explicit formula and are quite popular given their ability to capture a wider variety of joint dependence structures [20]. Table 1. Equations for the Archimedean Copula functions, where u and v are defined as two marginal distribution and θ is the dependence parameter.

Family C (u, v) Parameter
In our study, we use Akaike information criterion (AIC) and Root Mean Square Error (RMSE) to select best fitted Archimedean Copulas. The smaller the AIC and RMSE value, the better results for the Copula to interpret the correlation between vegetative growth and climatic factors. Moreover, we use diverse tail features of different Archimedean Copula to detect the bivariate correlation between NDVI and climatic factors (Figure 4). For each Archimedean Copula, the distribution shows unique features. If two variables can be described by the Clayton Copula, there is a strong correlation between the variables at the lower tail, while the variables are gradually independent at the upper tail. In contrast to the Clayton Copula, the Gumbel Copula shows a strong correlation between variables at the upper tail, while at the lower tail the variables are gradually independent. The Frank Copula has symmetrical tail characteristics, which means it cannot capture the asymmetrical tail characteristics between random variables. The correlation coefficient of the variables at the tail is zero, which indicates that the variables at both tails are gradually independent [19].

Composite Granger Causality Test for Legacy Effect and Casual Relationship
Originally, the stationary sequence has been required to operate the Granger causality test. To extend the Granger causality test to fit a more general data sequence, Engle-Granger cointegration test (EG test) has been regarded as an irreplaceable premise to broaden the scope of applicability in our study [21]. The Engle-Granger cointegration test can be used to prove two non-stationary data sequences may have stationary relation over long time period. In order to examine whether the data sequence has passed the EG test, the Durbin-Watson (DW) coefficient was introduced, when the DW approaches 2, the data sequence would be recognized to obey the normal distribution and can be recognized to have stationary relation. In addition, the vector autoregression (VAR) model [22] has been applied to investigate the legacy effect between NDVI and climatic factors, using the Akaike information criterion (AIC). Although the Granger causality test can be used on any lag length between these sequences, the significance of the Granger causality test is better based on the optimal lag length which was indicated by the minimum AIC.
The overall processes are listed in Figure 4: (1) The sequences have been first tested by the Augmented Dickey-Fuller Test (ADF) for stationarity, (2) the non-stationary sequence has to been further examined by Engle-Granger cointegration test, while the stationary sequence can be analyzed by the VAR model directly to disclose the legacy effect, and (3) the optimal leg length detected by the VAR model would be used to conduct the Granger causality test.

Nonlinear Bivariate Correlation between NDVI and Precipitation/Temperature
When the most appropriate marginal distribution has been determined, the Root Mean Square Error (RMSE) and Akaike information criterion (AIC) were used as two criteria to quantify the deviation between these data sequences. The most appropriate Copulas from Archimedean family (including the Clayton, Frank, and Gumbel Copulas) was selected to identify the relationship between seasonal NDVI-P and NDVI-T. The results of Copula selection are listed in Table 2. And the results of the premise process for identifying the marginal distribution has been mentioned in the supplementary Table S1.
Based on the AIC and RMSE values, the Clayton and Frank Copulas are the most appropriate Copulas in fitting the joint distribution of seasonal NDVI-P and NDVI-T. Particularly, the RMSE and AIC values of the Frank Copula in fitting the NDVI-T series during the late spring and summer are the lowest, indicating that NDVI-T can be well described by the Frank Copula during this period. Likewise, the Clayton Copula showed better agreement with the joint distribution of NDVI-T during early autumn in the upstream. Whereas, different from the homogeneity of the Copula for NDVI-T during late spring and summer, the Copulas varied by seasons and sub-watershed for NDVI-P. For NDVI-P, the Frank Copula was better to fit the joint distribution for all three sub-watersheds during the summer. In contract to the other sub-watersheds, the Frank and Clayton Copula was employed to describe the relationship between NDVI-P in upstream during the late spring and early autumn, respectively.

Time Series Sequence Identification
The premise of the composite Granger causality test is that the data sequence needs to be stationary or have two data sequences that are of the same order to meet the requirement of a cointegration relationship [23]. If the vegetation and climatic data sequence passes the ADF test, then the sequence has been recognized as stationary. However, if the sequences do not pass the ADF test, then the Engle-Granger cointegration test is introduced to determine whether there is a cointegration relationship between the NDVI and precipitation/temperature. The time sequence identification results for the NDVI, precipitation and temperature are listed in Table 3. The highlighted NDVI sequence in the midstream is non-stationary with a p-value greater than 0.05, while the other sequence for these three factors is stationary. To identify whether the non-stationary sequence of the NDVI has a cointegration relationship with precipitation or temperature, the Engle-Granger cointegration test was conducted ( Figure 5). The Durbin-Watson (DW) coefficient for NDVI-P and NDVI-T is 1.56 and 1.89, respectively. The specific statistical results are listed in supplementary Table S2. The criteria for whether there exists a cointegration relationship, depends on whether the DW coefficient is closer to 2 [24]. According to Figure 5, the NDVI-P has a relatively weak cointegration relationship in the midstream. Therefore, there would be no obvious causal relationship at the midstream between NDVI-P.

Optimal Lag Length and Granger Causality
After confirming the stationary and Engle-Granger cointegration relationship between the sequences of the variables, the Granger causality test can be used to discriminate the interactive relationship among the variables. In our study, the optimal lag length between the NDVI-T and NDVI-P was identified before the Granger causality test was carried out by the VAR model.
Notably, the minimum AIC value for detecting the optimal lag length varied between sub-watersheds for the NDVI-T/P ( Figure 6). The optimal lag length of the NDVI-T between the upstream, midstream and downstream is 2, 4, and 4 months, respectively. Because the result for the precipitation either shows a stationary or cointegration relationship with NDVI in the midstream, the lag length for NDVI-P has only been identified in the upstream and downstream with 2 and 3 months, respectively. The Granger causality test was conducted at the 95% confidence level. The null hypothesis for each NDVI-T/P in different sub-watershed is that "the climate factors/NDVI is not the Granger cause of the NDVI/climate factors". When p-value is less than 0.05, the null hypothesis would be rejected. The results in Table 4 indicate that for all three sub-watersheds with stationary or modified climate and NDVI data sequences, the temperature/precipitation has been identified as the Granger cause of the NDVI. and NDVI data sequences, the temperature/precipitation has been identified as the Granger cause of the NDVI.

NDVI Response to Temperature
In most cases, the correlation between temperature and NDVI obeys the Frank Copula, which has no significant features for tail distribution. This consistency among different sub-watersheds and seasons suggests that temperature within the specific range in the study catchment has a direct effect on the vegetation growth rate [25]. The features showing no obvious upper or lower tail indicates that different from humid regions where a strictly positive relationship between the NDVI and temperature was found [26], the vegetation growth in a semi-arid area is not strictly positive along

NDVI Response to Temperature
In most cases, the correlation between temperature and NDVI obeys the Frank Copula, which has no significant features for tail distribution. This consistency among different sub-watersheds and seasons suggests that temperature within the specific range in the study catchment has a direct effect on the vegetation growth rate [25]. The features showing no obvious upper or lower tail indicates that different from humid regions where a strictly positive relationship between the NDVI and temperature was found [26], the vegetation growth in a semi-arid area is not strictly positive along with the temperature, while a partially positive correlate relation has been noticed. The possibly reasons are listed below: first, due to the lack of rainfall, excessive warming may accelerate evaporation of soil water content and result in drought, which may lead to a negative effect on the biophysical process of plants.; second, given the excessive temperature, vegetation will prevent its own water loss by reducing leaf area and light saturation point, finally resulting in a decrease in vegetative cover [27]; third, an increase in temperature may promote autotrophic respiration and the transpiration rate, which may result in accelerating the consumption of organic matter and causing a decrease in the net productivity of vegetation [28]. Therefore, we conclude that vegetation growth in this semi-arid area present higher sensitivity to temperature than in humid areas.
Notably, different Copula patterns were found between the upstream (Clayton Copula) and midstream/downstream (Frank Copula) in early autumn. A larger and significant correlation coefficient was found at the upstream where the temperature and NDVI are both relatively low, which obey the Clayton Copula. This may result from the differences in land use type. The upstream region covers a large amount of cropland (>30%), where the NDVI has been largely reduced along with the temperature during early autumn by the agricultural activities, including harvesting and ploughing.
Furthermore, the two to four month lag length of NDVI in response to temperature for distinct sub-watersheds reveals that the vegetation may not have a synchronous response to temperature changes. As vegetation mainly accumulates nutrients and organic matter in spring and grows in summer and autumn [29], the distinct dominant and diversity of the vegetation types would be the main cause of different legacy effect for NDVI in response of the temperature. The agricultural area (upstream) presented a faster responce to temperature than the forest conservation area in the midstream and forest-urban mixed area in the downstream, which may be due to the short growth cycle and rapid growth rate of the crops [30]. The mature forest, protected by the national and local government, requires more energy to growth, which is produced from the photosynthetic processes [31] and finally result in the slowest response rate in the forest-dominant area (midstream and downstream).

NDVI Response to Precipitation
During late spring, the correlation between NDVI and precipitation in the upstream not only shows different fitting features (Frank Copula) to those of the midstream and downstream (Clayton Copula) but also has a short response time (2 months), which indicates that the crops are more sensitive to precipitation for growth in spring. When taking the longer lag length (3 months) for the NDVI in response to precipitation in the downstream into consideration, the forest and herbaceous in the forest-urban area (downstream) may retain a longer period for growth than the crops. These results imply that, in semi-arid area, the crops may have a better ability to take advantage of the restricted available water than forest to boost growth in spring [32]. However, the forest would present a longer time to accumulate nutrient in spring and exhibit stable growing rate during the whole growth season even in autumn [33].
In summer, the correlation between NDVI and precipitation was accumulated within a specific range, following the Frank Copula, in all three sub-watersheds. In contrast to previous studies in central Asia, which indicate a positive correlation between NDVI-P [34]. The fitted Frank Copula in the study catchment suggested that excessive rain in summer can suppress plant growth by limiting root respiration rates, and lead to the reduction on growth and coverage [35], in semi-arid area. Additionally, similar to NDVI-T, the harvest work could be the key reason for reducing the leaf cover and a lower tail relationship (Clayton Copula) between the NDVI and precipitation in the upstream, during the early autumn. Furthermore, the utilization of ground water resource in the upstream could accelerate the depletion of the available water for vegetation growth, especially in the autumn with less rainfall to replenish the ground water in semi-arid area [36]. The lack of a causal relation for the midstream sub-watershed, may due to the great amount of ground water conservation ability, which could help the vegetation to prevent drought or flood for a stable growth rate [30].

Implication
The multi-statistical methods, which expedite identification of the characteristics of the relationship between the NDVI and climatic factors (temperature and precipitation), can be used to aid in understanding the vegetation response to climate change, and thus be helpful for environmental management or investment decisions. To better fulfill the expectation of residents for the natural environment, national and local policy makers are anxious to find an effective way to balance economic development and ecosystem service provision [37]. However, as a result of rapid urbanization and population growth in China [38], it is difficult to establish a pure national park or conservation area without anthropogenic disturbance. Our mathematical work provides several potential clues for identifying a method for proposing cost-effective regulation, in order to better enhance or restore ecological conservation areas. According to the composite Granger causality test, temperature/precipitation has a considerable contribution to the vegetation compared to the impact of the vegetation on the regional temperature/precipitation. Therefore, the proposals for the conservation management are mainly based on the factors, which can be handled by humans to alleviate the negative impact of climate change on vegetation growth, such as the water supply and regulation criteria.
Our research suggests that, compared to the forest conservation and forest-urban mixed area, in the agricultural area in the upstream an important priority is to ensure the crops can acquire sufficient water in spring, especially when drought occurs. In addition, ground water resources play a critical role in vegetation growth, especially in the semi-arid area. Ground water conservation can well stimulate vegetation growth, when there is less rainfall. However, ground water deficits would largely affect agricultural efficiency [39]. Therefore, measures, such as utilization of groundwater and surface water for irrigation, as well as the development of a backup water source, should be regulated, in order to guarantee the growth of the crops. Different from forest and herbaceous, crops are more sensitive to water supply to promote growth, more water can motivate the dynamic processes in nutrient and mineral solution in the soil, which can be absorbed by vegetation under water-soluble conditions to yield more organic matter [40]. Then, physiological and biochemical processes such as photosynthesis and transpiration of crops can promote the growth rate [41]. Therefore, further researches and experiments regarding how to quantify the needed water for distinct plants should proceed for agricultural management in future.
Besides the climate factors, anthropogenic disturbance has been identified to be another key driver for vegetative growth [42]. It is acknowledged that over-grazing and over-cultivation were two major causes of reduced vegetative growth, especially in semi-arid areas, where the ecosystem is fragile [43]. Over-grazing would threaten vegetation species composition, variety and richness [44], while the over-cultivation may significantly deplete ground water for both human and crops [36]. In our study catchment, owing to the Grain for Green Program and environmental protection regulation since 1962, grazing and cultivation activities have been limited. In order to honor the great efforts for environmental protection, the LRB is internationally known for the Saihanba National Forest Park, who has won the award for Champions of the Earth in 2017 from the United Nations Environment Programme [45]. Some policymakers think that the ecological immigration, which aims at moving the local people to other regions would be a reliable measurement to regulate the nonnegligible anthropogenic disturbance in the National Park. However, the national and local policy makers and stakeholders can benefit from our findings to carry out more cost-effective regulations for better conserving the ecosystem services. Two indicators in our study have indicated that, for the policy makers, a more cost-efficient regulation to alleviate the negative effects of human activities in the conservation area or a national park would be to provide effective limitations to regulate human activities rather than migrate all residents outside of the conservation area. One indicator is that, as shown in Figure 2, the difference in the NDVI between the midstream and downstream is relatively narrow, which may be mainly due to the altitude gradient. With an increase in altitude, the temperature, atmospheric pressure and partial pressure of CO 2 would decrease [46], which result in the lower organic synthesis rate [31]. The other indicator is the consistency of the NDVI change trend, fitted Copula function and the same lag length identification between the midstream and downstream. These facts may be among the indicators that when appropriate regulation has been carried out, limited urbanization or anthropogenic disturbance would not have largely negative effects on the conservation area. Several researches have indicated that the limited and regulated urbanization processes, including low impact development, would strengthen the resistant ability for residential areas and reduce the negative impact on ecosystems [47]. Therefore, how to quantify and control the accepted human activities and urbanization should be further discussed.
To conclude, the multi-statistical methods proposed in this study provide a general measure to bridge the data gaps between stationary or non-stationary sequence and reject the pseudo regression. In addition, this method can be applied to not only identify the relationship between the NDVI and climatic factors but also detect the legacy effect and causal relationship. Our methods would be suitable for both researchers and policymakers with different intentions. For policy makers, our study implies that (1) the agriculture area should be ensured to have plenty of water prior than the forest-dominant area; and (2) the restricted regulation to limit human activities would be more cost-effective to protect the forest ecosystem in study catchment. For researchers, our research could be meaningful in following aspects: (1) The multi-statistical methods may improve the understanding of the relationship between two sequences, (2) the further studies on water demand for specific crops would help policy makers to carry out appropriate irrigation measures; and (3) how to quantify the accepted human activity in forest-urban mixed area would help local residents to achieve sustainable goals.

Limitation and Uncertainty
In our study, temperature and precipitation data have been analyzed as independent variables to affect the vegetation growth. However, in the natural ecological system, the temperature and precipitation-related climatic variables generally interact and determine plant growth. The independent indicator using in our study may therefor generate uncertainty regarding influences of precipitation and temperature on vegetation growth. Therefore, the interactive indices, such as combination of effect of relative air humidity, soil water content, evaporation rate, and (bio)climatic index, merit further studies, especially experimental studies, for vegetation growth.

Conclusions
This study provides a multi-statistical assessment method to address the characteristics of the NDVI and climatic factors (temperature and precipitation). This methodology transferring from economics to the environmental sciences field is promising and provides a comprehensive and simpler applicable approach to indicate a basis for evaluating the relationship and casualty between vegetative growth and climatic factors. The conclusions of this study are as follows: (1) Except for the NDVI-P in the midstream, which has no presented obvious causality relationship, both temperature and precipitation were identified as the Granger causes for the NDVI changes in most of the study area.
(2) The NDVI in semi-arid area is different with the other climatic zone that is not strictly positive to temperature and precipitation change. The excessive warmth and precipitation would inhibit the vegetative growth during the growing seasons.
(3) Crops show the two months lag length in response to temperature, while other forest-dominated area has presented four months lag in response to the temperature. Thus, crops respond to temperature faster than forest and herbaceous at the semi-arid area.
(4) Based on the NDVI change characteristics and assessment using these two statistical methods, we found that limited human activities in the ecological conservation area would not result in major negative effects on vegetative growth.