Dynamic Vegetation Responses to Climate and Land Use Changes over the Inner Mongolia Reach of the Yellow River Basin, China

Abstract

Ecological protection and high-quality development of the Yellow River Basin (YRB), China, aroused remarkable concerns from China’s Central Government, and has been a major national strategy. The Inner Mongolia reach of the Yellow River Basin (IM-YRB) is a typical dryland with pervasive vegetation restoration through the actions of the ecological projects that have been conducted in recent years. However, how climate changes and human activities, such as land use and land cover (LULC) changes, jointly impact vegetation variations in this region remains poorly understood. Here, using an explainable machine learning technique, we evaluated linkages between the kernel normalized difference vegetation index (kNDVI) and air temperature, precipitation, soil moisture, and LULC changes, and relevant marginal contributions of these four drivers to the observed vegetation changes. The grassland fraction on a pixel level was selected as the quantitative LULC variable due to its key role in regional LULC. We found that interannual kNDVI changes in most areas of this study region were negatively sensitive to temperature, but positively sensitive to precipitation and soil moisture, with area fractions of 71.74%, 96.93%, and 89.33%, respectively. The area fraction of negative kNDVI sensitivity to LULC was roughly equivalent to that of positive kNDVI sensitivity. The contributions of air temperature, precipitation, soil moisture, and LULC to overall kNDVI changes were 21.54%, 33.32%, 32.19%, and 12.95%, respectively. Moisture conditions also play a critical role in vegetation changes, which was reflected by the fluctuating growth of kNDVI as interannual changes in precipitation. Nonetheless, kNDVI changes are also affected by LULC, and LULC became the dominant factor behind the kNDVI anomalies over the grassland restoration regions from barren over the IM-YRB. This research provides theoretical support for dryland vegetation restoration under the influence of climate change.

1. Introduction

Vegetation plays a crucial role in modulating the land vs. atmosphere exchange in energy, water, carbon, and momentum, and provides pivotal ecosystem services [1,2,3,4]. Changes in the vegetation cover can impact the global carbon and water cycles, affecting the climate patterns as a result [5,6]. Vegetation changes are driven by concurrent climatic anomalies, environmental changes [4,7,8], and by human activities, such as land use and land cover changes (LULC) [9,10,11]. Specifically, land cover changes are closely related to shifts in vegetation conditions, which are viewed as one of the main drivers behind terrestrial ecosystem productivity [11]. Therefore, a deep understanding of the vegetation changes and related interactions with their climate and anthropogenic drivers is critical to describe and predict the vegetation dynamics under the changing climate and develop effective policies and strategies to promote environmental conservation and sustainability.

The increasing suite of remotely-sensed vegetation indices renders it possible to monitor vegetation dynamics, such as the NDVI (normalized difference vegetation index) [12], the LAI (leaf area index) [13], the EVI (enhanced vegetation index) [14], the NPP (net primary productivity) [15], the kNDVI (kernel normalized difference vegetation index) [16], and so forth. Based on the remotely-sensed measure of vegetation changes at different scales, a huge body of studies have since appeared aiming to enhance the human understanding of drivers behind vegetation changes and their responses to climate changes [2,4,6,17,18]. First of all, the distribution of global vegetation is strongly controlled by precipitation or rainfall patterns [18,19]. Climate warming also has an important role in vegetation change or plant phenology [2,5,20]. Meanwhile, increasing the temperature exacerbates drought as the increase in evapotranspiration [15]. Using autocorrelation multiple linear regression to fit the NPP and meteorological components, Zhang et al. [4] assessed vegetation dynamics under climate changes across the Qinghai–Tibet Plateau. With a process-based model of plant growth, Higgins et al. [20] showed vegetational responses to soil moisture changes in dry–warm regions and to temperature variations in low-temperature regions. However, the attribution of dynamic vegetation responses to environmental variables and human activities have been poorly understood due to compound impacts from complex drivers and human interferences, as well as the spatial heterogeneity of these drivers.

Vegetation changes in ecologically vulnerable regions have aroused pervasive concerns [4,21,22,23]. Vegetation changes over the drylands are highly sensitive to both climate changes and human disturbances [24,25]. Global warming is generally expected to amplify the existing spatial patterns of moisture conditions, leading to the drylands becoming even more drier [26]. Nonetheless, satellite observations have shown a significant positive trend in dryland vegetation cover in recent decades [7,23]. Except for the increasing precipitation beyond “dry get drier”, human activities also have the potential to drive the growing dryland vegetation cover [23,25,27]. For example, India, the North China Plain, the US Great Plains, and South-East Australia all experienced intensive agricultural expansion, while large-scale ecological restoration programs have been implemented in semi-arid or sub-humid North and West China [28,29]. Great efforts have been devoted to unravel the interactions between human and natural factors and their impacts on vegetation cover changes at regional scales [30,31,32]. However, the overwhelming effects of anthropogenic land use changes enable the human–environment nexus to become a great challenge in understanding the drivers behind these vegetation changes.

Variable methods are available addressing the attribution of dynamic vegetation responses to climate changes and LULC. Zhu et al. [8] used process-based ecosystem models and the optimal fingerprint method for the evaluation of the different effects of climate change and LULC on global LAI variations. Using a residual method, Li et al. [30] quantified human contributions to vegetation changes over the Loess Plateau, China. This residual method was further improved using nonlinear fitted vegetation dynamics based on binary regression or machine learning techniques [33,34]. The process-based method is difficult to manipulate with considerable uncertainty from the model itself, data, and scales [34]. Meanwhile, the residual method is easy to use. However, variable vegetation–climate relationships in space and time cause overestimations of the human contribution to vegetation changes by the residuals [35]. Comparatively, machine learning methods can well capture nonlinear relationships between the vegetation and drivers [13].

The Inner Mongolia reach of the Yellow River Basin (IM-YRB) is a typical dryland and has experienced successful vegetation restoration in recent decades [32]. Ecological conservation policies significantly modified land cover types. However, the contribution of LULC to vegetation dynamics was not directly evaluated considering the climatic influences. Here, considering kNDVI as the indicator of vegetation, we attempted to attribute vegetation changes to climatic and LULC variations during the 2000–2020 period using an explainable machine learning approach [36]. Thus, the objectives of this paper are to: (1) characterize LULC and vegetation changes over the IM-YRB; (2) quantify the sensitivity of vegetation changes to climatic factors and LULC; and (3) differentiate the marginal contributions of climatic and LULC changes to vegetation variations. This study highlights impacts of LULC on vegetation cover changes under the changing climate, providing theoretical grounds for ecological conservations over the IM-YRB.

2. Materials and Methods

2.1. Study Region

The IM-YRB (106°34′–112°78′E, 39°09′–41°84′N) is located at the northernmost region of the Yellow River Basin and the Loess Plateau, China, covering an area of approximately 96,000 km² (Figure 1). The IM-YRB is a typical agro-pastoral ecotone with high ecological vulnerability and has been undergoing vegetation restoration since the early 2000s. The IM-YRB is also located at the northwest end of the East Asian Monsoon. It is semi-arid and arid in climate with an average annual precipitation of 297.25 mm. The annual precipitation decreases from the southeast to northwest direction. The croplands in the central IM-YRB rely on water resources from the Yellow River and the ground water for irrigation with two irrigation districts, i.e., Hetao and Tumochuan. The mountains and plateau are dominated by the grasslands and barrens.

2.2. Data

2.2.1. Remotely Sensed Data

To describe the LULC and vegetation changes following the implementation of the ecological conservation policies, we selected data covering the 2000–2020 period in this study. The annual land cover product of China (CLCD) were sourced from https://0-doi-org.brum.beds.ac.uk/10.5281/zenodo.5816591 (Version 1.0.1) [37]. Eight land cover classes were available: cropland, forest, shrub, grassland, water, barren, impervious, and wetland, with a spatial resolution of 30 m. To reduce redundancy, forest and shrub were classified into one land cover category, and the same was conducted for the water and wetland classes (Figure 1a).

NDVI is a vegetation indicator that has been widely used for quantifying the vegetation biomass [18]. We used cloud-free spatial composites with temporal resolution of 16 days and spatial resolution of 500 m (MOD13A1.061) (https://lpdaac.usgs.gov/products/mod13a1v061/, accessed on 26 December 2022). Using the Google earth engine (GEE), the 16 day composite NDVI was processed as monthly average with a spatial resolution of 1 km.

In addition, we obtained SRTM (shuttle recovery topography mission) elevation data (with a spatial resolution of 90 m) from https://www.gscloud.cn/, accessed on 26 December 2022.

2.2.2. Climate Data

We obtained the monthly temperature data from https://0-doi-org.brum.beds.ac.uk/10.5281/zenodo.5111989, accessed on 26 December 2022 and the spatial resolution was determined as 1 km [38]. This dataset was interpolated spatially using the Gaussian process regression (GPR) method and in situ observations. The monthly precipitation data were sourced from https://0-doi-org.brum.beds.ac.uk/10.11922/sciencedb.01607, accessed on 26 December 2022, with a spatial resolution of 1 km [39]. This dataset was interpolated with in situ observed precipitation data collected from over 2400 meteorological stations across China. The spatial interpolation was performed using the anusplin4.4 [40]. The soil moisture data were sourced from http://data.tpdc.ac.cn/zh-hans/data/49b22de9-5d85-44f2-a7d5-a1ccd17086d2/, accessed on 26 December 2022, and the spatial resolution was 1 km [41]. The 10-layer soil moisture data were developed using the random forest model based on the data from the China Meteorological Administration. We processed the time scales of the soil moisture from daily to monthly and the mean soil moisture from 0 to 100 cm in depth. Temperature, precipitation, and soil moisture data covered the same period as the remotely sensed data, from 2000 to 2020, respectively.

2.3. Methods

2.3.1. Analyses for LULC

We evaluated LULC changes of the study region during the 2000–2020 period. Here we showed two typical LULC conversions covering the 2000–2020 period: (1) cropland abandonment/expansion; and (2) grassland restoration from the barren land.

To extract the cropland abandonment grids, we conducted spatial analysis of cropland data for $yea{r}_i$ and $yea{r}_{i+1}$ (i = 2000, …, 2019). We identified cropland abandonment grids where the croplands were converted to other land cover classes, and the abandonment year was $yea{r}_{i+1}$. Similarly, the grassland restoration/cropland expansion grids were extracted where other land cover classes were converted to grasslands/croplands, and the expansion year was $yea{r}_{i+1}$. It is worthy to note that some grids underwent multiple instances of cropland abandonment/expansion. In this study, we only displayed this spatial process based on the last cropland abandonment/expansion over the study region.

In the IM-YRB, conversions of grassland play an important role in regional LULC changes [42]. Hence, the fraction of grassland area in 1 km × 1 km grid was taken to quantify LULC changes during the 2000–2020 period. The grid-based fraction of grassland area can describe the land cover conversions between grassland and other land cover categories, e.g., increased grassland fraction in a grid cell either shows an abandoned cropland or vegetation restoration from barren. This process obtained quantitative time series LULC data with the same spatial resolution as the vegetation data and climate data.

2.3.2. Vegetation Changes

Vegetation dynamics were analyzed by utilizing kNDVI, a nonlinear generalization of the NDVI [16]:

$$\mathrm{kNDVI}=\tanh \left({\mathrm{NDVI}}^2\right)$$

(1)

The kNDVI has recently been proposed as a robust proxy for ecosystem productivity [16]. The kNDVI has been documented to be more closely related to primary productivity, to be resistant to saturation, bias, and complex phenological cycles, and has been evidenced as a strong proxy with an enhanced robustness against noise and stability when compared to alternative products, such as the NDVI [43].

Interannual vegetation changes were characterized for the growing-season-averaged kNDVI. According to the land surface phenology data obtained from the vegetation index and satellite-based phenology product (https://vip.arizona.edu/, accessed on 26 December 2022), the long-term growing season duration is nearly 150 days for crops and is specifically over 210 days for the natural vegetations over this study region. Therefore, the averaged kNDVI during the period from April to October was taken to characterize the interannual vegetation changes during the period of 2000–2020. The subsequent analysis was conducted for the same growing-season periods.

We evaluated the Mann–Kendall trends in interannual kNDVI. The Mann–Kendall trend test is a non-parametric technique that is typically used to detect monotonic trends in time series. In this study, trends at p  <  0.05 were considered to be statistically significant. Based on the slope of the linear model for each grid, we calculated the fitted variation of kNDVI from 2000 to 2020 as shown in Equation (2). These values were then compared to the kNDVI in 2000 to show vegetation changes over this period.

$${\mathrm{kNDVI}}_{p,ij}=\frac{{\mathrm{slope}}_{ij}\times 20}{{\mathrm{kNDVI}}_{2000,ij}}\times 100\%,$$

(2)

where ${\mathrm{kNDVI}}_{p,ij}$ indicates the percentage of kNDVI changes to 2000, ${\mathrm{slope}}_{ij}$ and ${\mathrm{kNDVI}}_{2000,ij}$ are the slope of the linear model from 2000 to 2020 and the kNDVI value in 2000 for each grid, respectively.

2.3.3. Attribution and Sensitivity of the Vegetation Changes

In this study, we attributed kNDVI changes to variations in temperature, precipitation, soil moisture, and LULC during the growing season of the 2000–2020 period. An approach of explainable machine learning [36] was used to evaluate the relationships between the kNDVI and these variables. To eliminate the available long-term trends, we subtracted the long-term mean signals from each individual variable. The growing-season-averaged kNDVI anomaly was treated as the target variable, while the corresponding growing-season-averaged temperature, precipitation, soil moisture, and fraction of grassland anomalies in each 1 km × 1 km grid were used as the predictors.

We trained random forests (RF) models and applied SHapley Additive exPlanations (SHAP) [13,44] to screen out the marginal contributions of each predictor for the target variable. To train the RF model for the core grid cell, we collected all predictor and target data from one grid cell and the surrounding grid cells (3 × 3 shape) during the 2000–2020 period.

We assessed the dynamic vegetation responses to climate and land use changes using vegetation sensitivity to the climate and LULC factors. We defined kNDVI sensitivity as the slope estimated from Theil–Sen regression between the SHAP dependence for kNDVI and the anomalies for each predictor, assuming linearity in grid cell-level interactions between the kNDVI and each predictor [45]. The overall SHAP values of each variable reflect its contribution to the kNDVI dynamics. The final step was to sort the absolute SHAP values to evaluate the overall importance of these four variables to the kNDVI anomalies observed at each grid. The models were implemented in Python 3.10 environment.

3. Results

3.1. LULC Changes of IM-YRB from 2000 to 2020

Based on the CLCD data of 2020, the dominant LULC classes in the study region were grassland, cropland, and barren, accounting for 58.50%, 24.60%, and 11.69% of the total land area, respectively. The impervious, forest, shrub, wetland, and water only covered 5.21% of the region. To demonstrate temporal changes in the LULC, we extracted the percentage of these classes during the 2000–2020 period as a LULC series (Figure 2).

Figure 2 demonstrated an increased grassland and decreased barren that occurred in the 2000s as the prominent LULC variations in the study region. The LULC data for 2000 and 2020 demonstrated that 3571.50 km² of land areas were subjected to conversion from barren to vegetated land, accounting for 24.47% of the barren in 2000 and 3.71% of the entire study area. Meanwhile, cropland, the prime land type in the study region, varied remarkably from 2000 to 2020. The “Grain for Green” strategy initiated from 2000 pushed land conversions from slope cropland to forest and grassland. However, cropland continued to expand slowly and gradually in low-lying regions. Meanwhile, urbanization drove increased impervious areas by >140% from 2000 to 2020 (1.17~2.81%). As a result, cropland accounted for 24.39% of the total area in 2000 and 24.60% in 2020, respectively.

Spatial patterns of grassland restoration and cropland abandonment/expansion during the 2000–2020 period showed land conversions from barren to grassland on the southern side of the Yellow River and the western part of the IM-YRB (Figure 3). Two typical areas were delineated as examples to display spatial details (Figure 3b,e). The hue of the color indicates the year that grassland restoration and cropland abandonment/expansion occurred. The dark color represents the earlier years while the light color represents the newer years. However, this increase in the grassland area was offset by other land conversions, such as cropland expansion, resulting in limited expansions of the grassland.

Variations in the cropland area (Figure 2) can be attributed to the abandonment and expansion of cropland in the study region (Figure 3d–f). Most of the cropland expansion occurred at the periphery of the irrigation districts. Meanwhile, a considerable amount of cropland far from the irrigation districts was abandoned. The difference in the area of cropland abandonment and expansion caused the net increase or decrease in the cropland area. Based on the SRTM data, the mean slope for the cropland expansion area was 2.06° during the 2000–2020 period, while the mean slope for the cropland abandonment area was 3.47°, respectively. This phenomenon suggests a trend of cropland centralization instead of further fragmentation. It reflects the direct result of the “Grain for Green” strategy and may also be related to the soil fertility decline of the slope cropland.

3.2. Spatiotemporal Pattern of Vegetation Changes

Based on the linear changes in the kNDVI during the 2000–2020 period, fitted vegetation changes compared to the kNDVI value in 2000 were presented to show the spatial distribution of temporal vegetation variations (Figure 4 and Figure S1).

We detected a significant overall increase in the kNDVI over the study region during the 2000–2020 period. High changing magnitudes in the kNDVI can be detected in the cropland expansion area around the Hetao irrigation district, the east of the Tumochuan irrigation district, and the cropland abandonment area in the southeast of the study region (Figure S1). Although with a relatively small linear slope (Figure S1), the kNDVI increased significantly compared to the kNDVI in 2000 (Figure 4) in the south part of the study region, showing land conversion from barren or sparse vegetated area to grassland due to vegetation restoration.

Notably, no significant changes in the kNDVI of stable cropland inside the irrigation districts was observed, and this was investigated using the monthly kNDVI changes (Figure S2). Different from natural vegetation, the growing season of the cropland here is from May to September. Although the cropland kNDVI markedly increased in July and August (Figure S2), the average kNDVI during the growing-season periods did not significant increase. As a result, these kNDVI changes imply a shortening growing season of these crops.

3.3. Sensitivity of the kNDVI Dynamics to the Climate and LULC

We identified a higher fraction of areas with a negative kNDVI sensitivity to temperature (71.74%) when compared to those with a positive kNDVI sensitivity to temperature (28.26%) (Figure 5a). Meanwhile, fractions of the area with a positive kNDVI sensitivity to precipitation (96.93%) and SM (89.33%) were much higher than that of negative (Figure 5b,c). Significant sensitivity (p < 0.05) indicates negative impacts of temperature increments on kNDVI dynamics and vice versa, accounting for 81.87%, 91.94%, and 83.45% of the total area in the study region, respectively. Grassland fraction was taken as the LULC variable, and positive and negative sensitivity depends on different LULC conversions. Moreover, 83.35% of the total area showed a significant sensitivity (p < 0.05) of the kNDVI dynamics to grassland changes. The percentage of negative kNDVI sensitivity to LULC (49.62%) was roughly equivalent to that of positive kNDVI sensitivity (50.38%).

The natural vegetation regions in the northern and western semi-arid and arid areas, as well as most croplands, exhibited a negative sensitivity of kNDVI change to temperature. Natural vegetation in dryland is often under drought stress, which is highly related to temperature variation in addition to moisture conditions. The negative sensitivity observed over the irrigated cropland was likely caused by the shortened growing season, as previously described in Section 3.2. However, there are still some regions that exhibited a positive kNDVI sensitivity to temperature, such as the less arid southeast area and some forest and shrubs in the northeast area.

Without incorporating the changing magnitudes, the spatial patterns of kNDVI sensitivity to precipitation and SM were deemed to be similar. A positive sensitivity to precipitation indicates that water resources play a significant role in supporting the kNDVI in this semi-arid and arid region. However, this sensitivity magnitude was relatively small in western arid natural vegetations and was even negative in the Kubuqi Desert (Figure 5b,c), which is a typical barren sand land that underwent grassland restoration. This feature is more prominent over this area and the northwestern area for the sensitivity to SM. The negative kNDVI sensitivity to precipitation and SM suggests that vegetation restoration did not result from the additional moisture supplements over this arid region. Excluding the climate drivers, it could be inferred that vegetation restoration was primarily driven by human behaviors. Meanwhile, human-induced vegetation restoration leading to the decrease in SM was also observed.

Different from temperature, precipitation, and SM, kNDVI sensitivity to LULC is spatially sporadic. Negative kNDVI sensitivity to grassland dynamics primarily corresponded to the land cover conversion from cropland, forest, and shrub to grassland. The kNDVI of these land covers is higher than that of grassland. In the northeastern mountainous area, the negative sensitivity of kNDVI to the increase in forest and shrub was stronger compared to that of the increase in cropland (Figure 5d). On the other hand, the positive kNDVI sensitivity to grassland dynamics was primarily associated with the grassland restoration, including natural vegetation areas in the west and abandoned cropland in the southeast IM-YRB. Additionally, there were also some areas that were observed with a positive sensitivity within the irrigation district, which may be attributed to the expansion of the construction lands.

In addition, we conducted a partial correlation analysis of these four variables with kNDVI dynamics (Figure S3). The distributions of the partial correlation coefficient were basically consistent with the overall kNDVI sensitivity based on the SHAP dependence method, ensuring the credibility of our findings.

3.4. Attribution of kNDVI Dynamics to Climate and LULC

Here, we quantified the fractional contributions of temperature, precipitation, SM, and LULC to kNDVI changes (Figure 6). Figure 6 highlights that precipitation, SM, temperature, and LULC are the first, second, third, and fourth important factors influencing the kNDVI dynamics in this study region, respectively. As the dominant variable, Figure 6a demonstrated that these four variables accounted for 55.94%, 33.14%, 6.44%, and 4.48% of the total area, respectively. Based on the overall SHAP values of these four variables, we found that temperature, precipitation, SM, and LULC contributed 21.54%, 33.32%, 32.19%, and 12.95% to the kNDVI dynamics observed during the 2000–2020 period, respectively.

Vegetation changes across most of the study region were affected by precipitation and SM variations, accounting for 89.08% and 76.14% of the total area as the first and second dominant drivers, respectively (Figure 6a,b). The influence of precipitation on kNDVI changes was widespread in most grassland and cropland. Meanwhile, the major contribution of SM to kNDVI changes was primarily concentrated in the north and east regions (Figure 6a). Interestingly, this northeast area corresponded to the forest and shrubs area in the mountains (Figure 1a), indicating a higher relevance of SM to kNDVI than precipitation over woody vegetation than grassland. The temperature and LULC were dominant in Figure 6c,d, indicating lower contributions of temperature and LULC to kNDVI compared to precipitation and SM over the entire study region.

Notably, even with the typical cropland abandonment/expansion, precipitation and SM still played dominant roles in the kNDVI dynamics over the southeast area. This was due to the small proportion and fragmentation of the cropland abandonment/expansion area each year, as well as the minimal differences in the kNDVI between the cropland and grassland just before and after abandonment. To further analyze this phenomenon, we extracted the kNDVI dynamics and precipitation changes during the 2000–2020 period based on the grids of abandoned cropland, new cropland, stable cropland, and grassland (with area proportion > 50% in 1 km × 1 km grid) in the area of Figure 3e, as shown in Figure 7.

Figure 7 reveals that the kNDVI dynamics of the abandoned cropland, new cropland, stable cropland, and grassland were temporally corresponded with the variation in precipitation, which is consistent with the dominant role of precipitation shown in Figure 6a. Despite the absence of an apparent increasing trend in precipitation, all kNDVIs exhibited an undulating increase trend during the 2000–2020 period. However, the slopes of the linear fits for this kNDVI series differed, indicating the effects of LULC under the primary influence of precipitation change. The slope for stable grassland and stable cropland was 0.0043 (adjusted R² = 0.7056) and 0.0035 (adjusted R² = 0.6131), respectively, suggesting a greening rate of natural vegetation and crop in this area without the influence of LULC. The higher greening rate for stable grassland signifies a gradual thicketization of the grassland. The slope for abandoned cropland was 0.0039 (adjusted R² = 0.6212), which was lower than that of the stable grassland and higher than that of the stable cropland. This was deemed to result from the shorter duration of grassland thicketization after cropland abandonment than stable grassland. The slope for new cropland was 0.0036 (adjusted R² = 0.6481), which was slightly higher than that for stable cropland. This finding was deemed to have been obtained due to both grassland greening before land cover change and subsequent reclamation.

In contrast to the coupled influences of precipitation and SM, evident impacts of the temperature and LULC on kNDVI dynamics concentrated across the croplands in the central study region and grasslands in the western study region (Figure 6a). This corresponded to 6.44% and 18.01% of the total area of kNDVI dynamics being dominantly or second dominantly influenced by temperature, respectively, and were mainly distributed over the irrigation district and Kubuqi Desert. Furthermore, 4.48% and 5.85% of the total area of kNDVI dynamics were dominantly or second dominantly influenced by LULC, respectively, being mainly distributed over the Kubuqi Desert and west areas of the study region. Regions with the dominant role of LULC in kNDVI dynamics were consistent with the distribution of grassland restoration. This was mainly ascribed to the significant increase in the kNDVI from barren to grassland. The role of LULC in kNDVI dynamics was deemed to be more important than precipitation, suggesting that human activities may have a more significant impact on vegetation restoration in these arid areas than natural revegetation.

4. Discussion

4.1. Quantification of LULC in the Attribution of Vegetation Changes

Since the early 2000s, vegetation coverage has significantly increased over the Loess Plateau, China [31,46]. This study observed grassland restoration from barren and cropland in the IM-YRB, inducing a widespread increase in the kNDVI. These changes imply a deep imprint of human activities in vegetation dynamics in the recent 20 years. Here we extracted changes in the grassland fraction at the grid level to evaluate the fractional contributions of the climatic and LULC factors to vegetation dynamics. We found that the kNDVI changes were significantly sensitive to LULC over 83.35% of the total area, and grassland change contributed 12.95% to the kNDVI dynamics in the study region. Therefore, using high-resolution spatial-temporal data to quantify LULC in attribution analysis can help to better evaluate the impacts of human activities on vegetation dynamics.

Different from the indirect residual analysis method for the impacts of human activities on vegetation coverage change [30,34,47,48], our approach directly determines the marginal contribution of the climate factors and LULC to vegetation changes. Meanwhile, we also differentiated the sensitivity of the kNDVI dynamics to LULC and other climatic variables. We found that the overall contribution of LULC to interannual kNDVI changes was 12.95% over this study area, which is lower than the contribution of human activities to vegetation variation in the Loess Plateau [30,48].

4.2. The Sensitivity of Vegetation Changes to Climate Factors

Here we used the SHAP dependence method to analyze the overall sensitivity of the kNDVI changes to climate factors and LULC and found that interannual kNDVI changes in most areas of this study region exhibited a negative sensitivity to temperature but were positively sensitive to precipitation and soil moisture. These findings are consistent with previous studies in drylands [49,50,51]. Meanwhile, a notable 23.14% of the total areas displayed a significantly positive kNDVI sensitivity to temperature. In general, temperature has a negative effect on surface SM under drought conditions [52], leading to negative effects on vegetation growth. However, this effect would be mitigated by additional deep SM for vegetation growth when the root zone is relatively deep, such as the southeast area and some forest and shrubs in the northeast of this study region.

We observed a negative sensitivity of kNDVI changes to SM in the western grassland and barren of the study region. Considering the overall kNDVI growth in this area (Figure 4), this negative sensitivity indicates a decreased SM. Previous studies have discussed the potential impact of vegetation restoration on regional SM deficits [21,53,54]. In the context of vegetation restoration and climate warming, areas with a positive sensitivity to precipitation and a negative sensitivity to SM suggest that additional precipitation may not be sufficient to meet the demands for increased evapotranspiration in the recent 20 years. In comparison, vegetations in areas with a negative sensitivity to both precipitation and SM may be more vulnerable.

4.3. Different kNDVI Increases under Moisture Conditions and LULC

Based on an explainable machine learning method, we evaluated the attribution of temperature, precipitation, SM, and LULC to regional kNDVI dynamics. We showed that the increase in the kNDVI was primarily driven by precipitation changes in this study region during the 2000–2020 period [46,48]. This influence of precipitation on the kNDVI is similar to the structural overshoot of vegetation [12,55]. A structural overshoot indicates that favorable climate conditions in the past could stimulate vegetation growth to surpass the ecosystem carrying capacity, leaving it vulnerable to climate stresses. Despite a transient decline, our study region exhibited an overall increase trend for the kNDVI. For instance, a sudden drought in 2011 and 2015 caused a dramatic decrease in the kNDVI, while the kNDVI increased even higher than that before in 2012 and 2016 (Figure 7). This phenomenon demonstrates the resilience of vegetation and the sufficient land carrying capacity for vegetation growth over the past 20 years. However, this resilience may decline under future climate change [43], or as continued vegetation growth. Studies based on CMIP6 indicate an increasing trend in the precipitation over this region in the decades to come [56,57]. The vegetation may not have reached a climax community in the grassland, forest, or shrubland of the eastern area over the past 20 years. Therefore, we anticipate that the growth of kNDVI will likely continue in the future.

The increase in the kNDVI is affected by LULC, which is typical in areas with cropland abandonment or expansion. In addition to the distinct slope of kNDVI increase in the stable grassland and stable cropland, these findings highlight the influences of cropland abandonment/expansion on the slope of kNDVI increase. Specifically, the slope of kNDVI increase was lower in the cropland abandonment area than that in the stable grassland and was higher than that in the stable cropland. However, the specific value depends on the mean year of cropland abandonment, which we did not quantitatively extract due to the small proportion and fragmentation of the cropland abandonment areas each year. This difficulty was further compounded by the spatial resolution of the kNDVI (1 km × 1 km). Additionally, different abandonment or expansion years introduce uncertainty into the mean slope of kNDVI increase for the corresponding grids in Figure 7.

Based on the results of the different effects of climate and LULC changes on vegetation in this study area, we suggest that human-induced vegetation restoration over the drier western area should pay close attention to changes in SM. Meanwhile, further study should focus on the resilience of vegetation following the sudden decline in precipitation over dryland.

5. Conclusions

With the implementation of ecological projects, vegetation coverage in the IM-YRB, China, has increased over recent decades. Based on the explainable machine learning technique and the SHapley Additive exPlanations (SHAP), we evaluated relationships between the kNDVI and temperature, precipitation, SM, and LULC anomalies. In the model, the fraction of grassland area on a pixel level was selected to quantitatively describe the LULC changes due to the key role of grassland in the LULC changes of this study area. We achieved the following interesting and important findings:

(1)

The overall sensitivity of the kNDVI dynamics to precipitation and SM is positive, while it is negative for temperature, with area fractions of 96.93%, 89.33%, and 71.74%, respectively. However, the opposite sensitivity is also observed in some west grassland and barren areas;

(2)

The fractional contributions of temperature, precipitation, soil moisture, and LULC to kNDVI anomalies are 21.54%, 33.32%, 32.19%, and 12.95%, respectively. As the dominant factor, increasing precipitation implies resilience of the vegetation and sufficient land carrying capacity to vegetation growth over the past 20 years;

(3)

The different slopes of kNDVI increase indicates that the dominant role of precipitation in vegetation dynamics is also affected by LULC, particularly in areas with cropland abandonment/expansion. In the drier western area, human land use became the dominant factor in affecting the kNDVI anomalies. However, we should pay close attention to changes in SM.