Introduction

Satellite observations have shown that since the 1980s, global vegetation has exhibited an overall greening. Despite this, some areas show browning trends [5, 28, 32, 50, 56]. After the mid-1990s, due to the intensification of drought condition, the intensity and extent of vegetation browning expanded, resulting in a stalling or reversal of vegetation greening in different parts of the world [28, 31, 50]. Although the trends of vegetation growth in some regions have declined, there are still regions that continue greening as a result of land-use management. In recent years, to protect and restore degraded forest ecosystems, many countries have implemented a series of forest restoration projects, such as the Three North Shelterbelt Development Program and the Grain to Green Program in China, the New York Forest Declaration and the Bonn Challenge [4, 20, 39, 52]. In addition, expanding the sown area through irrigation and fertilization increases crop yield [25], which not only enhances food security, but also promotes regional greenness. Studies have shown that since 2000, the increase in global vegetation growth is mainly caused by the greening of forests and cropland [4]. Many studies have attached great importance to the impact of afforestation on regional vegetation growth [4, 21, 33, 40]. However, the quantitative effects of expansion of irrigated cropland (ICE) on NDVI and GPP are still uncertain.

Vegetation growth dynamics are influenced by a combination of factors, including biogeochemical factors, as well as land-use-related drivers (e.g., land use changes and agricultural management measures) [32, 56]. Globally, elevated atmospheric CO2 concentration (eCO2) is considered to be the most important driver of vegetation greening [56]. In arid regions, climate change has had an important impact on vegetation growth [9, 26, 56]. Recently, as atmospheric CO2 concentration continues to increase, the warming-induced atmospheric aridity is intensifying, which generates a negative impact on vegetation growth [27, 31, 50]. Studies have shown that the physiological response of plants to eCO2 can enhance leaf photosynthesis, improve water use efficiency [7], and alleviate some of the impact of atmospheric drought on vegetation growth [2, 48]. However, some studies have reported that enhanced vegetation growth means more water and nutrients consumption, aggravating the water scarcity in arid regions [8, 22], and the CO2 fertilization effect will decrease due to water and nutrient limitations [17, 38, 42, 44]. These indicate that in drier conditions, the fertilization effect of CO2 may decrease. A critical gap in our understanding pertains to the synergy and antagonism between different drivers.

The inland arid region of Northwest China (IANC) is located in the hinterland of the Eurasian continent (Additional file 1: Fig. S1a). It is an important strategic place of the ancient and modern ‘Silk Road’. It covers an area of 2.115 million km2, equivalent to about half of the area of the European Union. Four of the eight deserts in China are completely contained in the IANC, and two deserts are adjacent to it. Among them, the Taklimakan Desert is the largest desert in China and the second largest mobile desert in the world. The agriculture in IANC is highly dependent on irrigation, and the oasis agriculture formed at the edge of the desert has become important ecological barrier. Since the mid-1990s, large-scale ICE has taken place to meet the surging food demand, increasing the area of irrigated cropland by about 55% (Additional file 1: Figs. S2b, S3). In addition, over the last few decades IANC has experienced remarkable climatic warming, and the relative humidity decreased significantly after 1990s (Additional file 1: Figs. S2c, S4c), which has had an impact on regional ecological vegetation. It is urgent to clarify the impacts of various factors on regional vegetation growth in such an ecologically fragile environment, and it is necessary to explore the interaction between human activities and ‘natural’ factors, as well as different ‘natural’ factors to ensure food security and ecological health in this and other arid regions. In this paper, satellite-based GIMMS3g NDVI and MODIS NDVI data sets were used to show the greening pattern in IANC from 1982 to 2018. We then quantified the effects of biogeochemical factors (including atmospheric CO2 concentration and climate) and land-use-related drivers (including ICE and change of nitrogen fertilization per unit area) on IANC’s NDVI variation in the growing season (from May to October). This was achieved through trend decomposition and three machine learning algorithms. In addition, the cumulative GPP in the growing season was calculated using a revised remote sensing-based light use efficiency (LUE) model, and model experiments were performed to quantify the effects of various factors on ecosystem productivity.

Methods

Data sets

The data used in this paper include satellite remote sensing data, meteorological data, reanalysis data, and statistical yearbook data. The time period, temporal and spatial resolutions, and sources of the data are shown in Additional file 1: Table S1.

Vegetation greenness

As an indicator of vegetation greenness, the satellite-based NDVI data was used in this study, including the GIMMS3g NDVI and the MOD13A3 NDVI (see details in Additional file 1: Table S1). The multi-year monthly average values of GIMMS3g NDVI in IANC were calculated, and the months with NDVI greater than 0.1 were designated as the growing season (from May to October). The growing season NDVI was integrated by taking the average value from May to October. The GIMMS3g NDVI is represented on a resolution of 0.0833°. The MOD13A3 NDVI was resampled to the same spatial resolution as the GIMMS3g NDVI by the nearest neighbour method, and the variation and magnitude of the GIMMS3g NDVI and the MOD13A3 NDVI were almost the same in the study area (Additional file 1: Fig. S5). Therefore, the MOD13A3 NDVI was used to supplement the GIMMS3g NDVI to obtain vegetation greenness information for a longer time period (1982–2018).

Human activities

Two land use data sets (see Additional file 1: Table S1 for details) were collected as indicators of spatial change in regional irrigated cropland area from 1982 to 2018, and the irrigated cropland area and its changes displayed in these two sets of data are highly consistent (Additional file 1: Fig. S3). By summing the cropland area of all 0.00833° pixels in each 0.0833° pixel, the data were aggregated to the same spatial resolution as the NDVI data.

We used data on nitrogen application per unit area in Xinjiang and the Hexi Corridor from the Xinjiang Statistical Yearbook and the Gansu Development Yearbook. In addition, spatial nitrogen application data were obtained from a global data set [23]. There was a significant linear correlation between the nitrogen fertilization data in the statistical yearbook and the global data set (Additional file 1: Fig. S6a), so the statistical yearbook data were used to correct and extend the data in the global data set.

Meteorological data and atmospheric CO2 concentration data

Meteorological data mainly came from the China National Meteorological Center. The data collected at 110 meteorological stations were interpolated to 0.0833° spatial grid data using the nearest neighbor method. The cumulative solar radiation for the growing season was calculated from the solar radiation data in the re-analyzed data set and the record of sunshine hours from the meteorological stations.

Atmospheric CO2 concentration data observed at the Global Atmospheric Background Station at Waliguan was obtained from the National Cryosphere Desert Data Center. Waliguan Station is the only atmospheric background observation station in the inland region of the Northern Hemisphere, and its observations reflect the atmospheric greenhouse gas concentrations and their changes in the mid-latitude inland regions of the Northern Hemisphere. In addition, longer observation from the Mauna Loa Observatory was obtained from the Global Monitoring Laboratory. The interannual variations in CO2 data from the two stations were highly consistent (Additional file 1: Fig. S6b), so the observations from Mauna Loa Observatory were used in calculations.

Trend change point detection

The Akaike Information Criterion (AIC) and segmented regression [1] (a linear regression that modifies its slope at a certain value of the predictor, or threshold) were used to determine the trend change points of interannual variations of observation elements. We used AIC to choose the model that best-fitted the data. When non-linear regression (e.g., quadratic function) was a better fit to the data, trend change point may be present, and in these cases, the trend change point is defined by segmented regression. The segmented regressions were achieved through the ‘SiZer’ package in R, and the bootstrap sampling on each database was set 1000 times, with a confidence interval of 95%.

After AIC judgement, we found that there were trend change points in the trends of ‘growing season’ mean GIMMS3g NDVI in the natural vegetation zone, the irrigated cropland area, and the relative humidity. The trend change points of the above three elements are 1990, 1995 and 1992 respectively (Additional file 1: Fig. S2). We selected 1995 as the dividing point between the early and later periods, because 1995 is not only the change point for the trend in irrigated cropland area, but also falls within the 95% confidence interval of the trend change points for natural vegetation NDVI and relative humidity (Additional file 1: Fig. S2).

Effect of ICE on NDVI calculated by trend decomposition

We used the trend decomposition method to quantify the effect of ICE on NDVI. We use a linear trend to represent the change in NDVI over a period of time. First, we calculated the growing season NDVI trend of each 0.0833° pixel (pel0.0833) from1982 to 2018. Secondly, the pels0.0833 with the proportion of cropland greater than 10% were identified as cropland pixels in 1982 and 2018, and then the cropland change pixels between 1982 and 2018 (hereinafter referred to as ICE pixels) were identified (Additional file 1: Fig. S7a). The pels0.0833 with the proportion of natural vegetation (uncultivated vegetation) greater than 60% coverage, both in 1982 and 2018, were defined as natural vegetation pixels. The proportions of cropland and natural vegetation were derived from the land use data products. We then decomposed the NDVI trends of ICE pixels by considering two cases:

  1. a)

    If there was natural vegetation in the ICE pixel before 1982, it was considered that the NDVI variation trend caused by ICE was the actual NDVI trend minus the reference trend value [Eq. (1)].

  2. b)

    If there was no natural vegetation in the ICE pixel before 1982, it was considered that the variation of NDVI was caused by ICE [Eq. (2)]:

    $${Trend}_{ICE}^{\mathrm{^{\prime}}}={Trend}_{actrual}-{Trend}_{ref}$$
    (1)
    $${Trend}_{ICE}^{\mathrm{^{\prime}}}={Trend}_{actrual}$$
    (2)

Where \({Trend}_{ICE}^{\prime}\) represents the NDVI trend caused by ICE in the ICE pixel, \({Trend}_{actrual}\) represents the actual NDVI trend in the corresponding pel0.0833, and \({Trend}_{ref}\) represents the reference trend value which is calculated at the grid with spatial resolution of 0.75° (grid0.75) (Additional file 1: Fig. S7c). Each grid0.75 contains 81 pels0.0833. The \({Trend}_{ref}\) is calculated by averaging the NDVI trends of all natural vegetation pixels in the grid0.75. It should be noted that we assume that at each grid0.75, all pels0.0833 have the same climate background and are taken to have experienced the same agricultural management measures. Through the above operation, we calculated the effects of ICE on NDVI in ICE pixels. However, the pels0.0833 that had cropland before 1982 and the pels0.0833 with small change in the irrigated cropland area were ignored. Where there was a significant linear correlation between cropland area change and NDVI variation (Additional file 1: Fig. S7b), the NDVI variation caused by ICE was spatially extrapolated through the actual irrigated cropland proportion of change of the pel0.0833 and the \({\overline{Trend}}_{ICE}^{0}\) of the grid0.75 where the pel0.0833 is located [Eqs. (3), (4)]:

$${Trend}_{ICE}^{0}={Trend}_{ICE}^{\prime}/{ICE}_{actual}^{\prime}$$
(3)
$${Trend}_{ICE}={\overline{Trend}}_{ICE}^{0}\times {ICE}_{actual}$$
(4)

where \({ICE}_{actual}^{\prime}\) is actual cropland proportion changes of the ICE pixels; \({Trend}_{ICE}^{0}\) is the NDVI trend caused by the change of unit cropland proportion in the ICE pixels, and the \({\overline{Trend}}_{ICE}^{0}\) is the average value of \({Trend}_{ICE}^{0}\) in a grid0.75; \({Trend}_{ICE}\) is the NDVI trend caused by the ICE (Additional file 1: Fig. S7d), and \({ICE}_{actual}\) is the actual irrigated cropland proportion changes in the pels0.0833.

Machine learning algorithms for estimating the effects of various factors on NDVI

Three machine learning algorithms (Additional file 1: Table S2) were used to model and estimate NDVI with variation in different biogeochemical factors (atmospheric CO2 concentration, photosynthetically active radiation, relative humidity and air temperature) and human activities (irrigated cropland area, nitrogen fertilization). Satellite-derived growing season mean NDVI data were used for model training and verification. If the numbers of target data are small, the machine learning algorithm may ‘over fit’. Therefore, we used grid0.75 as the calculation unit to expand the sample size. We had 81 pels0.0833 in each grid0.75, and each pel0.0833 contained 35 years of data, that is, there were 2835 samples at each grid0.75 for model construction. 70% of the samples were randomly selected for model training and the remaining 30% for testing. Through the random search and the cross validation, the optimal parameters were selected for each grid. We use the coefficient of determination (R2) between the simulated values and observed values in the test set as the evaluation metric for the model, and our results suggest that machine learning algorithms can effectively estimate NDVI (Additional file 1: Fig. S8). Among them, the simulation efficiency of the Random Forest and MLP Neural Network is superior to that of the Support Vector Machine. We synthesized the results of the three models by calculating their average.

We calculated the effects of various drivers on variation of NDVI through model experiments. Firstly, we made all factors change over time to calculate \({NDVI}_{ALL}\) (a normal model run). Next, we kept the variable x constant (e.g., ICE, eCO2) to calculate \({NDVI}_{x}\) with other factors changing over time. The difference between \({NDVI}_{ALL}\) and \({NDVI}_{x}\) is the NDVI variation driven by variable x:

$$\Delta {NDVI}_{xi}={NDVI}_{ALLi}-{NDVI}_{xi}$$
(5)

where \(\Delta {NDVI}_{xi}\) is the NDVI variation driven by variable x in the ith year. The trend of NDVI driven by x was calculated by the least square method. The sensitivity of NDVI to the ICE was analyzed through linear regression. The calculation equation is as follows:

$$\Delta {NDVI}_{ICEk}={\gamma }_{ICE}\times {\mathrm{\Delta ICE}}_{k}+\varepsilon$$
(6)

where \(\Delta {NDVI}_{ICEk}\) is the NDVI variation driven by ICE between 1982 and 2018 in the kth grid0.75, and \({\mathrm{\Delta ICE}}_{k}\) is ICE proportion (i.e., the proportion of irrigated cropland area to total vegetation area between 1982 and 2018) in the kth grid0.75; \({\gamma }_{ICE}\) is the sensitivity of NDVI to ICE; \(\varepsilon\) is the residual term.

Regional productivity estimation and attribution

A revised light use efficiency (LUE) model was used to estimate regional productivity of vegetation (GPP), and the calculation equation is as follows [10, 50]:

$$GPP=PAR\times FPAR\times {\varepsilon }_{max}\times {C}_{s}\times {T}_{s}\times {W}_{s}$$
(7)

where PAR is photosynthetically active radiation, calculated as 1/2 of the solar radiation (unit: MJ/m2/month); FPAR is the fraction of PAR absorbed by the vegetation canopy, which can be calculated by NDVI [Eq. (8)] [3]; \({\varepsilon }_{max}\) is the maximum LUE, and we take the fixed empirical value of 2.76 g C/MJ [12, 37, 50] to prevent the effects of some factors from being repeatedly considered; \({C}_{s}\), \({T}_{s}\) and \({W}_{s}\) represent the downward-regulation scalars for the respective effects of eCO2, air temperature and moisture conditions on the LUE [3, 10, 29, 50, 53]. The effect of atmospheric CO2 concentration on LUE (\({C}_{s}\)) was obtained from the work of Yuan et al. [50], and see details in the supplementary information. \({T}_{s}\) was calculated by mean air temperature (\({T}_{a}\), unit: \(^\circ{\rm C}\)) and the optimal air temperature which is the average temperature of the month with the largest NDVI in a year [10] (see details in the supplementary information). \({W}_{s}\) was calculated using regional evapotranspiration (ET) and potential evapotranspiration [PET, Eq. (9)]. The regional ET was obtained through the GLASS ET product [49] and PET was calculated using the Hargreaves equation [Eq. (10)]:

$$FPAR=1.24\times NDVI+0.618$$
(8)
$${W}_{s}=0.5+0.5ET/PET$$
(9)
$$PET={C}_{o}{R}_{a}({T}_{a}+17.8){({T}_{max}-{T}_{min})}^{0.5}$$
(10)

where \({T}_{max}\) and \({T}_{min}\) represent the maximum and minimum air temperature in the calculation period, respectively (unit: \(^\circ{\rm C}\)). \({R}_{a}\) is the theoretical solar radiation at the top of the atmosphere (unit: MJ/m2/month); the \({C}_{0}\) is the conversion coefficient with a value of 9.39 × 10–4.

In the cropland zone, we made some adjustments to the above LUE model. In addition to natural factors, we also considered the effects of ICE and nitrogen fertilization on GPP. The calculation is as follows:

$$GPP=PAR\times FPAR\times {\varepsilon }_{max}\times {C}_{s}\times {T}_{s}\times {W}_{sICE}\times {N}_{s}$$
(11)

We found that there was a significant linear correlation between the ICE and regional ET, so we substituted this relationship into Eq. (9) to obtain \({W}_{sICE}\):

$$ET=aICE+b$$
(12)
$${W}_{sICE}=0.5+0.5(aICE+b)/PET$$
(13)

here, a and b are the coefficients in the regression equation between ET and ICE, which were quantified by the ET of GLASS products and irrigated cropland area in the land use data set, and each pel0.0833 has specific values of a and b. Within a certain range, there is a stable proportional relationship between biomass and GPP [16, 24], so we collected nitrogen fertilization and biomass data for three main crops (corn, wheat and cotton) in IANC to quantify \({N}_{s}\). Wheat, corn, and cotton occupy 22.4%, 15.4%, and 31.8% of the planting area, respectively. Over the years, these three crops together have consistently comprised around 70% of the total cultivated area. The regression relationship between nitrogen fertilization and relative biomass (biomass after nitrogen application divided by biomass without nitrogen application) was established, and the AIC was used for judgment. The results showed that there was a significant quadratic function relationship between nitrogen fertilization and relative biomass (Additional file 1: Fig. S9a). We integrated data of different crops to establish a unified relationship (Additional file 1: Fig. S9b).

The revised LUE model was verified by the GPP observation data for the experimental sites, and the simulated results matched well with the observed values (Additional file 1: Fig. S10). GPP was then calculated by weighted averaging based on the proportions of cropland and natural vegetation in each of the pels0.0833:

$${GPP}_{j}={GPP}_{c}\times {\rho }_{c}+{GPP}_{nv}\times {\rho }_{nv}$$
(14)

where \({GPP}_{j}\) represents the GPP in the jth pel0.0833; \({GPP}_{c}\) and \({GPP}_{nv}\) represent the GPP of cropland and natural vegetation, respectively in the jth pel0.0833; \({\rho }_{c}\) and \({\rho }_{nv}\) represent the proportions of cropland and natural vegetation area in the jth pel0.0833.

We performed numerical simulations to evaluate the effects of various factors on GPP with the revised LUE model. In the LUE model, various factors affect GPP through two aspects: on one hand, drivers regulate GPP by affecting the FPAR which is calculated directly through NDVI, so we call this path the FPAR variable [Eq. (8)]; on the other hand, GPP is affected by various factors through controlling the downward-regulation scalars (e.g., \({C}_{s}\), \({T}_{s}\)) of the LUE, and we call this path LUE variable [Eq. (16)]:

$$GPP=PAR\times FPAR\times LUE$$
(15)
$$LUE=f(ICE,eCO2,NFER,RH,{T}_{a})$$
(16)

where NFER is the nitrogen fertilization and RH is the relative humidity. In the normal model run, all variables changed with time both in the FPAR path and LUE path, which we called \({GPP}_{ALL}\). Next, we calculated \({GPP}_{x}\) by allowing other factors to change over time while holding the variable x constant at its initial baseline level in the LUE path and substituting NDVI with the \({NDVI}_{x}\) in the FPAR path. \({NDVI}_{x}\) was calculated by keeping the variable x to maintain the initial value and other factors to change over time in the machine learning algorithm models mentioned above. The difference between \({GPP}_{ALL}\) and \({GPP}_{x}\) is the GPP variation driven by variable x [Eq. (17)]. In addition, we conducted another set of simulation experiments by substituting NDVI with the \({NDVI}_{x}\) in the FPAR path and making all variables change with time in the LUE path, which was named \({GPP}_{x}^{\prime}\). The difference between \({GPP}_{x}^{\prime}\) and \({GPP}_{x}\) is the GPP variation driven by variable x through the LUE path [Eq. (18)]:

$${\Delta GPP}_{xi}={GPP}_{ALLi}-{GPP}_{xi}$$
(17)
$$\Delta {GPP}_{xi}^{\prime}={GPP}_{xi}^{\prime}-{GPP}_{xi}$$
(18)

where \(\Delta {GPP}_{xi}\) is the GPP variation driven by variable x in the ith year, and \(\Delta {GPP}_{xi}^{\prime}\) is the GPP variation driven by variable x through the LUE path. The trend of GPP driven by x was calculated by the least square method. In the same way, the sensitivity of GPP to the change of ICE was analyzed through linear regression:

$$\Delta {GPP}_{ICEk}={\beta }_{ICE}\times {\mathrm{\Delta ICE}}_{k}+\varepsilon$$
(19)

where \(\Delta {GPP}_{ICEk}\) is the GPP variation driven by ICE between 1982 and 2018 in the kth grid0.75, and \({\mathrm{\Delta ICE}}_{k}\) is the ICE proportion (i.e., the proportion of the total vegetation area which is irrigated cropland area), a value which changes between 1982 and 2018 in the kth grid0.75; \({\beta }_{ICE}\) is the sensitivity of GPP to ICE; \(\varepsilon\) is the residual term.

Results

Vegetation growth changes in IANC

The GIMMS3g NDVI data showed the interannual variation of the average NDVI during the growing season in IANC, with the trends of NDVI calculated by the least square method, as shown in Fig. 1a, b, c. The trends were tested and the trend change points were determined by the AIC and segmented regressions. From 1982 to 2015, the NDVI increased significantly, especially in the cropland zone (Fig. 1a, b). The NDVI trends are 0.000572 yr−1 (p = 3.77e-8) and 0.00154 yr−1 (p = 1.5e-11) for the whole area and the cropland zone, respectively. The trend change point was found in the NDVI variation in the natural vegetation zone (mainly woodland and grassland) (Additional file 1: Fig. S2a). Since the 1990s, the growth rate of natural vegetation NDVI has slowed down. The study period was divided into two periods, before and after 1995. The NDVI in the natural vegetation zone increased significantly over the early period (from 1982 to 1995), but the growth trend in the later period (from 1995 to 2018) decreased by 84.4% compared with the early period, and the slope became statistically insignificant (Fig. 1c). Figure 2a, c, d shows the spatial distributions of GIMMS3g NDVI trends. After 1995, the intensity and extent of vegetation browning increased significantly, but the cropland zones in the central part of Northern Xinjiang, the Tarim River Basin, and the central and western Hexi Corridor continued greening (Fig. 2a, c, d, Additional file 1: Fig. S1). In addition, the MOD13A3 NDVI data also shows that from 2001 to 2018, the vegetation greenness in the cropland zone increased significantly during the growing season, and the zones with vegetation browning were mainly concentrated in the mid-range of the Tianshan Mountains and the northern mountainous area of Northern Xinjiang (Fig. 2b, Additional file 1: Fig. S1), which was consistent with the GIMMS3g NDVI result.

Fig. 1
figure 1

Variations of vegetation growth in IANC. Interannual variations of growing season (from May to October) mean GIMMS3g NDVI in a the whole region, b the cropland zone, and c the natural vegetation zone. Interannual variations of growing season accumulative GPP calculated by revised LUE model in d the whole region, e the cropland zone, and f the natural vegetation zone. The shaded areas in the figure represent the 95% confidence intervals of the fitted curves. In panels (c, f), we divided two periods, before and after 1995

Full size image
Fig. 2
figure 2

Spatial patterns of trends in GIMMS3g NDVI in inland arid region of Northwest China (IANC) (unit: yr.−1). Growing season mean GIMMS3g NDVI trends during a 1982–2015, c 1982–1995, and d 1995–2015. b Spatial patterns of trends in growing season mean MOD13A3 NDVI during 2001–2018

Full size image

We used the revised LUE model to estimate regional GPP [10, 50]. The results of the revised LUE model matched well with the measured GPP at the ground experimental sites, especially in the cropland zone (Additional file 1: Fig. S10). We also verified the result of revised LUE model by comparing it with other GPP products. The cumulative GPP in the growing season estimated by the revised LUE model is close to the GLASS GPP product, and is consistent with the variation of MODIS GPP, but the value is greater than that of MODIS GPP (Additional file 1: Fig. S11).

From 1982 to 2018, the cumulative GPP in the growing season in IANC increased significantly. The GPP trends are 2.77 g C m−2 yr−1 (p = 9.46e-13) and 9.29 g C m−2 yr−1 (p = 2.08e-16) for the whole area and the cropland zone, respectively (Fig. 1d, e). After 1995, the growth of GPP in the natural vegetation zone stalled. The increasing rate of GPP at the late period of the natural vegetation zone decreased by 79.1% compared with the early period (Fig. 1f).

Attribution of vegetation growth trend

The relative importance of various factors in the attribution calculated by the three machine learning algorithms are basically consistent (Additional file 1: Fig. S12a), and we integrated the results of the three models by averaging. The machine learning algorithms and trend decomposition method showed similar effects of ICE on NDVI (Additional file 1: Fig. S12b). Figure 3 shows the attribution of NDVI trend obtained from the model experiments using the machine learning algorithms. ICE was the dominant factor in NDVI variation in the cropland zone, explaining 61% (61.12 ± 9.78%) of NDVI variation (Fig. 3a). In the natural vegetation zone, NDVI variation was mainly due to eCO2 (Fig. 3b). Over the whole region, the effect of eCO2 also dominated, explaining 70% (70.25 ± 7.72%) of the NDVI variation (Fig. 3c). Among the climatic factors, change in relative humidity had the greatest effect on NDVI, and the relative humidity-driven NDVI trend is significantly negative in both cropland zone and natural vegetation zone (Fig. 3a, b). From 1982 to 1995, the regional relative humidity increased, and the NDVI trend driven by relative humidity was positive as well. In the late period (1995–2018), due to the significant decline of regional relative humidity, the NDVI trend driven by relative humidity also became negative (Additional files 1: Figs. S4c and Fig. 3d, e, f). After 1995, the effect of ICE on NDVI trend increased substantially, while the effect of eCO2 decreased (Fig. 3d, e, f). ICE showed an equally significant effect on NDVI variation as did eCO2 during the latter part of the study period (ICE: 0.356 ± 0.502, eCO2: 0.408 ± 0.992, unit: 10–3 yr−1) in the whole region (Fig. 3f).

Fig. 3
figure 3

Attribution of variation in NDVI calculated from model experiments with three machine learning algorithms. In panel (a), (b), (c), the effects of various factors on NDVI in the cropland zone, the natural vegetation zone and the whole region are shown respectively from 1982 to 2018, while panel (d), (e), (f) represent the attribution of NDVI variation in the two periods. The “early” notation indicates the period of 1982–1995, and “late” indicates 1995–2018. The drivers include irrigated cropland expansion (ES_ICE), nitrogen fertilization change (ES_NFER), elevated atmospheric CO2 concentration (ES_eCO2), photosynthetically active radiation change (ES_PAR), relative humidity change (ES_RH) and air temperature change (ES_Ta). ES_ALL notation indicates the result of model simulation with all factors changing over time. Error bars show the 95% confidence intervals of the trends. Two asterisks indicate that the trend is extremely significant (P < 0.01), and one asterisk indicates that the trend is significant (P < 0.05)

Full size image

Model experiments were conducted through the revised LUE model to quantify the effects of various factors on GPP variation. In the cropland zone, ICE is the dominant factor in GPP variation, explaining 58% (57.93 ± 11.66%) of the GPP trend (Fig. 4a). In the natural vegetation zone, as well as the whole region, the eCO2 had the greatest influence on GPP variation (Fig. 4b, c), which is consistent with the attribution of NDVI trend. However, over the whole region, the contribution of ICE to GPP trend reached 42% (41.99 ± 8.45%), which was close to that of eCO2 (49.79 ± 12.43%) (Fig. 4c). The effects of various factors on GPP trend were also analyzed in two periods (Fig. 4d, e, f), and the GPP trend driven by relative humidity change was positive at the early period and negative at the late period. The effect of eCO2 decreased during the late period, while the effect of ICE increased significantly, and the ICE-driven GPP trend over the whole zone exceeded that of eCO2 at the late period (ICE: 1.30 ± 0.229, eCO2: 0.800 ± 0.265, unit: g C m−2 yr−1) (Fig. 4f).

Fig. 4
figure 4

Attribution of trend in GPP calculated by model experiments of revised LUE model. In panel (a), (b), (c), the effects of various factors on GPP in the cropland zone, the natural vegetation zone and the whole region are shown respectively from 1982 to 2018, while panel (d), (e), (f) represent the attribution of GPP trend in the two periods. The “early” notation indicates the period of 1982–1995, and “late” notation indicates 1995–2018. The drivers include irrigated cropland expansion (ES_ICE), nitrogen fertilization change (ES_NFER), elevated atmospheric CO2 concentration (ES_eCO2), photosynthetically active radiation change (ES_PAR), relative humidity change (ES_RH) and temperature change (ES_Ta). ES_ALL notation indicates the result of model simulation with all factors changing over time. Error bars show the 95% confidence intervals of the trends. Two asterisks indicate that the trend is extremely significant (P < 0.01), and one asterisk indicates that the trend is significant (P < 0.05)

Full size image

ICE threshold

The percentage of the change of irrigated cropland area as a function of total vegetation area (ICE proportion) reached 5.25% by 2018. As the ICE proportion grows, the ICE-driven NDVI (or GPP) variation becomes greater, and eventually exceeds the effect of eCO2. The relationships of the ICE-driven NDVI and GPP variations against ICE proportion between 1982 and 2018 were calculated spatially [Eq. (6), (19)], and showed significant linear correlations (Fig. 5). This means that the sensitivities of NDVI and GPP variations to ICE are consistent in the spatial grids0.75, with values of 0.0022 for NDVI and 12.9 for GPP. During the whole period, only if the ICE proportion reaches 10.4%, ICE-driven NDVI variation can exceed the eCO2 effect in the whole region (Fig. 5a), while the threshold for ICE-driven GPP variation exceeded eCO2-driven GPP variation is 6.55% (Fig. 5b). Since regional irrigated cropland began to expand after 1995, the ICE-driven vegetation growth variation in the late period was close to that recorded during the entire period. The GPP variation driven by ICE outstripped that by eCO2 when the ICE proportion exceeded 4.37% in the late period. The ICE-driven NDVI variation reached a level equivalent to the eCO2 effect in the late period, when the ICE proportion increases to 6.38% (Fig. 5).

Fig. 5
figure 5

Threshold of the ICE-driven NDVI and GPP variation (\(\Delta {NDVI}_{ICE}\) and \(\Delta {GPP}_{ICE}\)) outpace that of eCO2. ICE proportion is the percentage of the change of irrigated cropland area to total vegetation area. a Relationship of the \(\Delta {NDVI}_{ICE}\) and ICE proportion. b Relationship of the \(\Delta {GPP}_{ICE}\) and ICE proportion. The shaded areas represent the 95% confidence intervals of the fitted curves. The red dots in the figure represent ICE proportion and \(\Delta {NDVI}_{ICE}\) (\(\Delta {GPP}_{ICE}\)) of different 0.75° grids in the region. The green dots represent ICE proportion and \(\Delta {NDVI}_{ICE}\) (\(\Delta {GPP}_{ICE}\)) between 1982 and 2018 in the whole region. The “ + ” represent ICE proportion and \(\Delta {NDVI}_{ICE}\) (\(\Delta {GPP}_{ICE}\)) between 1995 and 2018 in the whole region

Full size image

Discussion

Effects of ICE on NDVI and GPP

The term ‘Vegetation growth dynamics’ is often used with reference to changes in vegetation structures and functions, which are affected both by natural drivers and by human activities. As such, such growth dynamics can play a critical role in the material and energy flows of ecosystems [11, 32]. The growth dynamics of vegetation can be quantified by various indicators, such as vegetation greenness, vegetation cover, and ecosystem productivity [13, 21, 45]. Satellite-based NDVI represents vegetation greenness, which can reflect the capacity of vegetation to absorb PAR, as well as the extent of vegetation coverage. The ecosystem productivity quantified by GPP represents the carbon sequestration capacity of vegetation and is an important indicator of the amount of carbon cycling between ecosystems and the atmosphere [6, 36]. The annual precipitation in the desert oasis area in northwest China is less than 200 mm, and the natural vegetation cover is mainly desert steppe and low shrubs [47]. Since the mid-1990s, substantial areas of cropland have been developed and improved in the region through diverting water for irrigation to previously uncultivated areas and to grassy areas in IANC. This has resulted in enhanced vegetation growth regionally. The ICE in IANC supplements the water and nutrient conditions required for vegetation growth through irrigation and fertilization, which promotes regional vegetation greening, increases fixation of CO2 by plants and further increases ecosystem productivity [43, 51, 55].

After 1995, ICE in IANC had a significant positive effect on regional vegetation growth, and had a greater effect than NDVI on GPP (Figs. 3c, f and 4c, f). In the LUE model, GPP is calculated by multiplying the PAR, FPAR, and LUE [Eq. (15)], where FPAR is related to vegetation greenness and can be estimated by NDVI. The environmental factors that affect FPAR also affect the LUE. Figure S13 shows the GPP trends for the whole region from 1982 to 2018 driven by environmental variables through effects on LUE. During this period, the ICE contribution was the most significant of that of all driving variables (Additional file 1: Fig. S13c). As a result, ICE not only increases NDVI, which improves the ability of the vegetation canopy to absorb PAR, but more importantly, increases LUE, and further improves regional carbon sequestration.

ICE has slowed the decline of vegetation growth in the region

Plant photosynthesis is enhanced as by the climate change-driven enhancement of atmospheric CO2 concentration and this can significantly promote vegetation growth. In addition, the physiological response of plant leaves to eCO2 can improve water use efficiency [7]. Globally, eCO2 is considered the most important factor in vegetation greening [32, 56]. However, climate change is also increasing air temperatures, crop water losses and the development of drought on all continents of the world [15, 48], thereby adversely affecting plant growth [18, 35, 50]. Some studies suggest that the positive influence of eCO2 on vegetation growth can offset the negative impact of climatic drought, which can avoid further expansion of drylands [2, 48]. However, an increasing number of studies also show that the CO2 fertilization effect is reduced by the limitation of water and nutrient conditions in the ecosystem [17, 34, 42, 44]. The observed data from meteorological stations show that since the 1980s, the air temperature and precipitation in IANC have increased. The early precipitation trend is 2.00 mm/yr, and the later precipitation trend is 1.11 mm/yr, which decreased by 45% compared with the earlier period. However, the air temperature continued to rise, and the relative humidity increased significantly in the early period and decreased significantly in the later period. As a result, the regional climatic aridity intensified in more recent time (Additional file 1: Fig. S4). We found that as the relative humidity-driven NDVI and GPP trends changed from positive to negative, the fertilization effect of CO2 on NDVI and GPP also decreased (Figs. 3d, e, f and 4d, e, f), although the atmospheric CO2 concentration continued to increase (Fig. S6b). In the later period, the natural vegetation \(\Delta {NDVI}_{eCO2}\) is 0.338 \(\pm\) 0.094 (unit: 10–3 yr−1) and the \(\Delta {NDVI}_{RH}\) is − 0.15 \(\pm\) 0.119 (\(\Delta {GPP}_{eCO2}:0.644\pm 0.550,\Delta {GPP}_{RH}\): -0.298 \(\pm\) 0.568, unit: g C m−2 yr−1, Figs. 3e and 4e). The fertilization effect of CO2 alleviated the effect of drought, but its effect was not enough to alleviate the decline of NDVI and GPP trends. The trends of NDVI and GPP in cropland zone did not decrease as in natural vegetation zone, mainly due to the significant positive impact of ICE (Figs. 3d and 4d). The results of this paper show that in arid regions, it is ICE that mitigates the decline of vegetation growth trend, rather than the CO2 fertilization effect, in the case of intensified climatic drought. In addition, over the entire region, the effect of ICE on the NDVI or GPP variation, was close to or even has outpaced that of eCO2 at the late period (Figs. 3f and 4f). Our results indicate that the impact of ICE on vegetation growth has kept pace with the CO2 fertilization effect in IANC, which is due to the significant increase in irrigated cropland area and the reduced effect of eCO2 in the later years of the period under study.

In summary, ICE plays a significant role in the growth of vegetation in arid regions. We highlight the role of ICE, which on one hand confirms the importance of water availability in arid areas. Although the fertilization effect of CO2 has a considerable impact on vegetation growth, the limitation of water conditions still predominates in arid regions [19, 32, 54]. On the other hand, we emphasize the impact of human activities on fragile ecosystems, which continuously challenge natural regulatory capacities. This finding is consistent with many other research results [4, 30].

ICE may exacerbate regional water conflicts in the future

The fertilization effect of CO2 and human land use management promote regional vegetation growth, however, vegetation expansion, especially the increase of the afforestation and cropland, will exacerbate regional water conflicts in arid regions [8, 22]. In recent years, as the population grows rapidly and the demand for food increases, more and more land has been reclaimed. The data from the statistical yearbook show that the agricultural production in the northwest region has significantly increased, ensuring regional food security (Additional file 1: Fig. S14). More output means more water consumption, and agricultural water consumption also significantly increases (Additional file 1: Fig. S14). For the total amount of water resources in the region remains stable, the increase of agricultural water use will inevitably crowd out ecological water use, which will lead to the cutoff of downstream inland rivers and the formation of more groundwater table drop funnels. Terrestrial water storage observations from NASA’s Gravity Recovery and Climate Experiment (GRACE) satellite mission have revealed significant decline in regional water resources (Additional file 1: Fig. S14). The soil drought in some areas will intensify, and the natural vegetation will be affected and even face the threat of death, which will have an adverse impact on the regional ecology [14]. Facing the twofold pressure of water shortage and the adverse effects of climate change, maintaining an appropriate scale of irrigated cropland is critical for ensuring food security and ecological health in IANC in the future.

Assessment of data and methodological reliability

In this study, we employed multiple sources of data and a variety of methodologies to demonstrate the impact of natural factors and human activities on regional NDVI and GPP. To ensure the reliability and validity of our research findings, we conducted a comprehensive evaluation of the completeness and accuracy of the data, as well as the appropriateness of the methods. First, we utilized GIMMS3g NDVI and MOD13A3 NDVI data to study regional greenness variations and as input data for the LUE model. These data sets have been widely used and analyzed by many scholars and are broadly considered to have high reliability and representativeness [28, 41, 46, 50]. The driving data, such as climate data and atmospheric carbon dioxide concentration, are derived from authoritative observation organizations, ensuring high quality. Second, we integrated three machine learning algorithms to attribute NDVI, and subdivided larger grids into smaller ones to increase the data volume for training models, meeting the requirements of machine learning algorithms. The impacts of ICE on NDVI quantified by machine learning algorithms and trend decomposition methods are largely consistent, mutually validating each other. Finally, we used a revised LUE model to calculate GPP and attributed it. We collected real measurement data from four ground stations to validate the results of the improved LUE model, which showed that the LUE model performs well (Additional file 1: Fig. S10). Although there are certain limitations in the data and methods, we have taken measures to address them: (1) Some data have a short time range and can only be supplemented with other data to ensure the timeliness of the article. We have chosen data with consistent results for mutual verification and supplementation, and the supplementation period is short to ensure that it minimally impacts the main findings of our study. (2) The research objective of the text is to study the overall variation in regional vegetation growth indicators and their attribution. Using regional average values provides an overview of large-area trends, facilitating comparisons across the entire region, cropland zone, and natural vegetation zone. However, this approach lacks a more detailed and nuanced understanding of spatial variations, which should be addressed in future work through increased research on regional spatial changes.

Conclusions

Our results show that, from 1982 to 2018, vegetation growth was enhanced in IANC mainly under the influence of ICE and eCO2, with ICE being the dominant factor in cropland zones. Since the mid-1990s, the significantly decreased regional relative humidity has negatively affected vegetation growth, and the CO2 fertilization effects on NDVI and GPP also decreased. However, the irrigated cropland area has increased sharply, and its positive impact on vegetation growth approached or even exceeded that of eCO2 in the later period, which has delayed the slowing down of the growth trend. In the whole region, the change of irrigated cropland area accounts for 5.25% of the total vegetation area, and the variation of GPP driven by ICE reached a level equivalent to that by eCO2 at the late period when the ICE proportion reached 4.37%, while the ICE-driven NDVI variation will outpace the eCO2 effect when the ICE proportion raises to above 6.38%. Although our study highlights the positive effect of irrigation cropland expansion on vegetation growth in arid areas, further expansion of irrigated cropland in the future will exacerbate regional conflict between water demand and supply, which may adversely affect the growth of natural vegetation. In the future, the ecological impact of ICE needs to be in-depth studied.