1. Introduction

2. Materials and Methods

2.2. The MEDALUS model

The Environmentally Sensitive Area (ESA) index was based on the modified framework of the MEDALUS model [11], which uses four quality indices (QI): climate quality index (CQI), soil quality index (SQI), vegetation quality index (VQI), and management quality index (MQI). The framework has variations from that of the original method, such as the addition of information on Soil Carbon Content and Normalized Difference Vegetation Index (NDVI) (Table 1).

Driving factors of the desertification process were classified according to scores ranging from 1 (least sensitive to desertification) to 2 (most sensitive to desertification) for the determination of quality indices (Qis); values between 1 and 2 indicate intermediate susceptibility. Thus, obtaining the Environmentally Sensitive Areas (ESAs) took place in two stages: in the first stage, the QIs were obtained using the geometric mean of the variables that compose the indices using Equations 1 and 2 [9].

QI_x = (variable_1_ij * variable_2_ij * variable_3_ij … * variable_n_ij)^1/n

(1)

where: n = number of variables; x = quality indices (CQI, SQI, VQI, and MQI); iej = row and column of a pixel of the used rasters, respectively.

In the second stage, the previously obtained QIs were used to calculate the ESA:

ESA_ij = (CQI_ij * SQI_ij * VQI_ij * MQI_ij)^1/4

(2)

where: i and j = is the row and column, respectively, of a pixel of the rasters.

Thus, ESAs obtained from the listed QIs were classified into eight desertification sensitivity classes [22]: unaffected (N), potentially affected (P), fragile (F1, F2, F3), and critical (C1, C2 and C3) (Table 2).

Table 2. Types of Desertification Sensitivity Indices.

Table 2. Types of Desertification Sensitivity Indices.

Susceptibility	Weights	Subclasses
Not affected	< 1.170	No
potentially affected	1.170 - 1.225	P
Fragile	1.225 - 1.275	F1
	1.275 - 1.325	F2
	1.325 - 1.375	F3
Critical	1.375 - 1.425	C1
	1.425 - 1.530	C2
	> 1.530	C3

2.3. Quality scores

The respective QIs were obtained using digital databases in raster format. Part of these bases were already available in the GEE catalog. The unavailable ones in the platform were acquired and reclassified according to the weights (Table A1) established in the methodology to finally send them to the GEE platform for the creation of the ESAs map (Figure 1).

Figure 2. Flowchart showing the components associated with climate (CQI), soil (SQI), vegetation (VQI) and management (MQI) quality indexes to calculate the Environmental Sensitive Area (ESA). R – rainfall; Ar – aridity; As – basin orientation; St – soil type; Sl – basin slope; Tx – soil texture; Dd – drainage classes; OC – organic carbon content; Al – albedo; FR – fire risk; DR – drought resistance; EP – erosion protection; NDVI – normalized difference vegetation index; LU – land use; PD – population density. .

Figure 2. Flowchart showing the components associated with climate (CQI), soil (SQI), vegetation (VQI) and management (MQI) quality indexes to calculate the Environmental Sensitive Area (ESA). R – rainfall; Ar – aridity; As – basin orientation; St – soil type; Sl – basin slope; Tx – soil texture; Dd – drainage classes; OC – organic carbon content; Al – albedo; FR – fire risk; DR – drought resistance; EP – erosion protection; NDVI – normalized difference vegetation index; LU – land use; PD – population density. .

CQI was determined using the variables rainfall (R), aridity (Ar), and basin orientation (As) based on the GEE catalog bases. The mean annual total precipitation information (mm) came from the TerraClimate global dataset [23]. Aridity was calculated using the UNESCO aridity index [24] from the ratio between precipitation (mm) and potential evapotranspiration (mm). Potential evapotranspiration (ETP), for the period 2010 to 2019, necessary to obtain Ar, came from the Terra Net Evapotranspiration Digital Base, which uses MODIS satellite images by the Land Processes Distributed Active Archive Center (LPDAAC). Basin orientation was obtained using SRTM images from Hole-filled SRTM for the globe, version 4.

SQI was obtained using the variables soil type (St), basin slope (Sl), soil texture (Tx), drainage classes (Dd), organic carbon content (OC), and albedo (Al). Basin slope was extracted from the previously cited SRTM data. Information on texture, drainage classes and organic carbon content was extracted from the Harmonized World Soil Database [25], reclassified according to the defined weights, and sent to the GEE platform. Surface albedo was obtained from an annual mosaic for the year 2019, composed of pre-processed images from the Sentinel-2 (S2) satellite. Albedo calculation in S2 images was performed using coefficients defined in an approach proposed by [26].

VQI is composed of the sub-indices fire risk (FR), drought resistance capacity (DR), erosion protection capacity (EP) and the NDVI. The FR, DR and EP variables were obtained from the reclassification (Table 1) of a land use mosaic made available by MapBiomas [14] for the year 2019. The image, available on GEE, was then cut and reclassified according to the parameters established for each variable that composes the VQI. The NDVI for the study area was obtained using the mosaic composed of S2 images using the following equation:

NDVI = (IR – RED)/(IR – RED)

(3)

where IR is the infrared surface reflectance on band eight; RED is the red surface reflectance on band four of the S2 satellites.

MQI was composed of the variables intensity of land use (LU) and population density (PD) of the municipalities present in the basin under study. Land use intensity data come from the land cover classification image provided by MapBiomas and had scores assigned to each class, as previously determined in Table 1. The PD variable was acquired from the ratio between the number of inhabitants in a municipality and its respective area in km². As the last population census was carried out in Brazil in 2010, it was necessary to estimate the population in 2020.

Maps based on the Environmentally Critical Factors (ECF) index were created to complement the ESA map for the East Atlantic Basin. The composition of the ECF index was performed with the following equation [27]:

$${\mathrm{ECF}}_{\mathrm{ij}}={\sum }_{\mathrm{n}=1}^{15}\mathrm{f}({\mathrm{v}\mathrm{ariable}}_{\mathrm{nij}})$$

(4)

$$f(x)=\left\{\begin{array}{c}1, \mathrm{s}\mathrm{e} x \text{>} 1,425\\ 0, \mathrm{s}\mathrm{e} x \text{≤} 1,425\end{array}\right.$$

(5)

Pixels with values lower than the critical limit (C2) of 1.425 were not considered when processing. Thus, the composition of the ECF index was made from counting the number of variables that exceeded the critical value of 1.425 for each pixel.

3. Results

3.2. Environmentally Sensitive Area (ESA)

After applying Equation 2 using previously generated QIs, an ESAs map was obtained (Figure 4). The ESAs comprise the geometric mean of the four QIs generated in the first stage of the work.

In analyzing the ESA map, the north and central-west regions of the study area has areas with greater vulnerability to desertification. Conversely, areas located closer to the coastal zone and to the southwest of the basin were less susceptible. Only 5.0% of the basin area (Table 2) can be classified in the “Not affected” category (N), with emphasis on areas located in the region further to the coast and southwest of the basin.

The areas included in the “Potentially affected” (P) class comprise about 10.1% of the total area of the basin and are located, as well as the N class, in the coastal and southern regions. The three levels corresponding to the “Fragile” (F) class have a relatively uniform spatial distribution throughout the basin, occupying around 49.3% of the entire area. In general, subclasses F1 and F2 have similar spatial distribution to classes N and P, while subclass F3 is more present in the northern and central portion of the studied area, the latter being more associated with the “Critical” class (C). Class C corresponds to a portion of 35.7% of the basin and, as previously mentioned, the areas included in this classification (C1, C2, C3) are located to the north and center of the area under study.

Table 2. Desertification susceptibility classes and their respective areas within the East Atlantic Basin, Brazil.

Table 2. Desertification susceptibility classes and their respective areas within the East Atlantic Basin, Brazil.

Susceptibility	Weights	Subclasses	Area (km²)	Area (%)
Not affected	< 1.170	No	19177.49	5.0%
potentially affected	1.170 - 1.225	P	38883.42	10.1%
Fragile	1.225 - 1.275	F1	56033.22	14.5%
	1.275 - 1.325	F2	68675.25	17.8%
	1.325 - 1.375	F3	65941.47	17.1%
Critical	1.375 - 1.425	C1	54968.20	14.2%
	1.425 - 1.530	C2	70176.42	18.2%
	> 1.530	C3	12656.84	3.3%

When analyzing the areas characterized by high sensitivity to desertification, more specifically to the north of the basin, the presence of low-quality soils and high rates of aridity are predominant factors in these locations.

A joint analysis of the ECF index (Figure 5) and ESA allows a better understanding of the contribution of each model component to the desertification process. This analysis shows that certain areas with the same ESAs have variables with different sensitivity to desertification [27]. Thus, when analyzing the ECF maps, we can observe that variables related to the climate and soil of the study area are the most common factors for determining sensitivity to desertification in the East Atlantic Basin.

The PCA results (Table 3 and Figure 6) allowed a more objective assessment of the contribution of each edaphoclimatic and anthropic variable used for the composition of the ESAs maps for the East Atlantic Basin, enabling the optimization of the model in future research.

The first three PCs could explain 61.31% of the total variance, which can be considered a significant value due to the spatial extent of the study area. Principal component 1 (29.65%) was associated with variables related to VQI (Fire Risk, Resistance to Drought and Resistance to Erosion) and MQI (Land Use Intensity). PC1 was also associated with albedo, although it did not show a high correlation. Principal component 2 (18.92%) showed a significant correlation with the climatic variables that compose the CQI (Rainfall and Aridity), also correlating with the NDVI. Finally, PC3 (12.75%) was exclusively linked to the component variables of the SQI (Soil Texture, Drainage Classes, and Organic Carbon Content).

In analyzing the results in Figure 6, we can have a better understanding of the relationship between the selected variables, Qis, and ESAs, showing that NDVI, albedo, and VQI are positively correlated with ESAs.

4. Discussion

5. Conclusions

References

1. Liu, Z.; Si, J.; Deng, Y.; Jia, B.; Li, X.; He, X.; Zhou, D.; Wang, C.; Zhu, X.; Qin, J.; Ndayambaza, B.; Wang, B. Assessment of land desertification and its drivers in semi-arid alpine mountains: a case study of Qilian Mountains Region, northwest China. Remote Sens. 2023, 15. [Google Scholar] [CrossRef]
2. Wijitkosum, S. Reducing Vulnerability to Desertification by Using the Spatial Measures in a Degraded Area in Thailand. Land 2020, 9, 49. [Google Scholar] [CrossRef]
3. Ren, Y.; Liu, X.; Zhang, B.; Chen, X. Sensitivity assessment of land desertification in China based on multi-source remote sensing. Remote Sens. 2023, 15, 2674. [Google Scholar] [CrossRef]
4. Bakr, N.; Weindorf, D. C.; Bahnassy, M. H.; El-Badawi, M. M. Multi-temporal assessment of land sensitivity to desertification in a fragile agro-ecosystem: Environmental indicators. Ecol. Indic. 2012, 15, 271–280. [Google Scholar] [CrossRef]
5. Djeddaoui, F.; Chadli, M.; Gloaguen, R. Desertification susceptibility mapping using logistic regression analysis in the Djelfa area, Algeria. Remote Sens. 2017, 9, 1031. [Google Scholar] [CrossRef]
6. Wang, W.; Jiang, Y.; Wang, G.; Guo, F.; Li, Z.; Liu, B. Multi-Scale LBP texture feature learning network for remote sensing interpretation of land desertification. Remote Sens. 2022, 14, 348. [Google Scholar] [CrossRef]
7. Xu, D.; You, X.; Xia, C. Assessing the spatial-temporal pattern and evolution of areas sensitive to land desertification in North China. Ecol. Indic. 2019, 97, 150–158. [Google Scholar] [CrossRef]
8. Meza Mori, G.; Torres Guzmán, C.; Oliva-Cruz, M.; Salas López, R.; Marlo, G.; Barboza, E. Spatial Analysis of Environmentally Sensitive Areas to Soil Degradation Using MEDALUS Model and GIS in Amazonas (Peru): An Alternative for Ecological Restoration. Sustainability 2022, 14, 14866. [Google Scholar] [CrossRef]
9. Abuzaid, A. S.; Abdelatif, A. D. Assessment of desertification using modified MEDALUS model in north Nile Delta, Egypt. Geoderma 2022, 405, 115400. [Google Scholar] [CrossRef]
10. Lahlaoi, H.; Rhinane, Hilali, A. ; Lahssini, S.; Moukrim, S. Desertification assessment using MEDALUS model in watershed Oued El Maleh, Morocco. Geosciences 2017, 7, 50. [Google Scholar] [CrossRef]
11. Kosmas, C.; Ferrara, A.; Briassouli, H.; Imeson, A. Methodology for Mapping Environmentally Sensitive Areas (ESAs) to Desertification; European Commission: Luxembourg. 1999. [Google Scholar]
12. Wijitkosum, S. The impact of land use and spatial changes on desertification risk in degraded areas in Thailand. Sustainable Environment Research 2016, 26, 84–92. [Google Scholar] [CrossRef]
13. Vieira, R. M. S. P.; Tomasella, J.; Alvalá, R. C. S.; Sestini, M. F.; Affonso, A. G.; Rodriguez, D. A.; Barbosa, A. A.; Cunha, A. P. M. A.; Valles, G. F.; Crepani, E.; de Oliveira, S. B. P.; de Souza, M. S. B.; Calil, P. M.; de Carvalho, M. A.; Valeriano, D. M.; Campello, F. C. B.; Santana, M. O. Identifying areas susceptible to desertification in the Brazilian northeast. Solid Earth 2015, 6, 347–360. [Google Scholar] [CrossRef]
14. Tomasella, J.; Vieira, R. M. S. P.; Barbosa, A. A.; Rodriguez, D. A.; Santana, M. O.; Sestini, M. S. Desertification trends in the Northeast of Brazil over the period 2000-2016. Int J Appl Earth Obs Geoinformation 2018, 73, 197–206. [Google Scholar] [CrossRef]
15. Chen, R.; Yang, H.; Yang, G.; Liu, Y.; Zhang, C.; Long, H.; Xu, H.; Meng, Y.; Feng, H. Land-Use Mapping with Multi-Temporal Sentinel Images Based on Google Earth Engine in Southern Xinjiang Uygur Autonomous Region, China. Remote Sens. 2023, 15, 3958. [Google Scholar] [CrossRef]
16. Souza, C.M., Jr.; Z. Shimbo, J.; Rosa, M.R.; Parente, L.L.; A. Alencar, A.; Rudorff, B.F.T.; Hasenack, H.; Matsumoto, M.; G. Ferreira, L.; Souza-Filho, P.W.M.; et al. Reconstructing Three Decades of Land Use and Land Cover Changes in Brazilian Biomes with Landsat Archive and Earth Engine. Remote Sens. 2020, 12, 2735. [Google Scholar] [CrossRef]
17. Kadri, N.; Jebari, S.; Augusseau, X.; Mahdhi, N.; Lestrelin, G.; Berndtsson, R. Analysis of Four Decades of Land Use and Land Cover Change in Semiarid Tunisia Using Google Earth Engine. Remote Sens. 2023, 15, 3257. [Google Scholar] [CrossRef]
18. Velastegui-Montoya, A.; Montalván-Burbano, N.; Carrión-Mero, P.; Rivera-Torres, H.; Sadeck, L.; Adami, M. Google Earth Engine: A Global Analysis and Future Trends. Remote Sens. 2023, 15, 3675. [Google Scholar] [CrossRef]
19. Gorelick, N.; Hancher, M.; Dixon, M.; Ilyushchenko, S.; Thau, D.; Moore, R. Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote sensing of Environment 2017, 202, 18–27. [Google Scholar] [CrossRef]
20. Santos, L.S.; Mendes, L.A. Influence of the Equatorial Pacific Ocean Temperature Oscillation on the Precipitation Regime of the East Atlantic Hydrographic Region. Brazilian Journal of Physical Geography 2020, 13, 1502–1518. [Google Scholar] [CrossRef]
21. Brasil. Ministry of the Environment. East Atlantic Hydrographic Region Notebook. Brasília: Ministry of Health, 2006.
22. Karamesouti, M. , Panagos, P., Kosmas, C. Model-based spatio-temporal analysis of land desertification risk in Greece. Catena 2018, 167, 266–275. [Google Scholar] [CrossRef]
23. Abatzoglou, J.T.; Dobrowski, S. Z.; Parks, S. A.; Hegewisch, K. C. TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. Scientific data 2018, 5, 170191. [Google Scholar] [CrossRef]
24. United Nations Organization for Education, Science and Culture. Map of the world distribution of arid regions: Explanatory note. Paris: Unesco, 1979.
25. FAO/IIASA/ISRIC/ISSCAS/JRC. Harmonized world soil database (version 1.2). Rome, Italy and IIASA, Laxenburg, Austria: FAO, 2012.
26. Bonafoni, S.; Sekertekin, A. Albedo Retrieval from Sentinel-2 by New Narrow-to-Broadband Conversion Coefficients. IEEE Geoscience and Remote Sensing Letters 2020, 17, 618–1622. [Google Scholar] [CrossRef]
27. Ferrara, A.; Kosmas, C.; Salvati, L.; Padula, A.; Mancino, G.; Nolè, A. Updating the MEDALUS-ESA Framework for Worldwide Land Degradation and Desertification Assessment. Land Degradation & Development 2020, 31, 1593–1607. [Google Scholar] [CrossRef]
28. Lamqadem, A. A.; Pradhan, B.; Saber, H.; Rahimi, A. Desertification sensitivity analysis using Medalus model and GIS: a case study of the oases of middle draa valley, Marocco. Sensors 2018, 18, 2230. [Google Scholar] [CrossRef] [PubMed]
29. Lahlaoi, H. ; Rhinane, Hilali, A.; Lahssini, S.; Moukrim, S. Desertification assessment using MEDALUS model in watershed Oued El Maleh, Morocco. Geosciences 2017, 7, 50. [Google Scholar] [CrossRef]
30. Cunha, J.; Nóbrega, R. L. B.; Rufino, I.; Erasmi, S.; Galvão, C.; Valente, F. Surface albedo as a proxy for land-cover clearing in seasonally dry forests: evidence from the Brazilian Caatinga. Remote Sensing of Environment 2020, 238, 111250. [Google Scholar] [CrossRef]
31. Nabiollahi, K.; Golmohamadi, F.; Taghizadeh-Mehrjardi, R.; Kerry, R.; Davari, M. Assessing the effects of slope gradient and land use change on soil quality degradation through digital mapping of soil quality indices and soil loss rate. Geoderma 2018, 318, 16–28. [Google Scholar] [CrossRef]
32. Vieira, R.M.S.P.; Cunha, A. P. M. A.; Alvalá, R. C. S.; Carvalho, V. C.; Ferraz Neto, S. . Sestini, M. S. Land use and land cover map of a semiarid region of Brazil for meteorological and climatic models. Brazilian Journal of Meteorology 2013, 28, 129–138. [Google Scholar] [CrossRef]
33. Mallick, J.; AlMesfer, M.K.; Singh, V.P.; Falqi, I.I.; Singh, C.K.; Alsubih, M.; Kahla, N.B. Evaluating the NDVI–Rainfall Relationship in Bisha Watershed, Saudi Arabia Using Non-Stationary Modeling Technique. Atmosphere 2021, 12, 593. [Google Scholar] [CrossRef]
34. Silva, R.M.; Santos, C. A. G.; Maranhão, K. U. A.; Lima, R. F. P. Geospatial assessment of eco-environmental changes in desertification area of the Brazilian semi-arid region. Earth Sciences Research Journal 2018, 22, 175–186. [Google Scholar] [CrossRef]
35. Fokeng, R. M.; Fogwe, Z. N. Landsat NDVI-based vegetation degradation dynamics and its response to rainfall variability and anthropogenic stressors in Southern Bui Plateau, Cameron. Geosystem and Geoenvironment 2020, 1, 100075. [Google Scholar] [CrossRef]
36. Lee, E.; Piao, D.; Song, C.; Kim, J.; Lim, C. H.; Kim, E. Assessing environmentally sensitive land to desertification using MEDALUS method in Mongolia. Forest Science and Technology 2019, 15, 210–220. [Google Scholar] [CrossRef]
37. Budak, M.; Günal, H.; Çelik, I.; Yildiz, H.; Acir, N.; Acar, M. Environmental sensitivity to desertification in northern Mesopotamia; application of modified MEDALUS by using analytical hierarchy process. Arabian Journal of Geosciences 2018, 11, 481. [Google Scholar] [CrossRef]
38. Dutra, J. D.; Elmiro, M. A. T.; Garcia, R. A. Comparative analysis of methods applied in vegetation cover delimitation using Landsat 8 images. Soc. Nat. 2020, 32, 699–710. [Google Scholar] [CrossRef]
39. Liu, S.; Wang, T.; Kang, W.; David, M. Several challenges in monitoring and assessing desertification. Environmental Earth Sciences 2015, 73, 7561–7570. [Google Scholar] [CrossRef]
40. Chen, R.; Yang, H.; Yang, G.; Liu, Y.; Zhang, C.; Long, H.; Xu, H.; Meng, Y.; Feng, H. Land-Use Mapping with Multi-Temporal Sentinel Images Based on Google Earth Engine in Southern Xinjiang Uygur Autonomous Region, China. Remote Sens. 2023, 15, 3958. [Google Scholar] [CrossRef]