Applying Remote Sensing Methods to Estimate Alterations in Land Cover Change and Degradation in the Desert Regions of the Southeast Iberian Peninsula

: Numerous drylands worldwide have experienced degradation of both soil and vegetation in proximity to watering areas. Degradation can be observed in satellite imagery as fading radial brightness belts extending away from the water sources. The main objective of this study was to examine the spatio-temporal patterns of land degradation and rehabilitation in the drylands of the southeast Iberian Peninsula. The brightness index of tasseled cap was discovered to be the best form of spectral transformation for enhancing the contrast between the bright-degraded areas near the points and the darker surrounding areas far from and in between these areas. To comprehend the spatial structure present in spaceborne imagery of two desert sites and three key time periods, semi-variograms were created (mid-late 2000s, around 2015 and 2020). To assess spatio-temporal land-cover patterns, a kriging was used to smooth the brightness index values extracted from 30 m spatial resolution images. To assess the direction and intensity of changes between study periods, a change detection analysis based on kriging prediction maps was performed. These ﬁndings were linked to the socioeconomic situation prior to and following the EU economic crisis. The study discovered that degradation occurred in some areas as a result of the region’s agricultural activities being exploited.


Introduction
Overgrazing, according to the United Nations Environment Program (UNEP), is the practice of allowing a much greater number of animals to graze at a location than it can actually support. As a result, in terms of plant density, plant chemical content, community structure, and soil erosion, overgrazing by various types of livestock is perhaps the most significant anthropogenic activity that degrades rangelands and causes desertification [1].
Overgrazing degrades approximately 75 million ha of land globally, destroying its original biotic functions [2]. According to [3], land (soil and vegetation) degradation is particularly associated with areas surrounding natural or artificial water sources, such as wells or boreholes, in arid and semi-arid environments.
During grazing, domestic animals tend to concentrate near watering points, with the concentration gradually decreasing as the distance from water increases [4,5]. Typically, grazing occurs about 5 km from the watering point, but this can increase up to 20 km in extreme conditions [3,6]. It is also interesting to specify that the radial grazing pattern, known as the "piosphere", was named after the Greek root "pios", meaning "drink" [3]. On the other hand, and according to [7], "grazing gradient" refers to the spatial patterns in soil or vegetation resulting from grazing and can indicate land degradation. For this reason, most of the drylands used for pastoralism show degraded features caused by grazing activities.
Studies conducted using ground measurements and/or analyses of data obtained from remote sources have been conducted to observe the effects of changes in grazing on rangelands. A variety of biotic, abiotic, and environmental effects have been documented in the literature, such as those discussed in references [8][9][10][11]. Studies have looked at how vegetation cover, species richness and diversity, and the spread of biological soil crusts change depending on the distance from the watering point. Further, research has been conducted on soil chemistry and physical properties, as well as the effect of erosion and trampling. The majority of the ground-based measurements are taken from samples along transects, and these samples are taken from plots of varying sizes, ranging from 1 to 100 m 2 . It has been demonstrated that methods which assess only limited areas are heavily influenced by natural changes in vegetation composition and landscape features that cannot be confidently linked to the effects of grazing [12]. Traditional ground-based measurements are laborious and relatively costly, and distant sites are not available for regular repeatable sampling [7]. Ref. [4] postulated that grazing effects cannot be monitored when the measurement is conducted in a limited spatial scope. The challenges associated with semi-arid rangelands are exacerbated, yet remote sensing data can be utilized to counter them.
Since the radial pattern around watering points is clearly visible from satellite imagery, recent studies have been conducted by interpreting and modeling the remote sensing data [13,14], which can be processed in a semi-automated and repeatable manner over large and distant regions. For this reason, different remote-sensing models that utilize geographic information systems techniques have been created to calculate the spatial spread of various elements in the vicinity of watering points [5,[13][14][15].
Various studies have indicated that the circular pattern of grazing gradients around watering points results in similar patterns for various biotic, abiotic, and environmental variables. Vegetation cover, measured using indices such as SAVI (soil-adjusted vegetation index) or NDVI (normalized difference vegetation index) [5,[15][16][17], annuals and grass production [14,18], organic content, soil pH and bush encouragement, nitrate and phosphate [11], soil nutrient concentrations [19], and track density [20], are commonly used variables that follow this pattern. The improvement or decline in each variable does not change significantly a few kilometers from the watering point. However, in some cases, such as inversed, composite, and complex gradients [17,21], different responses were observed. The composite gradient is associated with non-native species replacing desirable ones, while the inverse gradient is observed in areas with dams and mostly woody vegetation. The complex gradient is seen in areas where vegetation growth is decreased in runoff and erosion areas, but increased in run-on and sediment deposits. According to [11], the grazing effects form a radial pattern with a sacrifice zone up to 50 m from the borehole; a nutritious grass zone up to 800 m, dominated by palatable and grass species; and a bush encroachment zone up to 2000 m. This distance is considered as the farthest point at which grazing has a significant impact.
Intensive grazing near water sources is a global occurrence. However, this research was carried out in the arid regions of Southeast Iberian Peninsula, focusing on the socioeconomic changes that took place in this region towards the end of the 20th century. In these desert areas, livestock farming is a significant sector of the economy [22].
In this regard, the importance of goats must be considered, since although they are important for the socioeconomic development of the study area, they are also important factors in the overgrazing process. This fact means that the canopy cannot regenerate in time, and therefore, erosive processes dominate the system [23].
A study conducted in 2006 [24] used remote-sensing methods to examine local desertification processes in the southeast Iberian Peninsula. In this current study, the authors present a different method for evaluating and mapping the impact of grazing around watering points. The main objective of the study was to investigate changes in the vegetation and soil patterns in the drylands of the Southeast Iberian Peninsula over time, in relation to socio-economic changes. Although similar studies have been carried out in other regions around the world, such as China [5,14,21], Australia [6], the South of Africa [11,14], America [14,15,20], the Mediterranean region [16], Eurasian steppes [17], and Tibet [19], this is novel, since it has not been carried out before in the drylands of the southeast Iberian Peninsula. In another vein, it is interesting to highlight that one of the primary features of the southeast Iberian Peninsula drylands is the extensive range of abiotic and biotic surroundings (this is the main difference from other studied areas). These areas are characterized by the prevalence of underdeveloped soils with minimal organic matter, limited aggregate stability, and nutrients, as well as a low water retention capacity, all of which contribute to the exacerbation of drought conditions for vegetation [23].
The study had four specific goals: (1) to analyze the spatial structure of imagery from two desert sites during three different time periods (mid-late 2000s, around 2015, and 2020) using semi-variograms; (2) to apply the kriging interpolation technique to smooth the brightness index values extracted from satellite images with a spatial resolution of 30-80 m, in order to assess spatial and temporal land-cover patterns; (3) to use kriging prediction maps for change detection analysis to determine the direction and intensity of changes between study periods; and (4) to relate these findings to the socio-economic situation before and after the EU economic crisis, which affected grazing intensity and, therefore, the land-use and land-cover state of the study sites.

Study Area
The region under investigation is situated in the southeastern part of Spain, in the Almeria province, as shown in Figure 1. It spans approximately 50 square kilometers and includes Tabernas, the largest town in the area, with a population of 4025 people. The nearest large urban center is the city of Almeria, which has a population of 199,237 and is located 37 km to the south of Tabernas.
The central region of the study area, where human activities are concentrated, is the Tabernas valley. It is surrounded by two mountain ranges running from east to west: the Sierra de los Filabres to the north and the Sierra Alhamilla to the south. In the western portion of the study area lies the Desert of Tabernas, which is regarded as the sole authentic desert in Europe meeting the desert criteria. Protection has been granted to both the Desert of Tabernas and a section of the Sierra Alhamilla as natural reserves [25].
The prevailing land cover types found in the region include vegetation consisting of bushes and grass, with or without trees [26]. Agriculture, primarily barley and some irrigated crops, is the primary form of land use, although it faces significant constraints due to the area's climatic and topographical conditions [27]. Nonetheless, within the past ten years, irrigated olive and almond plantations have been implemented. Currently, mining only occurs in a handful of gypsum quarries, and abandoned mines are prevalent in the area. While industry is not a major contributor to the local economy, the movie and entertainment industry has attracted numerous visitors in recent years. As a matter of fact, tourism is becoming increasingly important and could emerge as one of the primary economic activities in the future.
In another vein, according to the Andalusian Government [28], in the study area, the primary use is very scarce. This fact makes it possible for livestock activity to concentrate in areas where the canopy has reached a higher level of growth. It is evident that overgrazing, as well as livestock trampling and seasonal variation in rainfall regimes, leads to a progressive loss of existing vegetation, and as a result, a progressive loss of soil. Figure 1. Location of the study area (Desert of Tabernas "37°00′00″N; 02°27′00″W" and Sierra Alhamilla "36°59′20″N; 02°21′05″W") with land use maps obtained from the data in [25].
The central region of the study area, where human activities are concentrated, is the Tabernas valley. It is surrounded by two mountain ranges running from east to west: the Sierra de los Filabres to the north and the Sierra Alhamilla to the south. In the western portion of the study area lies the Desert of Tabernas, which is regarded as the sole authentic desert in Europe meeting the desert criteria. Protection has been granted to both the Desert of Tabernas and a section of the Sierra Alhamilla as natural reserves [25].
The prevailing land cover types found in the region include vegetation consisting of bushes and grass, with or without trees [26]. Agriculture, primarily barley and some irrigated crops, is the primary form of land use, although it faces significant constraints due to the area's climatic and topographical conditions [27]. Nonetheless, within the past ten years, irrigated olive and almond plantations have been implemented. Currently, mining only occurs in a handful of gypsum quarries, and abandoned mines are prevalent in the area. While industry is not a major contributor to the local economy, the movie and entertainment industry has attracted numerous visitors in recent years. As a matter of fact, tourism is becoming increasingly important and could emerge as one of the primary economic activities in the future.
In another vein, according to the Andalusian Government [28], in the study area, the primary use is very scarce. This fact makes it possible for livestock activity to concentrate in areas where the canopy has reached a higher level of growth. It is evident that overgrazing, as well as livestock trampling and seasonal variation in rainfall regimes, leads to a progressive loss of existing vegetation, and as a result, a progressive loss of soil.
In reference to the loss of sustainability of the study area, it is of great interest to highlight the difference between the erosive and desertification processes that occur. According to the Andalusian Government [28], the concept of desertification is more functional, and is now being considered as a disturbance that occurs in arid climates and that In reference to the loss of sustainability of the study area, it is of great interest to highlight the difference between the erosive and desertification processes that occur. According to the Andalusian Government [28], the concept of desertification is more functional, and is now being considered as a disturbance that occurs in arid climates and that leads the system (human beings-natural resources) to an irreversible loss of sustainability. On the other hand, erosion consists of the loss of soil by uprooting, transport, and subsequent accumulation, either by the action of wind or water. This phenomenon can be interpreted as a type of soil degradation resulting from desertification, but it is important to note that it is not the sole contributing factor. While wind erosion is the primary form of erosion associated with desertification globally, in the Iberian Peninsula, water erosion takes precedence, representing a prevalent environmental issue across much of Mediterranean Spain. This is particularly evident in the peninsular southeast, which encompasses the Desert of Tabernas and Sierra Alhamilla.
Regarding the climate of the study area, it should be noted, according to the Köppen and Geiger classification, that it is of the BSk type (local steppe).
The main climatic feature of the study area is its Mediterranean character, with mild temperatures and marked aridity. This is because the Betic mountains intersect the Atlantic fronts and leave this area in a rain shadow. The average annual precipitation is 239 mm and the number of rainy days per year ranges from 25 to 55, although only 6% of the rainy episodes exceed 20 mm. The average annual temperature is 17.9 • C and the average minimum of the coldest month is between 3 • C and 10 • C, with the maximum exceeding 40 • C in summer (sometimes reaching 48).
As is known, the study area is made up of loose sediments [29], which are highly salinized and easily washed away by rainwater. This fact, along with the sparse vegetation ( Figure 2c) and torrential rains, has caused the soil in these areas to erode almost completely ( Figure  slopes, erosion is intense and allows few plant species to take root, which is why they are often bare and give the landscape as a whole its characteristic desert appearance (Figure 2a

Materials and Methods
A total of 528 Landsat pictures were obtained from the Earth Explorer platform https://earthexplorer.usgs.gov/ (accessed on 1 October 2022) and utilized to cover the two research locations. These images were captured by the Thematic Mapper (TM) sensor of the Landsat 5 and Landsat 8 satellites at three different time intervals, as indicated in Table  1. The only land-use method in these areas is livestock grazing. The process of image processing commenced by transforming the digital numbers of

Materials and Methods
A total of 528 Landsat pictures were obtained from the Earth Explorer platform https://earthexplorer.usgs.gov/ (accessed on 1 October 2022) and utilized to cover the two research locations. These images were captured by the Thematic Mapper (TM) sensor of the Landsat 5 and Landsat 8 satellites at three different time intervals, as indicated in Table 1. The only land-use method in these areas is livestock grazing. The process of image processing commenced by transforming the digital numbers of the image into reflectance values [30,31]. From these values, various vegetation indices such as NDVI [32], SAVI [33], MSAVI [34], perpendicular vegetation index or PVI [35], and greenness and brightness indices [36] derived from tasseled cap were computed for each image subset. To determine the most appropriate index for distinguishing degraded from non-degraded land, a spectral separability analysis was carried out using the mean and standard deviation values of extreme classes in each scene for all the indices [37][38][39]. The results of this analysis showed that the brightness index (BI) had the highest separability value, and hence, it was selected for further analysis (see Equation (1)).
According to [36,37], BI can be inferred for different sensors using spectral band numbers (B i ) and appropriate BI coefficients (α i ). One of the benefits of using the BI index is its ability to compare different sensors that have different spectral bands. This is possible because the BI values can be normalized to create a single layer of data, making comparisons between sensors possible. Further, the BI index was originally designed to analyze soil properties, and there have been statistically significant results achieved using this method.
The following process involved averaging each window of 6 × 6 pixels in the TM (Landsat 5) and LC (Landsat 8) images, which reduced the resolution by a factor of 6. This resulted in a new pixel size of 171 m, which was small enough to assess changes within the 6 km range from the watering point, where degradation is expected. The subset size was reduced from around 1.5 million pixels to only 40,000 pixels, which made the data more manageable for processing. Note that even though the pixel size was reduced, imageto-image geometric correction was applied for better accuracy in the change detection process.
In another vein, geostatistics is based on the concept of a regionalized variable, which is a variable that can be characterized based on a number of spatial measurements. The fundamental idea behind geostatistics is that when spatial continuity is assumed, adjacent samples are expected to be more similar to each other than samples that are further apart. This spatial dependence can be statistically analyzed and described using parameters derived from a semi-variogram, which is a function that relates the semi-variance to the distance and direction between two samples. The semi-variance, defined as half the mean-squared difference between two samples that are a certain distance apart in a given direction, is used to quantify this spatial dependence. This approach provides a way to analyze and understand the spatial variation in a given variable (Equation (2)).
As is well known, Equation (2), referred to earlier, uses a vector h to define both the direction and the distance between two samples, which is commonly known as the lag.
In Equation (2), the semi-variance at lag h is denoted by γ(h), and N(h) represents the number of pairs of samples that are separated by a distance of h. The value of the regionalized variable at a given location i is represented by Z i . In addition to the lag, the variogram is characterized by three other parameters: the nugget, range, and sill. The nugget accounts for variability at zero distance and is attributable to errors in both sampling and analysis. The range represents the distance, often denoted as "a", beyond which spatial autocorrelation between sampling sites becomes negligible. The sill represents the variability of samples that are spatially independent. Empirical semi-variograms can be computed from a set of observations, and then a theoretical model can be fitted [40]. The exponential model is a popular theoretical model and is expressed in the form shown in Equation (3), where C 0 is the nugget, C 0 + C 1 is the sill, a is the range, and h is the lag.
Other important theoretical models are the Gaussian model (Equation (4)) and the spherical model (Equation (5)). Figure 3 shows the three aforementioned theoretical models. Other important theoretical models are the Gaussian model (Equation (4)) and the spherical model (Equation (5)). Figure 3 shows the three aforementioned theoretical models.
Regarding the Gaussian model, it exhibits a function with opposite curvature close to the origin. As the function progresses, it gradually approaches its maximum value. This model can be considered to have an effective range of around "(3) 1/2 · a", where it encompasses 95% of its maximum variability. However, a significant drawback of the Gaussian model is its tendency to approach the origin with a gradient of zero. According to [40], this indicates the limit of random variation and the point where the underlying variation becomes continuous and twice differentiable. As a result, this behavior may lead to unstable kriging equations, prompting us to discourage the utilization of this particular model.
On a different note, the spherical model appears to be the most apparent choice for characterizing variation in three-dimensional rock formations [40]. However, it may not be as well-suited for describing variation in one or two dimensions, which are typically more relevant for soil and land resource surveys like the one conducted in this study. This function exhibits a gentler curvature at various scales, which is another reason why this model was rejected. Due to what has been commented above, in relation to the Gaussian and spherical models, in the present work, the exponential model will be used.
To better understand the spatial structure of imagery for a given date and location, variogram analysis was used at the studied sites. The decision to use semi-variograms was based on the similarity in spatial structure of most of the variables, which gradually increased or decreased as the distance from the watering point increased, until grazing effects were no longer observed. Semi-variograms were introduced to remote sensing by [41], which found that the variogram parameters could be directly linked to image features. In this study, the presence or absence of a sill in the variogram could indicate the presence of a radial grazing pattern around watering points, while the level of the sill could be used to assess the homogeneity (i.e., lower variance) of an area. The lag distance in the variogram could also be related to the walking distance of livestock, with observations becoming increasingly independent beyond a certain distance (the range). Regarding the Gaussian model, it exhibits a function with opposite curvature close to the origin. As the function progresses, it gradually approaches its maximum value. This model can be considered to have an effective range of around "(3) 1/2 · a", where it encompasses 95% of its maximum variability. However, a significant drawback of the Gaussian model is its tendency to approach the origin with a gradient of zero. According to [40], this indicates the limit of random variation and the point where the underlying variation becomes continuous and twice differentiable. As a result, this behavior may lead to unstable kriging equations, prompting us to discourage the utilization of this particular model.
On a different note, the spherical model appears to be the most apparent choice for characterizing variation in three-dimensional rock formations [40]. However, it may not be as well-suited for describing variation in one or two dimensions, which are typically more relevant for soil and land resource surveys like the one conducted in this study. This function exhibits a gentler curvature at various scales, which is another reason why this model was rejected. Due to what has been commented above, in relation to the Gaussian and spherical models, in the present work, the exponential model will be used.
To better understand the spatial structure of imagery for a given date and location, variogram analysis was used at the studied sites. The decision to use semi-variograms was based on the similarity in spatial structure of most of the variables, which gradually increased or decreased as the distance from the watering point increased, until grazing effects were no longer observed. Semi-variograms were introduced to remote sensing by [41], which found that the variogram parameters could be directly linked to image features. In this study, the presence or absence of a sill in the variogram could indicate the presence of a radial grazing pattern around watering points, while the level of the sill could be used to assess the homogeneity (i.e., lower variance) of an area. The lag distance in the variogram could also be related to the walking distance of livestock, with observations becoming increasingly independent beyond a certain distance (the range).
In each subset of the data, the empirical semi-variogram was calculated using 40,000 pixels, and a theoretical model that best fit the semi-variogram was estimated. The model parameters were then determined by minimizing the squared differences between the empirical semi-variogram values and those of the theoretical model. Once the best-fitting model was chosen, several criteria were applied to assess its accuracy and to fine-tune its parameters: (1) a cross-validation scatter plot was carried out; (2) the mean estimation error was inferred using Equation (6); and (3) the mean standardized squared estimation error was obtained through Equation (7).
In order to minimize the influence of local effects and to obtain a more even representation of surface brightness values across the study sites, researchers employed a linear geostatistical method called ordinary kriging interpolation. This method, according to [42], is useful for revealing spatial phenomena. The approach estimates the mean of the values within a searching neighborhood as a constant.
The study employed the BI differencing method as a post-processing change detection technique, which is a modified version of the vegetation index differencing method [43]. Its purpose was to evaluate the primary patterns of degradation (in magenta color) or rehabilitation (in green color) in the study areas. To achieve this, two sets of kriging maps (i.e., 2016-2005 and 2020-2016) were subtracted for each site (desert of Tabernas and Sierra Alhamilla). Figure 4 displays the relevant BI products, which were computed using Equation (1) and the corresponding coefficients. The bright areas scattered across the mid-to late-2000s images, indicating watering areas, are much less visible in the 2015 and 2020 images. In the more recent 2020 image of the desert of Tabernas and Sierra Alhamilla, many watering areas are absent, but there is a large, bright area present. The brightness index (BI) values, which indicate bare, degraded soil, were calculated using the same method for all sensors, and the brightness levels were uniformly stretched.

Results
The geostatistical analysis utilized all of the BI images. Initially, empirical semivariograms were created for each of the two sites for the three time periods. Field observations did not reveal any anisotropic patterns that could potentially govern the direction of grazing, such as linear dunes or other barriers. Therefore, an isotropic distribution was assumed in all cases. The results of the cross-validation analyses are outlined in Table 2. The chosen slope coefficient and the intercept coefficient of the exponential model were very close to unity and zero, respectively, indicating that the model was able to accurately replicate the observed values.
After the examination of several theoretical models [40], the exponential model (Equation (3)) was chosen as it yielded the best results during cross-validation. This choice is due, as can be seen in Table 2, to the fact that the exponential model has a lower spatial variance than the Gaussian and spherical models. In addition, it should be noted that at smaller pixel sizes [40], the exponential model fit better than the Gaussian and spherical models, the latter being the one that would best fit for larger pixel sizes (since the sill value would be lower) as long as data from the MODIS satellite had been used, for example.
Later, kriging interpolation techniques were utilized, employing exponential models (see Figure 5) with the parameters from Table 2. The final outcomes illustrating the distribution of BI values for the three periods and both sites are displayed in Figure 6. In Figure 6,  The geostatistical analysis utilized all of the BI images. Initially, empirical semi-variograms were created for each of the two sites for the three time periods. Field observations did not reveal any anisotropic patterns that could potentially govern the direction of grazing, such as linear dunes or other barriers. Therefore, an isotropic distribution was assumed in all cases. The results of the cross-validation analyses are outlined in Table 2. The chosen slope coefficient and the intercept coefficient of the exponential model were very close to unity and zero, respectively, indicating that the model was able to accurately replicate the observed values.  After the examination of several theoretical models [40], the exponential model (Equation (3)) was chosen as it yielded the best results during cross-validation. This choice is due, as can be seen in Table 2, to the fact that the exponential model has a lower spatial variance than the Gaussian and spherical models. In addition, it should be noted that at smaller pixel sizes [40], the exponential model fit better than the Gaussian and spherical models, the latter being the one that would best fit for larger pixel sizes (since the sill value would be lower) as long as data from the MODIS satellite had been used, for example.
Later, kriging interpolation techniques were utilized, employing exponential models (see Figure 5) with the parameters from Table 2. The final outcomes illustrating the distribution of BI values for the three periods and both sites are displayed in Figure 6. In Figure  6, the desert of Tabernas has BI values ranging from 0.66 to 0.95, while the Sierra Alhamilla sites have BI values ranging from 0.20 to 0.60.  The BI differencing method, explained in the concluding paragraph of Section 3, was utilized to analyze changes in grazing patterns over two specific time periods.  The BI differencing method, explained in the concluding paragraph of Section 3, was utilized to analyze changes in grazing patterns over two specific time periods. Figure 7 displays the maps illustrating these changes.  Note that the gray color of the utilized palette was enlarged on the value scale only so that it could be properly appreciated, even though its value is 0.

Discussion
Regarding the BI maps (Figure 4), the brightness of the Sierra Alhamilla in the mid-tolate 2000s was significantly different (or opposite) from the other two maps of the same region. This fact is due to a twofold reason. The first one refers to the topography of those areas of Sierra Alhamilla in which the brightness is different. However, according to [25], at the beginning of 2015, a reforestation of the bright areas was carried out in the mid-to-late 2000s in order to minimize the soil loss due to surface runoff phenomena. As a result of this fact, for around 2015 and for 2020, said bright area appears with a dark hue.
The fit of the exponential model was assessed using the least-squares measure, and the variogram parameters can be found in Table 2. Each variogram (see Figure 5) was processed with 17 lags, each spanning a distance of 800 m. The lag values were determined through a trial and error process to optimize the aforementioned criteria. When observed visually, the variograms from the mid-late 2000s and around 2015 appear quite similar, exhibiting a typical variogram shape. Two notable features can be observed. Firstly, the sill of the variogram around 2015 was lower compared to the mid-to-late 2000s, suggesting a decrease in variance in the latter year. Further, the range in the mid-to-late 2000s was 4800 m, while in around 2015, it was 6400 m, which corresponds to the reported walking distance of livestock from their drinking source [44]. In contrast, the variogram for 2020 only reached a similar sill level after a distance of 12 km. This range does not appear to be indicative of the grazing pattern in the region.
In Figure 6, the initial maps, specifically during the mid-late 2000s, reveal the presence of bands encircling the watering points, signifying the progressive degradation of land extending from the wells. These bands denote the grazing gradient or the impact of grazing. Regions depicted in dark red (values of 0.95 for the Desert of Tabernas and 0.6 for Sierra Alhamilla) in the images indicate watering areas. The adjacent belts in light red (values of 0.92 in the Desert of Tabernas Desert, and 0.55 in Sierra Alhamilla) indicate the limit from which the impact caused by grazing pressure is significant around the watering areas, exerting a strong influence on spatial variation. Such regions are referred to as the "sacrifice zone" by [11]. On the other hand, the orange and yellow colors represent a mixed zone where the effects of grazing and natural variability coincide or establish a stable balance. This zone can be likened to the outskirts of the grazing impact and highlights the primary migration routes of livestock. The zone exhibiting green shades signifies an area where natural variability surpasses the impact of grazing, and was termed the "grazing reserve" by [11].
On the other hand, in Figure 7, and specifically in the change detection map for the period of 2016-2005, a large portion of the area appears grey, indicating that there were no noticeable alterations in grazing patterns. The remaining areas are represented in magenta, indicating degradation, and green tones, representing rehabilitation. Consequently, between the mid-to-late 2000s and around 2015, the overall trend was a positive one, with land-cover conditions improving or undergoing rehabilitation. It is worth noting that around 2015, the area became more uniform, as indicated by the reduced variation (lower sill) in Figure 5 and the narrower range of colors in Figure 6a,b. Conversely, during the second period (2020-2016), a significant portion of the area experienced degradation processes (magenta colors), while the area undergoing rehabilitation (green tones) diminished considerably (Figure 6a,b).
As a result of the 2007 economic crisis, there was a progressive decrease in the livestock census in the study area due to increased costs, which is why many farmers had to abandon economic activity due to lack of profitability. This fact had a positive influence on the environment, since it led to the recovery of those areas where livestock pressure was highest. This is one of the reasons why, from the mid-to-late 2000s to around 2015, a recovery of the study area can be observed (see Figure 7a). Despite the above, and although, in 2016, the Andalusian Government started the Recovery Plan for the Desert of Tabernas and Sierra Alhamilla [45], the increase in degradation (Figure 7b) is mainly due to climatic causes.
To explain the climatic effect in the study area, Table 3 shows the data corresponding to the sum of both the average monthly temperature (in Celsius) and the monthly precipitation (in mm) in each of the years subjected to study. All climatic data were downloaded from Version 4 of the CRUTS monthly high-resolution gridded multivariate climate dataset [46]. From Table 3, as a result of climatic peculiarities of the study area, the need to infer a new aridity index was seen, and it was obtained to help explain the results shown in Figure 6. We have called this index the RAMALL (RAmírez, MAdueño, López, and Leiva) Aridity Index. Equation (6) specifies how the RAMALL index is calculated.
where: i = months of each year; j = year for which to calculate RAMALL aridity index; MP j = monthly precipitation (mm) in the year "j"; AMT j = average monthly temperature ( • C) in the year "j"; AMT j+1 = average monthly temperature ( • C) in the year "j + 1". Values greater than 1 indicate a low aridity index, while values lower than 1 indicate a high aridity index. The reason for using temperature and precipitation data is based on the physical principle that evapotranspiration increases with temperature and, therefore, precipitation becomes less effective. However, unlike methods based on evapotranspiration (such as the Meigs aridity index [47] or the UNEP aridity index [48], among many others), the proposed index (RAMALL aridity index) takes into account the effect of the soil on the increase in environmental temperature close to the surface, as is specified by [39]. Figure 8 shows the evolution of this new RAMALL aridity index in the study area between 2005 and 2020.
After analyzing the data shown in Table 3, as well as Figure 8, the RAMALL aridity index presented mean values of 1.05 (low aridity index) between 2005 and 2015 and 0.94 (high aridity index) between 2016 and 2020.
As can be seen, the 2007 economic crisis had a positive effect on the average aridity during the 2016-2005 period (the RAMALL aridity index value was greater than 1). During this period, the water erosion phenomena [28] were minimized not only by the decrease in livestock pressure, but also as a consequence of the increase in plant variety that, being present in the soil seed bank, only colonized the area study as a result of said decreased grazing.
However, in the second period under study (2020-2016 period), for which RAMALL aridity index values were lower than 1, the progressive increase in the average annual temperature (Figure 9), along with the near-disappearance of grazing that occurred at the end of the first period (2016-2005), favored plant variety loss in the study area, and therefore, the shrub stratum was dominant. On the other hand, it should be noted that in 2014 (see Table 3), there was a decrease in rainfall levels, which, along with the progressive increase in temperature (see Figure 9), contributed even more to the desertification process in the study area. This fact has, as a serious consequence for the environment of the study area, led to an increase in erosive processes, and will also lead to both soil and habitat loss over time if the global warming trend continues [28]. After analyzing the data shown in Table 3, as well as Figure 8, the RAMALL aridity index presented mean values of 1.05 (low aridity index) between 2005 and 2015 and 0.94 (high aridity index) between 2016 and 2020.
As can be seen, the 2007 economic crisis had a positive effect on the average aridity during the 2016-2005 period (the RAMALL aridity index value was greater than 1). During this period, the water erosion phenomena [28] were minimized not only by the decrease in livestock pressure, but also as a consequence of the increase in plant variety that, being present in the soil seed bank, only colonized the area study as a result of said decreased grazing.
However, in the second period under study (2020-2016 period), for which RAMALL aridity index values were lower than 1, the progressive increase in the average annual temperature (Figure 9), along with the near-disappearance of grazing that occurred at the end of the first period (2016-2005), favored plant variety loss in the study area, and therefore, the shrub stratum was dominant. On the other hand, it should be noted that in 2014 (see Table 3), there was a decrease in rainfall levels, which, along with the progressive increase in temperature (see Figure 9), contributed even more to the desertification process in the study area. This fact has, as a serious consequence for the environment of the study area, led to an increase in erosive processes, and will also lead to both soil and habitat loss over time if the global warming trend continues [28]. Finally, it is necessary to emphasize that, both in the Desert of Tabernas and in Sierra Alhamilla, the influence of the arid climate has given rise, over time, to a landscape of mainly erosive morphology, whose main signs are slopes (of soft materials) in the form of gullies. With respect to the drainage network, it is necessary to indicate the action of the existing channels on the smoothing of the relief, resulting in a predominance of erosive Finally, it is necessary to emphasize that, both in the Desert of Tabernas and in Sierra Alhamilla, the influence of the arid climate has given rise, over time, to a landscape of mainly erosive morphology, whose main signs are slopes (of soft materials) in the form of gullies. With respect to the drainage network, it is necessary to indicate the action of the existing channels on the smoothing of the relief, resulting in a predominance of erosive action at the headwaters and, therefore, boxing in of the channels, as well as a predominance of sedimentation in the lower part. The impact caused by raindrops on bare soil gives rise to a wide variety of microforms in the landscape, a key characteristic of soil loss through surface runoff.

Conclusions
The aim of the present study was to analyze the changes in vegetation patterns over time in the Desert of Tabernas and Sierra Alhamilla (Almería, Spain), focusing on land degradation and rehabilitation. Specifically, our research aimed to understand how these changes were influenced by socio-economic factors following the economic crisis in 2007. To achieve this, the tasseled cap-derived brightness index (BI) was employed as a tool to depict the spatial distribution of land surface characteristics. The BI was chosen for its superior ability to differentiate among various spectral indices, its original purpose of examining soil properties, and its capability to compare different sensors with distinct spectral bands by normalizing them into a single BI layer.
The use of geostatistical analysis, specifically semi-variance analysis, was determined to be a suitable approach for understanding the spatial structure within imagery at a specific date and location. This method was chosen due to the similarity between the variograms' shapes and the directional changes observed in various biotic, abiotic, and environmental variables along the grazing gradient emanating from watering points in arid and semi-arid regions. By utilizing variogram models from the mid-to-late 2000s and around 2015, it becomes possible to quantitatively generalize the observed phenomena across the region of interest. The average range of 4800-6400 m aligns with the reported walking distance of livestock from their water sources. Comparing variograms from different years can provide insights into the temporal dynamics of the region, as evidenced by the failure to reach a sill in the variogram when grazing ceased at the end of 2015.
The kriging interpolation method was employed as a smoothing filter, whereby each pixel in the image was replaced with a solution derived from the variogram equation (specifically, an exponential model in this case), which was calculated using all other pixels. This approach aimed to reduce spatial errors and fine-scale variability, facilitating a more accurate delineation of degradation boundaries surrounding the watering points. The resulting maps of the study area revealed a radial pattern, indicating progressive land degradation emanating from the wells, commonly referred to as the grazing gradient. However, this pattern appeared blurred or absent in the 2020 maps. In comparison to other interpolation techniques, ordinary kriging is regarded as the superior linear unbiased estimator.
The index differencing technique was successfully utilized to analyze temporal changes. This study showcases the capability of satellite image analysis to track land-use and landcover modifications resulting from both the 2007 economic crisis and the significant decrease in grazing pressure.
In future work, it could be interesting to compare the existing similarity between the Desert of Tabernas and Sierra Alhamilla with the existing desertification in other areas of Mediterranean Europe. This study may be important in order to highlight the possible peculiarities that favor the preservation of native fauna and flora. Funding: This research received no external funding.

Data Availability Statement:
The data used to support the findings of this study can be made available by the corresponding author upon request.