The Resilient Recurrent Behavior of Mediterranean Semi-Arid Complex Adaptive Landscapes

Growing external pressures from human activities and climate change can exacerbate desertification, compromising the livelihoods of more than 25% of the world’s population. The dryland mosaic is defined by land covers that do not behave similarly, and the identification of their recurring or irregular changes over time is crucial, especially in areas susceptible to become desertified. To this aim, the methodological approach of this research is based on the integration of non-linear data analysis techniques, such as recurrence plots (RPs) and recurrence quantification analysis (RQA), applied to the Enhanced Vegetation Index (EVI), which is a functional ecological proxy of above ground net primary production. The research exploits the recurring change detected in vegetation cover over time to gauge the predictable (resilient) behavior of the EVI as well as its chaoticity in a semi-arid Mediterranean region (Apulia, Italy). Interestingly, the results have shown the spatial rendering of recurrence variables, confirming the well-known hot spots of soil degradation and desertification taking place in the region, which are characterized by greater EVI chaoticity, but they have also identified new potential candidate sites. As a result, the susceptibility to land degradation, as measured by the EVI-RQA approach, can help in measuring land desertification with evident operational benefits for landscape planning. The novelty of the research lies in the spatially explicit identification of resilient and less resilient areas to desertification that can support the definition of more targeted interventions and conservation priorities for better planning and sustainable management of Mediterranean drylands.


Introduction
Complex adaptive systems [1] such as socio-ecological landscapes (SELs) [2,3] show generally nonstationary and complex behavior and are usually assumed to respond smoothly to gradual change in climate, habitat fragmentation or exploitation. Such a condition can be sporadically interrupted by abrupt shifts towards a radically different regime [4]. Resilient SELs can adapt the behavior of the system (adaptability) to changes in external drivers and internal processes and remain in the current domain of stability or basin of attraction [2].
Research and practice in SELs are increasingly converging on a series of a few general lessons as they apply to sustainable development in drylands [5]: (1) the need to adopt an integrated approach, (2) a greater awareness of slowly changing conditions, (3) the recognition of non-linear processes since dryland systems are not in equilibrium, and Land 2021, 10, 296 2 of 18 therefore often show multiple ecological and social states, and (4) cross-scale interactions must be expected. In our approach, we have tried to address these lessons and to synergize them in the analytical techniques adopted.
Land degradation has been defined by the Millennium Ecosystem Assessment (MEA) as "a persistent net loss of capacity to yield provisioning, regulating, and supporting ecosystem services" [6]. The MEA refers to desertification as the same degradation process in drylands (i.e., semi-arid, arid and dry sub-humid zones). In this context, landscape services [7] directly or indirectly sustain the quality of human life and are provided by physical structures and processes called "supporting services" [6,8], among which the above-ground net primary production (ANPP) is a key-controlling factor determining most of the landscape service provisioning capacities [6]. Such essential supporting services influence the entire ecosystem functioning, governing most of provisioning and regulating services (e.g., food, water cycling and regulation, nutrient cycling) and even some cultural services, which tend to shift in concert with it [8][9][10]. Thus, a reduction in ANPP at a site is often viewed as land degradation [5] that, if exacerbated, can lead to a state of irreversible desertification [11]. As such, efforts to reduce the risk of undesirable state transitions such as land degradation and desertification should include an assessment of gradual or irregular changes or transitions in ANPP, rather than merely trying to control disturbances [4].
Vegetation cover is the primary determinant of SEL mosaics and is quite responsive to land cover transformations at all scales, as all climate changes and anthropogenic influences have a spatial outcome. Therefore, current integrated approaches to monitoring and assessing land degradation and desertification (LDD) often rely on two technological aspects: (1) the remote sensing of vegetation cover and (2) the analysis of satellite data (vegetation indices) through specific tools, such as geographic information systems (GIS) and time series analysis techniques [12].
However, while most of the integrated approaches to assess LDD focus mainly on the selection of one or few vegetation indices and try to link the generated spatial information with a variety of environmental and socio-economic data in a single framework, less emphasis is placed on the adequacy and effectiveness of time series analysis techniques to extract useful information from the analysis of time series behavior [12,13].
Since natural or human-dominated processes can have a distinct recurrent behavior, e.g., periodicities (as seasonal or Milankovich cycles), but also irregular cycles like ENSO (the El Niño-Southern Oscillation), the study of recurrence is becoming essential as it is a fundamental property of deterministic dynamical systems to be studied, and it is also typical of non-linear or chaotic systems [14]. It allows the measurement of resilience in terms of stability/predictability, which can be defined as the probability of a system's multiple states to persist and recur [15].
In this perspective, the purpose of this research is to demonstrate the effectiveness of the recurrence quantification analysis (RQA) [14] to provide for the region of interest (a) the estimate of the recurring resilience-related variables from the behavior of the Enhanced Vegetation Index (EVI) time series over the entire 14-year period, (b) the spatialization of the recurrence variables on the regional scale to identify hot spots with different stability/predictability as well as chaoticity of space-time functional dynamics because of climate change and/or water scarcity or human activities, and (c) the validation of this mapping with the confirmation of the well-known land degradation and desertification hot spots in the region as well as the identification of new sites as potential candidates to become desertified. To this aim, the available functional dynamics of vegetation covers in a Mediterranean semi-arid region have been analyzed in terms of possible recurrent and, therefore, predictable behaviors through: (a) time series of the Enhanced Vegetation Index (EVI), which is an ecological functional proxy of ANPP [16,17] and its changes can reflect the effects of changes associated with climate, land conversion, and farming practices as well as drought and fire, and (b) the use of non-linear data analysis techniques such as Land 2021, 10, 296 3 of 18 recurrence plots (RPs) and recurrence quantification analysis (RQA) [14] applied to the average EVI time series extracted from satellite multitemporal images.
It is interesting to note that these techniques, unlike those normally used, do not involve any manipulation of the original data and are fully independent from binding hypotheses and signal transformations. They are also applicable to relatively short time series and provide resilience-related recurrence variables that refer to stability/predictability and, conversely, chaoticity (the condition of being chaotic) [14]. We argue that the greater the stability/predictability of time series belonging to a land cover class, as shown by the recurrence of states within preferred bounds (adaptability), the higher the likelihood that such a class will be resilient to contagious and non-contagious stress, often showing, but not necessarily, a regular recurrent behavior [18]. Conversely, all the unfavorable biophysical factors that promote desertification in Mediterranean drylands, such as erratic precipitation (mainly during the winter), high summer temperatures with frequent drought events, poor and erodible soils, extensive deforestation with frequent wildfires and land abandonment [19], are expected to reverberate on the chaotic EVI temporal behavior, as highlighted by RP-RQA.

Study Area
Dryland landscapes are, in general, characterized by water scarcity and high climate variability, with strong weather extremes. In the Mediterranean basin, drylands are characterized by a strong multifunctionality since they coevolved with humans for millennia [20,21]. These landscapes are characterized by a progressive land degradation with the consequent loss of crucial services such as food production, soil quality regulation, water regulation, climate regulation, identity value, and people's sense of belonging to semi-arid landscapes and their associated recreational services [22,23].
The study area ( Figure 1) represents a typical Mediterranean semi-arid region [24], where more than 82% of the area is constituted by agroecosystems, critically dependent on water availability [25], where extensive forest clearances have occurred since the middle Holocene [26,27], slowly leading to desertification [28,29]. well as drought and fire, and (b) the use of non-linear data analysis techniques such as recurrence plots (RPs) and recurrence quantification analysis (RQA) [14] applied to the average EVI time series extracted from satellite multitemporal images. It is interesting to note that these techniques, unlike those normally used, do not involve any manipulation of the original data and are fully independent from binding hypotheses and signal transformations. They are also applicable to relatively short time series and provide resilience-related recurrence variables that refer to stability/predictability and, conversely, chaoticity (the condition of being chaotic) [14]. We argue that the greater the stability/predictability of time series belonging to a land cover class, as shown by the recurrence of states within preferred bounds (adaptability), the higher the likelihood that such a class will be resilient to contagious and non-contagious stress, often showing, but not necessarily, a regular recurrent behavior [18]. Conversely, all the unfavorable biophysical factors that promote desertification in Mediterranean drylands, such as erratic precipitation (mainly during the winter), high summer temperatures with frequent drought events, poor and erodible soils, extensive deforestation with frequent wildfires and land abandonment [19], are expected to reverberate on the chaotic EVI temporal behavior, as highlighted by RP-RQA.

Study Area
Dryland landscapes are, in general, characterized by water scarcity and high climate variability, with strong weather extremes. In the Mediterranean basin, drylands are characterized by a strong multifunctionality since they coevolved with humans for millennia [20,21]. These landscapes are characterized by a progressive land degradation with the consequent loss of crucial services such as food production, soil quality regulation, water regulation, climate regulation, identity value, and people's sense of belonging to semiarid landscapes and their associated recreational services [22,23].
The study area ( Figure 1) represents a typical Mediterranean semi-arid region [24], where more than 82% of the area is constituted by agroecosystems, critically dependent on water availability [25], where extensive forest clearances have occurred since the middle Holocene [26,27], slowly leading to desertification [28,29].  and a temperature increase, with its maximum in the summer leading to an increased frequency of high-temperature events [30,31]. The main landscape changes responsible for land desertification are given by the merging of arable lands [32], a contraction of vineyards, and the expansion of olive groves, plantations, and urban areas, with large increases in energy use, water consumption, and fertilizer applications, leading to extensive biodiversity loss ( Figure 1) [33]. As a result, the main vulnerability of the region relates to its persistent water scarcity and its water budget deficit (about 350 mm/year), requiring regular water imports from nearby regions and forcing the exploitation of aquifers [34,35].
The introduction of new farm machinery and strains of cereals and trees, along with the overexploitation of water resources for irrigated agriculture and the extensive application of fertilizers, represent the main socio-economic drivers of soil erosion and LDD in Mediterranean Europe [35] as well as in Apulia [29]. The net irrigation requirements are expected to increase in the next century, with a maximum of 65% for intensive olive groves due to the excessive pumping of groundwater [25]. This will increase water stress by the intensity, extent, timing, and duration of water availability [36] as well salinity intrusion in groundwater, one of the main drivers for the salinization of coastal soils.
The northern and central parts of the region are dominated by arable lands, whereas the central and southern parts of the region are dominated by extensive century-old as well as intensive olive groves, and other rural areas ( Figure 1). The few natural habitats are unevenly distributed as forested areas concentrated in the Gargano peninsula, as well as a few remnants of evergreen forests interspersed with olive groves. In addition, natural drylands are present with shrub and/or herbaceous vegetation and permanent pastures (Figure 1), including semi-arid grasslands, garigues, steppes and pastures, all collected here under the name of Grasslands, present in Daunia and Murgia Apennines.

Materials: Real World Time Series of Vegetation Indices
The normalized difference vegetation index (NDVI) has been widely used over the last few decades in the satellite-based modeling of ANPP for terrestrial vegetation [37][38][39]. However, several limitations have been observed when applying the NDVI index in dry land area research [23,40]; therefore, in this research, the enhanced vegetation index (EVI) [40] was preferred as it has been used in the satellite-based modeling of ANPP [41], is found more linearly correlated with the green leaf area index (LAI) in crop fields [42], is less prone to greenness saturation and has higher robustness to atmospheric conditions and soil background relative to NDVI [13,40,43]. Moreover, perennial and annual vegetation types that green up in deserts at rates and times characteristic for their distinct species can be easily captured by EVI data [44].
Integrated approaches based on time series of vegetation indices have mostly used linear techniques focusing on different aspects of temporal complexity. Recently, the "normalized spectral entropy", a non-linear analysis, has been developed to describe the degree of regularity (orderliness) and to characterize heterogeneity within an ecological time series based on its power spectrum [18,21,45], pointing to the system's self-organization strength [46]. However, because of non-stationarity, noise, relatively short time series [47] and frequently complex shapes of recurrent periodicities, more advanced methods are needed for quantifying non-linear dynamics [21,48].

Methods
The methodology is based on two steps: First step: satellite data harmonization and cross-correlation with temperature and precipitation.
The Moderate Resolution Imaging Spectroradiometer (MODIS, https://modis.gsfc. nasa.gov/ accessed on 25 June 2018), onboard the NASA Terra and Aqua satellites, provides both NDVI and EVI images with high temporal frequency and medium spatial resolution, allowing us to gauge the continuous flow of ANPP. MODIS EVI images (USGS, https: //glovis.usgs.gov/ accessed on 25 June 2018) are already filtered from daily data, although Land 2021, 10, 296 5 of 18 further filtering was necessary to keep the most reliable EVI. To this end, a simple filter (to avoid high drops of EVI between two high EVI values) was applied as in [49], and then, to build the most reliable EVI profile, a simple gap-fill approach has been used to implement the method of Xiao et al. [16,43].
Finally, cross-correlation analyses of EVI with temperature and precipitation were performed on a ten-day basis, acquired from the Italian National Environmental Information System (http://www.sinanet.isprambiente.it/it/sia-ispra accessed on 25 June 2018), and related to six weather stations distributed throughout the Apulia region inside or in proximity to the types of land cover analyzed in this study. This was performed to evaluate how the EVI and associated RAQ metrics were related to precipitation and temperature. Understanding this is central to being able to pull out signals of other perturbations, as well as responses to drought.
A brief overview of RP and RQA is reported below, but more details of the principles and algorithms can be found in [14,21].
A landscape state consists of state variables describing all the important properties of dynamical systems; therefore, it can be mathematically depicted by a vector → x i = u i , u i+τ , . . . , u i+(m−1)τ with embedding dimension m and time delay τ. The temporal evolution of the state is expressed by the phase space trajectory of the dynamical system. Formally, the RP is based on the recurrence matrix: where N is the number of measured states → x i , i.e., the length of the time series reduced by the number (m − 1) τ; ε is a threshold distance; and ||·|| is a norm, i.e., the spatial distance between two points x i and x j in the phase space trajectory. The choice of m is usually based on counting false nearest neighbors when increasing m and choosing the value of m where the number of false nearest neighbors goes to zero [14]. The delay τ must be selected to minimize the autocorrelations between the time series points, e.g., by using mutual information. The parameter ε is crucial as systems often do not recur exactly but just approximately to a formerly visited state. The appropriate threshold ε can be selected by following three methods, i.e., mean or maximum phase space diameter, recurrence point density, and standard deviation of the observational noise [14].
The RP is a graphical representation of this recurrence matrix (Figure 2), where black points represent those time points where the spatial distance between two states → x i and → x j is falling below the threshold ε and, therefore, the system recurs [14,54]. For illustrative purposes, the RPs of three typical synthetic time series (A, B, C) are presented ( Figure 2). In particular, the trajectory with regular oscillations ( Figure 2B) is the least common in nature, but the possibility given by RQA to analyze more complex but recurrent oscillation trends represents the strength of the RQA approach. Diagonal lines in the RP indicate that the temporal evolution of states is similar at different times and, therefore, can point to deterministic processes; vertical or horizontal lines in the RP indicate that some states do not change or change very slowly for some time and can suggest laminar or persistent states. just approximately to a formerly visited state. The appropriate threshold ε can be selected by following three methods, i. e., mean or maximum phase space diameter, recurrence point density, and standard deviation of the observational noise [14].
The RP is a graphical representation of this recurrence matrix (Figure 2), where black points represent those time points where the spatial distance between two states ⃗ and ⃗ is falling below the threshold ε and, therefore, the system recurs [54], [14]. For illustrative purposes, the RPs of three typical synthetic time series (A, B, C) are presented ( Figure  2). In particular, the trajectory with regular oscillations ( Figure 2B) is the least common in nature, but the possibility given by RQA to analyze more complex but recurrent oscillation trends represents the strength of the RQA approach. Diagonal lines in the RP indicate that the temporal evolution of states is similar at different times and, therefore, can point to deterministic processes; vertical or horizontal lines in the RP indicate that some states do not change or change very slowly for some time and can suggest laminar or persistent states. Periodic motion is reflected in the RP by uninterrupted diagonals, and the vertical distance between them corresponds to the period of the oscillation. These uninterrupted diagonals represent the "predictability time", i. e., the time that two trajectories are nearby before diverging out in state space, and the longest diagonal represents the maximum time.
We first estimate the value of τ from the first minimum of the average mutual information and space-time separation plots, and then m for ≤5% false-nearest neighbors. The estimated minimum embedding dimension is the one that initially gives a significant fraction of false nearest neighbors with the threshold ε, which in our case study varies from 0.08 to 0.13 among land covers.
RQA allows the estimation of recurrence variables [48], all concurring to provide a quantitative picture of spatiotemporal dynamics. In particular, the recurrence rate (REC) is the ratio of all recurrent states, namely recurrence points closer than ε to all possible states and, therefore, the probability of the recurrence of a certain state. The ratio of recurrence points forming diagonal structures to all recurrence points is called the determinism (DET) and gives indications about the predictability of service provisioning by the dynamical landscape. The laminarity (LAM) is the proportion of recurrence points forming vertical lines and is related to the number of laminar phases in the system (intermittency). LAM reveals laminar states, i. e., where some states do not change or change slowly and, hence, the state is stable for some time. The divergence (DIV) is simply 1/Lmax, where Lmax is the maximum predictability time given by the length of the longest diagonal line in the plot excluding the main diagonal line. DIV gauges the rate at which trajectories diverge and is the hallmark for dynamic chaos. The higher the DIV, the more chaotic (less stable) the signal is [14]. Periodic motion is reflected in the RP by uninterrupted diagonals, and the vertical distance between them corresponds to the period of the oscillation. These uninterrupted diagonals represent the "predictability time", i.e., the time that two trajectories are nearby before diverging out in state space, and the longest diagonal represents the maximum time.
We first estimate the value of τ from the first minimum of the average mutual information and space-time separation plots, and then m for ≤5% false-nearest neighbors. The estimated minimum embedding dimension is the one that initially gives a significant fraction of false nearest neighbors with the threshold ε, which in our case study varies from 0.08 to 0.13 among land covers.
RQA allows the estimation of recurrence variables [48], all concurring to provide a quantitative picture of spatiotemporal dynamics. In particular, the recurrence rate (REC) is the ratio of all recurrent states, namely recurrence points closer than ε to all possible states and, therefore, the probability of the recurrence of a certain state. The ratio of recurrence points forming diagonal structures to all recurrence points is called the determinism (DET) and gives indications about the predictability of service provisioning by the dynamical landscape. The laminarity (LAM) is the proportion of recurrence points forming vertical lines and is related to the number of laminar phases in the system (intermittency). LAM reveals laminar states, i.e., where some states do not change or change slowly and, hence, the state is stable for some time. The divergence (DIV) is simply 1/Lmax, where Lmax is the maximum predictability time given by the length of the longest diagonal line in the plot excluding the main diagonal line. DIV gauges the rate at which trajectories diverge and is the hallmark for dynamic chaos. The higher the DIV, the more chaotic (less stable) the signal is [14].
Recurrence variables can be interpreted in terms of landscape adaptability; a few of their corresponding theoretical descriptors are given in Table 1  Table 1. Relationships between some recurrence variables (stability metrics) and theoretical descriptors of landscape adaptability.

Recurrence Variable-Stability Metrics System's Adaptability Theoretical Descriptors
Determinism (DET): the percentage of recurrence points, which form diagonal lines. It gauges the stability and predictability of the dynamical system [14].
Stability: concentrates on the adaptive capacity of the system to remain within the same domain because of mutually reinforcing structures and processes [56].
Resilience: the probability that a system's multiple states will persist [15].
Laminarity (LAM): the percentage of recurrence points, which form vertical lines (laminar phases) in the system where some states do not change or change slowly (intermittency) and, hence, the state is trapped for some time [14].
Constancy: when the state of the system remains the same [58].

Divergence (DIV)
: the inverse of the length of the longest diagonal line (L max). It gauges the rate at which trajectories diverge and is the hallmark for critical transitions and dynamic chaos [14].
Precariousness: the proximity or trajectory of an ecosystem state to a threshold that if crossed will result in a regime shift or critical transition [4,59]. Chaoticity: The condition of being chaotic.

The longest diagonal line in the plot (L max):
the maximum predictability time; the smaller the divergence, the higher the system's trajectory stability and predictability [14].
Stability: qualities of a stability domain that resist changes in state or shifts in regime [60].
Resilience: the probability that a system's multiple states will persist [15].
Recurrence plots (RPs) and recurrence quantification analysis (RQA) were performed using the R Package "fNonlinear" [61]. Confidence intervals for the recurrence variables were derived from RQAs of 100 random Enhanced Vegetation Index (EVI) time series extracted from each land cover.

EVI Time Series and Cross-Correlations with Climate Data
Mean EVI time series of the Apulia region provide a picture of land cover types fluctuating around some trends (increasing or decreasing) or stable average, but with clear differences in the regularities of periodicities and amplitudes of the processes (Figure 3). Broadleaf Forests present the highest average EVI with clear periodicities, with a maximum EVI in late spring (May) and a minimum in winter (February). In contrast, Olive Groves show irregular periodicities with lower amplitude, with a maximum EVI in winter (January) and a minimum in summer (July). Arable Lands show periodicities that clearly reflect the cycle of agricultural practices such as sowing, irrigation and growing, and harvesting, with a maximum EVI in spring (April) and a minimum in July. Grasslands, including semiarid grasslands, garigues, steppes and pastures, show irregular periodicities with higher amplitude than Olive Groves but with a similar average EVI, with reduced photosynthetic activity in the summer.  Cross-correlations of the EVI with temperature and precipitation demonstrate that most of the time series are consistent with a climatic forcing and variability ( Figure 4). However, cross-correlations for Forests are antithetical to those of other types of land cover as forests' greater self-organization makes them almost independent of precipitation events. In contrast, there is an evident direct correlation between the EVI of Forests and temperature. For Arable Lands, the cross-correlations with both precipitation and temperature are negative, highlighting that water availability is not dependent on precipitation events, but on groundwater for irrigation. The cross-correlation between precipitation and Olive Groves behaves similarly to Grasslands, which shows a strong dependence on precipitation events. However, Olive Groves are negatively related to temperature, with the highest EVI in winter. Cross-correlations of the EVI with temperature and precipitation demonstrate that most of the time series are consistent with a climatic forcing and variability ( Figure 4). However, cross-correlations for Forests are antithetical to those of other types of land cover as forests' greater self-organization makes them almost independent of precipitation events. In contrast, there is an evident direct correlation between the EVI of Forests and temperature. For Arable Lands, the cross-correlations with both precipitation and temperature are negative, highlighting that water availability is not dependent on precipitation events, but on groundwater for irrigation. The cross-correlation between precipitation and Olive Groves behaves similarly to Grasslands, which shows a strong dependence on precipitation events. However, Olive Groves are negatively related to temperature, with the highest EVI in winter.

RP and RQA Analysis
The visual appearance of recurrence plots ( Figure 5) provides many hints a functional dynamics of land covers and their stability and predictability.

RP and RQA Analysis
The visual appearance of recurrence plots ( Figure 5) provides many hints about the functional dynamics of land covers and their stability and predictability.
Distances between diagonal lines, i.e., the period of the oscillation appear rather similar among the types of land cover. Higher predictability, determinism, recurrence rates and Lmax shown by Forests and Arable Lands mean that the dynamics of both land cover types are rather persistent and predictable, with strong periodic components ( Table 2). Relatively high laminarity means intermittency in the EVI of Forests whenever some states do not change or change slowly.
Forests and Arable Lands are more stable and resilient than Olive Groves and Grasslands, with effects on the provision of landscape services associated with the first two land cover classes. Forests and Arable Lands are rather quasi-periodic, but with determinism (DET) greater for Forests than for Arable Lands (Table 2). For Arable Lands, some states do not change or change slowly as an effect of human-managed balancing feedback loops at work (irrigation versus drought, soil fertilization versus soil impoverishment, etc.). Distances between diagonal lines, i. e., the period of the oscillation appear rather similar among the types of land cover. Higher predictability, determinism, recurrence rates and Lmax shown by Forests and Arable Lands mean that the dynamics of both land cover types are rather persistent and predictable, with strong periodic components ( Table 2). Relatively high laminarity means intermittency in the EVI of Forests whenever some states do not change or change slowly. Table 2. Results of RQA for EVI time series for the four land cover classes (Forests, Arable Lands, Olive Groves, Grasslands) in the Apulia region (95% confidence intervals in brackets).

Forests
Arable Lands Olive Groves Grasslands   The three-dimensional reconstruction of the EVI time delays in the phase space for Forests shows that the trajectories remained approximately in the same phase space in fourteen years, exhibiting the greatest stability and predictability with a minimum EVI in winter (Figure 6). Such a three-dimensional reconstruction supports the existence of much less temporal variation in CO 2 fluxes than commonly assumed for trees subjected to periodic photosynthetic down-regulation [62]. at work (irrigation versus drought, soil fertilization versus soil impoverishment, etc.).
The three-dimensional reconstruction of the EVI time delays in the phase space for Forests shows that the trajectories remained approximately in the same phase space in fourteen years, exhibiting the greatest stability and predictability with a minimum EVI in winter (Figure 6). Such a three-dimensional reconstruction supports the existence of much less temporal variation in CO 2 fluxes than commonly assumed for trees subjected to periodic photosynthetic down-regulation [62].

Maps of Recurrence Variables and Hot Spots Undergoing Degradation and Desertification (LDD)
The most significant recurrence variables related to predictability (Determinism, DET; Laminarity, LAM) and chaoticity (Divergence, DIV) for Apulia were mapped (Figure 7B-D) for the four most representative land cover types (Forests, Olive Groves, Arable Lands, Grasslands) together with recurrence estimates of other Corine land covers ( Figure  1) [63] derived from randomly selected MODIS pixels.

Maps of Recurrence Variables and Hot Spots Undergoing Degradation and Desertification (LDD)
The most significant recurrence variables related to predictability (Determinism, DET; Laminarity, LAM) and chaoticity (Divergence, DIV) for Apulia were mapped ( Figure 7B-D) for the four most representative land cover types (Forests, Olive Groves, Arable Lands, Grasslands) together with recurrence estimates of other Corine land covers (Figure 1) [63] derived from randomly selected MODIS pixels. The maps clearly show regions characterized by a greater recurrence (high predictability, DET, LAM, Figure 7C,D) mainly linked to the National Park of Gargano (1 in Figure  7A), Murgia of Trulli (7 in Figure 7A), the northern Daunia mountains (2 in Figure 7A) and extensive olive groves (9,10,11 in Figure 7A).
Conversely, evident and coherent space clusters, less predictable and with more chaoticity (higher divergence, DIV) ( Figure 7B), are clearly identifiable in central Alta Murgia (6 in Figure 7A), in the Gargano plateau (1 in Figure 7A) and in the southern Daunia mountains (2 in Figure 7A). In all these areas, the Grasslands class, including semi-arid Murgia of Trulli (7 in Figure 7A), the northern Daunia mountains (2 in Figure 7A) and extensive olive groves (9,10,11 in Figure 7A).
Conversely, evident and coherent space clusters, less predictable and with more chaoticity (higher divergence, DIV) ( Figure 7B), are clearly identifiable in central Alta Murgia (6 in Figure 7A), in the Gargano plateau (1 in Figure 7A) and in the southern Daunia mountains (2 in Figure 7A). In all these areas, the Grasslands class, including semiarid grasslands, garrigue, steppes and permanent pastures, has its highest concentration. These three areas correspond exactly to already known regional hot spots [19,29,31,64,65] with land desertification in progress.
More specifically, the central Alta Murgia area (6 in Figure 7A) is a limestone plateau primarily used as pastures with widespread surface features: karst and layers of outcropping rock with discontinuous and thin soils where xeric grasslands cover almost 24% of the total area. In this area, the Alta Murgia National Park is present (67,739 ha), which covers most of the territory and represents one of the most important conservation areas in Europe for the presence of habitats listed in the 92/43/EEC Habitat Directive with the highest priority for conservation, such as semi-natural dry grasslands and scrubland facies on calcareous substrates (Festuco-Brometalia), pseudo-steppe with grasses and annuals of the Thero-Brachypodietea* and eastern sub-Mediteranean dry grasslands (Scorzoneratalia villosae) [66].
In this area, the greatest decrease (−26%) in the extent of xeric grasslands from 1990 to 2018 is already known, and it has been mainly due to their conversion to Arable Lands [67] by extensive interventions of mechanical "stone crushing" (up to pulverization) of the soil-calcareous rock with the emergence of less productive shrublands [67]. This practice is banned nowadays.
As a result, in reclaimed arable lands, crops (mainly wheat) cover the thin soil for few months a year, and this increases the exposure to extreme weather events. The soil structural characteristics worsen as the percentage of organic matter decreases and the inert fine particles increase. Therefore, during heavy rains, even with moderate slopes, part of these low-quality soils can slide along the slopes [19]. Therefore, some pixels belonging to Arable Lands also contributed to the chaoticity in the area.
Another relevant cluster of high divergence with more chaotic greenness behavior is represented by the Gargano karst plateau (1 in Figure 7A), which occupies, from west to east, the central nucleus of the Gargano National Park. It is characterized by rangelands with bald and stony ridges alternating with woods. It is a karst environment, characterized by fields of sinkholes and tombs, with the alternation of semi-arid grasslands and pastures with trees, arable land and wooded areas. Here, the frequency of fires [21,68] has historically promoted the adaptation of evergreen trees such as Pinus halepensis and grasses such as Hyparrhenia hirta and Brachypodium ramosum, which can recover after fire disturbances [21]. The severity of the fire event influences the surface hydrological response and the consequent loss of soil. In the years following the fire, soil loss could be quite high and soil degradation processes much greater than before the event [19].
Finally, the third spatial cluster of high divergence is the South Daunia mountains (2 in Figure 7A), mainly associated with semi-arid grasslands, garrigue, steppes and permanent pastures, but also with land degradation for slope instability and subsequent landslides during seasonal rainfalls because of land use conversion and land abandonment [69,70]. Most of the unfavorable biophysical factors affecting Alta Murgia also occur substantially for southern Daunia.

Discussion
There are several approaches proposed for land degradation and desertification (LDD) analysis. One of the most popular spatially evaluates variables and indicators without considering their dynamics to map specific LDD issues through the recombination of thematic layers (e.g., [19,31,34,71]). However, all these studies are static and limited to the generation of desertification risk maps.
On the other hand, remote sensing-based analyses of vegetation cover decline rely on a wide range of change detection methods [72]. One approach is based on image classification when multitemporal Land Use Land Cover (LULC) maps are compared to identify changes (e.g., [73]). Such classification-based analyses, however, fail to evaluate the temporal evolution of LDD processes and can lead to high mapping inaccuracies due to error propagation between LULC classification results. Trend analyses of multiyear satellite images allow the capture of gradual LDD processes [74], but they are not independent from constraining assumptions like linearity. Trend analyses were routinely employed for LDD assessment using coarse-and multi-temporal imagery [75].
The present non-linear approach, based on the integration EVI-RQA, represents a clear novelty in the approaches for LDD assessment. It is a multi-temporal satellite analysis where non-linear time series analysis techniques are applied to extract information from time series of EVI pixels useful to analyze their complex behavior. It can also be useful in making the MEA definition of LDD operational, as poor and erodible soils, erratic precipitation (mainly during the winter), high summer temperatures with frequent drought events, land use conversion, extensive deforestation with frequent wildfires, and land abandonment can significantly reduce the ability of vegetation cover to provide, regulate and support landscape services. All these factors reverberate on the greater chaotic behavior of the EVI time series over time, mainly for the general Grasslands class (divergence, DIV; Table 2, Figure 7B). Therefore, greater chaoticity (less resilience) of the EVI over time can indicate LDD.
The results have shown that Forests and Arable Lands are more stable and resilient than Olive Groves and Grasslands, with effects on the provision of landscape services associated with the first two land cover classes. Forests and Arable Lands are rather quasiperiodic, but with determinism (DET) greater for Forests than for Arable Lands (Table 2). Specifically, for Arable Lands, some states do not change or change slowly, as highlighted by relatively high DET just with some intermittency, so that some states change slowly because of agricultural practices characterized by some periodic components. However, under human control, there are stress-responsive life-history strategies as for wheat (Triticum turgidum ssp. durum), whose cultivars exhibit traits adapted to both drought and heat [76].
The recurrent behavior of EVI trajectories as a proxy of ANPP indicates that services provided by Forests and Arable Lands are expected to continue with current balancing feedback loops at work. Some Mediterranean forests show, according to Garbulsky and colleagues [77], a low temporal variability in their EVI, being the coefficient of variation less than 10%. This is confirmed by the high DET, REC, Lmax and LAM values and the relatively very low DIV values for Forests (Table 2). Moreover, Mediterranean forests exhibit the highest determinism of the EVI (A N P P) ( Table 2), suggesting a strong self-organizing component resistant to the noise structure of the hydroclimatic forcing. These results are like those for European temperate forests, which are largely able to maintain relatively constant physiological settings and water retention even during severe drought and heat events [78]. Apulian evergreen forests are apparently in balance with the Mediterranean climate. This because their sclerophyllous oaks (Quercus ilex, Q. calliprinos and Q. coccifera) and Pinus halepensis represent the woody climax at the end of successional dynamics that can be achieved where biophysical factors remain favorable [79,80].
Permanent Olive Groves have seen their concentration and intensification increase like in other Mediterranean countries [81]. Olive Grove time series show shorter and more interrupted diagonals ( Figure 5); therefore, this land cover class is less periodic and less persistent than the century-old olive groves, which are more sparsely distributed probably because the management is based on traditional ecological knowledge and it is considered a best sustainable practice with benefits for the provision of landscape services. Hence, the recurrence rate of Olive Groves is lower than Forests and Arable Lands (Table 2). Their behavior appears more like Grasslands as satellite images also capture the grassy vegetation of the undergrowth. Olive Groves are moderately irrigated, and those laminarity states mainly occur because of balancing physiological feedback controls that determine under summer water stress less stomatal conductance and lower photosynthetic rates to avoid drought-induced hydraulic failure [82], and hence lower values of EVI. As a result of the adaptation of single olive trees to drought, traditional olive groves present an opposite periodicity with a maximum EVI in winter, and this can explain the lower amplitude of their periodic component.
In addition, the slight increasing trend of EVI time series of Olive Groves (Figure 3) may reflect the merging and enlarging of olive groves with the transition from ancient and traditionally managed farms, consisting of widely spaced century-old olive trees surrounded by a grassland matrix used for grazing animals, to modern and mechanical managed farms where olive trees are tightly packed with fewer grasses and more shrubs [81]. This replacement can be interpreted in terms of biodiversity loss as well as a decrease in the landscape services provided by century-old olive groves.
Natural Grasslands are also known for physiological balancing feedback loops inhibiting or reducing photosynthetic activity to adjust their responses to drought. Despite these adaptations, divergence is greatest for Grasslands, implying a more chaotic (less stable) signal (Table 2). This occurs because the structural dynamics and eco-physiological activity of xeric grasslands are mainly driven by discrete inputs or "pulses" of increasing seasonal precipitation [83] (Figure 4). Furthermore, there is an added spatial chaotic effect due to the unceasing conversion of natural Grasslands to Arable Lands through stone removal and crushing. From Table 2, it is evident that the divergence of Grasslands is six times greater than that of Forests. Finally, predictability decreases significantly from closed canopies to open drylands, i.e., from Forests to Grasslands, while chaoticity increases ( Table 2). If degradation and desertification worsen due to climate change and/or water shortages, xeric grasslands are likely the first system to degrade, and, therefore, this land cover type should be monitored as it represents the weakest state through which desertification can occur and spread.

Conclusions
The spatialization of recurrence variables from non-linear RP-RQA represents a clear novelty among the integrated approaches used for the monitoring and assessment of soil degradation and desertification (LDD). Its strength lies in the applicability to relatively short time-series, its complete independence from binding hypotheses and signal transformations, and in the provision of resilience-related recurrence variables that can be applied spatially to inform planning and management choices. This is fundamental to identify new sustainable landscape planning, intervention, and management choices, based on functional stability and predictability patterns. The approach could be used alone or alongside the traditional ones, based mainly on static structural knowledge.
The approach presented here measures persistent recurrence and chaoticity over time, and it appears to provide an efficient method to quantify and spatialize functional dynamic stability and predictability in dryland SELs, as well as to gauge the phase and function instabilities, which are otherwise undetectable. The mapping of the recurrence variables summarizes the (ir)regularity of the different land cover functional dynamics, providing a powerful insight to discriminate between areas with different stability/predictability in the landscape.
The effectiveness of this approach is validated by its ability to confirm already known hot spots with ongoing land desertification processes. However, its potentiality is to identify new areas of interest where to focus new more specific in situ investigations. In this perspective, the divergence (DIV) variable, represented by the chaoticity (lower resilience) of the EVI time series over time, could be used as a novel index, making the MEA definition of LDD operational.
The challenging goal of achieving the neutrality of soil degradation (LDN) by 2030 set by the United Nation Sustainable Development Goal (SDG) requires new conceptual strategies for analyzing and designing landscape sustainability models. In the case of drylands, they could include the strategic conservation of forests and xeric grasslands as well as the strategic planning of landscape conversion to reduce the intensity of water and soil stress.
Finally, it is well known that the future is not likely to be a simple extrapolation of the recent past because unexpected events can suddenly happen, as well as the recent exacerbated effects of climate change. However, having the knowledge of recent recurrence dynamics of vegetation cover from non-linear RP-RQA and of the spatial distribution of recurrence variables can assist in management and policies and foster the stability and short-term predictability (resilience) of dryland dynamics, guaranteeing the provision of landscape services.