Introduction

Soil erosion by wind and water is one of the critical causes degrading soils over the globe, which affects 84% of the total degraded lands area, with respective contributions of 28% and 56%1. Wind erosion causes many environmental problems with on-site and off-site consequences, while the total global cost for its economic losses is estimated about US$ 400 billion per year2. This phenomenon presents the highest rate in Asia, followed by Africa, South America, North and Central America, and lastly Europe3. Intensified wind erosion by climate change due to higher temperature, less precipitation and expanded desertification, is a real threat for the soil essential nutrients, soil isotopes, agriculture productivity, resilient survival and food security4,5. Therefore, prediction and mapping of global land susceptibility to soil erosion by wind is essential to mitigate its negative impacts, and especially important for resilience policies under climate change and a warmer world scenario.

To the best of our knowledge, there are several studies dealing with mapping of land susceptibility to wind erosion hazard, geochemical, mineralogy issues and related dust emissions using machine learning (ML) techniques, google earth engine, and artificial intelligence models at local and regional scales like in Central Asia6,7, northwest China8,9, Morocco10, Iran11,12, Iran-Iraq borders13, Middle East14, Europe4. However, up to date there is no study providing a wind erosion hazard map at a global scale based on ML new generation models such as deep learning (DL) and/or its integration with active learning (AL).

Nowadays (2021 to present), DL consisting of the individual and hybrid DL models, as new generation of ML models, have been successfully applied to map land susceptibility to soil erosion hazard by wind at local15 and regional scales16. ML and DL algorithms are helpful tools for analysis and mapping of large datasets with high spatial resolution and for classification purposes of land susceptibility over large spatial scales such as continents16,17.

Active learning (AL)—a unique variety of machine learning—focuses on the exploration of datasets, it assumes various samples in the same dataset and tries to select the samples with the highest value to build the training dataset18. AL has been well-motivated in many ML problems especially when labels are difficult, time-consuming and expensive to obtain. AL consists a class of techniques that are applied to train ML algorithms with limited labeled data by picking the most valuable data points to be labeled19.

One of the lesser-known aspects of predictions provided by the DL predictive models is their interpretability. This aspect plays a key role in the better understanding of predictive model’s output20,21,22,23. Interpretability is the degree to which a human can understand the cause of a decision. Due to the black box nature of ML models, it is difficult to understand how they arrive at their conclusions, because these models are trained on massive amounts of data and employ complex algorithms that are difficult to interpret24. A wide range of techniques has been employed to develop the interpretable models, and one of the interpretation classification techniques can categorize models into the local (Shapley Additive explanation (SHAP), local surrogate (LIME), etc.) and global techniques (e.g., partial dependence plot (PDP), permutation feature importance (PFI), etc.)25,26.

Therefore, such kind of research should aim for assessment and better understanding of the environmental conditions (topographic, atmospheric, climatic, pedological, vegetation, etc.) under which land degradation by wind erosion is most likely to occur. This study introduces an interpretable DL-AL model aiming to assess the important issue of global land degradation and wind erosion hazard. The main objectives are summarized as: (i) to identify the effective variables affecting soil erosion by wind using Harris Hawk Optimization (HHO) algorithm—a feature section algorithm for the identification of the important factors controlling wind erosion, (ii) to map land susceptibility to wind erosion over the global lands using DL models (e.g., recurrent neural network (RNN) and gated recurrent units (GRU)) with AL (RNN-AL, and GRU-AL) and without AL (GRU and RNN), and (iii) the development of an interpretable DL-AL model via three interpretation techniques (e.g., synergy matrix, ALE plot and SHAP decision plot). This is the first time that land susceptibility is mapped over a global scale via an interpretable DL-AL model. The suggested methodology has high capability of providing spatial accurate maps of wind erosion at different scales (from local to large scales), and on the other hand, this approach is capable to address the nature of the DL black-box and to interpret relationship between input variables and model’s output (or to estimate the contribution—or impact—of the input variables on the model’s output). The development and application of this approach at a global scale is the main advantage and innovation of the current work, while the results may be especially useful for geological, atmospheric and socio-economic sciences, regarding the deleterious consequences of wind-erosion hazard on soil nutrients, agriculture productivity, food security, dust storms, floods, etc.

Materials and methods

Datasets and wind erosion inventory map (WEIM)

Several variables may control, or directly and indirectly affect, the wind erosion, which exhibit large variability on space and time27. In this research, based on literature review, we initially mapped spatially thirteen likely variables that may control land susceptibility to wind erosion. These thirteen variables include 7 soil characteristics (topsoil carbon content, clay fraction, gravel fraction, organic carbon, sand fraction, silt fraction, and soil moisture), aspect, digital elevation model (DEM) or elevation, land cover, vegetation cover (normalized difference vegetation index (NDVI)), and climatic variables (precipitation and wind speed)28,29,30,31,32.

The data for all soil characteristics with exception of soil moisture was taken from the website https://daac.ornl.gov. This dataset, with a spatial resolution 50 m × 50 m, describes the global soil variables from the Harmonized World Soil Database. The soil moisture was downloaded from https://climate.northwestknowledge.net/TERRACLIMATE/index_directDownloads.php33.

The data for precipitation and wind speed were taken from websites as https://www.worldclim.org, and https://www.noaa.gov, respectively. Vegetation cover was mapped by Google Earth Engine (GEE), while NDVI products with a spatial resolution of 250m were taken from MODIS. Land cover was also mapped by GEE and MODIS land cover product. DEM with a spatial resolution of 30 m × 30 m was downloaded from the website https://earthexplorer.usgs.gov, and aspect was extracted from DEM in ArcMAP 10.5. After extraction of the data required for all variables (soil characteristics and climatic variables), the maps generated using inverse distance weighted (IDW) interpolation technique in ArcMAP 10.5. In comparison with another interpolation techniques (e.g., Kriging, Co-Kriging, Spline, and Natural Neighbor), IDW exhibited the lowest error, and was chosen as the most accurate technique. Overall, all input variables applied for the mapping wind erosion were converted to a consistent spatial resolution of 50 m × 50 m. The global spatial maps for all the thirteen variables impacting on wind erosion are presented in the Supplementary Figs. 1–13.

The spatial maps of the effective variables and WEIM are necessary to build the predictive model for mapping land susceptibility to wind erosion hazard. WEIM indicates the areas with and without the signs of wind erosion (Fig. 1) and was produced based on World Bank17,34,35 and Nobakht et al.36. In this study, the MODIS imagery was applied to map WEIM. Then, we randomly used 70% (258 points) and 30% (129 points) of total datasets as training and test (validation) datasets, respectively (see Fig. 1), in order to establish and test the model’s performance over the global scale. The cross-validation (CV) technique was used to evaluate the model and to test its performance. Finally, the predictive models were built for mapping land susceptibility to wind erosion. Indeed, from a technical point of view, interpolations with high accuracy can be performed with much lower count of points than 387.

Figure 1
figure 1

Location of the training (yellow dots) and test datasets (red dots) to build the predictive DL-AL model for mapping global land susceptibility to wind erosion. Emphasis is given in the global deserts and semi-arid lands. The training and test datasets are overlaid on the Esri world imagery (Earthstar Geographics) using ArcGIS 10.8 environment.

Full size image

HHO feature selection algorithm to identify important variables

Feature selection is a key step in the spatial modelling of the natural and environmental hazards by ML and DL models, and this step can improve the model’s performance significantly by removing the non-important variables and select and entrance important variables controlling wind erosion into the modelling phase. HHO is a population-based and nature-inspired technique37,38 and was applied here to identify effective variables conditioning on the land susceptibility to wind erosion hazard. The application of the HHO has been reported in fields such as predicting of temperature of motor stator39, shear tension40, and other fields, but its application in the wind erosion has not been reported yet. HHO has three phases consisting of exploration, transition and exploitation37,38,41.

Exploration phase

In this phase, it is planned to wait mathematically, seek and detect the required hunt. The iter + 1 (the Harris Hawks ___location) can be written as follows:

$$X\left( {t + 1} \right) = \left\{ {\begin{array}{*{20}c} {X_{rand} \left( {iter} \right) - \left. {r_{1} } \right|X_{rand} \left( {iter} \right) - 2 X_{r} \left( {iter} \right)\; \ldots \; \ldots \;if\;\; q \gg 0.5} \\ {X_{rabit} \left( {iter} \right) - X_{m} \left( {iter} \right) - r_{3} (LB + r_{4} \left( {UB - LB)} \right)\; \ldots \; \ldots \;if\;\; < 0.5} \\ \end{array} } \right.$$
(1)

where Xrabit and Xrand indicate the rabbit ___location and the elected hawk randomly among the obtainable population ri, respectively. i = 1, 2, 3, …, q are random numbers varying in [0,1], iter represents the local iteration and Xm as the mean ___location for hawks can calculated as follows:

$$X_{m} \left( {iter} \right) = \frac{1}{N}\mathop \sum \limits_{i = 1}^{N} X_{i} \left( {iter} \right)$$
(2)

where, Xi and N represent the ___location and size of the hawks, respectively.

Transition phase

This phase describes the transition from exploration to exploitation, and the energy of rabbit (E) can be calculated by:

$$E = 2 E_{0} \left( {1 - \frac{iter}{T}} \right)$$
(3)

where, T indicates the upper size of the repetitions E0 = ɛ (− 1, 1) as the premier energy during every step.

Exploitation phase

In this step, the Harris Hawks perform the surprise pounce by attacking the intend prey detected in the previous phase. When |E| ≥ 0.5 and |E| < 0.5 the soft and hard besiege happen, respectively37,38,39,40,41.

According to this technique, from the thirteen features possibly impacting on land susceptibility to wind erosion mentioned above (section “Datasets and wind erosion inventory map (WEIM)”), eight features consisting of wind speed, topsoil carbon content, topsoil clay content, elevation, topsoil gravel fragment, precipitation, topsoil sand content and soil moisture were selected as the most effective variables and were entered into the modelling phase. On the contrary, the rest five variables classified as non-effective (e.g., aspect, land cover, vegetation cover, topsoil organic content, and topsoil silt content) were excluded from further analysis.

Simple RNN and GRU deep learning models

RNN—a type of convolutional feedforward neural network—such as reinforcement learning and long short-term memory (LSTM) applied for unsupervised learning in many application domains42,43,44. RNN is consisting of input (x0, x1, …, xi, xi+1 …), hidden (h0, h1, …, hi, hi+1 …) and output (y0, y1, …, yi, yi+1 …)42.

The function of RNN can be written as follows:

$$h^{\left( t \right)} = \sigma \left( {Wh^{{\left( {t - 1} \right)}} + Ux^{\left( t \right)} } \right)$$
(4)
$$o^{\left( t \right)} = \sigma \left( {Vh^{\left( t \right)} } \right)$$
(5)

where, h(t) and o(t) indicate the hidden state and model’s output at time t, respectively. W, U and V represent the weight of state transfer, the weight of input and output, respectively.

GRU—a version of LSTM45—optimizes the LSTM network architecture, while it keeps LSTM performance. In comparison with LSTM, which has three gates (input, forget and output), GRU has two gates consisting of update and reset. Reset gate determines how to combine the new input with the previous memory, while update gate defines how much of the previous memory to retain. Forget gate controls range and remaining values in the cell. Mathematically, GRU can be written with the following equations:

$$rt = \sigma \left( {W\left( r \right)xt + U\left( r \right)ht - 1} \right)$$
(6)
$$zt = \sigma \left( {W\left( z \right)xt + U\left( z \right)ht - 1} \right)$$
(7)
$$hti = tanh\; \left( {rt\; o\; Uht - 1 + Wxt} \right)$$
(8)
$$ht = \left( {1 - zt} \right) \;o\; hti + zt \;o\; ht - 1$$
(9)

where, rt and zt indicate reset and update gates, respectively. hti and ht represent the new memories and hidden state, respectively. \(\sigma (.)\), W(r) and W(z) are the sigmoid function, and weight matrices, respectively.

GRU integrates the forget and input gates into a single update and merges the cell state and hidden state along with some other changes43,45,46.

AL-DL modeling to map wind erosion

AL—as an iterative training approach selecting the best samples to train—has been shown to reduce labeled data requirement, when training deep classification networks47,48. The critical constraint behind AL is that the performance of the ML model may increase with fewer training labels if it is allowed to choose the data from which it learns19. Therefore, the simple RNN and GRU models with (RNN-AL and GRU-AL) and without AL were applied to map land susceptibility to wind erosion hazard. The integrated framework based on AL and DL models for mapping land susceptibility to wind erosion is presented in Fig. 2.

Figure 2
figure 2

Schematic approach of the AL-DL integrated framework used to map land susceptibility to wind erosion.

Full size image

Predictive model’s performance assessment

There are many statistical indicators to assess the performance of the predictive models, based on regression and classification approaches. In this study, three curve and plots consisting of area under of receiver operating characteristic (AUROC) curve, cumulative gain (CG) and Kolmogorov Smirnov (KS) statistic plots were applied for assessing the performance of the four DL models consisting of simple RNN, RNN-AL, GRU and GRU-AL.

Development of an interpretable DL model to map wind erosion

Although ML and DL models performed very well and provided more accurate results, due to the nature of their “black boxes”, the application of these models among researchers and experts is minimal21,22. Therefore, to address the “black box” issue associated with ML and DL models, the application of the interpretation techniques to develop interpretable ML and DL models has been increasing remarkably35,49. Here, we applied three interpretation techniques consisting of synergy matrix, ALE and SHAP plots to interpret and explain the feature contributions and its impact on the model’s output, as well as how feature impacts on model prediction.

The synergy of feature xi relative to another feature xj quantifies the degree to which predictive contributions of xi rely on information from xj. Synergy as percentage of the feature importance is expressed as a percentage ranging from 0% (full autonomy) to 100% (full synergy).

ALE plots average the changes in the predictions and accumulate them over the grid:

$$\begin{aligned} \hat{f}_{S,ALE} \left( {x_{S} } \right) & = \mathop \int \limits_{{z_{0,S} }}^{{x_{S} }} E_{{X_{C} \left| {X_{S} = x_{S} } \right.}} \left[ {\hat{f}^{S} \left( {X_{s} ,X_{c} } \right)\left| {X_{S} = z_{S} } \right.} \right] dz_{s} - constant \\ & = \mathop \int \limits_{{z_{0,S} }}^{{x_{s} }} \left( {\smallint {_{XC} \hat{f}^{S} \left( {z_{s} ,X_{c} } \right)dP(X_{C} |X_{S} = z_{S} )d } } \right) dz_{S} - constant \\ \end{aligned}$$
(10)

where f is the prediction function at feature value (s) xS, averaged over all features in XC. E is the marginal expectation over the features in set C.

SHAP, introduced by Lundberg and Lee50, is based on game theoretically optimal Shapley values. The Shapley value for feature Xj in a predictive model can be written as follows:

$$Shapley\; \left( {X_{j} } \right) = \mathop \sum \limits_{{S \subseteq N\backslash \left\{ j \right\}}} \frac{{k!\left( {p - k - 1} \right)!}}{p!}\left( {f\left( {S \cup \left\{ j \right\}} \right) - f(S)} \right)$$
(11)

where, p represents the number of features (p = 8 in our case; the important variables), and N\{j} is a set of all possible combinations of features excluding Xj,. S indicates a feature set in N\{j} and f(S) is the model prediction with features in S, respectively. f(S \(\cup\) {j}) is the model prediction.

SHAP decision plot indicates how the sophisticated models arrive at their predictions or make decisions50. At SHAP decision plot, x-axis and y-axis represent the model’s output and model’s features, respectively. In this plot, the features are sorted by descending importance, and the importance is calculated over the observations plotted. All the critical steps for the spatial modelling of global land susceptibility to wind erosion by an interpretable DL-AL model are presented in the flow chart followed in this study (Fig. 3).

Figure 3
figure 3

Research diagram (flow chart) used for mapping of global land susceptibility to wind erosion hazard.

Full size image

Results

Validation of land susceptibility models to wind erosion

CG, KS plots and AUROC (for training and test datasets) were employed for the assessment performance of the four predictive models used for mapping land susceptibility to wind erosion hazard at the global scale (Figs. 4, 5, 6). CG plot visualizes the percentage of the target class members (Y-axis indicates cumulative gains) vs decile (X-axis). This plot has two reference lines or models (random model (with the black dotted lines in the Figs. 4 and 5) and wizard model (with the green dotted lines in the Figs. 4 and 5) that tell us how good or bad our model is performing. The random model line indicates what proportion of the actual target class would be expected to select when no model is used at all, whereas the wizard model is the upper bound of what the model can do. Overall, according to this plot, the model will always lie between wizard and random models, while closer to wizard model shows better performance and away from the wizard indicates worse performance. The CG plot for training and test datasets indicates that GRU-AL is the most accurate model for mapping land susceptibility to wind erosion, while the prediction accuracies of the models evaluated by the CG plot were classified as follows: GRU-AL > GRU > RNN-AL > RNN.

Figure 4
figure 4

Cumulative gain (CG) plots and Kolmogorov–Smirnov (KS) plots used to assess four predictive models (e.g., RNN, RNN-AL, GRU and GRU-AL) based on the training dataset.

Full size image
Figure 5
figure 5

Cumulative gain (CG) plots and Kolmogorov–Smirnov (KS) plots used to assess four predictive models performance (e.g., RNN, RNN-AL, GRU and GRU-AL) based on the test dataset.

Full size image
Figure 6
figure 6

AUROC curves for assessing the performance of predictive models (e.g., RNN, RNN-AL, GRU and GRU-AL) based on training datasets (a,b), and test datasets (c,d).

Full size image

Based on the training dataset, the accuracy of the model predictions assessed by KS plot was classified as follows: GRU-AL > GRU > RNN-AL > RNN, similar to previous case. Overall, according to test dataset, the accuracy classification of four predictive models evaluated by KS plot were as follows: GRU-AL > RNN-AL > GRU > RNN. Both plots (CG and KS) indicate that the performance of all four models for the training dataset is more accurate than the respective performances based on the test datasets.

The results of the performance validation for four predictive models consisting of RNN, RNN-AL, GRU and GRU-AL, based on both datasets (training and test), are presented in Fig. 6. The current findings represent that GRU-AL (with AUROC = 96 for training dataset, and AUROC = 76 for test dataset) performed better than GRU and RNN-AL (AUROC = 94 and 74 for training and test datasets, respectively), and RNN (AUROC = 91 for training dataset, and AUROC = 71 for test dataset). A difference of the AUROC by around 20 when comparing training and test data can indicate overfitting, and it occurs when the predicting algorithm ends up in describing too well the in-sample data, and incorporating not only the information related to the data-generating process51. According to Yesilnacar52, AUROC values of 90–100 and 70–80 indicate an excellent and good prediction performance, respectively. In the current study, the order of the accuracies of the model predictions assessed by AUROC was as follows: GRU-AL > GRU and RNN-AL > RNN.

Land susceptibility maps to wind erosion hazard

The values of the land susceptibility classes for each pixel (or the probability of wind erosion occurrence) predicted by the models (GRU, GRU-AL, RNN, and RNN-AL) vary between zero and unit. In this study, the susceptibility probability predicted by models for each pixel was classified into 5 wind erosion hazard classes consisting of: 0–0.2 (very low class), 0.2–0.4 (low class), 0.4–0.6 (moderate class), 0.6–0.8 (high class), and 0.8–1 (very high class)29,53. Finally, the predictive probability of wind erosion occurrence estimated by four models was converted into the developed maps of land susceptibility to wind erosion hazard in ArcMAP 10.5 software.

Figure 7 represents the land susceptibility maps to wind erosion hazard generated by RNN (Fig. 7a), RNN-AL (Fig. 7b), GRU (Fig. 7c), and GRU-AL (Fig. 7d) models. According to all three performance measures (e.g., CG plot, KS plot, and AUROC curve), GRU-AL outperformed the GRU, RNN-AL and RNN. Therefore, based on the results of the most accurate model—GRU-AL—very low, low and moderate susceptibility classes to wind erosion hazard cover 38.5%, 12.6% and 10.3% of global land area, respectively, whereas 12.5% and 26.1% of the total land belong to the high and very high susceptibility classes, respectively. Overall, the most vulnerable areas in terms of land susceptibility to wind erosion hazard around the world can be divided into fourteen regions consisting of: southern coast of Greenland, Central Australia, Central Asia and its southern regions (Thar Desert, India), East China, the Middle East, North Africa, South Africa (Namib, Kalahari Deserts), East Africa, Western USA, Central USA, Northern part of USA, parts in South America (mostly Atacama Desert), while some other regions in Europe and Russia may also consist of hotspots for land susceptibility to wind erosion.

Figure 7
figure 7

Land susceptibility map to wind erosion hazard predicted by RNN (a), RNN-AL (b), GRU (c), and GRU-AL model (d). The susceptibility probability predicted by models for each pixel was classified into 5 wind erosion hazard classes consisting of: 0–0.2 (very low class), 0.2–0.4 (low class), 0.4–0.6 (moderate class), 0.6–0.8 (high class), and 0.8–1 (very high class). The background map is the Esri world terrain base map (relief/ocean map) using ArcGIS 10.8 environment.

Full size image

Interpretability of predictive model’s output

Synergy matrix

Synergy is expressed as percentage of the feature importance. It is additive and sums up to 100% for any pair of features. Importantly, neither the relationship is necessarily symmetrical, while one feature may replicate or complement some or all of the information provided by another feature, while the reverse need not be the case. The synergy matrix of the features controlling land susceptibility to wind erosion is presented in Fig. 8. According to this matrix, for any feature pair (A, B), the first feature (A) is the row and the second feature (B) the column. The results present that wind speed and precipitation exhibited high synergy with carbon content, DEM and soil moisture. Furthermore, Carbon content presented high synergy with DEM and soil moisture on the predictions provided by the model. Overall, exploration of the synergy (correlation) between features is necessary because high synergy between pairs of features should be considered carefully when investigating the impact, as the values of both features jointly determine the outcome. Current findings show that the synergy of carbon content with soil moisture (26%) is the most contributed impact on the model’s predictions, followed by the synergy of carbon content with DEM (20%). Indeed, for the vast majority of the pairs of features, the synergy is very low, mostly below 10%.

Figure 8
figure 8

The synergy matrix for features controlling land susceptibility to wind erosion.

Full size image

ALE plots

The Accumulated Local Effects (ALE) plots for the prediction model of land susceptibility to wind erosion by eight features consisting of carbon content, clay content, DEM, gravel content, precipitation, sand content, soil moisture and wind speed are presented in Fig. 9. According to the results, soil carbon content and precipitation display negative strong effects on the prediction of land susceptibility to wind erosion. For the clay content feature, the ALE plot indicates that the predicted wind erosion probability is low, on average up to percentage of 7%, and no trend after that threshold exists, so this feature does not significantly affect the model predictions. About the DEM feature, the elevation exhibits a positive strong effect on the prediction. On the other hand, the gravel fragment feature does not present a significant effect on the predictions. For the sand content feature, the ALE plot shows that the probability of the predicted land susceptibility to wind erosion is low (or a negative effect is presented), on average, for less than 75% of the values, but this feature shifts to a positive strong effect on the model’s predictions for sand content > 75%. The ALE plot for the soil moisture feature highlights the strong negative effect on the predictions by the model, while the wind speed feature exhibited also negative effects on the model predictions.

Figure 9
figure 9

ALE plots of eight variables impacting on wind erosion.

Full size image

SHAP decision plot

The SHAP decision plot of the model’s performance is presented in Fig. 10. At this plot, each observation’s prediction is represented by a colored line. At the top of the plot, each line strikes the x-axis at its corresponding observation’s predicted value. This value determines the color of the line on a spectrum. Moving from the bottom to the top of the plot, SHAP values for each feature are added to the model’s base value. This indicates how each feature contributes to the overall prediction. According to SHAP decision plot (Fig. 10), soil moisture, DEM and precipitation are considered as the three most important features controlling land susceptibility to wind erosion, where sand content and gravel fragment exhibited the lowest impact on the model’s output. Overall, soil moisture is identified as the most important variable controlling wind erosion and this can be attributed to several potential mechanisms. For example, soil moisture can impact on the vegetation cover as a roughness element, which can affect the soil carbon content and wind speed. On the other hand, soil characteristics (e.g., clay, sand and gravel contents) were identified as the least important variables. This could be explained by the fact that soil properties are relatively stable properties and do not change frequently like soil moisture.

Figure 10
figure 10

SHAP decision plot to interpret model’s output.

Full size image

Discussion

Land susceptibility to wind erosion hazard worldwide

An important finding that is revealed from this research is that the southern coast of Greenland can be a new region, which suffered by wind erosion after snow and glaciers melting, nowadays. The motion of glaciers generates fine-grained sediments that is delivered via meltwater to proglacial floodplains. Therefore, when the glaciers retreat, a great land area with fine-grained sediments may be exposed to wind activity, while the phenomenon of Greenland glaciers retreat becomes especially important during the last decade due to climate change and global warming54,55,56. Glaciers watersheds located in this part of Greenland receive a large amount of fine- to coarse-grained fluvially-sediment57,58,59. Then, seasonal discharge flows deposit these sediments, especially sand and silt-grained particles, across outwash plains following desiccation, and in absence of dense vegetation, they entrained by strong downglacier or other kinds of winds60. The local dust emissions at the southern coastal region of Greenland are increasing after 2000, due to a decreasing seasonal snow-cover area arising from coastal Greenland warming, while this area is experiencing drastic changes nowadays61,62.

Australia’s land (with exception of the northern part, a narrow strip in Eastern coast and the spots in Central Australia) was characterized by high and very high susceptibility to wind erosion (Fig. 7). The arid lands of Simpson Desert and Lake Eyre Basin in Central Australia presented high and very high susceptibility to wind erosion, while some sedimentary environments (e.g., the dried beds of the ephemeral lakes), de-vegetated burnt areas by wildfires63,64, floodplains that are closely linked to sediment supply from flood events, and sand dunes are the prominent dust sources in the Australia’s interior65. Yang and Leys66 revealed that western and central regions (e.g., Lake Eyre Basin) of Australia are the most susceptible lands to wind erosion hazard, in close agreement with the current findings.

Current results also highlighted the dominance of very high and high susceptibility classes on lands over Africa. Due to high climatic erosivity, high soil erodibility and low vegetation cover in the lowland plains, and along the coastal areas of East Africa, some countries such as Sudan, Somalia, Eritrea and Kenya are nowadays highly susceptible to wind erosion29. Nearly, all countries located in North Africa (e.g., Algeria, Morocco, Tunisia, Sahara, Niger, Chad, Libya, Egypt, Mauritania) and South Africa (e.g., Southern African Republic, Namibia, Botswana, southern Angola, Zambia and Zimbabwe) exhibited high and very high susceptible lands to wind erosion. According to Fenta et al.29, more than 25% of lands in East and Southern Africa were classified at moderate to high-risk classes of wind erosion. The main reasons for high and very high land susceptibility to wind erosion in Eastern Cape Wild Coast, South Africa are the sparse vegetation cover, and the strong-wind exposed coastal regions covered by unconsolidated and erodible sandy soils67. One of the key sources of dust in South Africa is the Namib Desert, and the largest dust emissions are from source areas located on recently deposited fluvial surfaces, which are not active in the current age68. Due to sand soils and sparse vegetation cover, the most susceptible lands to wind erosion around the western part of Southern Africa are the Kalahari basin are the western coast of Namibia in the Namib Desert69. On the other hand, due to inappropriate human activities and land mismanagement associated with rained agriculture, other semi-arid lands have converted into dust sources in South Africa70.

The central Asian region consists of several deserts (e.g., Taklimakan, Balkhash dryland, Karakrum, Aralkum, Kyzylkum, Hexi Corridor, etc.) covering an area of about 7.75 M km2, and therefore, has been recognized as an important source for global dust emissions, contributing about 17–20%71. Current findings from all the ML models revealed that Central Asia and NW China, as a unit region, is the largest wind erosion region over the globe, while these drylands belong to high and very high susceptibility classes (Fig. 7). These results are in agreement with those of a previous work17, which extensively explored and identified the land susceptibility to wind erosion over Central Asia with high spatial resolution. So, extensive discussion about the main “hotspot” areas of land susceptibility to wind erosion can be found there17. Dust emissions originated from Central Asia can be transported thousands of kilometers downwind due to blowing of various wind regimes (westerlies, northerlies)58, and affect large parts of east and south Asia and countries like China, Japan, Korea, Iran, Afghanistan, Pakistan, etc.72.

According to current results, the Middle East was identified as a region with high and very high susceptibility to wind erosion hazard, and therefore, can play a key role in the regional and global dust cycle. These findings are consistent with several previous studies in this region, which characterized areas in central, east and south Iran as major hotspots for soil erosion and dust emissions14,35.

The lands of the southern great plains and Chihuahuan Desert in southwestern USA presented the highest susceptibility to wind erosion. In addition, the wind erosion in West Texas and eastern New Mexico States occurs in localized source areas, while the playas and cropped lands are the most vulnerable sources in this region73. The Chihuahuan Desert is a significant regional dust hotspot in North America74. Kandakji et al.75 found that cultivated croplands (enclosing 43% of the dust points), ephemeral lakes and dried playas in southwestern USA are the most important sources for production of dust in USA.

In Southern America, the most important wind erosion regions are located in Argentina (Pampa area) and in some parts of Peru and Bolivia. Atacama Desert76 and Patagonia are the two most important deserts in Southern America, and current results revealed that the southern parts of this continent, along with some hotspots in the eastern part of Brazil, exhibited high and very high land susceptibility to wind erosion.

Although limited compared to the regions discussed above, according to the EU Thematic Strategy for Soil Protection, about 42 million hectares are affected by wind erosion in Europe, with a serious threat in agricultural lands4. Therefore, land degradation due to wind erosion may affect the semi-arid Mediterranean areas77, as well as vast areas in northern Europe1,4. Our results reveal higher susceptibility to wind erosion in parts of northern Europe like Germany, the Netherlands, Denmark, as well as in areas of Spain, Alps (possibly due to extensive melting of the glaciers) and in Carpathian Mountains. Some hotspots are also predicted in the Scandinavian Alps in Norway (Fig. 7). However, these results are not consistent very well with a previous work in Europe4. The differences may originate from the different spatial resolutions associated to datasets used to map wind erosion, as well as to different approaches, models and techniques. At any case, the high susceptibility observed in Alps is a new finding of the current research.

According to the classification by78, 81.3% and 13.8% of the investigated lands in Europe were characterized by slight and moderate erodibility, respectively, whereas 4.9% were in the high erodibility class. This classification also revealed a division of Europe in three regions according to erodibility namely (i) the northern Europe (Poland, Denmark, the Netherlands and northern Germany) with high erodibility values due to sandy plains, (ii) the central-eastern Europe with moderate values and (iii) the southern Europe (Mediterranean area) with low wind-erodible fraction values4. Current results revealed higher susceptibility in the Mediterranean region, which are consistent with those presented by Borrelli et al.4.

Research advantages and limitations

Our modelling approach indicated that a new developed interpretable model based on DL-AL learning models and interpretation techniques is a promising tool to generate global spatial maps for soil erosion by wind with high accuracy. Same technique seems to be also applicable for other environmental hazards (e.g., soil erosion by water, floods, landslide, dust storms, etc.), at different scales from small (field and catchments) to large (global scale). In comparison with another modelling tools such as physical models and field-based techniques, interpretable models are low-cost techniques and can discover the nature of the model’s black-box by converting black-box to grey-box. On the other hand, the high spatial resolution of the input variables (e.g., DEM, satellite images used to calculate vegetation indices, suitable distribution of meteorological stations for providing data such as temperature, wind speed and precipitation) is necessary to provide spatial accurate maps, especially at small scales. Overall, the current research revealed optimistic results regarding the examination of global land susceptibility to wind erosion and we recommend that future research may employ the interpretable DL-AL model developed in this study for spatial modelling purposes of natural and environmental hazards to provide spatial accurate maps from local to global scales.

Conclusions

In this study, a new interpretable spatial modelling framework, integrating DL-AL models and SHAP interpretation technique, was developed to map global land susceptibility to wind erosion hazard. Among four deep learning models (e.g., RNN, RNN-AL, GRU and GRU-AL) used in current study, the GRU-AL model was found to be the most accurate in the spatial mapping of land susceptibility to wind erosion. Methodologically, the AL technique can improve the performance of the predictive DL models, and the interpretation techniques can clear the black box aspect of DL models.

Due to nature of DL-AL model’s black box, development of an interpretable model is necessary to analyze and interpret the contribution and importance of the input variables (features) on the predictive model’s output. Current findings revealed that soil moisture impacted on several variables (e.g., vegetation cover, wind speed, soil organic carbon) and finally, it can affect significantly the predictions provided by the model. Furthermore, this variable exhibited high synergy with soil carbon content regarding the global land susceptibility prediction to wind erosion hazard.

GRU-AL model results revealed that 38.5%, 12.6%, 10.3%, 12.5% and 26.1% of the global lands presented very low, low, moderate, high and very high susceptibility to wind erosion hazard, respectively. The modelling results indicated that apart from the desert, arid and semi-arid areas in the vicinity of the global dust belt, high latitudes such as the southern coast of Greenland, Alaska and some hotspots in other regions such as Siberia, Iceland, Canada displayed high and very high land susceptibility to wind erosion. This is a likely consequence of the global warming and retraction of ice sheets in the northern latitudes, thus contributing to land degradation and making these regions vulnerable to wind erosion. Overall, our suggested spatial modelling approach with some remote sensing techniques, can be a promising technique to contribute to future research, especially to map soil erosion by wind and another natural disasters such as soil erosion by water, floods, landslides, etc.