Distinguishing Trajectories and Drivers of Vegetated Ecosystems in China's Loess Plateau

Abstract

Terrestrial ecosystems can exhibit various behaviors in response to climate change and human activities. Nonlinear and abrupt shifts in ecosystems are particularly important as they indicate substantial modifications in ecosystem structure and function, posing a threat to the provision of ecosystem services. Here we distinguish between linear, curvilinear, and abrupt shifts in ecosystem productivity from 2000 to 2020 in China's Loess Plateau. We utilize spatial Random Forest models to analyze the driving factors behind these change patterns. Our findings indicate that 84.1% of the ecosystems experienced a positive change in plant productivity, while a small proportion (2.5%) exhibited a negative change. Abrupt shifts are prevalent in both positive and negative changes in ecosystem productivity, with negative changes often manifesting as abrupt shifts (79.3%). Negative changes in plant productivity, particularly the negative abrupt shifts, are primarily associated with human activities characterized by increased nighttime light and urbanization. Land conversion to forest is linked to a curvilinear trajectory in plant productivity, characterized by nonlinear changes with acceleration. Higher water availability and a wetter environment are more likely to promote positive changes in plant productivity. Moderate warming trends contribute to abrupt positive changes in plant productivity, while high warming trends are associated with increased negative abrupt and curvilinear changes. We highlight the importance of accounting for diverse change behaviors in ecosystems for the development of targeted conservation and restoration measures.

Introduction

Terrestrial ecosystems provide numerous benefits to humans, including soil conservation, carbon sequestration, and regulation of water flows (Díaz et al., 2018). Long-term satellite data have revealed a significant greening trend in global vegetated areas since the 1980s, particularly in China and India (C. Chen et al., 2019; Piao et al., 2020). However, human activities and global change are causing pervasive shifts in global vegetation dynamics, which are expected to escalate in the coming decades (McDowell et al., 2020). A recent study demonstrated the prevalence of abrupt vegetation shifts in global drylands, indicating that about 50% of all dryland ecosystems experienced abrupt positive or negative temporal dynamics in terms of gains or losses of Normalized Difference Vegetation Index (NDVI) (Berdugo et al., 2022). Moreover, research on ecosystem resilience has revealed emerging signals of reduced vegetation resilience worldwide, highlighting the vulnerability of global ecosystems to enduring negative shifts when exposed to external disturbances (Boulton et al., 2022; Feng et al., 2021; Forzieri et al., 2022; Smith et al., 2022). Consequently, the concurrent processes of vegetation greening, reduced vegetation resilience, and widespread abrupt shifts raise great uncertainties about the future sustainability of vegetated ecosystems.

Multiple factors contribute to abrupt shifts in vegetated ecosystems. Urbanization, for instance, can rapidly diminish plant productivity by converting vegetated areas into impervious surfaces (van Vliet, 2019; L. Zhang et al., 2022). Given that over half of the world's population lived in urban areas in 2018 (55%), and this number is projected to increase to 68% by 2050 (United Nations, 2019), the era of rapid urbanization is expected to witness a rise in abrupt negative changes in global vegetation ecosystems. Regarding climatic factors, research has demonstrated that higher interannual variability in precipitation is associated with increased vulnerability to abrupt shifts in dryland ecosystems, both in terms of positive and negative changes (Berdugo et al., 2022). Regions characterized by less precipitation, greater precipitation variability, low water availability, and closer to human activity often exhibit low vegetation resilience, increasing their susceptibility to abrupt transitions when faced with external disturbances (Boulton et al., 2022; Forzieri et al., 2022; Smith & Boers, 2023; Verbesselt et al., 2016). Nevertheless, it is important to note that abrupt positive changes in vegetation dynamics can also occur through ecological restoration, revegetation, and land abandonment (Berdugo et al., 2022).

Previous research has primarily focused on linear change patterns in vegetation metrics, which represent smooth and monotonic changes in ecosystems (C. Chen et al., 2019; Hua et al., 2017; Yan et al., 2020). However, simple linear models alone may not be sufficient to detect abrupt shifts in ecosystems (Berdugo et al., 2022). Recent research has begun to address this gap by developing and applying methods for detecting abrupt changes in ecosystems. For instance, the BFAST (Breaks For Additive Season and Trend) package in R (Verbesselt et al., 2010) has been used to detect abrupt changes in remote-sensing-based vegetation indices (e.g., Chen et al., 2022). This package uses iterative estimation to determine the timing and amount of abrupt changes within a time series (Verbesselt et al., 2010). While this approach can identify abrupt changes, it does not assess the persistence of the new state over time. Therefore, there is a need for methods that align more with the definition of abrupt shifts, which involve sudden changes followed by persistence (Berdugo et al., 2022; Ratajczak et al., 2018; Scheffer et al., 2001). To address this, a recent study proposed a multi-model-based trajectory diagnosis approach that distinguishes between linear (linear models), curvilinear (quadratic models), and abrupt change trajectories (step models) of NDVI (Berdugo et al., 2022). In this study, we adapt this approach to a unique Chinese ecosystem, considering the variability in change behavior.

The Loess Plateau, an ecologically vulnerable region in China, has historically experienced severe soil erosion and land degradation (Fu et al., 2022). In response, the Chinese government implemented the “Grain for Green” program in 1999. This program aimed to prevent soil erosion, mitigate flooding, and store carbon by increasing forest and grass cover on hillslopes previously used for agriculture (Bryan et al., 2018). Meanwhile, the Loess Plateau has also experienced rapid urbanization over the last two decades. The urban population fraction has risen from 34.9% in 2000 to 63.8% in 2020 (Fu et al., 2023). Consequently, the combination effects of ecological restoration, urbanization, and global climate change are expected to result in diverse ecosystem dynamics in the region. Despite the importance of these changes, previous studies have primarily focused on linear changes in vegetated ecosystems on the Loess Plateau. Nonlinear and abrupt shifts have received less attention. Here, we applied the trajectory diagnosis method proposed by Berdugo et al. (2022) to differentiate between linear, curvilinear, and abrupt change patterns in ecosystem kNDVI (kernel NDVI, a recently proposed vegetation index; Camps-Valls et al., 2021) from 2000 to 2020 on China's Loess Plateau. We quantified the relative occurrence of each change pattern to identify the spatial distribution of the most prevalent types of kNDVI dynamics across the Loess Plateau. To investigate the associations between kNDVI trajectories and potential social and environmental factors, we employed spatial Random Forest models. We aim to determine whether there is diversity in the change trajectories of vegetated ecosystems on the Loess Plateau, identify the dominant trajectory types, and explore the drivers that contribute to these different change trajectories.

Materials and Methods Data Sets and Preprocessing Vegetation Index

Monthly NDVI data were collected from the Moderate Resolution Imaging Spectroradiometer sensor (MOD13A3 Version 6). This data set covers the period from February 2000 to December 2020 and has a resolution of 1 km. We converted the NDVI values to kNDVI (kernel NDVI) using the equation: kNDVI = tanh(NDVI²) (Camps-Valls et al., 2021). Compared to alternative products like NDVI and near-infrared reflectance of vegetation, kNDVI has been recently proposed as a more reliable vegetation index that closely reflects ecosystem primary productivity (Camps-Valls et al., 2021; Forzieri et al., 2022). It also demonstrates greater resistance to saturation, bias, and phenological cycles, and exhibits enhanced robustness to noise and instability across various spatial and temporal scales (Camps-Valls et al., 2021). Considering these advantages, we selected kNDVI as the preferred metric to analyze vegetation dynamics in the Loess Plateau. Finally, we aggregated the monthly kNDVI values into annual values using the maximum value synthesis method (Liu et al., 2022).

Socioeconomic Data

Data on the total, urban, and rural population at the county level for the years 2000 and 2020 were collected from the fifth and seventh national population censuses in China (H. Wu, Shi, et al., 2021). To assess changes in population pressure within counties, we calculated the difference in total population between 2020 and 2000 for each county. Additionally, we calculated the difference in urban population fraction, which is the ratio of urban population to total population, between 2020 and 2000 at the county level. This difference reflects the changes in population pressure specifically between urban and rural areas, indicating instances of urban-rural population migration within each county. A positive value indicates an increase in urban population fraction, signifying a higher number of individuals moving to cities and thereby reducing population pressure in rural areas.

We obtained annual night light data with a 1 km resolution from HARVARD Dataverse (Y. Wu et al., 2022) for the years 2000–2020. This data set provides improved time-series data similar to Defense Meteorological Satellite Program-Operational Linescan System (DMSP-OLS) (1992–2021) by integrating DMSP-OLS and SNPP-VIIRS (the Suomi National Polar-orbiting Partnership satellite platform-Visible Infrared Imaging Radiometer Suite) (Y. Wu et al., 2022). To analyze the data on a county level, we calculated the average night light values for each county, transforming the grid-scale data into county-scale data. Next, we assessed changes in night light by calculating the difference between the night light values in 2020 and 2000.

Land Use and Afforestation Data

Land-use data for the years 2000 and 2020 was obtained from GlobeLand30 land-cover product (Chen et al., 2014). This data set has a 30 m resolution and provides land classification accuracy of approximately 85.72% (Kappa = 0.82). The land-use data was reclassified into five types: cultivated land, forest (including shrubland), grassland, water bodies (including wetland), artificial surfaces, and others (such as bare land and permanent snow and ice). To analyze land use conversions, we calculated the proportion of land transformed into forest, grassland, agriculture, and artificial surfaces from 2000 to 2020 relative to county areas.

To assess the ecological restoration efforts at the county level, we collected data on the average annual afforestation area percentage from 2002 to 2019 (Y. Chen et al., 2021). This data provides information on the extent of ecological restoration activities conducted in each county during the study period.

Climate Data

Monthly data on precipitation, temperature, and potential evapotranspiration (PET) were obtained from the National Earth System Science Data Center and the National Science & Technology Infrastructure of China (Peng et al., 2019). These data cover the period from January 2000 to December 2020 and have a resolution of 1 km. To examine inter-annual variability, the monthly-scale data was aggregated into yearly-scale data. Furthermore, data on palmer drought severity index (PDSI) and downward surface shortwave radiation (DSSR) were acquired from TerraClimate (Abatzoglou et al., 2018). These data cover the years 2000–2020 and have a resolution of approximately 4,500 m. The PDSI is a drought index that monitors long-term drought patterns using the supply-and-demand concept of the water balance equation (H. Yu et al., 2019). On the other hand, DSSR represents the radiative energy reaching the Earth's surface per unit of time and surface area in the wavelength range of [0.3 μm, 3.0 µm], which is essential for plant photosynthesis (Letu et al., 2020). To analyze climate trends at the county scale, we calculated the average values of climate indicators for each county, transforming the data from a grid-scale to a county-scale format. Linear regression slopes between the climate indicators and time (years) were then calculated. Additionally, we assessed climate variability by calculating the variation of precipitation and temperature using the coefficient of variation (Forzieri et al., 2022).

To analyze the occurrence of drought events, we collected a high-resolution data set of the Standardized Precipitation Evapotranspiration Index (SPEI) from the Science Data Bank of China (Xia et al., 2023). This data set covers the entire Chinese Mainland from 2001 to 2020 and provides SPEI data at multiple time scales (1 month, 3 months, 6 months, 12 months, 24 months) with a resolution of 1 km. For our study, we specifically selected the 3-month-scale monthly SPEI to examine drought characteristics in the Loess Plateau. To identify drought conditions using the SPEI drought index, a threshold value of −1 is commonly employed (Haile et al., 2020; Yao et al., 2023). We calculated the drought frequency by dividing the number of drought months by the total number of months (Haile et al., 2020). Additionally, we determined the drought intensity by summing the absolute values of the SPEI for months with SPEI values below −1, and then dividing this sum by the number of drought months (Haile et al., 2020).

Soil and Topography Data

Soil sand content, nitrogen content, and organic carbon content data were obtained from ISRIC (International Soil Reference and Information Centre)—world soil information, with a resolution of 1 km. These soil indicators consist of measurements at six standard depths. To simplify the analysis, we aggregated each soil indicator into a single raster by calculating the mean. Subsequently, we computed the average soil sand content, nitrogen content, and organic carbon content for each county, thereby transforming the grid-scale data into county-scale data. Annual soil moisture data from 2000 to 2020, with a resolution of approximately 4,500 m, were acquired from TerraClimate (Abatzoglou et al., 2018). We calculated the average values of soil moisture for each county and then calculated the linear regression slopes between the soil moisture and time (years).

We obtained digital elevation model (DEM) data with a resolution of 90 m from the National Earth System Science Data Center and the National Science & Technology Infrastructure of China. Using this DEM data, we calculated the slope. Similarly, we computed the average elevation and slope for each county, thereby converting the data from a grid-scale to a county-scale format.

Determining the Trends and the Shapes of kNDVI Trajectories

We chose kNDVI as a surrogate of ecosystem productivity. Only sites with complete kNDVI time series from 2000 to 2020 (i.e., no missing values) were retained for trajectory analysis. To determine the shape and type of kNDVI trajectories, we utilized the trajectory diagnosis method proposed by Berdugo et al. (2022). This classification categorized kNDVI trajectories into three shapes: linear (monotonic trends or no trends), curvilinear (nonlinear trajectories with acceleration or deceleration), and abrupt shift (discontinuous trends with sudden changes followed by persistence) (Berdugo et al., 2022). This classification process was performed in R environment (R Core Team, 2021). The main technical processes are briefly introduced below, see Berdugo et al. (2022) for more comprehensive methods.

To identify different change patterns of kNDVI, various models were employed: linear models for linear shapes, quadratic models for curvilinear shapes and smooth trend/slope changes, and step models for abrupt shapes (Berdugo et al., 2022). Prior to fitting these models to kNDVI and time (year), we normalized the kNDVI data by subtracting the mean value of the kNDVI time series at each location and dividing it by the standard deviation, enabling the comparison of coefficients across models. Linear and quadratic models were developed using the “lm” function in R, while step models were developed using the “chngpt” package (Fong et al., 2017). Additionally, we fitted a linear model without a slope parameter to examine the fitting of a no-trend typology, where the kNDVI trajectory follows a linear shape with only an intercept and no slope (Berdugo et al., 2022).

To select the best-fitted model for each trajectory, we compared the akaike information criteria values of each fit (Berdugo et al., 2022). To decrease uncertainty in trajectory classification, we employed the bootstrap method. Each regression was bootstrapped 100 times, and each bootstrapped iteration between models was compared (Berdugo et al., 2022). It is important to note that step fits with change points occurring within the first or last 3 years of the trajectory were excluded to avoid anomalous points at the extremes (Berdugo et al., 2022). Instead, they were fitted using other more suitable models (linear or curvilinear). Finally, we grouped trajectory shapes (linear, curvilinear, or abrupt) into positive (increasing trend), negative (decreasing trend), or neutral (no trend) types (Berdugo et al., 2022). Linear trajectories were considered neutral if the best-fit model was the no-trend model. Abrupt trajectories were grouped into positive or negative based on the step parameter. Curvilinear trajectories were grouped into positive or negative based on the combination of “a” and “b” parameters of their corresponding quadratic regression (the model form: y ∼ ax² + bx + c) (Berdugo et al., 2022).

Finally, we identified eight kNDVI trajectories: positive linear, positive curvilinear, positive step, negative linear, negative curvilinear, negative step, neutral linear, and neutral curvilinear. Using this information, we performed basic statistical analysis to assess the predominant types of kNDVI trajectories in the Loess Plateau, their distribution across different land-use types, and the occurrence frequency of positive and negative abrupt shifts in kNDVI throughout the study period.

Identifying Potential Drivers of kNDVI Trajectories

To analyze the factors influencing kNDVI trajectory trends and shapes, we employed spatial Random Forest models for six significant kNDVI trajectories: negative linear, negative curvilinear, negative abrupt, positive linear, positive curvilinear, and positive abrupt. The explained variables were the area proportions of each kNDVI trajectory relative to county areas. Additionally, we considered 21 potential social and environmental factors as explanatory variables (Table S1 in Supporting Information S1). The spatial Random Forest models were built using “spatialRF” and “ranger” packages in R (Benito, 2021; Wright & Ziegler, 2017).

Spatial Random Forest models are advantageous as they minimize spatial autocorrelation in model residuals and provide reliable variable importance scores by incorporating spatial predictors that capture the spatial structure of the training data (Benito, 2021). Because we built the model for explanatory rather than predictive purposes, we utilized all available samples for model training. The spatial predictors were generated based on the distance matrix of data points, which we obtained by extracting the latitude and longitude of each county area unit's centroid. Distances were calculated using the distm() function from the “geosphere” package (Hijmans, 2022). Distance thresholds of 0, 100, 200, 400, and 800 km were set based on the frequency distribution of the distance matrix (Figure S1 in Supporting Information S1).

Prior to model construction, several preparatory steps were performed. First, we assessed the spatial autocorrelation of the response variable and the predictors across various distance thresholds using the plot_training_df_moran() function (Benito, 2021). The results indicated a strong spatial autocorrelation, especially for climate factors within the 100–400 km range (Figure S2 in Supporting Information S1), underscoring the need for spatial Random Forest models. Second, we examined potential variable interactions among the most important predictors using the_feature_engineer() function (Benito, 2021). In this function, we set the parameter importance.threshold as 0.5 (uses 50% best predictors) and cor.threshold as 0.6 (max correlation between interactions and predictors). The tested variable interactions were then added to predictors. Third, we addressed multicollinearity in the predictors by utilizing the auto_cor() and auto_vif() functions (Benito, 2021), which help identify and reduce redundancy in the predictors based on criteria such as bivariate R-squared and variance inflation factor. Only predictors that passed the multicollinearity test were retained for subsequent model analysis.

The spatial Random Forest models were constructed using the rf_spatial() function, with the default “mem.moran.sequential” method applied for building, ranking, and selecting spatial predictors (Benito, 2021). As the initially set random forest hyperparameters may not be optimal for the data sets in this study, we performed hyperparameter tuning using the “rf_tuning()” function to select sensible values for three critical hyperparameters (Benito, 2021): the number of regression trees in the forest, the number of variables considered at each tree split, and the minimum number of cases required on a terminal node. To measure model uncertainty, we executed each model 30 times using the “rf_repeat()” function (Benito, 2021). Model performance was assessed through measures such as the R-squared and root mean squared error, obtained by predicting the out-of-bag data (fraction of data not used to train individual trees) (Benito, 2021). Variable importance, indicating the increase in mean error (computed on the out-of-bag data) when a predictor is permuted, was also evaluated (Benito, 2021). Moreover, response curves were analyzed to understand the relationships between dependent variables and predictors.

Results Trajectories of Ecosystem Productivity

The majority of areas in the Loess Plateau have exhibited positive changes in kNDVI, with 84.1% showing an increase in plant productivity (Figures 1a and 1b). Among these positive changes, linear trajectories accounted for nearly half (49.4%), suggesting a consistent and steady increase in vegetation productivity. It is worth noting that abrupt changes are also prevalent, accounting for 32.5% of the positive changes. In contrast, only a small proportion showed negative changes in kNDVI, accounting for 2.5% of the Loess Plateau. Negative changes were more frequently characterized by abrupt trajectories (discontinuous trends with sudden changes followed by persistence), accounting for 79.3% of these cases (Figure 1b). The remaining 13.4% of the Loess Plateau displayed neutral changes in kNDVI, indicating a relatively stable level of plant productivity.

Figure 1. Trajectories of ecosystem productivity (measured by kNDVI) in China's Loess Plateau. (a) Map of classified annual trajectories of kNDVI based on their trend and shape between 2000 and 2020. (b) A pie plot displaying the percentage of neutral, positive, and negative kNDVI trajectory trends. The ring surrounding the pie plot indicates the percentage of linear, curvilinear, and abrupt trajectory shapes for each trajectory trend. Panel (c) shows an example time series of each classified trajectory (y-axis: standardized kNDVI values; x-axis: years, 2000–2020, dashed line indicates fitted trajectory).

We conducted spatial analysis to examine the correspondence between kNDVI trajectories and land use patterns (Figures 2a and 2b). The findings revealed that forest ecosystems commonly exhibited positive curvilinear trajectories, characterized by nonlinear changes with acceleration. In comparison, grass and cultivated ecosystems showed relatively fewer instances of this trajectory pattern. Positive abrupt trajectories were observed across grass, forest, and cultivated ecosystems, with no significant differences among them. The occurrence of positive abrupt changes in kNDVI was mainly concentrated in the first decade, reaching its peak in 2011 (Figure 2c). Conversely, negative abrupt trajectories were predominantly found in artificial and cultivated ecosystems, suggesting a potential association with human activities. Notably, there was an annual increase in abrupt negative changes in kNDVI until 2010 (Figure 2d), which could be linked to the rapid urbanization and expansion of urban areas in the Loess Plateau over the past two decades.

Figure 2. Correspondence between kNDVI trajectories and land use patterns (a, b) and the frequency of abrupt changes in kNDVI by year (c, d). (a) Map of land use types in 2020 (Data source: Globalland30, http://www.globallandcover.com/). (b) Area proportions of kNDVI trajectories within each land use type. Histogram of the number of abrupt changes in a given year for positive (c) and negative (d) abrupt changes.

Drivers of Ecosystem Productivity Dynamics

To understand the factors that influence different change trajectories in kNDVI (Figure 3), we employed spatial Random Forest models. Our results showed that three models (negative linear, negative curvilinear, and positive abrupt) had R-squared values exceeding 40%, while the other three models (negative abrupt, positive linear, and positive curvilinear) had R-squared values surpassing 60% (Table S2 in Supporting Information S1). By examining the variable importance and response curves (Figures 4–7), we identified important social and environmental factors that explain the observed change trajectories in kNDVI.

Figure 3. Area proportions of kNDVI trajectories relative to county areas. Negative changes including negative linear (a), negative curvilinear (b), and negative abrupt (c). Positive changes including positive linear (d), positive curvilinear (e), and positive abrupt (f).

Figure 4. The distribution of importance scores of the predictors. Each spatial Random Forest model is repeated 30 times using the “rf_repeat()” function and the distribution of importance scores of the predictors across executions is plotted in violin plots using the “plot_importance()” function. Both functions are from R package spatiaRF (Benito, 2021). Panels (a) to (f) are the results of model for negative linear trajectory, negative curvilinear trajectory, negative abrupt trajectory, positive linear trajectory, positive curvilinear trajectory, and positive abrupt trajectory, respectively. “PET..pca..Land conversion to settlement” and “Drought frequency..x..Soil nitrogen content” in panel (d), digital elevation model..x..Drought frequency in panels (e) and (f) are the retained variable interactions through testing all possible interactions among the most important predictors using the_feature_engineer() function in R package spatiaRF (Benito, 2021). Interactions computed via multiplication are named a..x..b, while interactions computed via PCA are named a..pca..b, where a and b are variable names. PET, potential evapotranspiration; DEM, digital elevation model; PCA, principal component analysis. These explanations of variable interactions apply to Figures 5–7.

Figure 5. The response curves of models for negative linear (a–h) and positive linear (i–p). Each spatial Random Forest model is repeated 30 times using the “rf_repeat()” function and the response curves of models are plotted with “plot_response_curves()” function. The median prediction is shown with a thicker line. Both functions are from R package spatiaRF (Benito, 2021).

Figure 6. The response curves of models for negative curvilinear (a–h) and positive curvilinear (i–p). Each spatial Random Forest model is repeated 30 times using the “rf_repeat()” function and the response curves of models are plotted with “plot_response_curves()” function. The median prediction is shown with a thicker line. Both functions are from R package spatiaRF (Benito, 2021).

Figure 7. The response curves of models for negative abrupt (a–h) and positive abrupt (i–p). Each spatial Random Forest model is repeated 30 times using the “rf_repeat()” function and the response curves of models are plotted with “plot_response_curves()” function. The median prediction is shown with a thicker line. Both functions are from R package spatial RF (Benito, 2021).

Change in nighttime light was the most important variable in explaining three negative trajectories, particularly the negative abrupt trajectory (Figures 4a–4c). These negative trajectories demonstrated an increase in occurrence with higher levels of nighttime light (Figures 5a, 6a, and 7a). Additionally, we observed that counties with high urbanization rates exhibited a greater prevalence of negative abrupt changes in kNDVI (Figure 7h). These consistent findings highlight the prominent role of human activities in driving negative changes in plant productivity, particularly when such changes occur abruptly.

We found that the linear trajectory of vegetated ecosystems was linked to land conversion to grass. As land conversion to grass increased, there was a corresponding increase in positive linear changes in kNDVI (Figure 5j), accompanied by a decrease in negative linear changes (Figure 5d). In comparison, the curvilinear trajectory, characterized by nonlinearity with acceleration, was associated with land conversion to forest (Figures 6g and 6i). The land conversion to forest was the most important variable in explaining positive curvilinear changes in kNDVI (Figure 4e). This observation was supported by the spatial correspondence between kNDVI trajectories and land use patterns depicted in Figure 2b, where forest ecosystems exhibited relatively higher proportions of positive curvilinear trajectories compared to grass and cultivated ecosystems. However, it is important to note that the relationship between land conversion to forest and curvilinear trajectory is not simply linear. Increasing land conversion to forest initially led to a rapid rise in the occurrence of positive curvilinear trajectories, which stabilized soon thereafter (Figure 6i). On the other hand, beyond a certain point, higher levels of land conversion to forest were associated with an increase in negative curvilinear trajectories (Figure 6g).

A higher increasing trend in precipitation was associated with increased positive changes (Figures 5o, 6l, and 7p) and decreased negative changes in plant productivity (Figures 5e and 7c). We observed a nonlinear relationship between precipitation variability and the curvilinear and abrupt trajectories of kNDVI. When precipitation variability was low, there was an increase in positive curvilinear changes in kNDVI (Figure 6n). Conversely, high precipitation variability was associated with increased negative curvilinear (Figure 6h) and decreased positive abrupt changes in kNDVI (Figure 7m). Moreover, we found that a higher increasing trend in PDSI (indicating wetter conditions) tended to be associated with increased positive changes in plant productivity (Figures 5n, 6p, and 7o). As the trend in soil moisture increased, there is an increase in positive abrupt changes in kNDVI (Figure 7n), while negative linear and abrupt changes decreased (Figures 5c and 7f). In comparision, higher drought frequency and intensity were associated with an increase in negative linear changes in kNDVI (Figures 5g and 5h). Higher drought frequency was also related to increased negative curvilinear changes in kNDVI (Figure 6c), while reducing positive curvilinear changes (Figure 6o). These concise and consistent relationships indicate that higher water availability and wetter environments have a higher likelihood of promoting positive changes in plant productivity.

We identified nonlinear relationships between the trend in temperature and the curvilinear and abrupt change behavior of plant productivity. The trend in temperature emerged as the most important variable in explaining positive abrupt changes in kNDVI (Figure 4f). Notably, there was a rapid increase in positive abrupt changes within the moderate warming range of approximately 0.01°C per year (Figure 7i). However, beyond this threshold, a high increasing trend in temperature corresponded to an increase in negative abrupt changes in kNDVI (Figure 7g). Similarly, the negative curvilinear changes initially decreased, but after reaching the inflection point at around 0.01°C per year, they started to increase alongside an escalating trend in temperature (Figure 6f).

Discussion

Our study investigated patterns of change in vegetated ecosystems of the Loess Plateau and analyzed the factors driving these changes using spatial Random Forest models. The results indicate that plant productivity in the majority of areas in the Loess Plateau has showed positive changes over the past two decades (84.1% of Loess Plateau; Figures 1a and 1b), consistent with recent studies documenting a greening trend in the region (Naeem et al., 2021; X. Wu et al., 2019; Y. Yu et al., 2020). Furthermore, our analysis goes beyond the limitations of solely examining linear change trajectories, allowing us to identify diverse change patterns in the Loess Plateau's vegetated ecosystems (C. Chen et al., 2019; Hua et al., 2017; Yan et al., 2020). We find that abrupt trajectories are prevalent in both positive and negative changes in plant productivity, with negative changes often manifesting as abrupt shifts (Figure 1b). These non-linear and abrupt shifts are of particular interest to scientists and policymakers as they often indicate fundamental alterations in vegetation structure and function, posing a threat to the ecosystem services crucial for human well-being (Berdugo et al., 2022; Ratajczak et al., 2018; Scheffer et al., 2001). Therefore, our results emphasize the need to consider multiple change trajectories in vegetated ecosystems in future research.

We found that the change in nighttime light is the most important variable explaining three negative trajectories (Figures 4a–4c). As urbanization rates increase, negative abrupt changes in kNDVI also increase (Figure 7h). These results suggest that negative changes in vegetation productivity in the Loess Plateau are primarily related to urbanization and settlement expansion. Considering China's rapid urbanization (M. Chen et al., 2016), decision-makers and land managers need to account for the diverse impacts that urbanization poses on vegetated ecosystems (Ren et al., 2022; van Vliet, 2019). It is important to note that, in addition to direct negative effects such as vegetation removal, urbanization can also have indirect positive effects on vegetation growth. For example, a recent study demonstrated a widespread positive indirect effect of urbanization on vegetation growth in 672 cities worldwide, attributed to human management practices and higher air temperatures in urban areas compared to surrounding natural areas (L. Zhang et al., 2022). This implies that, through sound urban planning and management in the Loess Plateau, such as the construction of green infrastructure and spaces, the indirect positive effect of urbanization on vegetation growth can partially offset the direct loss of vegetation caused by land conversion to urban areas (L. Zhang et al., 2022).

Previous studies have reported the contribution of ecological restoration programs, like the Grain for Green project, to vegetation greening in the Loess Plateau (T. Chen et al., 2022; He et al., 2022; Song et al., 2022). However, our results indicate that the afforestation degree at the county level does not effectively predict the different trajectories of vegetated ecosystems (Figure 4). This result is contrary to the expectation that counties with higher afforestation degrees would correspond to greater gains in positive curvilinear and abrupt shifts in plant productivity. The discrepancy might be attributed to statistical biases associated with county-level afforestation data and the failure to consider tree survival rates (He et al., 2022; F. Wang et al., 2020).

In contrast, our research reveals a significant relationship between land use conversions and changes in vegetated ecosystems. Negative abrupt trajectories in kNDVI were mainly observed around artificial surfaces, corresponding to situations where vegetation was removed for urban landscapes (Ren et al., 2022; van Vliet, 2019; L. Zhang et al., 2022). Linear trajectories in kNDVI were associated with grassland conversions, while curvilinear trajectories primarily corresponded to forest conversions. This result aligns with expectations, as ecosystem productivity is more likely to exhibit a positive curvilinear trajectory, characterized by nonlinear changes with acceleration, when forests are planted on barren slopes and natural regrowth occurs on abandoned land through subsequent successional processes (Khorchani et al., 2021; F. X. Wang et al., 2007). However, we observed a nonlinear relationship between forest conversions and the curvilinear trajectory in kNDVI, which may be due to variations in tree survival rates, ecosystem carrying capacity, or tree management efforts (Fu et al., 2022; He et al., 2022; F. Wang et al., 2020). These findings carry important implications as they highlight the variability in the response of vegetated ecosystems to human-induced landscape alterations. Furthermore, they can serve as a management guide for optimizing land use patterns and enhancing the success of vegetation restoration efforts (Berdugo et al., 2022).

We found that higher water availability and a wetter environment promote positive changes in plant productivity. In contrast, increased precipitation variability and more frequent droughts are associated with negative curvilinear changes in kNDVI. This could result from reduced vegetation resilience due to climate variability (Forzieri et al., 2022; Smith & Boers, 2023). The negative feedback loop between low water availability, drought conditions, and vegetation degradation can accelerate the decline in vegetation productivity, hindering its capacity to recover from disturbances (Forzieri et al., 2022; Smith & Boers, 2023; Zemp et al., 2017). The question of whether Northwest China is undergoing a transition from a warm-dry climate to a warm-wet climate remains controversial and scientifically significant (Shi et al., 2007; L. Yang et al., 2022). Existing studies have demonstrated that Northwest China has undergone warming and wetter climate changes, with projections indicating that this trend will continue until the end of the 21st century under a medium emission scenario (J. Yang et al., 2021; L. Yang et al., 2022). Our results show that moderate levels of warming may contribute to positive abrupt shifts in ecosystem productivity. However, there is considerable uncertainty regarding the sustainability of vegetation growth in response to warming, particularly given the limited water resources of the Loess Plateau (Liang et al., 2018; S. Zhang et al., 2018). Notably, we observed that excessive warming trends corresponded to amplified negative abrupt shifts and accelerated adverse changes in ecosystem productivity. These results emphasize the necessity for vegetation management measures tailored to assess the challenges posed by climate change, contributing to the long-term preservation of ecosystems and the services they provide.

This study comprehensively quantified the diverse trajectories of ecosystem productivity in China's Loess Plateau and analyzed the potential driving forces behind these patterns of change. By doing so, it enhances our understanding of the complex dynamics within the region's ecosystems and provides valuable insights for future ecological conservation and sustainable development efforts. However, there are uncertainties that need to be considered. To examine the relationship between trajectories of ecosystem productivity and social and environmental factors, we utilized a spatial random forest model. This approach effectively controls for the effects of spatial autocorrelation, multicollinearity, and variable interactions on the results (Benito, 2021). Although we have largely controlled for spatial autocorrelation in the variables, it still exists in some models, particularly those displaying a positive curvilinear trajectory (Figure S3 in Supporting Information S1). Additionally, we aimed to explain the observed trajectories of ecosystem productivity by collecting as many relevant social and environmental factors as possible. Despite these efforts, the amount of explained variance in the three models (negative linear, negative curvilinear, and positive abrupt) remains relatively low. This suggests that these trajectories may be better explained by site-specific characteristics or the occurrence of specific events, rather than being driven by the overall dynamics of social and environmental factors (Berdugo et al., 2022).

Conclusions

In China's Loess Plateau, we investigate patterns of change in vegetated ecosystems from 2000 to 2020 and analyze the driving factors behind these change patterns using spatial Random Forest models. Results show that a significant majority of the ecosystems (84.1%) experienced positive changes in plant productivity, while a small fraction (2.5%) witnessed negative changes. Abrupt shifts are prevalent in both positive and negative changes in ecosystem productivity, with negative changes often manifesting as abrupt shifts (79.3%). Notably, negative changes, particularly negative abrupt shifts, are mainly associated with increased nighttime light and urbanization. Furthermore, we observed a strong association between kNDVI trajectories and land use conversions, wherein linear trajectories corresponded to grassland conversion, and curvilinear trajectories were linked to forest conversion. Increased water availability and a wetter environment were found to promote positive changes in kNDVI. Additionally, moderate warming trends contributed to abrupt positive changes in plant productivity, while high warming trends were associated with the amplification of negative abrupt and curvilinear changes. By considering a wider range of trajectory types, including curvilinear and abrupt changes, our findings provide a more accurate depiction of the diverse change behaviors exhibited by vegetation ecosystems within the Loess Plateau. Consequently, we stress the importance of accounting for these diverse change behaviors when devising conservation and restoration strategies.

Acknowledgments

This study was supported by the National Natural Science Foundation of China (41930649).

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

Data Availability Statement

Monthly NDVI (Normalized Difference Vegetation Index) data with a resolution of 1 km from February 2000 to December 2020 were downloaded from the MOD13A3 Version 6 (https://lpdaac.usgs.gov/products/mod13a3v006/). Population data at the county level can be obtained from Science Data Bank (H. Wu, Shi, et al., 2021). Land use/cover data with a 30 m resolution in the year 2000 and 2020 was obtained from GlobeLand30 land-cover product (https://www.webmap.cn/mapDataAction.do?method=globalLandCover). Annual nighttime light data from 1992 to 2022 was downloaded from Harvard Dataverse (Y. Wu, Shi, et al., 2021). Monthly precipitation, temperature, and potential evapotranspiration data with a resolution of 1 km from January 2000 to December 2020 were acquired from the National Earth System Science Data Center and the National Science & Technology Infrastructure of China: Precipitation (Peng, 2023b), temperature (Peng, 2023a), potential evapotranspiration (Peng, 2023c). Annual palmer drought severity index, downward surface shortwave radiation, and soil moisture with a resolution of about 4,500 m from 2000 to 2020 were obtained from TerraClimate (https://developers.google.com/earth-engine/datasets/catalog/IDAHO_EPSCOR_TERRACLIMATE#bands). High-resolution SPEI dataset in mainland China from 2001 to 2020 was downloaded from Science Data Bank (Xia et al., 2023). Soil sand content, soil nitrogen content, and soil organic carbon content with a resolution of 1 km were downloaded from ISRIC world soil information (https://files.isric.org/soilgrids/latest/data_aggregated/1000m/). DEM was obtained from Loess Plateau SubCenter, National Earth System Science Data Center, National Science & Technology Infrastructure of China (Q. Yang, 2017). Code for determining the trends and the shapes of kNDVI trajectories can be found in figshare (Berdugo & Dakos, 2022). The spatial Random Forest models were built using “spatialRF” and “ranger” packages in R (Benito, 2021; Wright & Ziegler, 2017). Code for data analysis and figures in this study can be found in zenodo (Z. Wang, 2023).