Trends, turning points, and driving forces of desertification in global arid land based on the segmental trend method and SHAP model

ABSTRACT

Desertification is a form of land degradation observed in arid, semiarid, and dry subhumid ecosystems. Assessing the global trends and drivers of desertification in arid land is crucial for developing effective land restoration policies and mitigating desertification. We aimed to evaluate global segmental trends in desertification in arid land using Moderate Resolution Imaging Spectroradiometer images from 2000 to 2022. By constructing a robust MSAVI-Albedo Desertification Distance Index (DDI), we assessed the segmented development characteristics of desertification on Google Earth Engine. Additionally, we employed the SHapley Additive exPlanations (SHAP) model and machine learning methods to analyze the individual and interactive driving mechanisms of desertification. The results indicated an overall reduction in desertification, with approximately 23.21 million km² (39% of the global arid region area) exhibiting negative trends in DDI. Approximately 31% of the land area showed a DDI of −0.002/y. Precipitation was consistently the primary factor influencing desertification, with an average SHAP value of 11.42. Secondary influencing factors included potential evapotranspiration, soil moisture, and vapor pressure differences. Notably, the coupling between precipitation and soil moisture exhibited the most significant impact on the desertification process, with SHAP coupling values of approximately 3.28 and 5.06 before and after the DDI turning point, respectively. These findings provide new insights into desertification on a global scale and offer valuable scientific support for promoting effective prevention and control measures.

KEYWORDS:

• Desertification
• Turning point
• SHAP
• Google Earth Engine
• Arid land

1. Introduction

Desertification is an intricate process of land degradation triggered by the interplay between climate change and human activity in arid lands. The influence of desertification on the sustainable development of human society is profound and represents one of the most pressing social, economic, and ecological predicaments confronting contemporary society (D’Odorico et al. 2013; Huang et al. 2020). Desertification has caused substantial environmental deterioration and economic damage worldwide. Desertification affects approximately one-fifth of the world’s population, with arid lands experiencing particularly pronounced effects (Meng et al. 2021). Intensified climate change and increasingly frequent human activities have amplified the global trajectory of desertification in arid land (Ma et al. 2021), posing imminent threats to human existence, productivity, and livelihoods. Ascertaining the global distribution characteristics and development laws of desertification is the foundation for scientifically formulating desertification control policies and promoting desertification control methods. However, the understanding of the global spatiotemporal evolutionary characteristics and underlying driving mechanisms of desertification in arid land are not comprehensive, impeding the effective advancement of desertification control measures.

Accurate identification of the precise location and extent of desertification and a systematic assessment of its driving mechanisms are prerequisites for effective desertification control. Regardless of the influence of technological conditions or the urgency of research, most of the literature has explored areas prone to desertification and those severely affected, for example, Western Asia (Akbari, Modarres, and Alizadeh Noughani 2020), Central Asia (Hu, Han, and Zhang 2020), North Africa (Yang et al. 2022), and the northern regions of East Asia (Meng et al. 2020). Global-scale assessments of desertification are crucial for understanding the global development process of desertification and formulating international policies and cooperation strategies for desertification prevention but few have been conducted. Huang et al. (2020) proposed a global desertification vulnerability index that combined climate change and human activities, providing a new perspective for the study of global desertification vulnerability. Despite Huang et al. (2020) offering a global perspective to assess the vulnerability to desertification, global assessments of desertification based on remote sensing image data remain limited. Remote sensing technology can reflect the status of surface desertification accurately and in real time, with minimal impact from desertification assessment models, which is especially critical for monitoring and assessing the global desertification process. Additionally, these desertification monitoring studies significantly enriched the field of global desertification monitoring by introducing various methodologies. Their insights have provided solid scientific groundwork essential for efforts to prevent and control desertification. The occurrence and development of land desertification typically follow a complex, non-monotonic pattern. Simply using linear trend methods to assess the evolution of desertification over a study period makes uncovering detailed changes difficult. Therefore, for subjects with long time series, inflection point detection methods must be used to determine the sudden change times of desertification and to conduct segmental linear trend analyses. However, some of these desertification monitoring studies had limited monitoring periods or primarily relied on linear analysis for trend assessment without accounting for segmental characteristics. Despite several studies on the analysis of the interactive effects of desertification driving factors (Lamqadem et al. 2018, 2018b; Sepehr et al. 2007), the interpretation of factor interactions has rarely been quantitatively described, which limits a comprehensive understanding of the driving mechanisms of desertification.

Data-driven machine learning models have been widely used to explore and construct nonlinear relationships between environmental variables and target variables (Huang, Duan, and Chen 2022; Meng et al. 2021). Among these models, eXtreme Gradient Boosting (XGBoost), LightGBM (LGBM), Random Forest (RF), and Multi-Layer Perceptron (MLP) are commonly used due to their efficiency and scalability in classification and regression tasks. Using these models to establish relationships between desertification indicators and environmental factors can successfully quantify the driving forces of desertification and predict future trends. However, machine learning models often act as “black boxes;” thus, understanding the connection between input variables and prediction outcomes is complex and difficult (Meng et al. 2021). SHapley Additive exPlanations (SHAP) uses the principles of game theory to calculate the independent contribution of each feature in machine learning model predictions and assign corresponding importance scores. Owing to its strong theoretical foundation and rich visualization tools, SHAP has been widely applied in the field of remote sensing. Hosseiny et al. (2022) provided crucial insights into the potential of nuanced importance metrics through explainable machine learning for urban land use and land cover classification. Descals et al. (2023) used SHAP to study the global climate’s limitation on photosynthesis phenology. Batunacun et al. (2021) employed SHAP to investigate severe grassland degradation in the Xilingol League, China, and its driving factors. These studies have demonstrated that SHAP not only increases the transparency of models in the field of remote sensing but also enhances the interpretability of the results, enabling researchers to improve the understanding of and trust in their models.

This study aimed to identify key areas and strategies for desertification prevention by revealing the spatiotemporal dynamics of global land desertification and its internal driving factors from 2000 to 2022, providing support for achieving the goal of zero growth in land degradation within the Sustainable Development Goals. The specific objectives were as follows: (1) achieve systematic monitoring of the global dynamic changes in desertification across arid land by using the GEE platform; (2) elucidate the long-term, change characteristics of desertification in arid lands and discern the segmental characteristics and rules of pixel-level desertification changes; and (3) uncover the mechanisms underlying pixel-scale transitions of desertification in arid land from 2000 to 2022. The analysis of the findings provides insights into the underlying causes of desertification dynamics.

2. Materials and methods

2.1. Remote sensing data

We used Moderate Resolution Imaging Spectroradiometer (MODIS) land products to comprehensively monitor global desertification in arid land. The data employed comprised two pivotal products: the MCD43A4 V6.1 Nadir Bidirectional Reflectance Distribution Function Adjusted Reflectance (NBAR) product (Zhen et al. 2023) and the MCD43A3 V6.1 Albedo Model dataset (T. He et al. 2019). The NBAR product facilitates the acquisition of 500 m reflectivity measurements for MODIS bands, which are subsequently refined through bidirectional reflectance distribution function adjustments. By amalgamating data from the Terra and Aqua spacecraft, this product effectively selects representative pixels over a 16 d period, ensuring optimal data quality and reliability.

The albedo data provides directional hemispherical reflectance (i.e. black-sky albedo) and bi-hemispherical reflectance (i.e. white-sky albedo) for each MODIS band. The spatial resolution and the synthesis method of the albedo data were consistent with those of the NBAR data.

The MODIS MCD43A4 version 6.1 NBAR dataset is generated daily using 16 days of Terra and Aqua MODIS data at a resolution of 500 meters (Huang et al. 2017). Nadir reflectance removes the influence of viewing angles, producing stable and consistent NBAR products (Geng et al. 2021). These data have undergone radiometric and geometric corrections, enabling them to provide high-quality, time series surface reflectance information. In the annual synthesis of the desertification monitoring index, clouds are detected and masked according to the “internal cloud algorithm flag” of the MODIS Quality band (Tian et al. 2018).

2.2. Environmental factors

The spatiotemporal variation and distribution of desertification are attributed to the combined influence of climate change and human activity. Through an extensive literature review and considering data availability, this study identified meteorological factors, topographic factors, and socioeconomic factors closely related to desertification or vegetation growth. Meteorological factors, namely, precipitation (Lamqadem et al. 2018), air temperature (Badreldin, Frankl, and Goossens 2014), vapor pressure difference (Li, Guan, and Zheng 2022), potential evapotranspiration (Meng et al. 2021), snow water equivalent (Liu et al. 2023), wind speed (Han et al. 2020; Sepehr et al. 2007), and soil moisture (Becerril and Mastachi 2021) directly affect the water cycle and energy balance in desertification areas. Topographic factors, namely, slope (Guo et al. 2023), aspect (Meng et al. 2020), and elevation (Yang et al. 2022), influence surface runoff, soil erosion, and microclimatic conditions, affecting the desertification process. Socioeconomic factors, namely, livestock quantity (Feng et al. 2015; Xie et al. 2023) and atmospheric carbon dioxide concentration (Piao et al. 2020), indirectly promote desertification by altering land cover and increasing greenhouse gas emissions.

Meteorological data were obtained primarily from TerraClimate (Abatzoglou et al. 2018), a global terrestrial monthly climate dataset with a spatial resolution of 0.04°. This dataset combines the high spatial resolution of the WorldClim dataset with a coarser temporal resolution through interpolation. TerraClimate provides monthly precipitation, temperature, Vapor pressure difference, wind speed, and potential evapotranspiration data and has been widely used in various research fields (Heinrich et al. 2023; Tiwari et al. 2023; Wu et al. 2023). The slope, aspect, and elevation were directly derived from the SRTM digital elevation model dataset (30 m spatial resolution) (Farr et al. 2007; Nasa 2013). The livestock quantity data were sourced from the Food and Agriculture Organization because this dataset provides comprehensive information and is widely recognized as a reliable source for livestock-related analyses. Additionally, the effects of elevated carbon dioxide concentration on vegetation recovery were considered (Park, Jeong, and Peñuelas 2020). Carbon dioxide concentration data were obtained from the National Oceanic and Atmospheric Administration. CarbonTracker data were obtained from the Atmospheric Administration’s Global Monitoring Laboratory Table 1. These datasets have been widely used to investigate the influence of carbon dioxide concentration on various ecological processes, including vegetation recovery (He et al. 2022; van der Laan-Luijkx et al. 2015).

Table 1. Environmental factors selected for modeling.

Category	Factor	Factor code	Resolution	Source/URL
Meteorology	Precipitation	PREC	4 km, interpolated to 500 m	TerraClimate dataset, https://www.climatologylab.org/terraclimate.html.
	Temperature	TEMP	4 km, interpolated to 500 m	TerraClimate dataset
	Potential evapotranspiration	PET	4 km, interpolated to 500 m	TerraClimate dataset
	Vapor pressure difference	VPD	4 km, interpolated to 500 m	TerraClimate dataset
Topography	Elevation	ELEV	30 m, interpolated to 500 m	Derived from digital elevation model (DEM), Shuttle Radar Topography Mission V3, https://www2.jpl.nasa.gov/srtm/.
	Slope	SLOPE	30 m, interpolated to 500 m	Derived from DEM
	Aspect	ASPECT	30 m, interpolated to 500 m	Derived from DEM
Soil	Soil moisture	SM	4 km, interpolated to 500 m	TerraClimate dataset
socioeconomic factors	livestock quantity	LIVS	5 arc minutes, interpolated to 500 m	Food and Agriculture Organization, https://www.fao.org/livestock-systems/en/.
	carbon dioxide concentration	CARB	3×2-degree, interpolated to 500 m	NOAA(National Oceanic and Atmospheric Administration), https://gml.noaa.gov/ccgg/carbontracker/download.php.

3. Methods

3.1. Arid land demarcation

Accurate identification and demarcation of arid land on a global scale form the fundamental basis for a comprehensive assessment of the dynamic trends of desertification (Liu et al. 2020). The aridity index (AI), introduced by the United Nations Environment Program, is a widely used measure for defining arid land (Ullah et al. 2022; Zarei and Mahmoudi 2021). The AI is the most rational and reliable indicator of climate change in such regions. AI = P/PET, where P denotes precipitation and PET represents potential evapotranspiration. Arid land can be divided into four categories based on the specific threshold values of AI: hyperarid (AI < 0.05), arid (0.05 ≤ AI < 0.20), semiarid (0.20 ≤ AI < 0.50), and dry subhumid (0.50 ≤ AI < 0.65) (Lian et al. 2021). This study encompassed all regions with an AI value below 0.65 (Figure 1).

Figure 1. Global arid index and arid land distribution map; (a) global arid index (b) arid land categories.

3.2. Developing an index for measuring desertification

Cloud cover in remote sensing images hampers the timely characterization of ground objects. Figure 2a displays the spatial distribution of the minimum number of effective observation times derived from MODIS satellite images between 1 January 2000, and 31 December 2022, elucidating the extent of cloud pollution in arid areas. Figure 2a also reveals that more than 100 effective observations are available from most areas, with West Asia, North Africa, and central and western Australia having over 250 high-quality observations. Hence, MODIS data are feasible for long-term land desertification monitoring.

Figure 2. Number of high-quality observations and desertification indicators in arid areas globally; (a) minimum number of yearly cloud-free observations during 2000–2022; (b) average annual MSAVI; (c) average annual albedo; (d) average annual DDI; (e) correlation MSAVI and albedo time series coefficient; and (f) annual average scatter plot of MSAVI and albedo.

The quantity of vegetation plays a crucial role in assessing the extent of desertification (Eskandari Dameneh et al. 2021). The vegetation index is widely used as a dependable measure to assess the conditions of vegetation growth. The MSAVI has demonstrated efficacy in desertification studies of arid land (Eskandari Dameneh et al. 2021; Meng et al. 2021; Wen et al. 2020). By leveraging the NBAR product data, the MSAVI was calculated using the following formula:

(1) $\mathit{MSAVI}=\frac{\left(2\mathit{NIR}+1-\sqrt{{\left(2\mathit{NIR}+1\right)}^2-8\left(\mathit{NIR}-\mathit{RED}\right)}\right)}{2},$(1)

where NIR represents the reflectivity of the near-infrared band, and RED corresponds to the reflectivity of the red band in the MODIS images. The global spatial distribution of the average MSAVI in arid land from 2000 to 2022 is characterized by predominantly high values in dry subhumid regions and low values mainly in hyperarid regions (Figure 2b).

Albedo is the ratio of the reflected surface flux to the incident solar radiation flux, which determines the extent of radiation energy absorption by the underlying surface (Guo and Wen 2020). Consequently, albedo is a significant parameter in the investigation of surface energy balance. With the exacerbation of desertification, there has been a simultaneous decline in surface water and surface roughness, increasing albedo. Thus, albedo is a crucial indicator for monitoring desertification (Lamaamri et al. 2023; Z. Wu et al. 2019). This study employed the MODIS surface albedo product MCD43A3 White-sky Albedo for the visible broadband spectrum to retrieve the albedo values. Figure 2c illustrates the spatial distribution characteristics of the average albedo from 2000 to 2022 within the global arid region. Notably, the hyperarid and arid regions of South America and Oceania exhibited the lowest albedo values, ranging from −2 to 0.15. The hyperarid land of North Africa and the arid and dry subhumid regions of northern Asia and North America displayed the highest albedo values.

The Desertification Distance Index (DDI) was calculated using a two-dimensional spatial feature model constructed by integrating the MSAVI and albedo (Guo et al. 2020). Different from the NDVI-albedo parameter estimation model (Lamaamri et al. 2023), the (1,0) non-desertification point serves as the reference point from which the distance to each point in all regions is evaluated, facilitating the characterization of the degree of desertification. This approach ensures that the independent evaluation results across regions remain comparable while maintaining the consistency of DDI, irrespective of the extent of the study area. Moreover, owing to DDI’s representation of continuous numerical data, it captures intricate details of land desertification, surpassing the capabilities of discrete data in denoting the degree of desertification. We employed DDI as the principal indicator for evaluating desertification in Mongolia. DDI was calculated using the following formula:

(2) $\mathit{DDI}=\sqrt{{\left(\mathit{MSAVI}-1\right)}^2-\mathit{Albed}{o}^2}.$(2)

Figure 2d illustrates the spatial distribution of DDI for arid areas. Regions such as North Africa, North Asia, Central Asia, and West Asia exhibited elevated DDI values. This phenomenon is particularly evident in North Africa and North Asia, where DDI reached a maximum of approximately 1. DDI was lowest near the 20° north and south latitudes, ranging from 0.3 to 0.65.

Figure 2e displays the correlation coefficient graph between MSAVI and albedo from 2000 to 2022, showing a negative correlation (correlation coefficient less than 0) within 90% of the areas, with 16% of the areas having a correlation coefficient less than −0.8. This result indicates that an increase in surface vegetation is often accompanied by a significant decrease in albedo. This strong negative correlation may be related to specific environmental conditions or vegetation types. The scatter plot in Figure 2f further supports this negative correlation, where the annual average correlation coefficient is approximately −0.6, indicating a moderate to strong negative correlation. This decrease suggests that on a spatial scale, as the vegetation index MSAVI increases, albedo decreases; the specific relationship between them is a notable foundation for constructing DDI.

3.3. Segmental trend analysis

Changes in the land desertification process are typically non-monotonous and exhibit significant spatial heterogeneity (Browning et al. 2017; Jiang et al. 2019; Yang et al. 2022). Therefore, conducting a segmental trend analysis of desertification at the pixel scale is necessary. Additionally, the development of desertification is a relatively stable long-term process. Although a reasonable assumption is the presence of multiple effective inflection points within a single pixel, DDI was calculated from annual MODIS image data and considered the length of the time series. Consequently, we adopted a binary linear regression method to detect changes in desertification trends (Meng et al. 2021). This method allows the identification of significant shifts in the desertification process.

(3) $\mathit{y}=\left\{\begin{array}{c}{\beta }_0+{\beta }_1\mathit{t}+\mathit{\epsilon },\mathit{t}\le \mathit{\alpha }\\ {\beta }_0+{\beta }_1\mathit{t}+{\beta }_2\left(\mathit{t}-\mathit{\alpha }\right)+\mathit{\epsilon },\mathit{t}>\mathit{\alpha }\end{array}\right.,$(3)

where y is the DDI trend, t is the year, α is the year at the turning point, β is the fitting parameter, and ε is the residual error. Because the slope obtained by linear regression is easily affected by outliers, this study used Sen’s slope analysis method (Gocic and Trajkovic 2013; Sen 1968) to calculate the DDI change trend before and after the mutation time and the Mann-Kendall method to test its significance (Silva et al. 2015).

3.4. Machine learning methods

We selected four commonly used machine learning models: Extreme Gradient Boosting (XGBoost), LightGBM (LGBM), Random Forest (RF), and Multi-Layer Perceptron (MLP). Each of these models has unique characteristics and can effectively handle both classification and regression tasks. The selection of these models was based on their outstanding performance and wide applicability in various fields such as finance, healthcare, and natural language processing. The following is a detailed introduction to each model.

XGBoost is a commonly used machine learning algorithm. It was developed by Chen and Guestrin (2016) and is commonly used owing to its effectiveness and efficiency in solving various machine learning problems. The primary use of XGBoost is for supervised learning tasks, such as classification and regression problems (Ma et al. 2021). It can manage numerical and categorical input features and be applied to various domains, including finance, health care, and natural language processing.

LGBM is a highly efficient gradient boosting framework developed by Ke et al. (2017). It was designed to be faster and more memory-efficient than traditional gradient boosting methods, making it suitable for large datasets with high-dimensional features. LightGBM supports both classification and regression tasks and can handle numerical and categorical features. It is widely used in various fields, including finance, healthcare, and recommendation systems, due to its high accuracy and scalability (Ke et al. 2017).

RF is an ensemble learning method introduced by Breiman (2001). It operates by constructing a multitude of decision trees during training and outputting the mode of the classes for classification tasks or the mean prediction for regression tasks. The method improves predictive accuracy and controls overfitting by averaging multiple deep decision trees. RF can handle a mix of numerical and categorical features and are robust to noise and overfitting, making them popular in various applications, including finance, biology, and engineering (Breiman 2001).

MLP is a type of artificial neural network that consists of multiple layers of nodes, each layer fully connected to the next. It was first conceptualized by Rumelhart et al. (1986). MLPs are used for both classification and regression tasks and can model complex relationships between inputs and outputs. They are particularly effective for problems where the relationship between inputs and outputs is nonlinear. MLPs have found applications in image recognition, speech recognition, and natural language processing (Rumelhart, Hinton, and Williams 1986).

We used the four machine learning models to construct nonlinear relationship models between environmental factors and the rate of change in DDI. Although machine learning models are widely popular for their excellent predictive performance, their decision-making processes are often not easy to understand due to their complexity. Thus, we used SHAP to provide detailed explanations for each prediction of the machine models, showing the specific contribution of each feature to the prediction outcome, increasing the transparency and credibility of the models.

3.5. Model interpretation

SHAP is a unified framework for interpreting the predictions of machine learning models (Lundberg and Lee 2017). It is based on the concept of Shapley values from cooperative game theory, which assigns each feature in a prediction model a unique contribution to the prediction outcome. The main goal of SHAP is to explain the output of a black-box model by attributing the prediction to different input features. This goal provides a means to understand the importance and impact of individual features on a model’s predictions (Zhang et al. 2023). SHAP has gained substantial attention in the field of interpretability and has been widely adopted in research and industrial applications (Feng et al. 2021; Mangalathu, Hwang, and Jeon 2020; Nazar et al. 2023). It provides a principled approach to understanding complex models and helps build trust and transparency in machine learning systems by explaining the reasoning behind their predictions. In cooperative game theory, Shapley values represent the average marginal contribution of each feature when considering all possible feature coalitions. In the context of machine learning, SHAP adapts Shapley values to compute the contribution of each feature in predicting a specific outcome. Shapely values can be calculated as follows:

(4) $\mathrm{g}\left(\mathit{x}{ }^{\mathrm{\prime }}\right)={\varphi }_0+\sum _{\mathit{i}=1}^M{\varphi }_i{x}_i^{\prime },$(4)

where M is input features’ quantities and ${\varphi }_0$ is the constant value when all $\mathit{x}{ }^{\mathrm{\prime }}$ features are missing. ${\varphi }_i\left(\mathit{f},\mathrm{X}\right)$ is defined as follows:

(5) ${\varphi }_i\left(\mathit{f},\mathrm{X}\right)=\sum _{\mathit{Z}{\mathit{ }}^{\mathrm{\prime }}\subseteq \mathit{X}{ }^{\mathrm{\prime }}}\frac{\left|\mathit{Z}{\mathit{ }}^{\mathrm{\prime }}\right|!\left(\mathit{M}-\left|\mathit{Z}{ }^{\mathrm{\prime }}\right|-1\right)!}{\mathit{M}!}\left[f_x\left(\mathit{Z}{ }^{\mathrm{\prime }}\right)-f_x\left(\mathit{Z}{ }^{\mathrm{\prime }}\mathrm{\setminus }\mathit{i}\right)\right],$(5)

where $\left|\mathit{Z}{ }^{\mathrm{\prime }}\right|$ is non-zero quantities in $\mathit{Z}{ }^{\mathrm{\prime }}$. ${\varphi }_i\left(\mathit{f},\mathrm{X}\right)$ is Shapely values (Lundberg and Lee 2017). SHAP facilitates the interpretation of feature contributions in model predictions on a local scale (local explanations) while enabling the computation of feature importance across the dataset (global explanations), offering insights into the relative significance of features for the overall model predictions. We applied the SHAP model to conduct an interpretability analysis of the results of the XGBoost model, quantifying the contribution rate of individual environmental factors and their interactions with the rate of change in DDI.

The technical flowchart in Figure 3 depicts the steps of this study. There are three main components: identification of desertification indicators, analysis of the segmental trend characteristics of desertification, and investigation of the driving mechanisms of desertification influenced by climate change and human activities.

Figure 3. Technical flow chart of the study.

4. Results

4.1. Fractional change characteristics of global desertification

Figure 4 presents a comprehensive analysis of DDI in arid regions, detailing segmental trends, change patterns, and turning points of time, as derived from the methodology outlined in Equation (3)(3) $\mathit{y}=\left\{\begin{array}{c}{\beta }_0+{\beta }_1\mathit{t}+\mathit{\epsilon },\mathit{t}\le \mathit{\alpha }\\ {\beta }_0+{\beta }_1\mathit{t}+{\beta }_2\left(\mathit{t}-\mathit{\alpha }\right)+\mathit{\epsilon },\mathit{t}>\mathit{\alpha }\end{array}\right.,$(3) . Before the turning point of DDI (denoted as T1 in Figure 4a), there was a predominant trend of decreasing DDI values, indicating a gradual reduction in desertification. The rate of DDI change was primarily within the range of −0.002/y and 0/y, constituting approximately 33% of total arid land. Conversely, the area exhibiting a DDI decrease rate surpassing −0.004/y, representing the most substantial reduction in desertification, accounted for 11% of total arid land. After the turning point (T2 in Figure 4b), the decreasing trend of DDI continued, indicating a further reduction in desertification. The rate of decrease after the turning point was higher than that during T1. The rate of decrease in DDI remained between −0.002/y and 0/y at T2, accounting for approximately 29%.

Figure 4. Comparison chart before and after DDI transition; (a) change trend of DDI before turning; (b) change trend of DDI after turning; (c) change type of DDI; and (d) time of DDI turning. Bar charts on the right show the proportion of each type on global arid land.

Figure 4c displays four main types of transitions in DDI: (1) decreases in the first and second stages (T1 < 0 and T2 < 0), (2) increases in the first and second stages (T1 > 0 and T2 > 0), (3) decreases in the first stage and increases in the second stage (T1 < 0 and T2 > 0), and (4) increases in the first stage and decreases in the second stage (T1 > 0 and T2 < 0). The most prevalent DDI transition types were T1 < 0 and T2 < 0, accounting for approximately 39%. TheT1 < 0 and T2 > 0 types accounted for 18%. Thus, the primary pattern of desertification was a decreasing trend that further decreased. Figure 4d provides insights into the distribution of the turning points in DDI. The majority of turning points were concentrated in 2005–2007, followed by 2016–2017, with a few occurrences in other years.

From a spatial perspective on trends, the majority of regions in Eurasia and North America exhibit consistent two-phase desertification development trends. However, this study observed a significant turning point in desertification in South America and Africa, potentially related to climate change and changes in land management practices. Such changes include fluctuations in precipitation, alterations in grazing management methods, and the implementation of ecological restoration projects (Meng et al. 2021). Additionally, with the demands of socioeconomic development, grazing pressure significantly increased in the later stages, elevating the risk of land desertification (Feng et al. 2015; Xie et al. 2023). Climate change also plays a substantial role in land degradation in these regions, and the turning phenomena of land desertification in Africa and South America could also be influenced by fluctuations in precipitation.

Figure 5 presents the statistical information regarding the trend of DDI across the arid land and the continents. In the first stage (before the turning point of DDI Figure 5a), the semiarid area exhibited the widest range of DDI values. The area of intensified desertification spanned approximately 1035.63 × 10⁴ km², the highest proportion of the total arid land (approximately 17.04%). The semiarid area had the greatest extent of desertification mitigation, approximately 997.71 × 10⁴ km² (16.76%). In dry, subhumid areas, the proportions of intensified desertification and desertification mitigation were relatively small, approximately 7%. In the second stage (T2, Figure 5b), the DDI trends occurred within the range from −0.002/y to 0/y. Notably, the semiarid area was the largest desertification area, approximately 13.67% of the total arid land area. Similarly, the semiarid area accounted for the largest proportion of the area experiencing desertification reduction, approximately 20.5% (1219.9 × 10⁴ km²). Dry subhumid areas exhibited the smallest extent of desertification aggravation or reduction among the four types.

Figure 5. Violin plot of DDI trends in arid areas and continents; (a) T1 DDI variation in arid areas; (b) T2 DDI variation in arid areas; (c) T1 DDI variation in continents; and (d) T2 DDI variation in continents.

Figure 5c and 5d shows an overview of the changing trends in DDI across all continents within arid land. Spatial variations in the degree of desertification were significant. Before DDI transition (Figure 5c), the Eurasian continent exhibited the largest area of DDI increase, approximately 1069.22 × 10⁴ km², 17.99% of the global arid area. Conversely, the smallest increase in DDI occurred in South America, 3.32% of the global arid area. Among these regions, the Eurasian continent had the greatest coverage of significantly intensified desertification (11.77%). In North America, the area of significant desertification intensification was the smallest, 89.33 × 10⁴ km². Regarding the areas where DDI decreased, Africa exhibited the widest distribution, approximately 23% of the global arid area, and South America exhibited the least distribution, approximately 4% of the global arid area. Among these regions, the areas with significantly reduced desertification were primarily in Africa, covering approximately 879.25 × 10⁴ km², and their distribution in North America was relatively small, accounting for 1.82% of global arid areas.

After the turning point of DDI (Figure 5d), desertification intensification was primarily observed in Eurasia, Africa and Australia, with relatively small distributions in South America. Among these regions, Africa and Australia exhibited the widest and narrowest range of significant desertification increases, respectively. Areas experiencing desertification reduction were mainly distributed in Africa and Eurasia, with the scope of desertification reduction in Africa being the most extensive, approximately 17.11% of global arid areas.

In T1 and T2, areas with semiarid and dry subhumid climates had the most significant increase in the proportion of land undergoing desertification. A likely explanation for this finding is that lands in semiarid and dry subhumid regions are ecologically fragile. Human activities are frequent in these areas, and ecological disturbances, primarily from agricultural cultivation, abandonment, and grazing, increase the probability of land degradation. These regions also become key focal points for later land science management and desertification control efforts.

4.2. Segmented characteristics of dominant explanatory factors

Figure 6 illustrates the segmental trends of the four representative explanatory factors: precipitation (PREC), potential evapotranspiration (PET), vapor pressure difference (VPD), and soil moisture (SM). During T1, the area of significantly increased PREC was approximately 2186.71 × 10⁴ km², and the spatial extent exhibiting a significant decrease in PREC was approximately 1456.39 × 10⁴ km² (Figure 6a). During T2, the increase area in PREC was considerably less than the corresponding decrease area (Figure 6b). Moreover, the predominant transformation types for PREC were T1 > 0 and T2 < 0, encompassing an area of approximately 4629.40 × 10⁴km². In the case of TEMP, PET, VPD, SWE, and VS, the areas exhibiting a significant increase in these five explanatory factors surpassed the areas displaying a decrease, and the prevailing transformation types for these factors were predominantly T1 > 0 and T2 > 0 (Figure 6c,f). Additionally, the decreased area of SM exceeded the increased area, with the primary transformation types being T1 < 0 and T2 < 0 (Figure 6g,h). For CARB, during T1 and T2, CARB had an increasing trend (T1 > 0 and T2 > 0).

Figure 6. Segmental trends of main explanatory variables; (a) PREC change trend during T1; (b)PREC change trend during T2; (c) PET change trend during T1; (d) PET change trend during T2; (e) VPD change trend during T1; (f)VPD change trend during T2; (g) SM change trend during T1; and (h) SM change trend during T2.

Global precipitation showed a significant trend of first increasing and then decreasing, with the Sahel region of North Africa as an example, where precipitation has a significant turning point and occupies a high proportion of the global precipitation turning area. This result is highly likely to influence the development trend of land desertification. PET and VPD, the main negative variables affecting vegetation growth, showed an increasing trends in both phases, indicating that global warming poses a risk of exacerbating global land desertification.

4.3. Influence of single factors on desertification

The relationship between environmental factors and DDI in arid land was established using XGBoost, RF, LGBM, and MLP. The precision of each model was calculated using 5-fold cross-validation and the evaluation metrics included Root Mean Squared Error (RMSE) and the Coefficient of Determination (R²). The results showed that RF performed the best across all metrics, with the lowest RMSE (20.5) and highest R² (0.89). XGBoost followed, with an RMSE of 22.5 and R² of 0.87. MLP model had moderate performance, with an RMSE of 36.5 and an R² of 0.66. LGBM performed the worst, with an RMSE of 42 and an R² of 0.57. Therefore, RF was selected to construct the relationship model between environmental factors and DDI.

Figure 7 shows the global and local SHAP values for the factors influencing desertification. Figure 7a depicts the global importance ranking of all factors during T1. Among these factors, PREC emerged as the most important factor, with the average magnitude of its SHAP value reaching + 9.92, followed by SM and CARB, ranked according to importance. Factors, namely, VS, SWE, LIVS, and TEMP exhibited relatively low importance, with average SHAP values below + 2. Among all the factors, SLOPE, VPD, and ASPE demonstrated the least importance; the total magnitudes of their SHAP values were lower than + 3.5. Figure 7b shows the ranking of global importance for all factors during T2. Similar to that during T1, PREC remained the most important factor, with the average magnitude of its SHAP value increasing to + 12.93. VPD and SM followed PREC in importance. Moreover, their importance in T2 was higher than that in T1. The remaining explanatory variables were less important than PREC, VPD, and SM.

Figure 7. Driving force diagram of explanatory factors of desertification process; (a) global importance of factors during T1; (b) global importance of factors during T2; (c) local importance of factors during T1; and (d) local importance of factors during T2.

A local explanations study using SHAP values was conducted to determine the significance and direction of the factors. Insights into the local explanatory information of all factors are provided in Figure 7c,d, where the x-axis represents the SHAP value, and the color represents the value of the factor. During T1, PREC, CARB, and SM were negatively correlated with global desertification, and PET was positively correlated with desertification. In addition, elevation did not significantly influence global desertification. However, the lower the elevation, the higher the likelihood of desertification. Therefore, based on Figure 7c, a focus should be on the possibility of desertification in low-elevation areas to prevent its occurrence and expansion. During T2, PREC and SM continued to exhibit a negative correlation with global desertification, and VPD showed a positive correlation. This result implies that as global VPD increases, the likelihood of it aggravating desertification increases.

4.4. Interactive effects of explanatory factors on desertification

Figure 8 illustrates the coupling effects between the factors before and after DDI transformation. During T1, the most significant couplings were observed between PREC ∩ PET, PREC ∩ SM, and PET ∩ SM in relation to the desertification process. The SHAP values for these couplings were approximately 3.2 (1.6 × 2). The couplings between PREC ∩ ELEV, PREC ∩ CARB, CARB ∩ SM, and PET ∩ CARB were relatively weak. There was almost no coupling between the explanatory variable groups SWE ∩ ASPE and SWE ∩ SLOPE, with the samples randomly distributed.

Figure 8. SHAP values of coupling effect between factors. (a) Coupling effect between factors during T1; (b) Coupling effect among factors during T2.

During T2, couplings between the factors were significantly enhanced. Notably, couplings between PREC ∩ SM, PREC ∩ VPD, VPD ∩ SM, and PREC ∩ CARB were stronger than other combinations of explanatory variables. Specifically, the SHAP coupling values for PREC ∩ SM and PREC ∩ VPD were as high as 5.06 (2.53 × 2) and 4.84 (2.42 × 2), respectively. Additionally, strong couplings were observed between VPD ∩ SM, PREC ∩ SWE, PREC ∩ CARB, PREC ∩ ELE, and VPD ∩ ELEV, with substantial SHAP coupling values. There was almost no coupling between the explanatory variable groups ASPE ∩ SLOPE and ASPE ∩ TEMP.

To facilitate a detailed discussion of the coupling effects between two explanatory factors, we generated dependence scatter plots (Figure 9). In these plots, the x-axis represents the principal factor value, the y-axis represents the SHAP value, and the color indicates the secondary factor value. The coupling dependence between the pivotal factors of desertification during T1 is shown in Figure 9a–c. Figure 9a demonstrates that as PREC increased, PET decreased, and their coupling effect mitigated the degree of desertification, although the influence was relatively limited (mostly less than −0.002/y). Conversely, as the PREC trend decreased, the PET trend gradually increased, which contributed to the exacerbation of desertification, exhibiting a relatively extensive impact ranging from 0 to 0.006/y. As shown in Figure 9b, SM increased with the increase in PREC, and their combined effects led to a mitigation of desertification. This result indicates that PREC and SM are positively correlated, and both variables have a negative correlation with the rate of change in DDI. In other words, as PREC and SM increase, desertification is alleviated. As shown in Figure 9c, PET and SM exhibited a significant negative correlation, with SM values decreasing as PET increased, while the SHAP value for PET increased significantly, indicating an intensification of desertification.

Figure 9. Coupling dependence scatter plots of main explanatory factors; (a) PREC ∩ PET, (b) PREC ∩ SM, (c) PET ∩ SM, (d) PREC ∩ SM, (e) VPD ∩ SM, and (f) PREC ∩ VPD.

The coupling dependence of the pivotal factors of desertification during T2 is shown in Figure 9d–f. When the PREC trend was less than 0, the SM exhibited a downward trend, indicating the intensification of desertification (indicated by a positive SHAP value). As PREC and SM increased, the degree of desertification demonstrated a decreasing trend, as low as −0.008/y. As shown in Figure 9e, SM decreased with increasing VPD, leading to the continued intensification of desertification. When the VPD trend was negative, SM rapidly increased, further intensifying desertification. Regarding the interaction between PREC and VPD (Figure 9f), it was exactly the opposite of that between PREC and SM. As PREC increased, VPD significantly decreased, leading to a rapid intensification of desertification, with the maximum annual increase rate of DDI being around 0.01/y. These dependence scatter plots provide valuable insights into the coupling effects between explanatory factors and their impact on desertification during stages T1 and T2.

Globally, natural factors, especially precipitation and SM, remain the primary influence on the process of desertification, however, the impact of human activities in specific regions must be considered. The reversal of desertification in countries such as China and India is mainly attributed to ecological restoration projects and the expansion of farmland (Piao et al. 2020). Moreover, in the second phase, the contribution of the factor of livestock numbers increased, second to SM. Thus, with socioeconomic development, rational grazing and scientific land management have become increasingly important measures for preventing and mitigating land desertification.

4.5. Accuracy assessment of desertification results

Vegetation cover can serve as an important indicator for characterizing the degree of desertification, and vegetation cover has been used as the criterion for assessing desertification (Hein and De Ridder 2006; Kundu et al. 2015; Symeonakis and Drake 2004; Wei et al. 2018). This section demonstrates the credibility of the DDI results by evaluating the correlation between DDI and vegetation cover. Vegetation cover was calculated using the pixel dichotomy model (Wei et al. 2018). Figure 10 shows the distribution of the correlation coefficients between DDI and vegetation cover (VC) globally for 2000–2022. Regarding the spatial distribution of correlation, 86% of the global area had a correlation coefficient less than −0.95, with the remaining 14% (correlation coefficient greater than −0.95) of the area primarily concentrated in extremely arid desert regions such as the Sahara Desert, the Arabian Peninsula, and the Taklamakan Desert. This finding may be due to the extremely low VC and DDI values in these regions, where the change in both values during the study period is less than the sensor noise value, causing spurious correlations. Figure 10 shows six sites from various continents to show specific correlations. All six scatter plots show a clear negative correlation (correlation coefficient greater than 0.99); as VC decreased, DDI correspondingly increased (i.e. the degree of desertification intensified).

Figure 10. Correlation graph between Global Desertification Distance Index (DDI) and Vegetation Cover (VC).

5. Discussion

5.1. Global desertification trends and their driving factors

According to statistical results from different arid regions or continents, there has been a significant reversal in desertification globally since 2000. The abundance of vegetation is an important indicator for characterizing desertification in arid lands and if vegetation is used as a reference, the results of this study are consistent with those of Piao et al. (2020), indicating that global vegetation greenness has been increasing during this study. The reversal of desertification is the mainstream trend, and regional differences are pronounced. China and India are the main areas of desertification reversal for opposite reasons. China has implemented the Grain for Green program, resulting in a significant increase in vegetation in the Loess Plateau and achieving zero growth in land degradation (He, Shang, and Zhang 2021; Wang et al. 2021). India has implemented agricultural intensification, converting arid land into farmland, which also achieved the goal of increasing vegetation from the satellite perspective but many scholars believe this method is unsustainable. Moreover, the amount of dust released during the fallow period of crop rotation is also large (Pi et al. 2018). As shown in Figure 4, the majority of areas in northern Russia on the Eurasian continent have experienced an intensification of desertification in both phases, which is directly related to the increase in atmospheric saturated VPD, increased evapotranspiration, and reduced SM. This research result refines the belief that rising temperatures are the cause of intensified desertification in Russia (Zolotokrylin et al. 2016). There is a distinct semiarid region in northeastern Brazil in South America. The present study showed that before 2016, most areas in this region showed a weakening trend in desertification, with few cities in the central part showing an intensifying trend. Our results are highly consistent with those of Bezerra et al. (2020). Areas where desertification intensified are characterized by tropical-arid or tropical-semiarid climates, with the Caatinga biome dominating (Bezerra et al. 2020). Our observations indicate a general decline in desertification trends in the Sahel region, located to the south of the Sahara in North Africa, aligning with the findings reported by Yang et al. (2022). However, different from Yang et al. (2022), we did not find that precipitation was the main driving force affecting the reversal of desertification in the region; instead, VPD, PET, and SM were the main influencing factors. This difference may be because Yang et al. (2022) based the driving model on spatial structure. Spatially, desertification is related to precipitation but land desertification is a dynamically changing geographical phenomenon, and analyzing the temporal development trend of desertification from fixed locations is an important method for dissecting the causes of desertification. Therefore, the attribution results we obtained can realistically reflect the causes of changes in desertification. This type of research, using the variability of geographical variables as the dependent variable in the attribution model, has been widely applied in other fields of geography (Huang et al. 2021; Meng et al. 2020).

5.2. Key areas for desertification control and suggested measures

As shown in Figure 5, lands experiencing an intensification of desertification are primarily concentrated in arid and dry subhumid regions, becoming globally vulnerable areas for desertification (He et al. 2023; Xue et al. 2019). The reasons for this finding include but are not limited to human activities such as excessive grazing, farming and abandonment of cultivated lands, and mineral development, coupled with the variability in precipitation caused by unstable monsoons in these regions, which further increases the frequency of surface vegetation fluctuations (D’Odorico et al. 2013). Therefore, areas for desertification control should first focus on regions heavily affected by climate change and human activities, such as the agro-pastoral ecotones and grazing areas (Yue, Geng, and Li 2023; Zhou et al. 2019). Regarding specific spatial distribution (Figure 4), regions where desertification intensified in both phases mainly included central and western Australia, eastern Russia, the north-central United States, south-central Canada, and a small part of eastern Mongolia. These regions have specific characteristics of land degradation; thus, desertification control should be tailored to local conditions. Central and western Australia and eastern Mongolia should protect and restore native vegetation and scientifically develop agriculture and animal husbandry based on water resources and livestock carrying capacity. North America should develop sustainable land use and management methods, improve agricultural practices, and enhance drought response capabilities. Desertification in eastern Russia is mainly attributed to climate change, deforestation, and improper use of water resources; thus, its control policies should strengthen forest management, limit deforestation, and enhance the capabilities of forest fire monitoring and firefighting.

This study systematically analyzed land desertification development information across global regions and clarified patterns of desertification development and driving forces, providing important information to support the formulation of environmental policies and land management.

5.3. Limitations and further research

This study has three notable limitations. First, as shown in Figure 2a, certain regions, particularly North Asia, exhibited a scarcity of observations, with the minimum number of observations per year below 100. Fewer cloud-free observations led to a significant reduction in the quality of regional satellite monitoring, which directly increased the uncertainty of the MSAVI observations, affecting the assessment results of global desertification trends. Therefore, further research could combine multi-source remote sensing data to reduce the influence of clouds, cloud shadows, and other factors on monitoring results and to improve research accuracy (Balha et al. 2021).

Second, because desertification is a slow-changing process, using the dichotomous trend analysis method in this study to assess the changing characteristics of desertification was reasonable, considering the relatively short study period of 23 years (Meng et al. 2020). However, certain arid lands worldwide may undergo desertification development processes that involve three or more stages of change. Therefore, further research could refine the segmental analysis approach to comprehensively capture and account for the multistage change characteristics of global desertification, deepening the understanding of desertification dynamics and increasing the accuracy of assessments of desertification patterns and trends.

Third, desertification can be characterized by factors other than changes in vegetation. For example, the transformation from grassland to shrubland landscapes is another form of land degradation or desertification (Bestelmeyer et al. 2018). This definition is because an increase in shrub cover, leading to a decrease in grassland area, increases the grazing concentration, accelerating grassland degradation (van de Koppel et al. 2002). Additionally, in areas dominated by shrubs, the low grass cover increases the likelihood of wind erosion, leading to the loss of soil organic matter and creating a vicious cycle of land degradation. However, such degradation is difficult to capture with medium-resolution satellite imagery such as MODIS and AVHRR. Thus, for such types of land degradation, further research could combine sub-meter satellite data and drone imagery to develop identification and assessment methods.

6. Conclusion

We utilized MODIS satellite images and the machine learning and SHAP model to systematically evaluate the segmental change characteristics of desertification in arid lands globally. The impact of climate change and human activities on global desertification was quantitatively analyzed, providing new insights into its prevention and control. The main conclusions are as follows:

From 2000 to 2022, desertification of arid land exhibited a decreasing trend globally, with a slower reduction rate before the turning point of DDI than after the turning point. Continuous desertification mitigation was the dominant turning type, and the transition time of DDI was primarily concentrated in 2005–2007 and 2016–2017. During T1, PREC was the most important factor, followed by SM and CARB. During T2, PREC remained the most important factor, followed by VPD and SM. Additionally, the coupling effect analysis revealed strong couplings between PREC and PET, PREC and SM, and PET and SM during T1, with a SHAP value of 3.2. During T2, the couplings among explanatory factors were significantly enhanced, particularly between PREC and SM, PREC, and VPD, with a SHAP value of 4.9 for these two groups.

In summary, fluctuations in PREC and SM are significant bimodal characteristic driving factors of global land desertification. In the later stages of the study, the impact of livestock numbers on desertification increased. Future desertification prevention and control should focus on dry subhumid and semiarid regions. The research revealed that the development trend of land desertification in some areas has multiple fluctuation characteristics, necessitating segmented methods for precise identification.

These findings highlight the complex dynamics of desertification and underscore the importance of considering both climatic and anthropogenic factors to understand its changes. This study provides valuable insights to guide effective measures to prevent and control global desertification.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The raw data required to reproduce the above findings are available to download from https://modis.gsfc.nasa.gov/data/(MODIS. images) and https://www.climatologylab.org/terraclimate.html (Terraclimate data). The processed data required to reproduce the above findings are available to download from 10.6084/m9.figshare.25989064.v1 and 10.6084/m9.figshare.25990534.v1.

Additional information

Funding

This work was supported by the Natural Science Foundation of China [grant number 42301494], the China Postdoctoral Science Foundation [grant number 2023M730944] and the Third Xinjiang Scientific Expedition Program [grant number 2021xjkk0305].