An Improved Conceptual Model Quantifying the E ﬀ ect of Climate Change and Anthropogenic Activities on Vegetation Change in Arid Regions

: Vegetation shows a greening trend on the global scale in the past decades, which has an important e ﬀ ect on the hydrological cycle, and thus quantitative interpretation of the causes for vegetation change is of great beneﬁt to understanding changes in ecology, climate, and hydrology. Although the Donohue13 model, a simple conceptual model based on gas exchange theory, provides an e ﬀ ective tool to interpret the greening trend, it cannot be used to evaluate the impact from land use and land cover change (LULCC) on the regional scale, whose importance to vegetation change has been demonstrated in a large number of studies. Hence, we have improved the Donohue13 model by taking into account the change in vegetation cover ratio due to LULCC, and applied this model to the Yarkand Oasis in the arid region of northwest China. The estimated change trend in leaf area index (LAI) is 1.20% / year from 2001 to 2017, which accounts for approximately half of the observed (2.31% / year) by the moderate resolution imaging spectroradiometer (MODIS). Regarding the causes for vegetation greening, the contributions of: (1) LULCC; (2) atmospheric CO 2 concentration; and (3) vapor pressure deﬁcit were: (1) 88.3%; (2) 40.0%; and (3) − 28.3%, respectively, which reveals that the largest contribution was from LULCC, which is probably driven by increased total water availability in whole oasis with a constant transpiration in vegetation area. The improved Donohue13 model, a simple but physics-based model, can partially explain the impact of factors related to climate change and anthropogenic activity on vegetation change in arid regions. It can be further combined with the Budyko hypothesis to establish a framework for quantifying the changes in coupled response of vegetation and hydrological processes to environment changes.


Introduction
Vegetation is one of the most active components of the Earth system and is often considered to be an important indicator of ecosystem change [1][2][3][4]. At the same time, vegetation also has impacts on hydrological processes, including direct effects such as root water uptake and transpiration, as well as indirect effects such as rainfall interception, soil infiltration, and hillslope runoff [5]. Using satellite observations, a number of studies have shown that global greening has become an indisputable fact in past decades [6][7][8][9][10][11]. There are many factors leading to vegetation change, including meteorological variables such as precipitation, temperature, and water vapor pressure [12][13][14], anthropogenic activity variables such as land use change and grazing [12,[15][16][17], and nutritional variables such as carbon and nitrogen [18][19][20]. However, in the context of global change, we still do not have a clear understanding on the attribution of vegetation greening, which is not conducive for predicting future vegetation changes or for understanding the feedback of vegetation changes on hydrology and climate [21].
Recently, numerous studies have shown that vegetation changes in arid regions are more sensitive to global change and anthropogenic activities than in other climatic regions [22][23][24][25], so arid regions have become natural experimental sites for understanding vegetation changes. Fensholt et al. [26] reported an increase trend in vegetation change together with their drivers in global semi-arid regions. Wang et al. [17] used three statistical methods to conclude that the vegetation change in the mountain of the Tarim River Basin was limited by temperature, while that in the oasis was controlled by runoff. Dardel et al. [27] used a linear regression method to analyze the causes of vegetation changes in the Sahel region in the past 30 years. To quantitatively explain vegetation changes, Wang et al. [17] used a simple conceptual model to quantitatively analyze the impact of change in the irrigation area and desert area on vegetation cover change in the Heihe River basin. Later, Shen et al. [15] improved that model by incorporating factors such as the gross domestic product (GDP) and precipitation, thus quantitatively explaining the impact of anthropogenic activities and climate change factors on vegetation change in the middle reaches of the Heihe River. However, the model used in Shen et al. [15] Based on a linear regression method, which still lacks a reliable physical basis. Differently, Donohue et al. [12] developed a physically based conceptual model from the perspective of water use efficiency (WUE) and explored the effect of CO 2 fertilization on the maximum vegetation coverage in global arid regions. However, the model is based on the underlying assumption that the vegetation coverage is uniform in space, and in another word, it is a point scale model. Therefore, it cannot take into account the effect of change in vegetation cover ratio due to land use and land cover change (LULCC). In fact, many studies have shown that LULCC plays an important role in vegetation change. For example, Bégué et al. [28] showed that land use change was one of the main drivers of NDVI change in the Bani catchment from 1982 to 2006. Brandt et al. [29] showed that vegetation changes in the Sahel were mainly caused by tree-and land-cover change. Wang et al. [17] showed that LULCC was the main contributor to the vegetation change in the Heihe River basin. Therefore, our objectives here are threefold: (1) To identify changes in the leaf area index (LAI) and its influencing factors, such as water availability, temperature, atmospheric CO 2 concentration, vapor pressure deficit, and LULCC, in the Yarkand Oasis; (2) to improve the Donohue13 model so that it can quantitatively explain the relationship of vegetation changes with the change in LULCC; and (3) to apply the improved Donohue13 model to quantitatively attribute the change in LAI of Yarkand Oasis.

Donohue13 Model
To understand the impact of CO 2 fertilization on vegetation greening, Donohue et al. [12] proposed a concept physical model, which relates vegetation change to atmospheric CO 2 concentration and other factors. The specific method is described as follows: The water use efficiency (WUE, g C/mm) at leaf scale, i.e., the ratio of CO 2 assimilation rate per unit leaf area (A l , g C/year) to transpiration rate per unit leaf area (T l , mm/year), is defined as: where C a (ppm) and C i (ppm) represent atmospheric and intercellular CO 2 concentrations, respectively, and v (kPa) refers to the vapor pressure deficit between the atmosphere and intercellular space. The relative effects of the variable changes in the equation on WUE are given by the following: Wong et al. [30] showed that at the leaf scale and under the condition of a given v, C i /C a was relatively constant for a specific mode of photosynthesis (C 3 or C 4 , whose photosynthesis dark reaction are different). However, (1 − C i /C a ) slowly increase with increasing v and was proportional to √ v, and this relationship was confirmed by both simulation and observation results [31][32][33]. Therefore, Equation (2) can be rewritten as: In arid regions, the LAI is usually very small, so the following reasonable assumption is made [34]: where T (mm/year) is the total regional transpiration rate and L (m 2 /m 2 ) is the LAI. Thus, the relative change is: Donohue et al. [12] assumed that the water supply is constant (i.e., dT T = 0), but reliability of the assumption was not tested. Therefore, dT T is retained and Equations (3) and (5) are combined to give: 2.1.2. Improved Donohue13 Model to Consider Changes in LULCC On the regional scale, the land cover can be simply classified into two types, namely vegetation cover land and non-vegetation cover land. The mean LAI of the whole region L region (m 2 /m 2 ) can be estimated as: where L is the mean LAI of the vegetation cover area and F veg is the ratio of the vegetation cover area to the total area. Thus, the relative change is: Combining Equations (6) and (8) leads to: According to the PETA hypothesis [35], the fractional response of leaf assimilation, A l , was proportional to that of water use efficiency, WUE, i.e., where α is resource availability index, i.e., where k is assumed as a coefficient related to species and its default value is 0.5. Generally speaking, when one or more key resources (such as water, light, and nutrient) are scarce, L region is usually low, whereas when resources are sufficient, L region is usually high [36]. Therefore, L region reflects actual region conditions, and α, which increases as the increasing of L region and varies between 0 and 1, can simply represent the actual resource availability. thus, Equation (9) can be simplified as: Finally, Equation (10) quantitatively expresses the relationship between the changes of F veg , T, C a , and v with the change of L region in arid area, making the changes in variables of different units and orders of magnitude comparable.

Study Area
The Yarkand Oasis (37 • 39 18"-40 • 23 31"N, 76 • 31 30"-80 • 23 38"E) is located in the extremely arid region of northwest China, with a total area of approximately 3.16 × 10 4 km 2 ( Figure 1) and is the largest oasis in Xinjiang, China [37], as well as the fourth largest irrigation area in China [38]. It is also a typical representative of the oasis-desert ecosystem (ODES) [39]. It is in the deep Eurasian continent and is under the temperate continental arid climate, with mean annual temperature of 12.0-12.6 • C and relative humidity of 46.6-54.9% from 1988 to 2017. The oasis region has low precipitation, with mean annual precipitation of 62.5-77.6 mm/year, and very high evaporative capacity, where the mean annual evaporation from a 20-cm evaporation pan is 1488-2499 mm/year from 1988 to 2017. The water availability of the oasis mainly comes from the Yarkand and the Tiznafu rivers, which originate from the upper mountainous regions. Both the Yarkand and the Tiznafu rivers are typical rivers fed by glacier and snow melting [40], with mean annual runoff of 69.
where k is assumed as a coefficient related to species and its default value is 0.5. Generally speaking, when one or more key resources (such as water, light, and nutrient) are scarce, Lregion is usually low, whereas when resources are sufficient, Lregion is usually high [36]. Therefore, Lregion reflects actual region conditions, and α, which increases as the increasing of Lregion and varies between 0 and 1, can simply represent the actual resource availability. thus, Equation (9) can be simplified as: Finally, Equation (10) quantitatively expresses the relationship between the changes of Fveg, T, Ca, and v with the change of Lregion in arid area, making the changes in variables of different units and orders of magnitude comparable.

Study Area
The Yarkand Oasis (37°39′18″-40°23′31″N, 76°31′30″-80°23′38″E) is located in the extremely arid region of northwest China, with a total area of approximately 3.16 × 10 4 km 2 ( Figure 1) and is the largest oasis in Xinjiang, China [37], as well as the fourth largest irrigation area in China [38]. It is also a typical representative of the oasis-desert ecosystem (ODES) [39]. It is in the deep Eurasian continent and is under the temperate continental arid climate, with mean annual temperature of 12.0-12.6 °C and relative humidity of 46.6-54.9% from 1988 to 2017. The oasis region has low precipitation, with mean annual precipitation of 62.5-77.6 mm/year, and very high evaporative capacity, where the mean annual evaporation from a 20-cm evaporation pan is 1488-2499 mm/year from 1988 to 2017. The water availability of the oasis mainly comes from the Yarkand and the Tiznafu rivers, which originate from the upper mountainous regions. Both the Yarkand and the Tiznafu rivers are typical rivers fed by glacier and snow melting [40], with mean annual runoff of 69.1 × 10 8 m 3 and 13.3 × 10 8 m 3 , respectively.

Data
The data used in this study included LAI, land cover, meteorological data (including precipitation, temperature, and relative humidity), runoff, and CO 2 concentration, ranging from 2001 to 2017 (Table 1). The moderate resolution imaging spectroradiometer (MODIS) LAI Collection 6 product MOD12A2H.006 [41] from 2001 to 2017 are used to quantitatively analyze the changes in vegetation LAI in the study area. Its spatial resolution is 500 m and temporal resolution is 8 days. The data of the MODIS land cover product MCD12Q2.006_LC_Type1 [42] from 2001 to 2017 are used to quantitatively analyze the land cover changes in the study area, which are classified according to the classification criteria of the International Geosphere-Biosphere Program (IGBP), with a temporal resolution of 1 year and a spatial resolution of 500 m. The precipitation, temperature, relative humidity, and wind speed data are from the China Meteorological Data Network (https://data.cma.cn), and they are the daily data from 2001 to 2017. The study area has a total of five meteorological stations, their locations were shown in Figure 1. The oasis is recharged by the Yarkand River and the Tiznafu River. The inflow stations are the Kaqun and Jiangka stations, and the downstream outflow is controlled by the Heiniyazi station (see Figure 1).

Data Processing
Regarding land cover, the ratio of the number of pixels representing vegetation types (i.e., open shrublands, grassland, croplands, and cropland/natural vegetation mosaics) to the total number of pixels in the region is denoted as F veg . It is assumed that there is no seasonal change in land cover. Table 2 lists the average ratio of various land cover types in the region from 2001 to 2017. The regional average LAI (L region ) could be estimated as the accumulated LAI value of vegetation pixels dividing by the total pixels number covering the study area. Then, the growing season was defined as the period from April to October [43], and the growing season mean LAI was calculated using the arithmetic average method according to MOD12A2H.006.
The vapor pressure deficit (v) was calculated using the following Equation [44]: where AT ( • C) is the atmospheric temperature at 2 m and RH (1%) is the relative humidity at 2 m. The Thiessen polygon method [45], which is a simple method to reflect the spatial variability, was used to calculate the regional average daily value from the daily v of each station in the Yarkand Oasis, and the arithmetic average method was used to calculate the growing season mean v according to daily v. The total vegetation transpiration in vegetation area T (mm/year) could be calculated as: where W (mm/year) represents the total water consumption in vegetation area and R ((mm/year)/(mm/year)) is the ratio of transpiration to evapotranspiration. In the study, the inter-annual variation of surface water storage is ignored since there is no over-year regulation reservoir in the Yarkand Oasis. The inter-annual variation of soil moisture was also ignored. Therefore, W was calculated as follows: where I (mm/year) is the upstream surface inflow into the oasis, O (mm/year) is the downstream surface outflow from the oasis. I and O are obtained by dividing inflow and outflow runoff (m 3 /year), respectively, by vegetation area (m 2 ) every year with an assumption that runoff is only consumed in vegetation region. P (mm/year) is the precipitation and ∆G (mm/year) is the groundwater variation.
According to the research report [46], net ground flow only accounted for~0.24% of the surface flow in the Yarkand Oasis. Therefore, net ground flow can be ignored in this analysis. The monthly precipitation, surface inflow, and surface outflow from 2001 to 2017 are given in Figure A1 in the Appendix A. Thus, the last part of W is groundwater variation (∆G), which is calculated as follows: where µ (mm/mm) represents the specific yield of groundwater, and as the study area is mainly covered by sandy soil, 0.1 m/m is used [47]. ∆h (mm/year) represents the variation in groundwater depth. Chen et al. [48] estimated a -18 m/year change in groundwater level from 1993 to 2011 based on 55 gauges across the oasis. Therefore, the value of ∆h is set to be 0.18 m/year for this study. The ratio R in Equation (12) is calculated based on Wei et al. [49] and Wang et al. [50]. According to Wei et al. [49], R sg = 0.69L sg 0.28 (18) where R crops and R sg represent the T/ET of crops and shrubs/grasses, respectively, and then L crops and L sg represent the LAI of crops and shrubs/grasses, respectively. According to Wang et al. [50]: R ans = 0.91L ans 0.08 (21) where R as , R ns , and R ans represent the T/ET of the agricultural systems, natural system, and overall systems, respectively, and L as , L ns , and L ans represent the LAI of the agricultural system, natural system, and overall system, respectively. In the Yarkand Oasis, the vegetation types can be categorized into crops and shrubs/grasses, and they also can be categorized into agricultural system and natural system. Therefore, we took a simple way to estimate R during every growing season, i.e., calculating the five different ratios by using Equations (15)- (19) and replacing L crops , L sg , L as , L ns , and L ans by L during every growing season, and then taking the arithmetic average of the five ratios as R. In addition, maximum and minimum of the five ratios are used to calculate maximum and minimum of T, which constitute the error limit of T. The CO 2 concentration (C a ) in this study area is assumed to be represented by that measured at the Waliguan Observatory, since both the Waliguan Observatory and Yarkand Oasis are located in the arid region of the Eurasian continent, and the distance between them is only approximately 2000 km. We could not obtain C a of 2017, so it was obtained by interpolation using the C a of the Hawaii Observatory, since there was a high correlation of C a at the Waliguan Observatory and Hawaii Observatory, with the R 2 of the linear regression up to 0.989 from 2001 to 2016 ( Figure A2). Then, the growing season mean C a were calculated as the arithmetic mean of the daily CO 2 concentration during the growing season.

Calculation of Change Rate
To estimate the changes of variables, the Mann-Kendall method was used to test the trend of a given time series [51,52], and then the linear regression method was used to obtain the trend line of the time series [53]. Finally, the interannual change of variables was calculated as: where dy y represents the interannual change of a variable (%/year), y n , y 1 , and y represent the end value, initial value, and average value, respectively, of the trend line, and n represents the time series length. Finally, T , dC a C a , and dv v were calculated using Equation (20).  Table 3 shows the annual average rate of change in variables in the Yarkand Oasis. The simulated LAI change dL region L region by using Equation (10) was 1.20%/year. The relative contribution rate of each factor is shown in Figure 3. Among all the factors, change in vegetation cover ratio due to LULCC had the largest contribution (88.33%). The unchanged transpiration had no impact on vegetation change (0.0%), which indirectly reflects the impact of water availability per vegetation area. In addition, the contributions from enhancing CO 2 concentration and increasing vapor pressure deficit were 40.0% and -28.3%, respectively. They had approximate values in magnitude but opposite signs. From the perspective of vegetation WUE, the facilitation of elevated-CO 2 was greater than the inhibition of elevated-v, thus improving the WUE. T was estimated from water consumption per vegetation area by multiplying the ratio of vegetation transpiration to regional evaporation, which was calculated according to Equations (15)- (19) and the maximum and minimum of estimated T were shown as the error bars in subplot (d).  Figure 3. Among all the factors, change in vegetation cover ratio due to LULCC had the largest contribution (88.33%). The unchanged transpiration had no impact on vegetation change (0.0%), which indirectly reflects the impact of water availability per vegetation area. In addition, the contributions from enhancing CO2 concentration and increasing vapor pressure deficit were 40.0% and -28.3%, respectively. They had approximate values in magnitude but opposite signs. From the perspective of vegetation WUE, the facilitation of elevated-CO2 was greater than the inhibition of elevated-v, thus improving the WUE.   T was estimated from water consumption per vegetation area by multiplying the ratio of vegetation transpiration to regional evaporation, which was calculated according to Equations (15)- (19) and the maximum and minimum of estimated T were shown as the error bars in subplot (d).

Applicability of the Improved Donohue13 Model
This study detected a 2.31%/year greening trend in the growing-season of the Yarkand Oasis from 2001 to 2017, which was consistent with the previous studies [24,54,55]. However, the simulated greening trend was 1.21%/year and approximately accounted for half of the observed trend. The difference between the simulated and observed trend might be caused by the following two reasons: 1) Data error: According to the evaluation of Yan et al. [56], MODIS LAI Product Collection 6 was underestimated in grassland and cropland, which causes an overestimation of trend by Equation (22) due to underestimation of y and constant of 1 () n yy − in Equation (22).
Furthermore, the change in vapor pressure deficit (v) might be overestimated. On one hand, these meteorological stations are located in urban areas, resulting in a stronger upward trend than vegetation areas in v [57]. On the other hand, Bachu and Yecheng, which are closed to the desert and represent inadequate water supply condition, had weaker upward trends of

Applicability of the Improved Donohue13 Model
This study detected a 2.31%/year greening trend in the growing-season of the Yarkand Oasis from 2001 to 2017, which was consistent with the previous studies [24,54,55]. However, the simulated greening trend was 1.21%/year and approximately accounted for half of the observed trend. The difference between the simulated and observed trend might be caused by the following two reasons:

1)
Data error: According to the evaluation of Yan et al. [56], MODIS LAI Product Collection 6 was underestimated in grassland and cropland, which causes an overestimation of trend by Equation (22) due to underestimation of y and constant of (y n − y 1 ) in Equation (22). Furthermore, the change in vapor pressure deficit (v) might be overestimated. On one hand, these meteorological stations are located in urban areas, resulting in a stronger upward trend than vegetation areas in v [57]. On the other hand, Bachu and Yecheng, which are closed to the desert and represent inadequate water supply condition, had weaker upward trends of v than regional mean v. While Maigaiti, Shache, and Zepu, located inside the oasis with an abundant water supply condition, had stronger increase trends of v (Table 4). Actually, Yarkand Oasis is located on the edge of Taklimakan Desert, and vegetation might be under inadequate water supply condition. Therefore, change in v of Bachu and Yecheng might be closer to represent the actual change in v of Yarkand Oasis. In a word, the change in v might be overestimated, which amplifies the inhibition of v on vegetation growth. 2) Seed optimization: Cropland was one of the main land cover type in Yarkand Oasis, accounting for 20.4% of the total oasis area in 2001-2017 ( Table 2). As agricultural science advances, crop seeds were continuously improved. Optimized seeds may increase water use efficiency (WUE) by themselves, which causes a lower C i /C a and furthermore a greater d(1 − C i /C a )/(1 − C i /C a ) [58]. Removing the simplification that (1 − C i /C a ) was proportional to √ v, Equation (12) can be rewritten as: Therefore, a greater d(1 − C i /C a )/(1 − C i /C a ) due to seed optimization might cause a greater value of simulated dL region L region . In Equation (10), the transpiration T implies the effect of water availability in the vegetation area, which had no impact on vegetation change in the Yarkand Oasis (0.0%). To calculate the vegetation transpiration per vegetation area (T), five empirical formulas, i.e., Equations (15)- (19), were used in this study. As shown in Table 5, the significance levels of change in T based on different T estimated Equations (15)- (19) were greater than 0.2, which indicates unchanged T. Therefore, it means that the uncertainty due to different formula for T could be acceptable. Moreover, Equations (15)- (19) focus on the partition of evapotranspiration between soil evaporation and transpiration. Besides, the water consumption includes water surface evaporation. However, in the Yarkand Oasis, water surface covers 0.6% area, which possibly led to a relatively small error in the calculated T. In addition, the impacts of anthropogenic activities, such as the optimization of irrigation methods, were not considered in our study. The optimization irrigation methods might decrease soil evaporation (E) of cropland and increase estimated transpiration (T). The increase in the T of cultivated land might affect the change of total T of the oasis, which needs further studying in the future. Table 5. The significant levels of change in transpiration (T) based on different T estimated Equations (15)- (19).

Equations
Equation ( The results show that the facilitation of the increase in WUE to LAI was 11.7% (0.14%/year), as calculated as (1 − α) dC a C a − 1 2 (1 − α) dv v = 40.0% (0.48%/year) − 28.3% (0.34%/year). It was similar to the FACE experiments, which show that an elevated CO 2 concentration may increase the WUE of vegetation and further promotes vegetation [59][60][61][62]. In addition, coefficient of k in Equation (11) was assumed to relate to species and was assigned 0.5 by default. In order to analyze the impact of different species, the values of 0.4, 0.5, and 0.6 were set, which reflected the effectiveness of different species on resource availability. As shown in Table 6, there were small differences in the contribution analysis caused by different coefficients of k. It indicates that the results were not sensitive to species. Table 6. The contribution analysis based on a different coefficient of k in Equation (11) v * denote the contribution rate of CO 2 concentration, vapor pressure deficit, and the difference between the two, respectively. Furthermore, to evaluate the applicability of the improved model to different land cover type, the results in grassland and cropland, respectively, were presented with the assumption that transpiration per grassland or cropland area was constant. As shown in Table 7, there were small differences between observed and simulated changes in L grass , indicating a good performance in grassland. In contrast, as shown in Table 8, the result in cropland was relatively poor with a difference of 0.97%/year between observed and simulated L crop , which might be caused by data error and seed optimization. Table 7. The observed and simulated LAI in the grassland of Yarkand Oasis and its causes (unit: %/year). dL grass L grass is the change rate of LAI in grassland, together with the observed and simulated indicating the observed from MODIS and the simulated using the improved Donohue13 model, respectively. dT T , dC a C a , dv v , and dF grass F grass denote the impacts from changes in transpiration, CO 2 concentration, vapor pressure deficit, and ratio of grassland area, respectively.

Impact of Water Availability
The observed unchanged transpiration in vegetation area verifies the assumption that water consumption per vegetation area is constant in arid regions [12]. For a specific species, transpiration consumption for growth could be considered constant. In such arid regions, the water availability per vegetation area almost only makes up for the lack of vegetation transpiration for growth and accompanying soil evaporation during the growing season, and does not have too much excess during the growing season or water storage in the non-growing season. Therefore, water determines vegetation, but vegetation does not determine how much water to consume.
Although the water availability per vegetation area was constant, the total water availability was actually increasing. As shown in results, the vegetation area ratio (F veg ) was increasing, which was probably driven by increased total water availability. In other words, more F veg could be fed by more total water availability with constant transpiration consumption per vegetation area. Furthermore, the water availability mainly originated from local precipitation and inflow from outside. In the region that precipitation dominates the water availability, precipitation has a determining influence on vegetation [63][64][65]. In the region that most water availability comes from outside, such as the Yarkand Oasis, the water is mainly supplied by the streamflow from glacier and snow melting and plays a dominant role in vegetation change, and this pattern of water supply and consumption is widespread in Central Asia [66,67]. Due to global warming, the water from glaciers and snowmelt in Central Asia has increased significantly [68][69][70][71], and it increases the total water availability in the region [72,73]. Thus it plays an important role in the regional vegetation greening [74]. However, a decrease in total water availability may be expected by the end of the twenty-first century as a result of depleted glaciers without refreezing [75,76], which means that the vegetation greening due to F veg driven by increased total water availability is unsustainable.

Impact of LULCC Change
Using the improved Donohue13 model, we found that the change of LULCC contributes the most (88.3%) to the LAI change in the Yarkand Oasis, which is consistent with the results of Wang et al. [17] in the arid region of the Heihe River basin. Similarly, a number of studies on vegetation change in the Sahel region showed that LULCC plays an important role [28,29,77]. In recent years, some large-scale studies on the attribution analysis of vegetation change have begun to consider the impact of LULCC after recognizing its possible important role [9,10]. However, due to factors such as the uncertainty of the model and coarse spatial resolution, those studies may have large uncertainty in assessing the role of LULCC. For example, Zhu et al. [10] estimated the contribution of LULCC to global greening (4%) but the results from different ecological models had significant differences in magnitude and even in signs ranging from -11% to 18.4%. Furthermore, in their study, the spatial resolution was coarse (0.5 • ), which made it difficult to fully reflect the change in LULCC. Therefore, it is possible that they underestimated the contribution of LULCC to global greening. The improved Donohue13 model provides a simple but feasible tool to re-evaluate the impact role of LULCC on global greening based on finer land use map.

Implication of Improved Donohue13 Model
It is a hot topic to study the vegetation response to changing environments, including climate change, elevated atmospheric CO 2 concentration, and anthropogenic activities. In previous studies, the commonly used methods are statistical analysis [54,78,79] and process-based model simulation [9,10,21,80]. Generally speaking, the statistical analysis is difficult to quantify; while the process-based model has uncertainties in input data, structure, and parameters, and at the same time, it cannot provide an intuitive way to reveal the relationship of vegetation change with environmental factors. The improved Donohue13 model is physical-based and its variable relationship is explicit, which may be used to explore the relationship of vegetation change with environmental factors.
The interaction of vegetation and hydrology forms a coupled vegetation-hydrological system, changed by climate change, elevated CO 2 concentration, and anthropogenic activities, which causes changes in regional hydrological processes and water resource [5,81]. In previous studies, the watershed coupled water-energy balance equation (i.e., Budyko framework) relates evapotranspiration (E) with precipitation (P) and potential evapotranspiration (E 0 ) [82], and the equation includes a unique empirical parameter (n) related to vegetation [83]. Under the Budyko framework, many studies have analyzed the response of runoff to changing environment [84,85], yet based on the assumption that n is constant. In fact, the parameter n is changing with changing environment, and one of its main drivers is vegetation for a particular catchment [86], with a relationship of n = f L region , · · · between n and LAI [87]. Moreover, the improved Donohue13 model can explicitly express the relationship between vegetation change dL region L region and environmental variables on the catchment scale and provide a method to estimate the change in n. Therefore, it is a possible way for investigating the response of coupled vegetation-hydrology systems under environment change to couple the improved Donohue13 model and the Budyko framework.

Conclusions
Vegetation shows a greening trend on the global scale in the past decades, which has an important effect on hydrological cycle. At the same time, many studies have shown that LULCC plays an important role in vegetation change. To quantitatively interpret causes for vegetation greening, this study improved the Donohue13 model by taking into account the change in ratio of vegetation coverage area to total regional area due to LULCC, furthermore applies this model to the Yarkand Oasis in the arid region of northwest China. The results show that the regional mean LAI during the growing season significantly increased at a rate of 2.31%/year from 2001 to 2017. Due to data error and seed optimization, there is not good performance of the improved Donohue13 model with a 1.20%/yr upward trend. On the vegetation greening in the Yarkand Oasis, the contributions of: (1) LULCC; (2) atmospheric CO 2 concentration; and (3) vapor pressure deficit were: (1) 88.3%; (2) 40.0%; and (3) -28.3%. It means that the greatest contribution was from LULCC, which is probably driven by increased total water availability in the whole oasis with constant transpiration in the vegetation area. The improved Donohue13 model, a simple but physics-based model, can help explain the impact of climate change and anthropogenic activity variables on vegetation change in arid regions. It can be further combined with the Budyko hypothesis to establish a framework for quantifying the changes in coupled response of vegetation and hydrological processes to environment changes.

Acknowledgments:
We would like to thank the three anonymous reviewers for their constructive comments and language editing.

Conflicts of Interest:
The authors declare no conflict of interest. Figure A1 shows that the monthly precipitation (P), surface inflow (I), and surface outflow (O) in vegetation area of Yarkand Oasis from 2001 to 2017. The monthly P can be estimated as the sum of the daily P, which is calculated from daily P of meteorological stations in the Yarkand Oasis using the Thiessen polygon method, during each month. The monthly I and O are obtained by the monthly inflow at the Kaqun and Jiangka stations and the monthly outflow at the Heiniyazi station, respectively.