This section of CASCADiS reports on our experimental findings regarding the role of a gradual increase in environmental pressure as trigger of sudden shifts in drylands towards degraded states.

Following a manipulative field experimental approach, we applied a gradually increasing pressure of vegetation removal, which aimed to mimic a gradient of accumulated impact of pressure (grazing and wood gathering). We then monitored resource loss (runoff and sediment yield), vegetation dynamics, bare-soil connectivity and soil-surface condition for ∼3.5 years. 

Before treatment  
After treatment

The main objectives of this research were to experimentally assess the occurrence of tipping points towards a degraded state in response to decreasing plant cover, and to gain insights about the mechanisms underlying such dynamics. Specifically we aimed to test the following hypotheses:

1. For a range of low cover values (close to suggested thresholds for irreversible degradation), a gradual decrease in cover results in a tipping-point response regarding resource loss (i.e., large, non-linear increase in runoff and sediment yield, as well as in rill development and erosion potential) once certain low cover value is surpassed.
2. Resulting from the above response, soil functioning sharply decreases with decreasing plant cover, contributing to the net loss of resources and eventually to an ecosystem shift to a degraded state.
3. Once the pressure is ceased, and depending on the plant cover that remains, the extant vegetation could either stabilize at its degraded state (if plant cover values were lower than a threshold value for degradation) or recover back to pre-treatment cover values and healthy state (if plant cover values were above that threshold).
4. Increased inputs of resources from upslope bare-soil inter-patches enhance the performance and growth of the individual downslope patches, reducing in turn overall bare-soil connectivity, and thus contributing to prevent or delay a shift towards a bare-soil state.

Study area

The study area (38°23'02.0"N; 0°34'55.4"W) is an experimental slope located in the semiarid area of the province of Alicante, Southeast Spain. It is a south-facing slope, slightly concave, and with an altitude of 170 m.a.s.l. (top of the plot). Slope angle is approximately 30%. The climate is semi-arid Mediterranean. Mean annual rainfall is 311 mm, which mainly falls in autumn and winter, and mean annual temperature is 18°C (Alicante weather station, 81 m.a.s.l., 38°22'21''N - 0°29'39'W; period 1981-2010). The geological substrate consists of Upper Cretaceous Albian blue-green marls and limestones, with an upper limit of Plio-Quaternary caliche and clay materials (Martín-Algarra and Vera 2004). Soil is loamy, with a high cover of rock fragments; the surface soil layer (0-20 cm) contains 21% clay, 40% silt, 39% sand, and 1.4% organic carbon.

Overview and detail of slope

The vegetation on the slope is a semiarid grassland-shrubland mosaic, with patchy distribution of perennial grasses such as Brachypodium retusum, Lygeum spartum and Hyparrhenia hirta, scattered shrubs of Thymelaea hirsuta, Rhamnus lycioides and Salsola genistoides, and a variety of chamaephytes, such as Fumana ericoides and Asparagus horridus. Total vegetation cover, before treatment application in the area, was 61% ± 7%. The only recent land use of the area was marginal grazing.

Experimental design

On a set of three large-size (∼ 300 m² each) experimental plots, we mimicked a gradual gradient of pressure (grazing and wood gathering) by removing manually part of the vegetation cover to intended final values of approximately 45, 30, and 15%. Each cover level was randomly assigned to each plot. Plots were then labelled as IC45, IC30 and IC15, in agreement with their respective initial cover values at the beginning of the experiment. We targeted these values because they are close to the threshold of ∼30% proposed by Thornes (1985) and Francis and Thornes (1990) as critical value below which runoff and sediment yield would substantially increase and drive the system to irreversible degradation.

Using aerial images of the experimental plots, we mapped the spatial distribution of the vegetation of each plot as binary maps (vegetation patches and bare-soil interpatches) with a resolution of 0.5 x 0.5 m pixel size. For each plot, and taking into account the original percentage of vegetation pixels, we calculated the number of vegetation pixels to be removed (changed to bare-soil pixels)  so that the resulting vegetation cover were the intended values of 45, 30, and 15%. Then, we randomly selected the vegetation pixels to be removed. The resulting maps (Fig. 4) were used as templates to guide the field implementation of the vegetation-removal treatment. With the help of parallel strings placed on the plots, we selected the 0.5 x 0.5 vegetated areas from which vegetation had to be removed in agreement with the guiding maps. Vegetation was removed by hand, uprooting the plants, or cutting them below the root collar when uprooting was not possible without digging. The soil surface was then compacted by stepping on it. The removal treatment was applied in November 2012.

Plant cover map

For monitoring purposes (see »Measurements and monitoring), we randomly selected 4 groups of 10 bare-soil pixels each: 10 pixels located immediately upslope a vegetation patch; 10 pixels located immediately downslope a vegetation patch, 10 pixels lateral to a vegetation patch, and 10 isolated interpatch pixels with no contact with vegetation.

Before the vegetation removal treatment, the experimental plots were inter-calibrated to assess natural variation between them in total runoff and sediment yield. For two years prior treatment application, we measured runoff and sediment yield (see »Measurements and monitoring) from the three plots. This long inter-calibration period allowed considering a wide range of variation in rainfall amount and intensity values for the comparative analyses of the hydrological response of the three experimental plots. The inter-calibration of the experimental units for the target processes studied is the standard procedure to control treatment effect when there are not replicates for the treatments applied. This approach, which is an extension of the “paired catchments” approach, has been widely used to assess the impact of land use or vegetation change on catchment hydrology (Bosch and Hewlett 1982). Such studies involve the use of two or more catchments with similar characteristics in terms of slope, aspect, soils, area, precipitation and vegetation located adjacent to each other. Following a calibration period, where all catchments are monitored, the catchments are subjected to different treatments. The change in the hydrological behaviour can then be attributed to changes in vegetation due to the treatment (Bosch and Hewlett 1982, Lane and Mackay 2001, Best et al. 2003).

Note: For full references to papers quoted in this article see

» References

Hydrological measurements

We measured runoff and sediment yield after each natural rainfall event for a period of 41 months (November 2012–March 2016). Runoff from each plot was collected in a gutter (Gerlach trough) connected to a large (1000 litres) water tank, where it was measured. For IC15 plot, we add a second water tank to increase the total capacity for storing runoff. The concentration of sediments in the runoff was estimated by desiccation of the runoff samples taken from the collection tanks after extensive stirring of the water-sediment mix and resuspension of sediments that had settled. Total sediment yield per rainfall event was calculated by adding up the amount of sediments in the tanks, estimated from the sediment concentration data, and the bed-load sediment accumulated in the gutters.

Rainfall was measured with a tipping-bucket rain gauge installed in the experimental site with a temporal resolution of 5 min. For each rainfall event we calculated total rainfall amount and maximum intensity in 15 minutes (I15). We then analysed the relationships between these rainfall properties and runoff and sediment yield.

Sediment and runoff monitoring

To assess the effect of the vegetation-removal treatment on the redistribution of water by runoff within the experimental plots, we installed TDR (Topp and Davis 1985) soil-moisture probes (0-15 cm depth) upslope of nine vegetation patches (Lygeum tussocks) per plot. These patches varied in the size of the bare-soil interpatch upslope each of them, with patches in plot IC15 having the largest upslope bare-soil interpatches, and patches in plot IC45 having the smallest ones. For selected well-forecasted large rainfall events, we measured soil water content (SWC) the day just before the event and the day after, estimating from the difference in soil moisture between these two days the water gains for each vegetation patch.

Differences in runoff and sediment yield between the experimental plots (treatment effect) were analysed by ANCOVA, using rainfall amount or rainfall intensity (I15) as covariable. The relationships between the size of the upslope interpatch area and the water input gained after a rainfall event were analysed by regression analysis. These analyses and all statistical analyses described hereafter were performed using R (R Core Team, 2013).

Microtopography and erosion potential

Starting in May 2014, we monitored the profile dynamics of main rills developed on the experimental plots after the vegetation-removal treatment. We identified in the field and marked the most conspicuous rills developed on each plot. We then installed a number of regularly spaced pairs of fixed rods transversal to the main axis of each rill. Periodically, we assessed the variation in the profile of the rills by using a simple two-dimensional microprofilemeter mounted on these pairs of rods.

Microprofile meter

On May 15, 2015, using a Remote Piloted Aircraft (RPA) equipped with a GPS system, we took 203 aerial images of the experimental area. From these images, we generated a high resolution orthophoto (5.7 mm per pixel) and a DEM (11.5 mm per pixel) of the area. Further treatment of the images we performed with QGIS and GRASS software. Using the high resolution DEM, we estimated a modified LS factor of RUSLE (Renard et al. 1997) for each pixel in a raster map with a resolution of 5 mm per pixel. The S-factor measures the effect of slope steepness, and the L-factor defines the impact of slope length. The combined LS-factor describes the effect of topography on soil erosion. Values of LS factor are dimensionless and represent relative values compared to soil loss from a field under the standard conditions of 22.1 m length and 9 percent steepness. Our calculation of the modified slope length (LS) factor was based on slope and (specific) catchment area, the latter as substitute for slope length (see, for example, Panagos et al. 2015). We used one of the hydrology modules (Module LS factor) available in SAGA (System for Automated Geoscientific Analyses) software, which provides numerous applications for data analyses focusing on DEMs. This module takes a Digital Elevation Model (DEM) as input and derives catchment areas according to Freeman (1991). Using this module we generated maps for each plot representing the modified LS for each pixel, so that larger values indicate higher erosion potential. Rills concentrate runoff from microcatchments and have local high slope angles; thereby they typically yield high LS values, which helps to easily identify them on a LS map.

Soil surface condition and functioning

We assessed soil surface condition and functioning following the LFA methodology (Tongway and Hindley, 2004). We estimated the LFA soil functional indices (stability, infiltration and nutrient cycling indices), by combining field measurements of eleven soil features (Table). Using 50 cm × 50 cm sampling quadrats, and following a rating scale (Tongway and Hindley, 2004), these soil features were measured for the four types of bare-soil microsites: upslope, lateral and downslope a patch of Lygeum spartium, and isolated bare-soil interpatch. Ten sampling cuadrats were assessed per type of microsite and plot in March and October 2013, September 2014 and July 2015.

Table: LFA functional indices and soil surface features used to estimate each of them (Tongway and Hindley 2004).

Vegetation performance, pattern, and dynamics

Plant cover was monitored by the point intercept method along three (25-m length) permanent transects per plot (125 points per transect). Plant cover and species number were recorded in November 2012, June 2013, August 2014, June 2015, and April 2016.

From the orthophoto obtained with the RPA on May 2015, we created a binary map and estimated total plant cover on each plot. Combining the orthophoto and the DEM, we calculated the index Flowlength (FL) for each plot. Flowlength quantifies the connectivity of bare-soil areas considering both vegetation pattern and topography. It is calculated as the average of the runoff pathway lengths from all the cells in a raster-based map of the area (Mayor et al. 2008).

Finally, to improve our understanding of the role of source–sink dynamics in plant productivity in drylands, and gain insights on the local (patch-scale) effect of pattern-driven resource redistribution, we monitored 10 individual Lygeum tussocks per plot and estimated the tussock basal growth (change in basal diameter with time), and reproductive effort (number of spikes per tussock basal area unit).  Since tussock growth may depend on the tussock size (in turn depending on age), we represented tussock growth relative to the initial size. We studied these relationships using regression analyses. We compared growth and reproductive effort between plot treatments using repeated measures ANOVA.

Note: For full references to papers quoted in this article see

» References

Inter-calibration of experimental plots

During the 2-years inter-calibration period, the hydrological response of the three experimental plots was very similar. There were no statistical differences between plots regarding the relationships between rainfall and runoff (Figure 1). Although we measured some between-plots differences in runoff production in response to small rainfall events, they were minor differences, probably related to differences in the runoff produced very close to the plot outlet. Conversely, the runoff produced by large rainfall events was almost identical in the three plots. We obtained similar results regarding sediment yield, as well as using other rainfall metrics, such as rainfall intensity (data not shown).

Figure 1

Given the lack of differences in the hydrological response of the three experimental plots, the differences found after treatment application were attributed to the treatment effect.

Hydrological response

During the study period reported here (Nov 2012-March 2016), we recorded 40 rainfall events that produced runoff on one or more plots (Figure 2). These events totalized 501 mm, 60% of total rainfall in the area in that period. The temporal variability of runoff-producing rainfall events was very high, with many medium-size events concentrated in autumn of 2012 and 2014, few very large events in spring 2013 and autumn 2015, and long dry intervals between these periods. It is worth noting the extraordinary dry period occurred between July 2013 and October 2014, only interrupted by six events that produced runoff. The year 2015 was also extremely dry, with only five events that produced runoff, yet one of them was particularly large.

Figure 2

Total runoff produced over the study period was 3.7 mm, 11.1 mm. and 13.2 mm, for plots IC45, IC30, and IC15, respectively, which means that 0,4% to 1,6% of total rainfall was lost as runoff from the experimental plots. At the event scale, runoff coefficients for runoff-producing events ranged between 0.01% and 4.22% (average 0.33%) for IC45; 0.33%); 0.01% and 13.94% (average 0.99%) for IC30; and 0.01% and 9.20% (average 1.30%) for IC15. Most of this runoff was produced during the first six months of the experiment, before the severe 2013-2014 drought, while the longer but drier period of July 2013-March 2016 yielded approximately half of the runoff produced over the first six months (Figure 3). The production of sediment yield followed the same pattern, In general, plots IC30 and IC15 produced much higher runoff and, particularly, sediments than IC45, with IC15 showing only slightly higher yields than IC30 (Figure 3).

 
Figure 3

Rainfall-runoff and rainfall-sediments relationships and ANCOVA analysis (Figures 4 and 5; Table 1) showed significant variation in the hydrological and erosional response of the three plots (i.e., significantly different regression slopes, as expressed by significant interaction effect treatment x rainfall; Table 1), with almost identical response to rainfall of plots IC30 and IC15, and much higher production of runoff and sediments in these two plots than in IC45 for any given rainfall intensity higher than ∼20 mm h-1. The best explanatory rainfall variable was the AI15 index, which captures the influence of both rainfall amount and intensity (Figure 4). Results for only rainfall amount or intensity were similar, yet rainfall amount explained less of the variation in runoff and sediment yield (Figure 5; note R2) and showed weaker treatment and rainfall-treatment interaction effects than AI15 (Table 1).

Figure 4
Figure 5

Table 1. ANCOVA results (F and p values) for treatment and rainfall (covariable) main effects and interaction effect on runoff and sediment yield. Treatments: IC45, IC30, IC15. Covariables: maximum rainfall intensity in 15 minutes (I15; mm h^-1), rainfall amount (A, mm), and AI15 index (mm² h^-1)

The residuals from the relationships shown above between rainfall and runoff or sediment yield (using the best explanatory variable: Maximum intensity in 15 minutes, I15) can be plotted against time since treatment implementation. This type of graphical representation (Figure 6) provides information on the trends in the hydrological behaviour once the effect of the magnitude and/or intensity of the rainfall have been removed from the analysis. For the three experimental plots, the residuals from the regression between runoff and AI15 index showed a constant trend, which indicates that the capacity of the plots for producing runoff (i.e. losing resources) did not significantly vary with time since treatment application. This means that, for any given type (amount-intensity) of rainfall, the hydrological response of the plots was more or less the same either soon or long time after treatment application.

Figure 6

The relationships between sediment yield and runoff were very similar for all plots (Figure 7, left), with IC30 showing a higher slope and therefore indicating a higher average sediment concentration in runoff than the other two plots. The residuals from these relationships plotted against the time elapsed since treatment application showed that sediment concentration in runoff did not vary with time, except for a hardly perceptible decreasing trend in IC15 (Figure 7, right).

Figure 7

Regarding water redistribution at the patch scale, soil water gain upslope Lygeum tussocks after a runoff-producing rainfall event, varied between plots (F= 3,76; p=0,037; ANOVA), with higher water gains upslope tussocks in CI15 and lower in CI45 (Figure 8, left). The variation in soil water gains depended on the size (length) of the upslope interpatch (Figure 8, right).

Figure 8

Microtopography

The assessment of rill-profile dynamics showed that the largest rills developed on plot IC15 (Figure 9), yet the average maximum depth for the assessment period was similar for rills on plots IC30 (5.1 ± 0.6 cm) and IC15 (5.0 ± 0.4 cm) and smaller for rills on IC45 (4.2 ± 0.4 cm). However, while rills on IC45 hardly varied their profile with time, rills on IC30 and IC15 showed a further increase in rill incision from May 14 to Sep 14 (Figure 9) followed by a partial refilling of the rill (from Sep 14 to Dec 14).

Figure 9
Figure 10

From the high resolution orthophoto taken in May 2015, we estimated the spatial variation and average value per plot of the LS factor, which captures the contribution of topography to erosion potential. Average LS-factor values were 2.23, 2.77, and 2.75, for IC45, IC30 and IC15, respectively. The spatial variation of LS factor highlighted the rill network of each plot, which showed longer and denser rills in IC30 and IC15 plots than in IC45 (Figure 10). While IC45 concentrated high LS values and a rill area on the top of the plot that did not connect with the plot outlet, IC30 and, particularly, IC15 showed highly connected network of rills and high LS areas.

Soil surface condition and functioning

In general, bare-soil interpatch areas (microsites) located next to a vegetation patch (e.g. lateral interpatch) showed slightly higher values for the LFA soil functional indices than interpatch areas located at certain distance from vegetation patches (isolated interpatch), particularly for the Infiltration and Nutrient cycling indices (Figure 11).

Figure 11

Isolated bare-soil interpatches on plots IC30 and IC15 showed no variation with time for the Infiltration and Nutrient cycling indices, whereas IC45 showed higher values for these two indices at the beginning of the monitoring period (March 2013) that rapidly decreased and then stabilized at slightly lower vales than IC30 and IC15 (Figure 11). Conversely, the stability index for isolated interpatches showed an increasing trend during the first 1.5 year of the monitoring period (March 13 - Oct 14), particularly for  IC45 plot, which changed from initial (March 2013) lower values to final (June 1015) higher values than IC30 and IC15.

Bare-soil interpatch areas located next to vegetation patches (lateral interpatches) showed a very similar pattern for Infiltration and Nutrient cycling indices: lower values and no variation with time for IC15; a clear decreasing trend for IC45; and a slight decreasing trend for IC30, which all together resulted in no differences between plots at the end of the monitoring period (Figure 11). The Stability index for lateral interpatches showed no variation or slight decreasing trend for IC15 and IC30, respectively. However, IC45 showed an increasing trend with time, which resulted in higher values of this index for IC45 than for the other two plots at the end of the monitoring period.

Likewise lateral interpatch areas, interpatch areas located upslope or downslope a vegetation patch also showed higher values than isolated interpatch areas. However, these two types of interpatch area showed less variation with time and between experimental plots (data not shown).

Vegetation performance and dynamics

The dynamics of plant cover and number of species reflected the impact of both the plant-removal treatment and the severe 2013-14 drought (Figure 12). During the first rainy period and growing season after the treatment application, plant cover increased in all plots. Then, in response to the drought, plant cover decreased in all plots, particularly in IC30. After the drought, plant cover stabilized around 20-25% cover in IC15 and IC30 plots, and slowed down the decreasing trend in IC45. Only during the last year of the monitoring period reported here, plant cover has showed certain degree of recovery, resulting in value around 30% for IC15 and IC30 plots, and close to 40% for IC45 plot (Figure 12, left), still very far from the original pre-treatment values. The number of species followed a similar dynamics (Figure 12, right).

Figure 12
Figure 13

In May 2015, total plant cover values estimated from a high-resolution aerial picture taken from a Remote Piloted Aircraft (Figure 13) fully matched the values estimated from the field transects for that period. The between plot variation in bare-soil connectivity (estimated by the Flowlength index; Mayor et al. 2008), did not match the variation in plant cover, as IC15 showed higher bare-soil connectivity than IC30 despite having a similar, slightly higher plant cover value. Average patch size for IC45, IC30 and IC15 were 292 ± 33 cm², 239± 22 cm², and 423± 85 cm², respectively.

The growth of Lygeum spartum tussocks (relative to their initial size) showed different trends between the different plots and the dates on which they were measured (Figure 12). During the first part of the assessment period (May 2013 - May 2014) individuals of IC15 and IC30 plots showed an increasing trend in growth with increasing initial size, while individuals of IC45 plot showed the opposite trend (Figure 14,left). During the second part of the assessment period (May 2014 - May 2015), the observed trends were almost opposite to the ones observed the previous period, so that individuals of plots IC30 and IC15 showed a decreasing trend in growth with increasing initial size, especially clear in individuals of IC30. Individuals of IC45 did not show any effect of the initial size on their growth during this period (Figure 14, right). However, these differences in trend between plots were not statistically significant (ANCOVA Analyses).

Figure 14
Figure 15

Production of flowers by Lygeum individuals showed a clear impact of the drought period, which reduced production to almost zero, but not did not show any significant effect of the treatment, except for a slight, non-significant decreasing trend from IC45 to IC15 during the first assessment period (Figure 15).

Note: For full references to papers quoted in this article see

» References

Discussion

Our results demonstrate that the gradual decrease of (already) low vegetation cover in a dryland system causes a non-linear increase in the production of runoff and sediments, with the change from 45 to 30% being critical for that increase, whereas no further significant increase was produced by a change from 30% to 15% cover. Similarly, rill development and the contribution of microtopography to erosion potential (LS-factor) were much lower in IC45 than in IC30 and IC15 plots, which did not present major differences between them.

These results support our two first hypotheses and are consistent with a number of conceptual frameworks and models that point to a sharp increase in bare-soil connectivity at low vegetation cover values as the main factor driving to the substantial change in runoff and erosion (Ludwig and Tongway 1995, Davenport et al. 1998, Wilcox et al. 2003, Mayor et al. 2013, Saco and Moreno-de las Heras 2013). Furthermore, the fact that plots IC30 and IC15 have a similar hydrological behaviour, with much higher resource loss than IC45 suggests that both IC30 and IC15 resulted in a particular degraded state, significantly different from IC45.

For our experimental system, a small reduction in vegetation cover within the range of 45%-30% of vegetation cover caused a critical change in the hydrological behaviour of the system. This sensitive range or threshold area could vary among dryland systems depending on a number of environmental factors such as soil type, topography, or rainfall regime, that have proven to be critical control factors of dryland hydrology (e.g., Francis and Thornes 1990, Boorman et al. 1995, Devito et al. 2005, Moody and Martin 2009).

A critical question is if the sharp increase in resource loss resulting from a decrease in vegetation cover triggers a feedback loop that further decreases plant cover, eventually leading to full degradation (ecosystem shift to a desert/unproductive state). Within CASCADE, we have incorporated this type of feedback to a model of dryland vegetation dynamics (Mayor et al. 2013), and we found that for a gradual increase of environmental pressure (e.g., aridity, grazing), these feedbacks decrease the amount of pressure required to cause a critical shift to a degraded state. However, the analysis of the runoff dynamics over the study period reported here showed no significant increasing trend for any of the experimental plots, which suggest that no degradation loop has been triggered so far. Instead, our results suggest that the disturbance and stress experienced by the experimental plots changed their state to a higher level of degradation from which they do not seem to be recovering, neither further degrading. Similarly, the rills initially developed on the plots did not show a consistent increasing trend in their incision level, but periods of incision and deposition that kept their profiles without major net changes.

A previous long-term study that assessed erosion losses after disturbances of varying intensities found irreversible loss of soil productivity below vegetation cover values around 20%, with successful recovery over time for cover values above this value (Gao et al., 2011). In agreement with this study and with the conceptual framework of sudden shifts in ecosystems (e.g., Scheffer et al 2001, Rietkerk et al. 2004), we hypothesized that once the pressure is ceased, the extant vegetation could either stabilize at its degraded state or recover pre-disturbance cover values and its healthy state, depending on the plant cover value that resulted from the disturbance. However, as discussed above, we found no evidence that pointed to the recovery of pre-disturbance vegetation cover for any of the plots. Only at the end of the monitoring period, 3.5 years after treatment application, we found a slight increasing trend in vegetation cover for all the three plots. However, the final values measured (∼ 30% cover for IC30 and IC15 plots and ∼ 40% for IC45) are still very far from pre-disturbance cover values in the area (61% ± 7%). The variation in the number of species matched the variation in cover, which suggests that decreasing cover resulted in the local extinction or reduction of the least suitable species for the new conditions created.

Sediment concentration in runoff showed a slight and nonsignificant decrease with time, which could be related to the final slight recovery in vegetation cover, but it could also reflect some exhaustion in the amount of easily detachable soil material with time, a rather common process in closed plots (Boix-Fayos et al. 2007).

It could be argued that the vegetation-removal gradient applied was a relatively mild pressure gradient, as it was applied only once and the soil was not excessively altered, as per triggering the hypothesized degradation loop. It must be stressed though that the combined effect of the experimental treatment and the severe drought that occurred in the area in 2013-2014 reduced vegetation cover to very low values, particularly in IC15 and IC30 plots (∼20%), low enough to represent the impact of high environmental pressure over the system. However, despite the clear and nonlinear increase in resource losses (runoff, soil) that resulted from the decrease in cover, the total amount of water lost through runoff only represented between 0.4 and 1.6% of total water input from rainfall, which probably is not enough relevant for triggering a degradation loop in the timeframe of our experiment.

It is important to note that while the vegetation-removal treatment implied a clear stress gradient between plots, the natural drought implied a homogeneous stress for all three plots, yet the impact could have varied as a function of plant biomass/cover (see below). Therefore, the effect of the drought might have masked to some extent the differences in stress due to the vegetation-removal levels.

Initially, the new degradation level of the experimental plots reflected the degree of vegetation removal. Thus, initial values (i.e., monitoring values during the first year after treatment application) of vegetation cover and bare-soil LFA infiltration and nutrient cycling indices linearly decreased from IC45 to IC30 and IC15. However, the natural drought occurred in the area during the second year of the monitoring period appeared to have larger adverse impacts in plots with higher vegetation cover, which seems to have driven both vegetation cover and bare-soil functional status to certain convergent trend between plots. The LFA stability index did not show as clear dynamics and between-plots differences as the other two indices, which can be explained by a lower sensitivity of this index in Mediterranean conditions (Mayor and Bautista 2012).

It is worth mentioning that the differences found in the functional conditions of the bare-soil interpatches as a function of the distance to vegetation patches (i.e., isolated bare-soil areas versus bare-soil areas next to vegetation patches) support the assumptions of dryland vegetation models that incorporate local facilitation as one of the critical processes that drive dryland vegetation dynamics and response to environmental change (e.g., Kéfi et al. 2007a), yet the functional differences found in our work do not seem to be as contrasting as assumed by these models. Our results also highlight the interacting role of drought and vegetation cover as control factors of soil functioning.

At the patch scale, larger bare-soil connectivity also implies larger inter-patch areas, which is beneficial for the performance of the downslope patch. It is known that the transfer of resources from bare-soil interpatches to downslope vegetation patches contribute to plant productivity (Aguiar and Sala 1999, Yu et al. 2008, Turnbull et al. 2012, Urgeghe and Bautista 2015, see also »Potential for shifts). Enhanced patch growth would then contribute to reduce the size of the bare-soil areas, which in turn would reduce the amount of resources transferred to downslope patches. In agreement with this rationale, it has been proposed that the relative amount of source and sink areas in functional landscapes remains within a certain range of variation around a hypothesized optimum source:sink ratio that maximizes both the availability of resources and the growth of vegetation patches (Puigdefábregas et al. 1999, Urgeghe at al. 2010).

We hypothesized that increased inputs of resources from upslope bare-soil inter-patches enhance the performance and growth of the individual downslope patches, which in turn reduce overall bare-soil connectivity, contributing to prevent or delay a shift towards a bare-soil state. Our study proved that runon inputs to individual patches increased with increased size of the upslope interpatch area, which resulted in larger runon input for patches in plots IC15, followed IC30 and IC45. The relationship between soil water gains and upslope interpatch area tended to a plateau, so that larger upslope areas beyond certain value did not further increase water gains. This could be due to a negative relationship between drainage-catchment size and runoff yield per unit area (e.g., Mayor et al., 2011), which is explained by the area-dependent opportunity for runoff to infiltrate before it reaches a downslope vegetation patch (Wilcox et al., 2003). It could also be that large water inputs from large upslope interpatches exceed the maximum capacity of plant patches to capture runon.

Regarding plant performance, we found that before the drought the larger the individual patch in IC30 and IC15 plots the greater the patch growth, suggesting a larger capacity for capturing and benefiting from runon. However, during and after the drought, these relationships shifted to decreasing trends (lower growth with increasing initial size), particularly for patches in IC30 plot, which could be explained by the expected lower growth rate with increasing plant age, as size co-varies with age, but could also be explained by a higher impact of drought and/or enhanced global resource loss on the largest patches or on patches that grew more before the drought.

These results highlight the importance of the redistribution of resources and feedbacks that operate at the patch scale, which contribute to maintain patch productivity under very low plant cover, and could counterbalance net losses from the system, preventing or delaying the shift to a bare-soil state and promoting instead different communities and patterns with different overall resource availability and redistribution (Puigdefábregas et al. 1999, Rietkerk et al. 2004). However, our results also highlight the complexity of the interactions between endogenous feedbacks and external pressures as control factors of dryland dynamics (D’Odorico et al., 2012). Thus, for example, the influence of patch size in plant performance depended on the level of drought stress.

Conclusions and final remarks

Overall, the results reported here support the conjecture that a reduction in plant cover beyond certain threshold values nonlinearly increases resource (water, soil) losses from the system and may trigger a change to a degraded state that functions under less resource-wise conservative conditions. For our case study, there was an area of particular sensitivity within the transition from 40% to 30 % vegetation cover, in which small changes in the cover percentage resulted in a turning point in the hydrological response of the system.  The pressure exerted over the system resulted from a gradient of vegetation removal combined with the impact of a severe drought naturally occurring in the area. All over the study period, and once the pressure ceased, there were no evidences pointing to any recovery trend towards pre-disturbance levels. However, our results did also not provide evidence for the onset of any degradation loop that could further degrade the systems. The three experimental plots shifted to degraded states with vegetation cover around 25-35% at the end of the study period (half of the pre-disturbance vegetation cover value). The assessment of the stability of these new degraded states falls beyond the timeframe of this work.

The major strengths of the research reported here include

1. the use of a manipulative approach to establish a pressure gradient, which allowed a large degree of control in the level of stress exerted over the systems, and
2. the comprehensive set of processes and responses addressed at the global (plot) and local (patch) scales, which provided a rather complete picture of the system behaviour and dynamics.

However, a clear limiting factor for this work is the relatively short duration of the experiment, as the systems have proven to be very dynamic and our research would benefit from a longer monitoring period that provides insights about longer-term dynamics. Acknowledging this, we extended the monitoring period beyond the original window scheduled, and we will continue the monitoring till the end of CASCADE project duration, incorporating any potential new finding to the final dissemination products of CASCADE. Nevertheless, the study period (∼3.5 years) reported here has provided a solid temporal framework for addressing most of the questions formulated in this part of the CASCADE research, and has yielded critical findings and insights that are relevant for improving modelling and theoretical developments on dryland dynamics. Thus, for example, our results on differential bare-soil functioning dynamics depending on the distance to a plant patch; the positive power-type relationship between upslope interpatch area and resource inputs from runon; or the changes in resource loss with decreasing cover are essential inputs for parametrizing and calibrating dryland dynamics models. From the management point of view, our results highlight the need for being conservative regarding the minimum vegetation cover that should be maintained through proper resource exploitation.
Finally, by incorporating the use of RPA (drone) systems to the monitoring of our experimental plots, we were able to very precisely assess vegetation patterns and surface conditions. This tool will largely enhance our capacity for assessing further dynamics in the system and eventually test the degree of stability of the new conditions generated.

Note: For full references to papers quoted in this article see

» References

References cited in this section of CASCADiS

• Aguiar MR, Sala OE. 1999. Patch structure, dynamics and implications for the functioning of arid ecosystems, Trends Ecology Evolution 14: 273–277.
• Bautista S, Mayor AG, Bourakhouadar J, Bellot J. 2007. Plant spatial pattern predicts hillslope runoff and erosion in a semiarid Mediterranean landscape. Ecosystems 10: 987-998.
• Best A, Zhang L, McMahon T, Western A, Vertessy R. 2003. A critical review of paired catchment studies with reference to seasonal flows and climatic variability. CSIRO Land and Water Technical Report 25/03. Murray-Darling Basin Commission, Canberra, AU.
• Bestelmeyer BT, Duniway MC, James DK, Burkett LM, Havstad KM. 2013. A test of critical thresholds and their indicators in a desertification-prone ecosystem: more resilience than we thought. Ecology Letters 16: 339-345.
• Bestelmeyer BT, Ellison AM, Fraser WR, Gorman KB, Holbrook SJ, Laney CM, Ohman MD, Peters DPC, Pillsbury FC, Rassweiler A, Schmitt RJ, Sharma S. 2011. Analysis of abrupt transitions in ecological systems. Ecosphere 2: art129. http://dx.doi.org/10.1890/ES11-00216.1
• Boettiger, C, Hastings  A. 2013. Tipping points: From patterns to predictions. Nature 493: 157-158. DOI.10.1038/493157a.
• Boix-Fayos C, Martínez-Mena M, Calvo-Cases A, Arnau-Rosalén E, Albaladejo J, Castillo V. 2007. Causes and underlying processes of measurement variability in field erosion plots in Mediterranean conditions. Earth Surface Processes and Landforms 32: 85-101.
• Boorman DB, Hollist JM, Lilly A. 1995. Hydrology of soil types: a hydrologically-based classification of the soils of the United Kingdom. Report no. 126. Institute of Hydrology, Wallingford, UK.
• Bosch JM, Hewlett JD. 1982. A review of catchment experiments to determine the effect of vegetation changes on water yield and evapotranspiration. Journal of Hydrology 55:  3-23.
• Clements CF, Ozgul A. 2016. Including trait-based early warning signals helps predict population collapse. Nature Communications 7, 10984. doi:10.1038/ncomms10984
• Dakos V, Carpenter SR, Brock WA, Ellison AM, Guttal V, Ives AR, Kéfi S, Livina V, Seekell DA, van Nes EH, Scheffer M. 2012. Methods for detecting early warnings of critical transitions in time series illustrated using simulated ecological data. PLoS ONE 7(7): e41010.
• Davenport DW, Breshears DD, Wilcox BP, Allen CD. 1998. Viewpoint: sustainability of piñon-juniper ecosystems—a unifying perspective of soil erosion thresholds. Journal of Rangeland Management 51:231–240.
• K. Devito K, Creed I, Gan T, Mendoza C, Petrone R, Silins U, Smerdon B. 2005. A framework for broad-scale classification of hydrologic response units on the Boreal Plain: is topography the last thing to consider?. Hydrological Processes, 19: 1705–1714. DOI: 10.1002/hyp.5881
• D'Odorico P, Okin GS, Bestelmeyer B. 2012. A synthetic review of feedbacks and drivers of shrub encroachment in arid grasslands. Ecohydrology 5: 520-530.
• Duarte CM, Agustí S, Wassmann P, Arrieta JM, Alcaraz M, Coello A, Marbà N, Hendriks IE, Holding J, García-Zarandona I, Kritzberg E, Vaqué D. 2012. Tipping elements in the Arctic marine ecosystem. Ambio 41, 44–55.
• Gao Y, Zhong B, Yue H, Wu B, Cao S. 2011. A degradation threshold for irreversible loss of soil productivity: a long-term case study in China. Journal of Applied Ecology 48: 1145-1154.
• García-Fayos P, Bochet E, Cerda A. 2010. Seed removal susceptibility through soil erosion shapes vegetation composition. Plant Soil 334:289–297
• Kéfi S, Guttal V, Brock WA, Carpenter SR, Ellison AM, Livina VN, Seekell DA, Scheffer M, van Nes EH, Dakos V. 2014. Early warning signals of ecological transitions: methods for spatial patterns. PLoS ONE 9(3): e92097.
• Kéfi S, Rietkerk M, van Baalen M, Loreau M. 2007a. Local facilitation, bistability and transitions in arid ecosystems. Theoretical Population Biology 71(3):367–379.
• Kéfi S, Rietkerk M, Alados CL, Pueyo Y, Papanastasis VP, ElAich A, de Ruiter PC. 2007b. Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems. Nature 449:213–217.
• Lane PNJ, Mackay SM. 2001. Streamflow response of mixed-species eucalypt forest to patch cutting and thinning treatments. Forest Ecology and Management, 143: 131-142.
• Ludwig JA, Tongway DJ. 1995. Spatial organisation of landscapes and its function in semi-arid woodlands, Australia. Landscape Ecology 10:51–63.
• Ludwig JA, Wilcox BP, Breshears DD, Tongway DJ, Imeson AC. 2005. Vegetation patches and runoff-erosion as interacting ecohydrological processes in semiarid landscapes. Ecology 86(2): 288–297.
• Martín-Algarra A, Vera JA. 2004. La Cordillera Bética y las Baleares en el contexto del Mediterráneo Occidental. In JA Vera (ed.), Geología de España. SGE-IGME/Madrid.
• Mayor AG, Bautista S, Bellot J. 2011. Scale-dependent variation in runoff and sediment yield in a semiarid Mediterranean catchment. Journal of Hydrology 397: 128–135.
• Mayor AG, Bautista S, Small EE, Dixon M, Bellot J. 2008. Measurement of the connectivity of runoff source areas as determined by vegetation pattern and topography: a tool for assessing potential water and soil losses in drylands. Water Resources Research 44: W10423.
• Mayor AG, Kefi S, Bautista S, Rodríguez F, Cartení F, Rietkerk M. 2013. Feedbacks between vegetation pattern and resource loss dramatically decrease ecosystem resilience and restoration potential in a simple dryland model. Landscape Ecology 28(5): 931–942. DOI: 10.1007/s10980-013-9870-4.
• Millennium Ecosystem Assessment. 2005. Ecosystems and Human Well-being: Desertification Synthesis. World Resources Institute, Washington, DC.
• Moody JA, Martin DA. 2009. Synthesis of sediment yields after wildland fire in different rainfall regimes in the western United States. International Journal of Wildland Fire 18, 96–115.
• Moreno-de las Heras M, Saco PM, Willgoose GR, Tongway DJ. 2011. Assessing landscape structure and pattern fragmentation in semiarid ecosystems using patch-size distributions. Ecological Applications 21(7): 2793–2805.
• Peters DPC, Havstad KM. 2006. Nonlinear dynamics in arid and semi-arid systems: Interactions among drivers and processes across scales. Journal of Arid Environments 65 (2006) 196–206.
• Puigdefábregas J. 2005. The role of vegetation patterns in structuring runoff and sediment fluxes in drylands. Earth Surface Processes and Landforms 30: 133–147. DOI: 10.1002/esp.1181.
• Puigdefábregas J, Solé A, Gutierrez L, del Barrio G, Boer M. 1999. Scales and processes of water and sediment redistribution in drylands: results from the Rambla Honda field site in southeast Spain. Earth Science Review 48: 39–70. DOI: 10.1016/S0012-8252(99)00046-X.
• R Core Team, 2013. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/.
• Renard KG, Foster GR, Weesies GA, McCool DK, Yoder DC. (coords.) 1997. Predicting soil erosion by water: A guide to conservation planning with the Revised Universal Soil Loss Equation (RUSLE). USDA-ARS Agriculture Handbook 703.
• Reynolds J, Stafford-Smith DM, Lambin EF, Turner II BL, Mortimore M, Batterbury SPJ, Downing TE, Dowlatabadi H, Fernández RJ, Herrick JE, Huber-Sannwald E, Jiang H, Leemans R, Lynam T,
• Maestre FT, Ayarza M, Walker B. 2007. Global desertification: Building a science for dryland development. Science 316, 847
• Rietkerk M, Dekker SC, de Ruiter PC, van de Koppel J. 2004. Self-organized patchiness and catastrophic shifts in ecosystems. Science 305: 1926–1929.
• Saco PM, Moreno-de las Heras M. 2013. Ecogeomorphic coevolution of semiarid hillslopes: emergence of banded and striped vegetation patterns through interaction of biotic and abiotic processes. Water Resources Research 49(1): 115–126. DOI: 10.1029/2012WR012001.
• Scheffer M, Carpenter S, Foley JA, Folke C, Walker B. 2001. Catastrophic shifts in ecosystems. Nature 413: 591-596.
• Scheffer M, Carpenter, SR. 2003. Catatrophic regime shifts in ecosystems: linking theory to observation. Trends Ecol. Evol., 18, 648–656.
• Schlesinger WH, Pilmanis AM. 1998. Plant-soil interactions in desert. Biogeochemistry, 42: 169–187.
• Solé R. 2007. Ecology: Scaling laws in the drier. Nature 449(7159):151–153.
• Thornes JB. 1985. The ecology of erosion. Geography 70: 222-236.
• Francis CF, Thornes JB. 1990. Runoff hydrographs from three Mediterranean vegetation cover types. In JB Thornes (Ed.), Vegetation and erosion. John Wiley & Sons, Ltd, pp. 363-384.
• Tongway DJ, Hindley N. 2004. Landscape Function Analysis: Procedures for Monitoring and Assessing Landscapes. CSIRO Publishing, Brisbane, Australia.
• Topp GC, Davis JL. 1985. Measurement of soil water content using timedomain reflectrometry (TDR): a field evaluation. Soil Science Society of America Journal 49(1): 19–24.
• Turnbull L, Wilcox BP, Belnap J, Ravi S, D’Odorico P, Childers, Gwenzi W, Okin G, Wainwright J, Caylor KK, Sankey T. 2012. Understanding the role of ecohydrological feedbacks in ecosystem state change in drylands. Ecohydrology 5: 174–183.
• Turnbull L, Wainwright J, Brazier RE. 2010. Changes in hydrology and erosion over a transition from grassland to shrubland. Hydrological Processes 24, 393–414.
• Urgeghe AM. Bautista S. 2015. Size and connectivity of upslope runoff-source areas modulate the performance of woody plants in Mediterranean drylands. Ecohydrology 8, 1292–1303. DOI: 10.1002/eco.1582.
• Urgeghe AM, Breshears DD, Martens SN, Beeson PC. 2010. Redistribution of runoff among vegetation patch types: on ecohydrological optimality of herbaceous capture of runon. Rangeland Ecology and Management 63(5): 497–504.
• van Auken, O.W. 2000. Shrub invasions of North American semiarid grasslands. Annual Review of Ecology and Systematics, 31: 197–215.
• van de Koppel JPM, Herman J, Thoolen P, Heip CHR. 2001. Do alternate stable states occur in natural ecosystems? Evidence from a tidal flat. Ecology 82: 3449–3461.
• von Hardenberg J, Meron E, Shachak M, Zarmi Y. 2001. Diversity of vegetation patterns and desertification. Physical Review Letters 87(19):198101.
• Wilcox BP, Breshears DD, Allen CD. 2003. Ecohydrology of a resource conserving semiarid woodland: effects of scale and disturbance. Ecological Monographs 73: 223–239.
• Yu M, Gao Q, Epstein HE, Zhang X. 2008. An ecohydrological analysis for optimal use of redistributed water among vegetation patches. Ecological Applications 18: 1679–1688.