Investigating and predicting spatiotemporal variations in vegetation cover in transitional climate zone: a case study of Gansu (China)

Vegetation ecosystems are sensitive to large-scale climate variability in climate transition zones. As a representative transitional climate zone in Northwest China, Gansu is characterized by a sharp climate and vegetation gradient. In this study, the spatiotemporal variations of vegetation over Gansu are characterized using the satellite-based normalized difference vegetation index (NDVI) observations during 2000–2020. Results demonstrate that a significant greening trend in vegetation over Gansu is positively linked with large-scale climate factors through modulating the water and energy dynamics. As a climate transition zone, the northern water-limited and southern energy-limited regions of Gansu are affected by water and energy dynamics, differently. In the water-limited region, a weakening Asian monsoon along with colder Central Pacific (CP) and warmer North Pacific (NP) Oceans enhances prevailing westerlies which bring more atmospheric moisture. The enhanced atmospheric moisture and rising temperature promote the local vegetation growth. In contrast, large-scale climate variations suppress the southwest monsoon moisture fluxes and reduce precipitation in southern energy-limited regions. In these energy-limited regions, temperature has more effects on vegetation growth than precipitation. Therefore, the greenness of vegetation is because of more available energy from higher temperatures despite overall drying conditions in the region. Based on the above mechanism, future scenarios for climate impacts on vegetation cover over Gansu region are developed based on the two latest generation from coupled climate models (Coupled Model Intercomparison Project Phase 5 and Phase 6; CMIP5 and CMIP6). In the near-term future (2021–2039), the vegetation is likely to increase due to rising temperature. However, the vegetation is expected to decrease in a long-term future (2080–2099) when the energy-limited regions become water-limited due to increasing regional temperatures and lowering atmospheric moisture flux. This study reveals an increasing desertification risk over Gansu. Similar investigations will be valuable in climate transition regions worldwide to explore how large-scale climate variability affects local ecological services under different future climate scenarios.


Introduction
Vegetation is a natural interface between soil, hydrology, ecosystem, and climate, and it is a sensitive indicator of regional environmental change (Cui and Shi 2010). Vegetation variability in different parts of the world varies greatly over the past decades. The vegetation has been increasing in the north of extratropical latitudes (Mao et al. 2016), and South Asia (Wang et al. 2017b), but there are opposite trends in boreal Eurasia (Piao et al. 2011) and Inner Asia (Mohammat et al. 1 3 2013). Shifts in vegetation are mainly attributed to global and regional climate changes (Cui and Shi 2010;Li et al. 2015;Xu et al. 2016), land-use changes (Dirnböck et al. 2003;Fernandes et al. 2011;Tasser and Tappeiner 2002), and the carbon dioxide fertilization (Los 2013;Schimel et al. 2000;Yang et al. 2016). Among these factors, climate variability has been recognized as the most direct and important driver for vegetation variations (Cui and Shi 2010;Yang et al. 2019).
Shifting climate patterns play an important role in vegetation spatiotemporal variability, especially in climate transition zones (Xia et al. 2019), including Qinling Mountains in China (Xia et al. 2019) and central Queensland in Australia (Krull et al. 2005). In such climate transition zones, ecosystems are unstable and highly sensitive to regional climate fluctuations (Hou et al. 2019). As a transitional climate zone between humid and arid regions in Northwest China, Gansu is characterized by a sharp climate and vegetation gradient . Due to global warming, the dry Northwest China has been transformed into a wet region since the last century (Dai et al. 2011), while Northeast China became drier (He et al. 2020). The combinations of wetting and drying trends in the different parts of Gansu make the entire ecosystem more vulnerable. The wetting and drying patterns might be due to shifting regional precipitation and temperature patterns over the region (Dai et al. 2011;Wang et al. 2017a). The regional precipitation and temperature distributions are controlled by large-scale climate variability by modulating regional the water and energy cycles (Ouyang et al. 2014;Xiao et al. 2015). Therefore, detecting changes in vegetation dynamics and identifying their linkages with regional and large-scale climate variability is crucial to regional ecological health assessments and regional economic development coordination under changing climate scenarios.
Over Asian regions, vegetation covers are closely related to the Asian monsoons, i.e., the East Asian monsoon (EAM; Jiang et al. 2006;Zhao and Yu 2012) and the Indian monsoon (IM; Chen et al. 2014;Lee et al. 2009). Vegetation dynamics are also related to sea-surface temperature (SST) modes in the Pacific (Erasmi et al. 2009;Jiang et al. 2011;Lü et al. 2012) and Indian Oceans (Li et al. 2017). For example, climatological anomalies in vegetation cover over Indonesia were associated with increases in extreme events (especially droughts) in response to El Nino Southern Oscillation (ENSO; Erasmi et al. 2009). Similarly, ENSO was demonstrated to affect the vegetation cover in China (Jiang et al. 2011;Lü et al. 2012).
All previous studies suggested that Asian monsoons and SST oscillations affect regional or local vegetation growths through the modulation of local climate variables. However, the relationships between vegetation and local climate variables including precipitation, evaporation, and temperature are not consistent in previous studies (Li et al. 2009;Xu et al. 2016;Zhao et al. 2011). When Li et al. (2009) found that precipitation and temperature are both significantly related to vegetation, Zhao et al. (2011) suggested that the role of temperature is insignificant. Moreover, Xu et al. (2016) found that grassland and cropland have different responses to precipitation, evaporation, and temperature. In this study, we establish hydroclimate mechanisms over Gansu based on both water and energy dynamics to link the large-scale climate variability with local vegetation. Based on this mechanism, we will produce the future projections of vegetation based on climate model outputs for 2021-2039 and 2080-2099. Overall, this study aims to (i) investigate the spatiotemporal changes in vegetation over Gansu; (ii) establish water and energy mechanisms between climate drivers and vegetation variation; and (iii) develop future scenarios for spatiotemporal changes in vegetation based on climate model outputs.
The paper is structured as follows. In Sects. 2 and 3, we introduce the materials and methods. In Sect. 4, the spatiotemporal vegetation variations and its linkages with local and large-scale climate variability over Gansu are investigated. In the final section, the implications and possible future applications related to desertification risks are summarized and discussed for climate transition zones.

Study area
In the inner land of Northwest China, Gansu has plateau terrain inclined from South-west to North-east, with an elevation from 598 to 5602 m above the sea level (a.s.l.) . It lies between 93°-110° E and 32°-44° N, with a total Fig. 1 The elevation with 16 locations (blue dots) over Gansu. The magenta (elevation contour at 3000 m a.s.l.) and purple lines (annual precipitation contour at 400 mm) divide the Gansu into three graphically regions: the Hexi Corridor arid region (HCAR), the Qinghai-Tibet alpine region (QTAR), the Loess Plateau semi-arid region (LPSR) area of 455,000 km 2 (Fig. 1). Based on Zhao (1983), Gansu can be divided into three natural geographical regions by 3000 m a.s.l. elevation contour and 400 mm contour of annual precipitation ( Fig. 1): the Hexi Corridor arid region (HCAR), the Qinghai-Tibet alpine region (QTAR), and the Loess Plateau semi-arid region (LPSR). The Gansu region straddles the alpine, semi-arid, and arid climatic zones, which make it vulnerable to climate change (Li et al. 2013). The annual average temperature is around 0-14 °C, and it varies greatly from the cold QTAR (western part) to the warm HCAR and LPSR (eastern part) (Wang et al. 2014b). Annual average precipitation varies from 50 in the northwest to 500 mm year −1 in the southeast region (Cheng and Falkenheim 2016). Over the Gansu region, precipitation concentrates between June and September during the Asia summer monsoon season, and it shows strong interannual variations (Li et al. 2013). The main land and vegetation types are deserts, grassland, and forest (Wang et al. 2014b). It has been suggested that intensive climate variations have placed a heavy pressure on the local fragile ecological and hydroclimate systems and they have impeded the sustainable agricultural and economic development in the region (Han et al. 2015;Wang et al. 2003).

Data
Vegetation indicators are extracted from satellite products, and meteorological variables and climate indices are derived from reanalysis datasets. The details of the datasets are summarized in Table 1.

Normalized difference vegetation index
Vegetation cover is widely and continuously monitored by satellite remote sensing (Yang et al. 2019). The normalized difference vegetation index (NDVI) is one of the most widely used remote sensing measures, and many NDVI products have continuous records over decades (e.g., Li et al. 2010b;Mkhabela et al. 2011;Zhang et al. 2003). Here, the NDVI dataset is derived from the Terra Moderate Resolution Imaging Spectroradiometer (MODIS) product MOD13C2 in a spatial resolution of 0.05 × 0.05° at a monthly timestep (Didan 2015). The NDVI dataset is extracted from the National Aeronautics and Space Administration (NASA) Land Processes Distributed Active Archive Center (LP DAAC; https:// e4ftl 01. cr. usgs. gov/ MOLT/ MOD13 C2. 006/). The MODIS NDVI has been widely applied in large-scale vegetation studies (e.g., Badreldin et al. 2014;Li et al. 2015;Xu et al. 2016). In this study, the NDVI is used to explore the spatiotemporal variability in vegetation cover between February 2000 and January 2020 (total of 20 years), a period that has not been largely explored in previous studies.

Moisture budgets
The vegetation growths are mainly controlled by water and energy balances. To investigate regional changes in vegetation and their link to water and energy variability, we examine moisture dynamics that associate with water and energy cycles (Peng and Zhou 2017). The atmospheric moisture conservation equation in flux form of vertical integration is written as follows (Trenberth et al. 2011): where W is the total column water vapor, AET is the actual evapotranspiration, P is the precipitation, and ∇ • Q represents the vertically integrated atmospheric moisture flux divergence (hereafter called MFD). The tendency term W t is small for long-term means. Equation (1) can be written as follows: The primary balance of moisture is thus between P + ∇ • Q (gaining moisture) and AET (losing moisture). For a region, water mainly comes from precipitation brought by horizontal moisture movements and vertically convective activities. To represent the horizontal water movement, vertically integrated moisture fluxes and MFD are used. AET is a key factor in the water cycle, but also an important part in the energy cycle in the form of latent heat (Trenberth et al. 2011). To study thermally vertical motion, the convective available potential energy (CAPE), a measure of energy available for lifting air parcels from the lower to the upper atmosphere, is used. Therefore, both horizontal and vertical  Huang et al. (2017) dynamics are used to examine the effects of water and energy on vegetation. From the above framework, AET acts as a bridge or fluxes between the water and energy cycles. The precipitation, AET, and temperature data are obtained from the ERA5-Land monthly averaged datasets between 1981 and 2020 (https:// cds. clima te. coper nicus. eu/ cdsapp# !/ datas et/ reana lysis-era5-land-month ly-means? tab= form) (Muñoz 2019). By including improved land surface processes, the ERA5-Land reanalysis datasets provide higher spatial resolution data (0.1° × 0.1°) than its driven climate reanalysis data (0.25° × 0.25°) (Muñoz 2019). Many studies recommended the use of Tropical Rainfall Measuring Mission (TRMM) data to estimate precipitation over China (e.g., Cao et al. 2018;Ferreira et al. 2013;He et al. 2020). In Fig. 14 of the Appendix, ERA5-Land and TRMM precipitation data are very comparable, with high-correlation levels and high significances. It suggests the reliability of ERA5-Land datasets over Gansu. MFD and CAPE, which are not available in ERA5-Land, are extracted from ERA5 data over the same period, but with a spatial resolution of 0.25° × 0.25° (https:// cds. clima te. coper nicus. eu/ cdsapp# !/ datas et/ reana lysis-era5-press ure-levels-month ly-means? tab= form) (Hersbach et al. 2019).

Water-limited and energy-limited environments
Vegetation growths are affected by both water and energy factors. In this study, to determine the dominating factors, we employ the concept of water-limited and energy-limited environments (Parsons and Abrahams 1994). Based on the Parsons and Abrahams (1994), a water-limited (energy-limited) environment is defined as the areas having a P PET −1 ( potential evapotranspiration) ratio lower (greater) than 0.75. The PET data is also derived from the ERA5-Land dataset. The majority of the QTAR area is energy-limited because the P and PET ratio is higher than 0.75, while the HCAR is water-limited with the P and PET ratio lower than 0.75, and the LPSR is a mixed energy-and water-limited region (Fig. 2). Therefore, vegetation growth in the HCAR is limited by the availability of precipitation but not temperature (Javadian et al. 2020;Parsons and Abrahams 1994), whereas in the energy-limited QTAR and some parts of LPSR, the growth of vegetation is mainly restricted by the temperature (Gokmen et al. 2013; Parsons and Abrahams 1994).
Among different kinds of Asian monsoons, and the SST indices in the Pacific and Indian Oceans, the Webster and Yang Monsoon (WYM), the Central Pacific El Nino oscillation (CP), and the North Pacific SST anomalies (NP) are found to be the more significant contributors to NDVI-based vegetation variability to Gansu (Figs. 15-16 of the Appendix). The WYM is computed using the difference between zonal winds at 850 and 200 hPa over the Indian region (0-20° N, 40°-110° E; Webster and Yang 1992). The SST indices for the CP and the NP are respectively calculated according to the definitions provided in Kao ad Yu (2009) and Mantua et al. (1997). The CP pattern is different from the eastern type of ENSO (EP; the traditionally defined ENSO type), and these Pacific SST patterns affect precipitation over China differently (Lv et al. 2019). The western North Pacific subtropical high (WNPSH) was suggested to be more strongly related to the CP than to the EP (Weng et al. 2011). The WNPSH Fig. 2 The distribution of the ratio of precipitation and PET. Warm color denotes the water-limited regions (i.e., a ratio less than 0.75), while cool color indicates the energy-limited regions (i.e., ratio larger than 0.75). The magenta and purple lines divide the Gansu into three graphically regions (c.f. Figure 1) plays significant role in regulating the hydroclimate system in China (Gao et al. 2020). Therefore, the CP is expected to be responsible for vegetation changes through its effects on local climate systems over China. In addition, the NP has been shown to contribute to precipitation variability over China (Ao and Sun 2016;Li and Li 2000;Zhou and Xia 2012), and it thus might affect vegetation growth in the transition regions of China.

Climate change scenarios
The Coupled Model Intercomparison Project Phase 5 (CMIP5) and Phase 6 (CMIP6) model outputs have been widely used for evaluating future conditions for vegetation (Zhao et al. 2020;Zhou et al. 2020). Compared to previous phases, CMIP5 models include more carbon processes and feedback mechanisms of climate systems, while CMIP6 have finer resolution with improved dynamical processes (Eyring et al. 2016;Taylor et al. 2012). In this study, the SST and wind fields from CMIP5 (Table2) and CMIP6 models (Table 3) are used to derive the atmosphere-ocean oscillation indices. Regarding CMIP5 models, we used historical scenarios between 1850 and 2005 and three future scenarios (RCP 2.6, RCP 4.5, and RCP 8.5) between 2006 and 2100. The RCP 2.6 scenario corresponds to a strongly declining emission scenario, leading to warming of well below 2 °C, which is compatible with the Paris Agreement (Wang et al. 2021). The RCP 4.5 scenario corresponds to an approximate doubling (medium emission scenario) in carbon dioxide relative to the pre-industrial level, whereas the RCP 8.5 scenario represents a more than the threefold increase (high emission scenario; Swain and Hayhoe 2015). It is worth noting that the SST and wind field data under RCP 2.6 scenario are unavailable for most CMIP5 models in Table 2. Therefore, only four models, i.e., bcc-csm1-1-m, IPSL-CM6A-MR, MPI-ESM-LR, and NorESM1-M, are used for RCP 2.6 scenario. For CMIP6, similarly, atmosphere-ocean oscillation indices were estimated based on the historical scenarios (1850-2014) and three Shared Socioeconomic Pathways (SSP), including SSP1-2.6, SSP2-4.5, and SSP5-8.5. For four CMIP6 models, total 51 simulations were used: 10 members for ACCESS, 25 members for CanESM5, 6 members for IPSL, and 10 members for MIROC. In this study, two future

Mann-Kendall test
The Mann-Kendall (MK) test is a nonparametric method to quantify the significance of linear temporal trends (Kendall 1975;Mann 1945). Previous work argued that the results of the MK trend test could be misleading if serial correlations and outliers are ignored (Hamed 2008;Hamed and Ramachandra Rao 1998;Khaliq et al. 2009). In this study, the MK test is modified, according to Hamed and Ramachandra Rao (1998), to examine the trend significance of NDVI and climate variability. Trend intensity is estimated based on Thiel-Sen's slope, which is robust to outliers (Sen 1968).

Generalized least square regression
Vegetation covers over Gansu are hypothesized to be related to monsoons and SST anomalies in Pacific Oceans. The generalized least square (GLS) regression models for NDVI values with the adjustments of serial correlations are expressed as follows: where NDVI m is the modeled NDVI, and WYM , CP , and NP are the regression coefficients for their corresponding monsoon and SST indices for the Central and North Pacific Oceans. The GLS model is only based on large-scale processes, as CMIP5 models are more likely to perform better to reproduce these processes than local precipitation and temperature (Wang et al. 2014a). In this study, for future vegetation projection, near term (2030s; 2021-2039) and long term (2090s; 2080-2099) will be used.

Bias correction
Before developing the empirical regression models, bias-corrections have been applied to winds and SST model outputs, to reduce differences between climate model outputs and reference data sets (ERA5 and ERSSTv.5). Based on cumulative distribution function (CDF), a quantile mapping method by Panofsky and Brier (1958) is used to reduce biases due to (3) NDVI m = 0 + WYM WYM + CP CP + NP NP scale gaps between the numerical model grid and the scale of investigated processes. This quantile mapping method has been widely applied in hydrological impact studies (Boé et al. 2007;Li et al. 2010a;Shukla et al. 2019) and regional climate change investigations (Fowler et al. 2007;Grillakis et al. 2013;Miao et al. 2016). For a variable x , the method can be expressed as follows: where F −1 r is the inverse CDF of the reference data set, i.e., ERA5 and ERSSTv.5, and F m is the CDF of modeled climate indices from CMIP5 and CMIP6. x m are the modeled variables and x BC are the bias-corrected outputs.

Evaluation of GLS model skills and bias-correction performances
To evaluate the performance of the GLS models for the NDVI estimation and the bias-correction procedures, the leave-one-out cross-validation method (LOOCV) is used. For each validation, n samples are randomly divided into a training set with n-1 samples and a test set with one sample. All the cross-validations are done based on 200-simulation ensembles. The cross-validation error (CVE) is defined as follows: where f (•) is the GLS model and bias-correction procedure developed from the training sets using Eq. (3) or (4). X(test) and Y(test) are the corresponding test sets. The value is closer to 0, and the performances of the GLS model and the bias-correction procedures are better. We then further examine the GLS model performance in simulating the historical NDVI variability using the mean square error (MSE) and the ratio of standard deviation (RSD). The MSE, one of the most common estimates of errors, is written as follows: where x r is the reference climate index of x i .
The RSD is defined as follows: where m and r are the standard deviations of modeled and referenced datasets. (4)

Vegetation patterns
Over the Gansu region, the mean annual NDVI shows a North-South gradient, with more vegetation in the North than in the South (Fig. 3a). In different geographical areas of Gansu, the averaged NDVI value of the LPSR is the highest (0.2 ≤ NDVI ≤ 0.7), whereas the averaged NDVI value of QTAR (0 ≤ NDVI ≤ 0.5) is higher than that of the HCAR (0 ≤ NDVI ≤ 0.2; Fig. 3a). All three regions show significant increasing trends in the NDVI (Fig. 3b). The increasing trend in the NDVI is however more pronounced in regions with greater mean vegetation cover (Fig. 3a-b): the LPSR area has the highest increasing NDVI trend, while HCAR has the lowest.

Local water and energy patterns and vegetation covers
The spatiotemporal patterns of vegetation growth are related to the variations of local water-energy dynamics. The local water-energy dynamics are usually represented by variations in local precipitation, AET, and temperature. In Fig. 4a, the average annual precipitation amount shows a similar spatial pattern to the averaged NDVI (Fig. 3a): The LPSR and the QTAR receive more than 480 mm year −1 , but most of the HCAR receive less than 200 mm year −1 . It suggests that vegetation is denser where precipitation is relatively higher. It is also noted that a decreasing trend in precipitation over the southern regions and an increasing trend in northern regions are in between 2000 and 2020 ( Fig. 4b). Regarding AET, mean spatial patterns and trends are consistent with precipitation ( Fig. 4a-d): (i) AET is very low in the HCAR, as there is little water for local evaporation; (ii) higher AET is found in the LPSR and the QTAR. Consistently with precipitation trend patterns, AET is thus decreasing in southern regions, but increasing in northern areas (Fig. 4d). Therefore, in terms of water dynamics, AET appears to balance the long-term changes in precipitations.
In the southern (northern) region with less (more) precipitation, the AET is less (more). Over the QTAR, the average temperature is below 0 °C, while the LPSR and the HCAR have a temperature ranging from 5 to 15 °C (Fig. 4e). As shown in previous global study results (Turkington et al. 2019), most of the Gansu region experiences warming (Fig. 4f). Figure 5 shows the regression maps of how the vegetation is related to local water and energy variations. The precipitation, AET, and temperature show significant positive relationships with vegetation covers (Fig. 5). Precipitation and temperature provide the water and energy to sustain the vegetation growth. More AET suggests more latent heat flux and water vapor in the atmosphere, helping the formation of precipitation (Yang et al. 2018). Therefore, AET favors vegetation growth by both providing more energy and promoting the local water cycle.

Large-scale atmosphere-ocean variability and vegetation covers
Large-scale atmosphere-ocean variability modulates regional energy and water circulations which affect local vegetation variations. In Sect. 4.2, the precipitation, AET, and temperature show significant positive relationships with vegetation covers. Figure 6 shows the impacts of monsoon Fig. 3 The mean NDVI (a) and NDVI trend (b) between February 2000 and 2020 January. Black dots indicate that the trends are statistically significant at p < 0.1 according to the MK trend test. The magenta and purple lines divide the Gansu into three graphically regions (c.f. Figure 1) (i.e., WYM) and Pacific SST variability (i.e., CP and NP) on regional water, energy, and vegetation. The WYM and CP indices show significantly positive relationships with precipitation, AET, and temperature ( Fig. 6a-b, d-e, g-h). According to the positive relationships between vegetation and precipitation, AET, and temperature (Fig. 5), it is  Figure 1) suggested that the WYM and CP could promote local vegetation growth by providing more water and energy (i.e., more precipitation and AET, and higher temperature). The positive contributions of WYM and CP to vegetation are also validated in Fig. 6j-k. As opposed to WYM and CP, NP mainly shows negative, but non-significant, relationships with precipitation, AET, and temperature over Gansu (Fig. 6c, f, i). Interestingly, the NP is significantly and positively related to vegetation (Fig. 6l), but it has non-significant relationships (even at p < 0.1) with regional water and energy variables over most regions (Fig. 6c, f, i). Therefore, the accumulated weak energy and water effects of NP appear to have significant impacts on vegetation growth. To investigate the mechanisms driving these positive relationships between vegetation variations and large-scale climate variability, horizontal (i.e., MFD) and vertical water and energy dynamics (i.e., CAPE) are examined in Figs. 7 and 8. Figure 7 shows how the monsoon and the Pacific SST oscillations are related to vertically integrated horizontal moisture movement. Over Gansu, climatological moisture patterns are controlled by prevailing westerly and Asian monsoons ( Fig. 7a; Ren et al. 2016). Prevailing westerly brings moisture from Euro-Atlantic regions to the north part of Gansu, while the south-westerly moisture fluxes from the Indian Ocean and the south-easterly fluxes from the Pacific bring the atmospheric moisture to South Gansu (Fig. 7a). Westerly winds and monsoon moisture fluxes generally converge in the Gansu midlands (Fig. 7a). Both WYM and CP are mainly positively related to south-westerly moisture fluxes to China, even though the CP pattern would be weaker (Fig. 7b-c). On the contrary, WYM and CP both suppress the westerly  Figure 1)  Fig. 6 The precipitation (a-c), AET (d-f), temperature (g-i), and NDVI (j-l) regressed map with WYM, CP, and NP between February 2000 and 2020 January. Black dots indicate significant values at the 0.1 significance level according to the t-test. The magenta and purple lines divide the Gansu into three graphically regions (c.f. Figure 1) moisture fluxes from the Euro-Atlantic regions (Fig. 7ac). Different from WYM and CP, NP is negatively related to south-westerly moisture fluxes but positively related to south-easterly fluxes (Fig. 7d). This suggests that a warm North Pacific would lead to less moisture from the Indian Ocean but more water vapor from the Pacific Ocean to China. Such opposite impacts of NP on the different moisture fluxes contributing to the Gansu water balance may explain its weak impacts on local water and energy variables (Fig. 6c, f, i). Different circulation effects may interact with each other, thus masking the NP impacts on local climate variables.
In Fig. 8b-d, WYM and CP have significantly positive relationships with CAPE, but NP suppresses CAPE over Gansu. Generally, CAPE is very small over Gansu (smaller than 200 J kg −1 ) in Fig. 8s. Low CAPE suggests that the atmosphere is not convective in the region. Moreover, there is no trend in CAPE during 2000 and 2020 (Fig. 17 of the Appendix).
Therefore, the weak CAPE in the region has a limited role in local thermodynamic effects on vegetation although largescale climate variability has impacts on the CAPE strength.

Impact scenarios for vegetation cover in Gansu
After investigating how WYM, CP, and NP are related to vegetation variations over the Gansu region, CMIP5 and CMIP6 outputs are used to explore how vegetation cover is likely to change over the twenty-first century over the Gansu region under different emission scenarios.

Bias-correction and cross-validation
WYM, CP, and NP indices are computed for CMIP5 and CMIP6 models, and are bias-corrected, before developing future scenarios for vegetation covers using the GLS regression models. Figure 9 shows the bias-corrected CDF results   Fig. 9 The comparison of the empirical CDF of referenced (i.e., ERA5 and ERSSTv.5 datasets, here called observation for simplicity), modeled, and BC-modeled WYM (a), CP (b), and NP (c) during the historical period from CMIP5. d-f is the same as a-c but for the CMIP6 of historical simulations from CMIP5 (Fig. 9a-c) and CMIP6 ( Fig. 9d-f) against the original CDF from the reference datasets. The WYM index derived from the CMIP5 outputs overestimates the minimums of WYM values (between − 20 and − 5) and underestimates the maxima of monsoon values (between − 5 and 20; Fig. 9a). The CP index derived from the CMIP5 model overestimates CP-Nina conditions and underestimates CP-Nino conditions (Fig. 9b). The modeled NP values from the CMIP5 outputs are all underestimated (Fig. 9c). The WYM, CP, and NP indices from CMIP6 outputs show similar results to CMIP5 (Fig. 9d-f). Specifically, the CP index from CMIP6 matches better with referenced data compared to that from CMIP5 (Fig. 9b, e). After the bias corrections, the CDF of the climate indices derived from CMIP5 and CMIP6 match well with their reference distributions (Fig. 9). Between simulated and referenced climate indices, the MSE of the WYM, the CP, and the NP are reduced by 41% and 40%, 95% (76%), and 96% (97%) for CMIP5 (CMIP6) outputs, respectively. For the cross-validation results of bias-correction procedures, the CVE values are 0.537 and 0.385 for WYM, 0.045 and 0.052 for CP, and 0.004 and 0.003 for NP from CMIP5 and CMIP6, respectively, and they are in the same magnitude as the MSE. Before projecting NDVI, the GLS model is also cross-validated (Fig. 10). The CVE values of the NDVI over Gansu are lower than 0.01, which is in the order magnitude of the MODIS NDVI accuracy, suggesting that the GLS regression model performs adequately for predicting NDVI variations over the study period (Fig. 10).

NDVI future scenarios
Using cross-validated GLS regression models, the biascorrected WYM, CP, and NP indices are used to project NDVI values for three RCP scenarios (RCP 2.6,RCP 4.5,and RCP 8.5) and three SSP pathways (SSP1-2.6, SSP2-4.5, and SSP5-8.5). For 16 locations (cf. Figure 1), the performances of simulated NDVI over the historical period (from February 2000 to January 2020) are evaluated for 13 CMIP5 (4 models for RCP 2.6) and four CMIP6 climate models using the MSE and the RSD. The MSE values are generally smaller than 0.01 and the RSD values are between 0.77 and 1.1 (close to 1) for all models from CMIP5 and CMIP6, all locations, and all scenarios (Fig. 11). These results indicate that the GLS regression models have an adequate accuracy for reconstructing NDVI variations based on large-scale climate indices. Figure 12 shows the median changes in the NDVI after 2020, based on 13 CMIP5 (4 models for RCP 2.6) and four CMIP6 models. In CMIP5 models, in the near-term future (2030s), NDVI values increase for all locations under RCP 2.6 ( Fig. 12a). NDVI values increase in the southern locations but decrease in the northern locations for both the RCP4.5 and RCP 8.5 scenarios (Fig. 12c, e). Moreover, under the RCP8.5 scenario, there are fewer locations where NDVI values increase, compared to the RCP4.5 scenario (Fig. 12c, e). It suggests that excess greenhouse gas (GHG) emissions may harm the vegetation growth even though the GHG like CO 2 is beneficial for the vegetation growth. Looking the long-term future (2090s), almost all locations of the Gansu region show a decrease in the NDVI for three scenarios (Fig. 12b, d, f). In CMIP6 models, in the near-term future (2030s), the NDVI values increase for most locations and under three SSP scenarios, compared to the study period (i.e., 2000Fig. 10g, i, k). It is also noted that the increase rate is larger over southern locations (Fig. 10g, i, k). In the long-term future (2090s), we also found an increase in the NDVI values for most locations under the SSP1-2.6 scenario (Fig. 12h). Meanwhile, NDVI values decrease for almost all locations under SSP2-4.5 (Fig. 12j). Interestingly, under SSP5-8.5 scenario, regional contrasts are found in future NDVI values (Fig. 12l): decreasing (increasing) in northern (southern) locations of Gansu, compared to 2030s.
Generally, both CMIP5 and CMIP6 show a tendency toward an increase in vegetation cover to increase in the near-term future (2030s). For the long-term future (2090s), CMIP5 and CMIP6 models show a decrease tendency in vegetation cover under most scenarios, but not in the SSP1-2.6 scenario. The projection changes suggest that current climate patterns will promote the vegetation growth over Gansu in the near-term future, but will eventually lead to the overall vegetation reduction in the long-term future. Moreover, the identified increasing vegetation under SSP1-2.6 scenario suggests that declining emissions could help to alleviate the vegetation reduction in the future. Fig. 10 The CVE for the estimated NDVI between February 2000 and January 2020 is based on the GLS models using the LOOCV method. The magenta and purple lines divide the Gansu into three graphically regions (c.f. Figure 1)  Fig. 11 The MSE (a) and RSD (b) between satellite-based and the modeled NDVI from the CMIP5 models under RCP 2.6. c-d and e-f are the same as a-b but for RCP 4.5 and RCP 8.5, respectively. The MSE (g) and RSD (h) between satellite-based and the modeled NDVI from the CMIP6 models under SSP1-2.6. i-j and k-l are the same as g-h but for SSP2-4.5 and SSP5-8.5, respectively

Discussion and conclusion
This study investigates spatiotemporal changes in vegetation cover over the transition zone of Gansu between 2000 and 2020, using a satellite-based NDVI product. Since 2000, Gansu has been increasingly greener, especially in the southern regions. Such vegetation recovery could be explained through the combined impacts of local water and energy dynamics associated with weakening WYM strength, a colder central Pacific Ocean, and a warmer North Pacific Ocean (Fig. 18 of the Appendix). The local water and energy budgets are controlled by horizontal and vertical dynamics. The vertical thermodynamic (i.e., unchanged CAPE values over the period; Fig. 17 of the Appendix) is found to have limited impacts on water and energy budgets over the Gansu region (Fig. 13). Therefore, the horizontal atmospheric dynamics play a major role in local vegetation variations of Gansu.
As a climate transition region, the water balance in Gansu is controlled by both Asian monsoons and prevailing westerly (Ren et al. 2016). The prevailing westerly brings moisture into North Gansu, which is mainly a water-limited region (Fig. 2). Asian monsoons carry moisture into South Gansu, which is mainly an energylimited region (Fig. 2). Then, the mechanisms of largescale climate variability through horizontal atmospheric dynamics on vegetation are separated into two types: energy-limited regions for South Gansu and water-limited regions for North Gansu (Fig. 13).
For energy-limited regions, the weakening WYM, cold phase of CP, and warm phase of NP (Fig. 18 of the Appendix) weaken the monsoon moisture fluxes to inland China, leading to lower precipitation than normal conditions (Fig. 13). Less precipitation means that less water would be evaporated, thus lower AET in the region (Fig. 13). Due to an increasing global temperature trend, the rising local temperature would continue to promote vegetation growths in the energy-limited region, despite locally drying conditions. In the waterlimited region, the weakening WYM, cold phase of CP, and warm phase of NP (Fig. 18 of the Appendix) enhance the prevailing westerly through weakening the Fig. 12 The median NDVI difference between the 2030s and study period (SP), and between 2090 and 2030s across all models from CMIP5 under RCP2.6 (a-b). c-d and e-f are the same as a-b but for RCP 4.5 and RCP 8.5, respectively. The median NDVI difference between the 2030s and study period (SP), and between 2090 and 2030s across all models from CMIP6 under SSP1-2.6 (g-h). i-j and k-l are the same as g-h but for SSP2-4.5 and SSP5-8.5, respectively. The green and red dots represent positive and negative changes, respectively. The magenta and purple lines divide the Gansu into three graphically regions (c.f. Figure 1) 1 3 The blue symbols are for positive or increasing relationships, and the red symbols are for negative or decreasing variations. The dotted box indicates no much changes for the variables (i.e., CAPE) southwest monsoon moisture fluxes (Fig. 13). Westerly and monsoon winds converge in the middle part of Gansu. When the monsoon becomes weaker, prevailing westerly becomes stronger and brings more moisture fluxes to North Gansu (Fig. 19 of the Appendix). More moisture fluxes over northwest China brought by prevailing westerly are consistent with previous studies (Peng and Zhou 2017;Ren et al. 2016). The warmer temperature and more precipitation promote AET. The increasing water and energy promote vegetation growth in the water-limited region.
In addition to climate effects, human activities also affect local vegetation variations. Figure 20 of the Appendix shows possible human effects based on the NDVI residuals, which are the differences between the observed and the GLS modeled NDVI. Both annual variations and residuals of NDVI between 2000 and 2020 can roughly be divided into three stages (Fig. 20 of the Appendix): Stage 1 is a steady increasing period between 2000 and 2013; Stage 2 is a plateau condition between 2014 and 2016, and stage 3 is another NDVI increasing period after 2016. The three-stage patterns of the residuals may indicate that the vegetation variations of the 2010s could partly be related to human activities. Gansu is one of the earliest pilot provinces which implement the Grain to Green Program (GTGP) (Li et al. 2015). In the last two decades, Gansu had two rounds of the program. The first round of the GTGP was between 1999 and 2013 and the second round started in 2014 (Gansu Forestry and Grassland Bureau, 2020; available in http:// lycy. gansu. gov. cn/ conte nts/ 79149. html). The increase in the NDVI values during 2000-2013 and after 2016 could thus be partly related to the first and second rounds of GTGP (Li et al. 2015). Between 2014 and 2016, the central government stopped the original GTGP subsidy which may be related to a break in the vegetation increasing trend in stage 2 of Fig. 20 of the Appendix.
For the future scenarios of vegetation cover, various studies have used the downscaling methods to improve the robustness (Maraun et al. 2010;Sunyer et al. 2015;Thomas et al. 2007). Downscaling procedures are indeed recommended in mountainous regions, especially where using local climate variables (i.e., precipitation and temperature) that which are strongly depending on the regional orography (Ding et al. 2018;Ma et al. 2018). However, in this study, the NDVI values are estimated based on large-scale climate teleconnections, which are unlikely to be affected by regional orography. Therefore, a downscaling procedure was not included in this study.
Based on the GLS model projections driven by largescale climate modes of variability derived from CMIP5 and CMIP6 models, the vegetation in Gansu, especially in the southern energy-limited regions, will keep increasing in the 2030s, as a response to climate variability and change. However, in the 2090s, Gansu will be more likely to experience a decline in vegetation cover based on most of the CMIP5 and CMIP6 projections. The continuous decreases in precipitation will thus lead to a transition from energy-limited regions toward more water-limited regions. Therefore, an increasing desertification risk should be considered for regional development planning and management, and more novel afforestation strategies based on changing monsoons and the Pacific SST conditions needed to be proposed. Overall, this study provides a framework to study possible increasing desertification risk, using climate scenarios from climate models, for the water-and energy-limited transition regions in the world.

Data availability
The NDVI data that support the findings of this study are available in the National Aeronautics and Space Administration (NASA) Land Processes Distributed Active Archive Center (LP DAAC; https:// e4ftl 01. cr. usgs. gov/ MOLT/ MOD13 C2. 006/). The climate data that support the findings of this study are openly available in ERA5-Land monthly averaged data from 1981 to present at http:// doi. org/ 10. 24381/ cds. 68d2b b30, ERA5 monthly averaged data on single levels from 1979 to present at http:// doi. org/ 10. 24381/ cds. f1705 0d7, and NOAA Extended Reconstructed Sea Surface Temperature (ERSST), Version 5 at http:// doi. org/ 10. 7289/ V5T72 FNM. The CMIP5 and CMIP6 model data are openly available in the Earth System Grid Federation (ESGF; https:// esgf-node. llnl. gov/).

Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.