Vegetation clump size and number as indicators for alternative stable states in semi-arid ecosystems

Vegetation clump size and number as indicators for alternative stable states in semi-arid ecosystems | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Vegetation clump size and number as indicators for alternative stable states in semi-arid ecosystems Z. Chen , K. Halder , J. A. Bonachela doi: https://doi.org/10.1101/2025.10.03.680294 Z. Chen 1 Department of Ecology, Evolution, and Natural Resources, Rutgers University , USA 2 Department of Biological Sciences, University of Calgary , Canada Find this author on Google Scholar Find this author on PubMed Search for this author on this site K. Halder 1 Department of Ecology, Evolution, and Natural Resources, Rutgers University , USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site J. A. Bonachela 1 Department of Ecology, Evolution, and Natural Resources, Rutgers University , USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: juan.bonachela{at}rutgers.edu Abstract Full Text Info/History Metrics Preview PDF Abstract Dryland ecosystems are vulnerable to desertification, a pressing issue in the face of global climate change. In these ecosystems, vegetation often grows in spatially periodic patterns that differ as aridity increases (gaps, labyrinths and clumps), which has been widely studied theoretically, aiming to assess the proximity of the system to desertification. While some theoretical models predict smooth desertification transitions, most typically predict an abrupt transition linked to the possibility of two alternative stable states (desert and vegetated states), predictions that are yet to be confirmed empirically. If this bistability of alternative stable states occurs, however, environmental fluctuations and the history of the ecosystem determine which state ultimately materializes. This uncertainty makes it harder to predict desertification, compounding the challenges posed by it. Here, we combine empirical data and theoretical methods to investigate the links between bistability and vegetation spatial organization, which can help identify the presence of alternative stable states. We found that, although the presence of vegetation clumps is not indicative of bistability, changes in clump morphology can provide reliable indicators of bistability and an impending desertification transition. Thus, our methodology indirectly identifies whether desertification will occur abruptly, and whether restoration efforts should consider a potential ecosystem history-dependence. Introduction In arid and semi-arid ecosystems, where the average precipitation only meets less than 50% of potential plant water demand, the temporal and spatial variation of vegetation is highly influenced by water availability ( Snyder and Tartowski, 2006 ). In these environments, water becomes the main growth-limiting resource, which induces plant competition. Locally, however, the self-organization of vegetation into patches can also facilitate plant growth as the root systems break soil biocrust thus improving water infiltration, and plant shading contributes to preventing water evaporation. These interacting mechanisms (short-range facilitation and long-range competition, typically referred to as scale-dependent feedbacks, SDF; Rietkerk and van de Koppel (2008) ) ultimately determine the spatial distribution of vegetation, which in these ecosystems results in periodic vegetation patterns found globally such as gaps, labyrinths and clumps ( Deblauwe et al., 2008 ; Borgogno et al., 2009 ). By invoking SDF, very different mathematical models have successfully reproduced the sequence of vegetation patterns observed in arid and semi-arid ecosystems as rainfall decreases (gaps, labyrinths, and clumps, see insets in Fig. 1 ) ( Klausmeier, 1999 ; HilleRisLambers et al., 2001 ; von Hardenberg et al., 2001 ; Lejeune, 2002; Okayasu and Aizawa, 2001 ; Rietkerk et al., 2002 ; Gilad, 2004 ; Ursino, 2005 , 2007 ; Martinez-Garcia et al., 2023 ). In these models, clumped vegetation patterns emerge for the driest condition, ultimately leading to desertification if rainfall keeps decreasing ( Deblauwe et al., 2011 ), evidencing a link between vegetation self-organization and the desertification transition. Download figure Open in new tab Figure 1: Schematic diagram illustrating a discontinuous vegetation transition in response to changes in rainfall. The transition features a bistability regime with simulated patterns: the blue dots mark the start of the bistability regime (shaded area), while the red dots indicate where the system ultimately transitions to the desert state. Because of bistability, however, the transition may occur in the rainfall range between the two dashed lines. While some models predict a smooth, progressive transition, most models predict an abrupt desertification and thus a potentially catastrophic transition ( Fig.1 ) ( Martinez-Garcia et al., 2023 ). Moreover, although recent mathematical one-dimensional approaches suggest that vegetation reorganizes spatially to increase resilience in a changing environment, which could ultimately lead to either smooth or abrupt transitions depending on whether the reorganization is successful ( Siteur et al., 2014 ; Baastianen et al., 2018; Vanderveken et al., 2023 ), the lack of data has prevented empirical corroboration of any of the predictions above. Determining the shape of the transition is nonetheless important because, theoretically, discontinuous transitions are accompanied by bistability and hysteresis, which would make desertification even more difficult to predict and revert. Hysteresis (tendency for the system to follow different paths depending on its history i.e., its previous states) means that the rainfall levels needed to restore the system after desertification (blue points in Fig.1 ) would be higher from the rainfall levels for which the system became a desert (red points, same figure). In addition, bistability means that the system can potentially reach either one of two alternative stable states (or even switch dynamically from one to the other) for the same value of the environmental variables: in Figure 1 , the system could be either vegetated or desertic within the rainfall range between the red and blue point, with the realized state decided by the initial states or environmental fluctuations ( Lejeune et al., 2002 ; Lejeune et al., 2004 , Meron et al., 2004 , Sherratt and Lord 2007 ; Kéfi et al., 2010 ). Therefore, within the wide bistability range that precedes desertification, rainfall levels do not unequivocally predict ecosystem state, and restoring past environmental conditions does not necessarily ensure ecosystem regeneration. Thus, given the uncertainty and ecological consequences associated with bistability, determining whether an ecosystem is within a bistability regime is crucial for environmental management and conservation. Existing empirical data and theoretical studies have traditionally focused on developing methods to identify and monitor ecological transitions. To our knowledge, however, no existing method can account for or utilize spatial self-organization to determine whether the system is within a bistability regime. The sequence of patterns as rainfall decreases, and clumps more specifically, have generally been proposed as early-warning signals for the desertification transition ( Rietkerk et al., 2004 ; Scheffer et al., 2009 ; Dakos et al., 2011 ). However, clumps can be found within and out of the bistability range ( Fig. 1 ). Moreover, the analysis of pattern morphology using Fourier transforms ( Deblauwe et al., 2011 ) has been used in theoretical studies to anticipate desertification: the distance between clumps (and therefore the wavelength where the fundamental peak of the Fourier transform occurs, see Supplementary Information) has been shown to increase as the transition is approached ( Siteur et al., 2014 ). However, whether such wavelength changes can be related to bistability is unknown. EVI time-series and other remote-sensing products have been used to measure temporal and spatial indicators such as correlation, skewness, and variance that behave in characteristic ways as the transition approaches ( Veldhuis et al., 2022 ). However, some of these indicators (e.g., critical slowing down) require local perturbation analyses, which are labor-intensive, and/or require long time series that capture the transition, which are difficult to find ( Dakos et al., 2011 ). Importantly, the technical and conceptual challenges above have also prevented confirmation of an aspect as fundamental as to whether the desertification transition is abrupt (i.e., catastrophic) or smooth. This information, and the associated ecological consequences (e.g., the potential presence of hysteresis), can be helpful for and should necessarily be considered by managers and decision-makers ( Deblauwe et al., 2011 ; Martinez-Garcia et al., 2023 ). The presence of bistability could corroborate indirectly the presence of an impending catastrophic transition. Here, we aim to fill this gap by studying the relationship between vegetation patterns and bistability. We analyzed pattern morphology from satellite images using a study region in Sudan, and a model tailored to that region, to better understand the differences across resulting vegetation patterns and how those differences relate to bistability. We parametrized the model using available data for our focal system, calculating key parameters such as plant and water diffusivity, typically estimated qualitatively in the literature. We used Fourier transforms to monitor pattern regularity, and found a previously overlooked feature — a secondary fundamental peak — that can capture the changes in main properties of the vegetation spatial distribution as well as anticipate the transition. The changes in this peak as rainfall decreased can ultimately be linked to changes in clumps that are unique to the bistability regime, therefore pointing to simple metrics that can be used to identify the bistability regime (if it exists) and thus assess ecological uncertainty. We tested our predictions using satellite images from our focal region, and observed similar pattern behavior to those observed within the bistability regime of the model. Thus, according to our predictions the focal area would be undergoing bistability and shows behavior indicative of an impending abrupt desertification transition. Our methodology can help assess the risk of transition and ecosystem health for ecosystems where spatial self-organization occurs, and therefore can potentially be used to design and test management policies for such ecosystems. Methods Study area Our study area is located in the western part of South Kordofan in Sudan ( Fig. 2 ; 11° 42’ 0’’-11° 0’ 0’’ N, 27° 54’ 0’’-28° 0’ 0’’E; ( Esri, 2009 )), part of a larger region where vegetation patterns have been reported in the past ( Deblauwe et al., 2011 ). Topographically, the area has an overall slope of 0° (ArcGIS Online at a 25 m spatial resolution ( Esri, 2013 )). At the effective plant root zone depth, the soil texture is loam and the root zone plant-available water holding capacity ranges between 55 and 83 mm (Africa SoilGrids online resources at a 1 km spatial resolution ( Leenaars and Hengl, 2015 ; Leenaars and Hengl, 2017 )). Thus, there is no significant difference across the study area in terms of slope and soil properties ( Supplementary Figure S1.1 ). Download figure Open in new tab Figure 2: Optical satellite image of study area and different patterns across rainfall levels. The blue squares delineate the rainfall pixels from CHIRPS Precipitation Data (1990-2020). The selected 1km×1km vegetation patterns on the right are illustrative examples of clumps, labyrinths, and gaps (from top to bottom). The mean annual precipitation in our study area ranged between 220 and 590 mm over the last three decades (CHIRPS Precipitation Data: 1990-2020), with a monthly rainfall level that decreases from south to north. As rainfall decreases, vegetation patterns change from gaps to labyrinths to clumps ( Fig. 2 ). We selected 14 rainfall pixels on the same longitude (27° 54’ 0’’-28° 0’ 0“E), which ranged these different types of patterns. The information above justified focusing on rainfall as the main environmental feature driving the observed changes in patterns. Focal vegetation species We focused on Acacia mellifera , a native shrub in semi-arid regions in Africa and western Asia including Sudan (CABI, 2022). A. mellifera tends to release its seeds in a small radius around the parent shrub and this aggregated dispersal can result in clump-like aggregated growth consisting of different age groups ( Hagos, 2001 ). Various A. mellifera patterns have been reported in Kordofan, Sudan (Wickens and Collier, 1971; Deblauwe et al., 2011 ). Modified Klausmeier model To represent our system, we chose the generalized version of one of the simplest and earliest models for self-organized vegetation patterns in semi-arid ecosystems, the Klausmeier model ( Klausmeier, 1999 ; Ursino, 2005 ). The model describes the dynamics of water density ( w ; kg H 2 O m −2 ) and plant density ( p ; kg drymass m −2 ): The first term in Eq. (1) represents precipitation (only water input in this model), which occurs at rate r ( kg H 2 O m −2 year −1 ). The second and third terms represent, respectively, water loss due to evaporation and transpiration and plant growth water uptake, where l ( year −1 ) is the water loss rate and v (( kg drymass m −2 ) −2 year −1 ) is related to water maximum uptake rate. The last term represents water spatial movement, here encoded as diffusion at a rate D w ( m 2 year −1 ). This diffusive term differentiates this from the original Klausmeier model, which instead considered directional water movement (e.g., due to slope terrain, represented with an advection term, that is, a first spatial derivative). The first term in Eq. (2) represents plant biomass growth, where y represents the yield of plant biomass per unit water consumed, in kg drymass ( kg H 2 O ) −1 . The second term represents plant natural mortality (at a rate m ; year −1 ), the only cause of plant biomass loss considered here. The last term in the equation represents plant dispersal, which occurs at a rate D p ( m 2 year −1 ). The model considers that water uptake and growth benefit from the presence of other plants nearby (increased infiltration and reduced evaporation), reason why those terms are quadratic in p . Without dispersal of water and plants, the model results in homogeneous (i.e. non-patterned) systems, with desert being the only possible outcome (i.e. only stable state) for rainfall below , and bistability of vegetation and desert for larger rainfall values ( Klausmeier, 1999 ). The fact that our study area showed a diversity of patterns emphasized the need to find tailored values for the parameters and, in particular, for the dispersal (i.e. diffusivity) parameters. Model parametrization Values for the parameters of the model were obtained directly from the literature or through calculations tailored to our focal region. Water-related parametrization To replicate the gradient of rainfall within the study area, we explored fixed values of the rainfall parameter r between a realistic range based on the monthly rainfall data obtained from CHIRPS Precipitation Data. We alternatively implemented a stochastic, dynamic rainfall function generated from the monthly rainfall data, aiming to mimic the seasonality within the focal semi-arid climate (see SI Section 2 for details). The water loss rate ( l ) was obtained by dividing the theoretical water evapotranspiration for one year by the yearly water availability (precipitation). To calculate the theoretical water evapotranspiration, we used the Thornthwaite Potential Evapotranspiration equation (see SI Section 3a , Eq. S3.2 ) for Al Nahud and Muglad, two cities with data collection that overlapped with the study area ( Thornthwaite, 1948 ; Alam et al., 2013 ). The estimated water loss rate range was from 3.7 to 5.4 year −1 (see SI Section 3a for detailed calculation). On the other hand, models for vegetation in arid/semi-arid ecosystems often follow a qualitative approach to parameterizing or modeling water diffusivity ( D w ), or reuse values from classic work ( Sherratt, 2016 ). Here, we calculated water diffusivity aiming to keep the same thoroughness as with the rest of ecological parameters. To this end, and under the assumption that w represents surface water, we calculated the estimated range of water diffusivity by using the classic Manning equation for the velocity of a liquid ( Manning, 1891 ), considering flow velocity in all directions. As a result, we estimated a maximum theoretical water diffusivity and a more conservative value within the same order of magnitude, D w = 170 m 2 year −1 , which is also in line with the few estimates available in the literature. See SI Section 3a for details. Vegetation parametrization For the yield of plant biomass per unit water consumed ( y ), we used the standard assumption in the literature for encroacher species ( Terminalia sericea and A. mellifera ) of 1500 liters of water consumed to produce 1 kg of dry mass ( NAU, 2010 ). For the mortality rate ( m ), the existing literature provided a range of values for A. mellifera with 0.03 and 0.54 year −1 as extremes (see SI section 3b ). On the other hand, we interpreted diffusivity parameters for plants ( D p ) as the successful dispersal event for the focal species while considering germination rate and first-year seedling survival rate. Our calculations provided a plant diffusivity between 1.15-1.88 m 2 per year (see SI Section 3b for detailed calculation). To acknowledge that plant dispersal also results from shrub growth (i.e. clonal reproduction), which happens more slowly, we used a more conservative value within that order of magnitude ( D p = 1 m 2 year −1 ). Finally, the parameter accompanying the water uptake term ( v ) had units that are not typical for an uptake rate (( kg drymass m −2 ) −2 year −1 ), and therefore we could not find a value following ecological arguments as above. Instead, we used a value within the range in the original model, which was calibrated for semi-arid ecosystems ( Klausmeier, 1999 ). Our choice ( v = 57 ) resulted in patterns within rainfall ranges comparable to those observed in the focal area; other values of the parameter did not qualitatively alter our conclusions. See section S3 for further details and Supplementary Table S3.2 for final parameterization. Simulations We ran our simulations in a 250×250 spatial mesh, using periodic boundary conditions to reduce finite-size and border effects. The spatial resolution ( dx ) was 4 meters (i.e., simulated area of 1 km×1 km), and the time resolution ( dt ) was 0.001 year. The value of dx was small enough to resolve pattern motifs with several pixels, and large enough to cover a big area that provided a statistically adequate number of those motifs. We used the Jacobi method to discretize the two-dimensional Laplacian that represents diffusion, and the Euler method to integrate the partial differential equations of our model due to its simplicity and ease of implementation (especially under a stochastic environment, see above). Following the standard protocol, we started simulations at the homogeneous stable state for the given rainfall value; we further perturbed randomly chosen locations representing a small percentage of the system by adding a certain amount of plant biomass (Table S2). We then let the system run 100,000 years to allow stable patterns to emerge. We ran 25 replicates and stored final snapshots of the plant density matrix for further analysis. Pattern analysis Image preparation We analyzed the patterns observed in snapshots from our simulations, as well as 1km×1km optical satellite images with 1.2 m spatial resolution for the focal study area extracted from Google Earth Pro. Because we lacked information about biomass density for the focal area, to ensure a meaningful comparison with the model we first converted all the images to black-and-white scale, and thus analyzed images only based on the presence or absence of vegetation at a given pixel. For this binarization, we used a standard unsupervised classification algorithm (Kmeans, see details in SI Section 4 ). Power spectrum analysis To study changes in spatial periodicity, we calculated the radially averaged power spectrum associated with each snapshot collected at each rainfall level. We first obtained the two-dimensional (discrete) Fast Fourier Transform (2D FFT) using the well-established Python algorithm in the NumPy package ( Harris et al., 2020 ). We further computed the radially averaged power spectrum by averaging across all radii from the center of the 2D power spectrum in all directions, normalized by their respective bin-width areas. Peaks in this radially averaged power spectrum denote periodicity in the original pattern (see SI Section 5 ). Clump size and count distributions As justified below, measuring clump size and number can be informative when the rainfall levels are close to catastrophic levels. In order to measure the size and number of clumps from snapshots, we used the two-pass Hoshen–Kopelman algorithm from the SciPy Python library to identify and label each clump. For the binarized versions of both the simulation and satellite images, we calculated the average number, size, and maximum size over replicates for each rainfall level. Non-bistable vegetation model We also assessed which features observed with the tailored model (if any) are specific to bistability. To this end, we qualitatively compared the results we obtained with the Klausmeier model to those obtained with another model for regular vegetation patterns that does not show bistability, and that instead predicts a smooth desertification transition (Martinez-Garcia et al. 2013). See Supplementary Information for further details. Results In our study area, different levels of annual rainfall led to different types of patterns ( Fig. 3a - 3c ). The clump pattern occurred in areas with rainfall levels of 430.8 to 463.1 mm/year, while the labyrinth pattern was observed in areas with rainfall levels of 463.6 to 481.4 mm/year. We found dense labyrinths that also showed features of the gap pattern in areas with rainfall levels ranging from 496.0 to 508.9 mm/year. We measured the average percentage of each snapshot that is covered by vegetation, as a proxy for vegetation density. The representation of this metric against rainfall levels can be used to indicate whether the system is in a vegetated or desertic state, thus providing the so-called phase diagram of the system (see schematic example in Fig.1 ). For our study area, however, the phase diagram did not clearly indicate the character of the expected desertification transition nor whether there is bistability ( Fig. S1.2 ), which justified the use of a model tailored to our focal system to provide a more complete picture of the transition. Download figure Open in new tab Figure 3: Binarized versions of the real vegetation patterns from our study area and simulated patterns. Panels (a), (b) and (c) are satellite images derived from Google Earth pro after binarization (i.e., black-and-white transformation, see Methods), with annual precipitation 430.8, 478.0 and 508.9 mm/year, respectively. Panels (d), (e) and (f) correspond to the simulation with annual precipitation r = 430, 475 and 510 mm/year, respectively. For all snapshots, black indicates the presence of our focal vegetation and white means bare soil. Simulated vegetation patterns and phase diagram With the parametrization described in Methods and exploring different (fixed) precipitations, our simulations of the Klausmeier model were able to generate clump, labyrinth, and gap patterns at similar rainfall levels to those observed in the study system ( Fig. 3d - 3f ). The patterns typically required over 50,000 simulated years to become spatially stationary, after which the system showed only minor changes. The patterns obtained with a data-informed stochastic rainfall function (see SI Section 2 and Fig. S2.1 ) were very similar to those obtained with a fixed rainfall level ( Fig. S2.2 ). Thus, for simplicity we used fixed rainfall levels for our subsequent simulations and analyses. We first studied the desertification transition by constructing the phase diagram. To that end, we measured the stationary value of the vegetation density (spatial average of the density across pixels for a given snapshot, averaged across replicates) for a given rainfall level ( Fig. 4 ). Given the sensitivity to initial conditions expected for systems with bistability, we used two different initial perturbation values (1.0 kgm −2 and 0.1 kgm −2 ). Both perturbations resulted in similar patterns for rainfall levels under 375 and over 442 mm/year: the desert state was observed for rainfall levels under 375 mm/year; clumps, when rainfall ranged from 443 to 460 mm/year; labyrinths were observed when rainfall ranged from 460 to 510 mm/year; and gaps were found when rainfall ranged from 510 to 525 mm/year. However, the two perturbations led to very different results from 375 to 442 mm/year: the simulation with an initial perturbation value of 1.0 kgm −2 generated clump patterns, while the simulation with an initial perturbation value of 0.1 kgm −2 resulted in extinction. The system, thus, showed bistability between 375 and 442 mm/year (bistability regime). Download figure Open in new tab Figure 4: Phase diagram obtained by representing the values of mean vegetation density as a function of rainfall levels. The bistability regime is indicated in pink. Radially averaged power spectrum To understand if the properties of the patterns (e.g., periodicity) changed across the rainfall gradient, we calculated the power spectrum for each simulated replica’s final snapshot and studied the changes occurring across different rainfall levels. As expected, the heatmap of the power spectrum obtained across rainfall levels after averaging over replicates showed a first fundamental peak for all patterns, with spatial frequencies ranging between 15 and 19 cycles/km (red dashed line, Fig. 5 ). The first fundamental peak represents the wavelength of the periodic motif in the pattern. For the clump pattern, for example, the increase of spatial frequency with increasing rainfall availability signifies shorter distances between clumps. Also as expected, the power spectrum exhibited harmonic peaks at twice (light pink dashed line, Fig. 5 ) and thrice the spatial frequency of the first fundamental peak for all simulated rainfall levels. Download figure Open in new tab Figure 5: Heatmap of the semi-log transformed power spectra for simulated patterns. Coordinates are the spatial frequency of the pattern and rainfall levels (ranging from 375 to 525 mm/year). Color represents the value of the radially averaged power spectrum for the given spatial frequency (see definition in SI; yellow means highest and dark blue lowest). The first and second fundamental peaks are highlighted as dashed red and pink lines, respectively. Light blue lines were used to separate different types of patterns (clump: 375-460 mm/year; labyrinth: 460-510 mm/year; gap: 510-525 mm/year). To the left of the purple dashed line is the bistability regime (375-442 mm/year). Surprisingly, the semi-log plot also revealed a second fundamental peak between the first fundamental peak and its first harmonic peak, located around 26 cycles/km. This second fundamental peak occurred when the underlying pattern was clumps or gaps, but not labyrinths. Inside the bistability regime, the height of this second fundamental peak increased as rainfall decreased (see colors in heatmap), growth that was accompanied by a decrease in height of the first harmonic of the fundamental peak. See SI Section 5 for details. We tried to identify this changing secondary peak (and therefore the bistability regime) in images from our study area by calculating the power spectrum from each image. Unfortunately, the heterogeneity and stochasticity associated with real vegetation patterns and environmental conditions prevented the second fundamental peak or any of the harmonics from being discernible ( Fig. S5.6 ). To test whether this behavior was specific of the Klausmeier model or a signature of bistability, we repeated the analysis with the model in (Martinez-Garcia et al 2013), which produces similar vegetation patterns but predicts a smooth desertification transition ( Fig. S6.1 ). As shown in Fig. S6.2 , the radial spectrum showed in this model a secondary fundamental peak that followed a similar qualitative behavior as rainfall decreases. Differently from the Klausmeier model, however, the peak was noticeable for all the clump-pattern region, and the location of the first fundamental peak (and therefore the distance between clumps in the pattern) did not change with growth conditions in the model with a smooth desertification transition. Clump size and number change with rainfall in simulated and real snapshots As justified in SI, all the radial spectrum phenomenology above resulted from changes in clump size and number. Importantly, the link between the second fundamental peak of the power spectrum and the size and count of clumps within the pattern across examples (see SI Sections 5 and 6 ) indicates a potential more direct pathway to identify bistability. Thus, we measured and compared clump size and count across rainfall levels for simulations and study area. For the tailored Klausmeier model ( Fig. 6a-b and S6.3d-f), outside the bistability region the average clump size decreased as rainfall decreased, whereas the number of clumps remained constant. In the bistability range, however, the net change in clump size and maximum clump size remained small, whereas the number of clumps decreased. For the model with a smooth transition, however, clump size decreased with declining growth conditions whereas the number of clumps remained constant ( Fig. S6.3g,i ). This behavior therefore indicates differences between the changes that the clump pattern experiences in the presence and absence of bistability, which is confirmed by the observation that clumps outside of the bistability region in the Klausmeier model show similar behavior than in the non-bistable model (Martinez-Garcia et al 2013). Download figure Open in new tab Figure 6: Typical size (left column) and number of clumps (right column) as a function of rainfall level. Panels (a)-(b) show simulation results, whereas panels (c)-(d) correspond to satellite images. Note that “clump-labyrinth transition” corresponds to patterns that converge (very slowly) to labyrinths but still showed some clumps when the simulation was stopped. The end of the bistability region in simulations is marked with a small purple mark on the horizontal axis. For simulations, we show here a rainfall range equivalent to that observed in the study area; see Fig. S6.3 for a wider range of rainfall values, as well as a comparison with the results from the model in (Martinez-Garcia et al. 2013) and the measurement of maximum size. For satellite images ( Fig. 6c-d , and S6.3a-c), clump size slightly increased as rainfall increased, then jumped to higher values after clumps transitioned into labyrinths. Maximum size remained constant then jumped for labyrinthine patterns. Lastly, the number of clumps increased with rainfall levels but decreased for labyrinthine patterns. Discussion A tailored, ecologically grounded parameterization for the modified Klausmeier model We used here a modified version of the original Klausmeier model that replaces a directional advection term for water movement with a diffusion term. The latter was needed as the data we compiled for our study area indicated a region with little to no slope, indicating no directional water movement. Existing work has utilized this modification to the Klausmeier model, or similar models, to study vegetation patterns, typically in one spatial dimension. Here, the comparison with remote-sensing images of our study area required a focus on two dimensions. We informed the model with a bottom-up parametrization approach by which we obtained tailored, ecologically relevant parameter values, then studied whether patterns emerged for realistic rainfall levels. We identified a native plant species ( A. mellifera ) that is highly likely to underlie the observed vegetation patterns, and derived ranges for the parameters in Eq. (S2)-(S4) that corresponded to this species. For plant diffusivity, for example, we considered successful plant dispersal events by taking into account seed dispersal rate, outside canopy ratio, germination rate, clonal reproduction, and seedling survival rate. Regarding water diffusion, there is surprisingly very little information in the vegetation pattern modeling literature about the physical processes involved, and thus very little justification for the values typically used for water diffusivity ( Sherratt, 2016 ). The water variable in the Klausmeier model, for example, does not discern between superficial or subterranean water. This distinction is important because, ultimately, only underground water can be utilized by the root system. Here, we chose to interpret the water variable as superficial water and assume perfect water infiltration (i.e., all superficial water is available for the root system), which facilitated the calculation of the diffusion coefficient. Thus, we drew inspiration from the Manning equation for the velocity of a liquid and calculated a plausible range of diffusion rates. Despite the importance of the diffusivity parameters for vegetation pattern models, there is no definition or standardized way to calculate them, with further difficulties stemming from a lack of clarity regarding the definition of plant or water diffusion in these models. The calculations we used for our parameterization provide an ecological justification and a range that can be used to constrain parameter values in future theoretical work. Nonetheless, interdisciplinary efforts in which theoreticians and empiricists agree on definitions for the different components of these models have the potential to provide clear instructions to measure these important parameters ( Martinez-Garcia et al., 2023 ). Our parameterization for the model resulted in emerging patterns that were comparable to those of the study area, which served as a qualitative validation. Clumps and bistability can be characterized through previously unidentified features of the power spectrum The radially averaged power spectrum is typically used to differentiate between types of patterns based on the existence or location of a primary fundamental peak ( Deblauwe et al., 2011 ). This was also true in our study, as the location of the first fundamental peak jumped whenever patterns changed from clumps to labyrinths or from labyrinths to gaps (more subtly for real-world images). In addition, we further identified a second fundamental (i.e., non-harmonic) peak. Regardless of whether this is a real fundamental peak or one of the Fourier modes of the pattern, this peak and the primary fundamental peak showed a distinctive behavior when the system is bistable: an increase in rainfall caused an increase in the location of the first fundamental peak, and the decrease in height (and eventual disappearance) of the second fundamental peak. The latter resulted from different degrees of overlap between wells and peaks of the power spectrum as rainfall decreased (see SI Section 5 ), in turn reflecting changes in clump size and number. Thus, the phenomenology above should be expected from clump patterns and their changes within the bistability region. Indeed, the model with a smooth desertification transition showed instead a constant primary peak and a declining but non-negligible secondary peak as growth conditions improved. Thus, the different behavior of the primary fundamental peak observed in these two models identifies it as an indicator for bistability. Moreover, although the changes in this peak with rainfall have been suggested as an indicator of ecosystem health worsening as aridity increases ( Rietkerk et al. 2004 ; Deblauwe et al. 2011 ), our analyses rule out this possibility. Instead, our results reveal that the emergence and increase in height of the secondary fundamental peak is conserved across models. Therefore, this secondary peak can be used as an alternative early-warning indicator of impending desertification. The stochasticity inherent of real-world patterns, however, prevented the secondary peak from being discernible in the study-site images we analyzed. Importantly, the link between power spectrum and pattern features provides an understanding of the mechanisms underlying the behavior of the peak that enables bypassing the power spectrum (which is typically noisy for real systems, as shown here) and shifts the focus to properties that are easily measurable with satellite images. Clump size and count behavior are a signature of bistability and early-warning signals for desertification As changes in clump number and size are the primary factors contributing to the changes in the fundamental power-spectrum peaks, and these are general properties that can be easily captured from snapshots of vegetation patterns, we examined these properties for both the simulations and satellite images. Our analysis of simulated patterns with the tailored Klausmeier model shows that a decrease in rainfall leads to a reduction in clump spatial frequency (i.e., increase of wavelength) associated with a decrease in the number of clumps. Existing work has reported increasing pattern wavelength with aridity for a larger region containing our focal area in Sudan ( Deblauwe et al., 2011 ), and occurring in the modified Klausmeier model as the system approached the desertification transition ( Siteur et al., 2014 ; Bastiaansen et al., 2020 ). These observations refer to the location of the first fundamental peak of the power spectrum, which we also observe moves to lower frequencies (i.e., higher wavelengths) as rainfall decreases, indicating more distance between clumps and therefore fewer clumps within the observation window. As explained in those references, as aridity increases fewer clumps can survive, with a decline in number of clumps that becomes more subtle if rainfall changes are slow and gradual and/or for less-regular patterns (like the ones expected in real-world scenarios). The slight reduction of clump size initially observed as rainfall decreased within the bistability region is in line with an overall decline in biomass expected from increased water scarcity. Nonetheless, the net change in average clump size observed when considering the whole bistability region was almost negligible, as it was the case for the maximum clump size. This, and in particular the constant maximum size observed before the transition, may be indicative of a size above which the advantages of local facilitation cannot outweigh the limitations of both local and long-range competition. Additional empirical work, however, is needed to test this hypothesis. To the best of our knowledge, our study is the first to show a correlation between those changes and the existence of a changeable second fundamental peak. Moreover, changes in height for this second fundamental peak (and the first harmonic of the first fundamental peak) seem to be a more noticeable reflection of the changes in clump size and number than the changes in location of the first fundamental peak. Importantly, the behavior of these peaks and how it reflects on clump size and number can be used to characterize the bistability regime (i.e., the rainfall range in which both the vegetated and desertic state are stable). In this regime of the Klausmeier model, the number of clumps declined with declining rainfall, but average and maximum clump size barely did so. In contrast, the model without bistability showed markedly declining clump sizes but constant clump number. We hypothesize that, as rainfall decreases within the bistability region, the expected increase in stability of the desertic state and decrease of stability of the vegetated state lead to an increase in competition range that results in the fewer clumps observed in the Klausmeier model. This is initially supported by the fact that, in the model with a smooth desertification transition, competition range is fixed and thus worsening growth conditions barely affect clump numbers. Nonetheless, we are currently investigating the mechanisms underlying the differences between vegetation patterns in the presence or absence of bistability, which will be reported elsewhere. The behavior of average and maximum clump size as well as that of the number of clumps can therefore potentially be used to assess ecosystem health in semi-arid systems. After applying the method above to the satellite images from our focal area, we observed a behavior resembling that obtained with the Klausmeier model. The hump-shaped clump count, the near constant maximum clump size and the (still) slightly declining average size observed would indicate that the system is at a state similar to that of the model mid-way into the bistability region, when the first harmonic of the primary fundamental peak becomes negligible. Thus, with caution due to the apparent outlier observed in the number of clumps of regions of the focal system with low rainfall, the conclusion would be that the study area would be within a bistability regime. Consequently, the system would be at risk of desertification, as environmental changes and stochasticity may trigger the transition (i.e. a jump from the vegetated state to a desertic state) even if aridity levels are not unusually low or extreme. Management action would be required to avoid such outcome. The phenomenology observed here thus provides guidance to manage systems at risk of desertification by, for example, monitoring clump number and size over time ( in situ or through satellite images) as rainfall levels change. If, as expected from climate change, rainfall decreases over the years, and clump number decrease but average and maximum size remain approximately constant, the system is within a bistability regime and requires attention or management efforts to control large external disturbances like grazing or human activities, which can tilt the system towards desertification ( Kéfi et al., 2007 ). Thus, our study provides a complete workflow that combines empirical and theoretical methods to investigate pattern morphology and its relationship with bistability. Importantly, because this workflow identifies whether the system is within a bistability regime, it provides information about the character of the transition and path to recovery: if bistability is detected, then an abrupt transition and hysteresis are to be expected. Moreover, because the bistability regime is typically wide, our methodology offers a more conservative approach to identify that the health of the ecosystem is deteriorating. Therefore, the information provided by our approach can be used to develop policies that encourage sustainable land use practices and protect the livelihoods of local communities from land degradation, or to adapt agricultural practices to better suit the changing conditions of the land. Limitations Certain limitations of this study also call for caution when using our methodology, thus suggesting further research directions. Such limitations stem from missing information regarding the study area, as well as technical limitations. First, the lack of high-resolution historical satellite images (including for our focal area) made temporal analysis impossible due to the lack of such images before the 1990s. Also, despite our efforts to ensure that all site replicates had similar environmental conditions (other than rainfall), other factors for which we did not have information such as grazing and soil nutrients could also influence the pattern. Additionally, we assumed that the main plant species underlying the patterns was A. mellifera, but other species are present in the ecosystem, and calculated ranges for plant and water parameters based on proxy observations and data. Finally, we made important simplifying assumptions about both water and vegetation transport. All these aspects could be improved with ground-truthing observations and dedicated experiments in situ . Regarding the model, the phenomenological simplified depiction of the water-vegetation interactions, although mathematically convenient, leads to ecological conundrums. For example, the uptake term depends on a parameter v that does not have ecologically meaningful units and, although to some extent can be interpreted as an uptake rate, this lack of meaning prevented us from using bottom-up arguments to tailor its value to our specific system. Importantly, our main results were not altered qualitatively when testing other values for this parameter, which only moved the desertification and pattern transitions along the rainfall axis. Other, more detailed models could be used to represent vegetation in semi-arid ecosystems, although they are more difficult to parametrize due to the additional terms that they incorporate. One important aspect to keep in mind when using such models is that, although they may show similar vegetation patterns, the dynamical properties of the model can differ significantly (see SI Section 7 for an example). Future work will focus on how dynamical properties affect pattern behavior within the bistability region. Finally, a common feature of vegetation models with bistability is that the rainfall level at which patterns first emerge in the upper branch of the phase diagram (red point in Fig.1 ) is determined by the associated perturbation strength (here, 1.0 kg m −2 , see Methods). Our focal region did not cover a range of rainfall levels that included desertification, and thus we lacked information to parameterize this perturbation. Nonetheless, this parameter only affects the width of the bistability region, not the patterns or pattern transitions. Future empirical work is needed to directly link the pattern and associated power spectrum changes reported here to the actual desertification transition. Supplementary Information Section 1. Study area information Our focal area is located in the western part of South Kordofan in Sudan, shown in the World Topographic Map ( Figure S1.1a ; 11° 42’ 0’’-11° 0’ 0’’ N, 27° 54’ 0’’-28° 0’ 0’’E; Esri, 2013 ). This area is part of a larger region where vegetation patterns have been reported in the past ( Deblauwe et al., 2011 ). For the topographic information of our study area, the slope map ( Fig. S1.1b ) was derived from ArcGIS Online at a 25 m spatial resolution ( Esri, 2013 ). We derived soil water holding capacity ( Fig. S1.1c ) and soil texture data ( Fig. S1.1d ) at a 1 km spatial resolution from the Africa SoilGrids online resources ( Leenaars and Hengl, 2015 ; Leenaars and Hengl, 2017 ). The vegetation cover percentage was calculated for each rainfall pixel that showed clump-like patterns, which corresponded to rainfall levels ranging from 430 to 465 mm/year ( Fig. S1.2 ). Within each rainfall pixel, we took five replicates at a 1 km spatial resolution and calculated the mean vegetation cover percentage across the replicates. Error bars show the corresponding standard deviation. The clump-labyrinth transition was labeled with triangles, indicating that the pattern is a mixture of clump and labyrinth patterns. This phase diagram neither indicated nor ruled out the presence of bistability between vegetated and desertic states for any rainfall levels. Download figure Open in new tab Figure S1.1: Focal area and main topographic properties. The study area (shown in panel (a)) is flat (slope between 0° and 1°, panel (b)). Panels (c) and (d) show that soil water-holding capacity and soil texture are similar at the effective plant-root zone depth across the region. Download figure Open in new tab Figure S1.2: Phase diagram obtained representing the values of the mean vegetation cover percentage as a function of rainfall levels for our focal area. Section 2. Rainfall functions To check whether simulated patterns are robust under stochastic and seasonal changes in rainfall, we replaced the constant value used for the rainfall parameter with a stochastic rainfall function generated from the monthly rainfall data. This function allowed us to resolve seasonality within the focal semi-arid climate, as well as the inherent climatic fluctuations. To implement the seasonality part of the function, we first classified pixels according to the observed pattern (clump, labyrinth, or gap). We then estimated the medians and interquartile ranges of monthly precipitation over a 30-year period, which we fitted to a skewed-Gaussian function that provided us with a seasonal resolution of rainfall. The resulting stochastic rainfall function was (see Fig. S2.1 ): Where r is annual precipitation rate (mm/year), M is the specific month of the year, ω = 1.8, ψ = 8.8 and λ = −1.5. To represent the stochasticity inherent to rainfall, we added a noise component to Eq. (S2.1) . Download figure Open in new tab Figure S2.1: Rainfall function R(M) obtained with data collected at sites that show different types of patterns (clump, labyrinth, and gap). The dots represent medians across years, and error bars represent interquartile ranges. We fitted a skewed-gaussian function to represent seasonal precipitation. Fig. S2.2a-c shows the patterns observed in the regions where the rain data were collected; panels g-i show the simulated patterns obtained with these rainfall curves and added stochasticity. Thus, we added a random number generated from a normal distribution with mean 0 and standard deviation 0.3 R(M), the latter representing the observation that drier months show a lower rainfall variability (see Fig. S2.1 ). The resulting random rainfall values were distributed normally with mean R(M) and standard deviation 30% of the given rainfall. We used these simulations with stochastic rainfall to check whether vegetation patterns changed seasonally or remained stable across seasons. As shown in Figure S2.2 , both the stochastic rainfall simulations and fixed rainfall simulations exhibited similar patterns at similar annual rainfall levels, indicating that the vegetation patterns generated by this model are robust under stochastic and seasonal changes in rainfall. Therefore, for simplicity we used fixed rainfall values for our subsequent simulations and analyses. Download figure Open in new tab Figure S2.2: Simulated vegetation patterns and real patterns from our study area. Panels (a), (b) and (c) are satellite images derived from Google Earth pro, with annual precipitation 430.8, 478.0 and 508.9 mm/year, respectively. Panels (d), (e) and (f) correspond to simulations with fixed annual precipitation r = 430, 475 and 510 mm/year, respectively. Panels (g), (h) and (i) correspond to the stochastic rainfall simulations with annual precipitation r = 430, 475 and 510 mm/year, respectively. For all snapshots, black indicates the presence of vegetation and white means bare soil, with panels resulting from the binarization protocol explained in Section 4 . Section 3. Model parametrization For all model parameters, we used available information to constrain the range of possible values. Our main goal was to identify a plausible order of magnitude for each parameter. Whenever data were not available, we used theoretical calculations. In the absence of either option, which was the case for one parameter only, we resorted to sensitivity analyses to select a value that helped replicate the spatial patterns of our focal region. 3a. Water-related parameters Water diffusivity (D w ) To calculate the diffusivity parameter, we first calculated the cross-sectional average velocity of water flow using the Manning equation ( Manning, 1891 ): where V is the cross-sectional average velocity (m/s), n is Manning’s Roughness Coefficient for different surface material, R is the hydraulic radius (m, see Table S3.1), and S is the slope (elevation difference/section length) View this table: View inline View popup Download powerpoint Table S3.1: Chanel section geometry and associated parameters. To calculate the hydraulic radius, we assumed that the intense rainfall during the rainy season creates a very shallow sheet flow with width B = 1 m and depth h = 0.0001 m in a horizontal area with very small local slope (0-0.1%). These two values allowed us to calculate the hydraulic radius for a channel with a given rectangular cross-section (Table S3.1). For Manning’s roughness coefficient, we assumed that the surface material for the shallow sheet flow was similar to a clean and straight natural stream with no rifts or deep pools, which led to a value of 0.03 ( Arcement and Schneider, 1989 ). With this information and a local slope of 0.05%, Eq. (S3.1) resulted in a flow velocity of 0.00161 m/s. Water diffusivity was determined by taking into account the flow velocity in all directions. Naively, because we assumed that w represents surface water (see main text), if water can travel 0.00161 m per second, then water moving in all directions will initially cover a surface π (0.00161) 2 per second. The calculated value was thus 255 m 2 per year (similar to the one used in Siteur et al. (2014) ), which we considered as a maximum possible diffusivity since the area added per unit time is that of an annulus instead of a disc. With this in mind, we chose the more conservative D w = 170 m 2 per year, a value within the same order of magnitude that also ensured numerical stability, is close to the one used in Klausmeier (1999) repeatedly adopted in theoretical work, and is in line with the value calculated in Ursino (2005) . Water loss rate (l) To calculate the loss rate, we first calculated the potential evapotranspiration (PET) using the following equation ( Thornthwaite, 1948 ): PET is the estimated potential evapotranspiration (mm/month), T d is the average daily temperature (degrees Celsius) of the month, N is the number of days in the month, L is the average day length (hours) of the month, I is a heat index which depends on the mean temperature, and the exponent a is given by the formula: We set the averaged day length of the month ( L ) to 12 hours and the number of days in the month ( N ) to 30 days. For the rest of parameters, we used two cities within the study area as a reference, Al Nahud and Muglad, for which data were available. The average daily temperature ( T d ; Al Nahud: 27.8 ℃; Muglad: 28.5 ℃) and precipitation data (1960-1990) for these two cities (Al Nahud: 315 mm; Muglad: 471 mm) were derived from Alam et al. (2012). Based on the average daily temperature and the heat index table ( Seed and Malik, 1968 ), the heat index I for Al Nahud was 13.49 and for Muglad was 13.94. Therefore, we could calculate the PET for each city and get the estimated water loss rate by dividing the annual PET to the mean annual precipitation. The estimated water loss rate range was from 3.7 to 5.4 year −1 . 3b. Plant-related parameters Plant diffusivity (D p ) We interpreted the diffusivity parameter for plants ( D p ) as resulting from successful dispersal events for the focal species, considering germination rate and first-year seedling survival rate. There is no literature indicating the specific seed dispersal strategy of A. mellifera , but it has been pointed out that the seeds of A. mellifera disperse over a relatively small distance near their parent plant ( Donaldson, 1969 ). This means that the primary seed dispersal strategy can be barochory, and thus the seed dispersal of A. mellifera would primarily rely on gravity and assisted by wind. The average monthly wind speed in the last two decades for our area of interest was under 5 m/s (FNMOC 10 m Surface Winds, erdlasFnWind10). Moreover, according to the Beaufort Wind Scale, the sample area typically experiences light to gentle breeze, which is unlikely to make a long-distance seed dispersal. Thus, we assumed that the seed dispersal distance outside the canopy equals the radius of the canopy (2.04 m, Figure S3.1 ) and calculated the circle of diffusion area based on the schematic diagram for the study of the effect of wind direction on seed dispersal (Hagos, 2002). Download figure Open in new tab Figure S3.1: Schematic diagram showing A. mellifera seed dispersal. The annulus in light gray represents the area of successful seed dispersal. Therefore, we considered the seeds that disperse outside the canopy as a potential successful dispersal event and assumed the dispersal radius to be twice that of the canopy, which leads to a dispersal area of 39.22 m 2 . We further weighted the success of dispersal by the survival probability after a year (and thus establishment probability). On one hand, 14% of A. mellifera seeds disperse outside the canopy (Hagos, 2002); on the other hand, the germination rate of A. mellifera seeds is around 50-80% ( Bein et al., 1996 ), and the calculated first-year survival rate of seedlings outside the canopy is around 42% ( Joubert et al., 2013 ). Altogether, the plant diffusivity calculated from the data above ranged 1.15-1.88 m 2 per year, resulting from the following equation: In addition, clonal reproduction of the shrub also contributes to its spatial spreading even in the absence of wind. To also account for this possibility, we used the more conservative value D p = 1 m 2 per year, in line with the one used in Klausmeier (1999) . Other plant-related parameters For the yield of plant biomass per unit water consumed ( y ) we used the assumptions made in the literature for encroacher species ( Terminalia sericea and A. mellifera ) of 1500 liters of water consumed to produce 1 kg dry mass ( NAU, 2010 ), i.e. y = 0.00067 kg drymass per kg water . In addition, two studies have mentioned the mortality rate of A. mellifera : Joubert et al. (2013) reported that around 1% of labeled seedlings survived after 6 years, which means the probability for A. mellifera seedlings to survive for 6 years is 0.01. Assuming that the survival probability does not change over the years and life stages, A. mellifera has a probability 0.46 to survive for one year, and therefore the annual mortality probability is 0.54. A later study showed that the mortality rate of shrubs over a 35-year period is around 61.3% ( Joubert et al., 2017 ), which means that the probability to survive for one year is 0.97, and thus the annual mortality probability equals 0.03. These data provided a range of values for A. mellifera mortality probability ( m ) with 0.03 and 0.54 per year as extremes. The parameter accompanying the water uptake term ( v ), on the other hand, has units that are not typical for an uptake rate (( kg drymass m −2 ) −2 year −1 ), and therefore we could not find a value from available data or by following ecological arguments as above. A sensitivity analysis revealed that the main role of this parameter is to shift the rainfall values where patterns and desertification occur, but does not change the patterns. Due to the lack of information, we chose this parameter to ensure that (once the rest of parameters were fixed to realistic values, see above and Table S3.2), the range of rainfall values for which patterns emerged was comparable to those observed in the focal system. Thus, we used a value of v = 57 , but other values of v that we tried ( v = 25, 50, 60, 75, 100 ) did not alter our conclusions qualitatively. View this table: View inline View popup Table S3.2: List of model parameters, values and units. For parameter values or ranges directly obtained from the literature, we included the corresponding reference; for values that we calculated based on ecological and/or physical arguments, we included the section with the corresponding explanation of the steps. Section 4. Processing of vegetation images through unsupervised classification and clump detection Given the lack of information to map image colors into vegetation density, it is standard practice to binarize images before analyzing them. Here, in order to detect and discriminate vegetation from the surrounding bare soil, we used an unsupervised classification algorithm known as KMeans ( Arthur and Vassilvitskii, 2007 ) that is widely used in image analysis, particularly in remote sensing applications for land cover classification, change detection, and vegetation analysis ( Celik, 2009 ; Olaode et al., 2014 ; Ahmed and Akter, 2017 ; Wan et al., 2018 ). KMeans clustering involves the partition of a dataset into k optimal distinct clusters, through the iterative minimization of a function ψ . More specifically, when provided an integral input k on a dataset χ ∈ ℝ d , the KMeans algorithm calculates k cluster centers C to minimize the sum of Euclidean distances from each point in the dataset to its nearest cluster center. Mathematically, this is achieved by minimizing the following function ( Arthur and Vassilvitskii, 2007 ): Because our aim is to binarize images, in our case we applied a binary KMeans classification ( k = 2) to our satellite and simulated images, which resulted in a classification of image pixels into two clusters, as either “low” intensity pixels (representing bare soil) or “high” intensity pixels (representing vegetation) (see example in Fig S4.1). For our analyses, we used the KMeans implementation from the scikit-learn Python library (sklearn.cluster.KMeans), which features several enhancements over naïve KMeans. Specifically, we used the KMeans++ initialization strategy, which improves clustering performance and convergence stability by selecting initial cluster centers that are probabilistically distant from one another ( Arthur and Vassilvitskii, 2007 ). Once images were binarized, we proceeded to identify and quantify the motifs, i.e. vegetation clumps. Thus, we used the two-pass Hoshen–Kopelman algorithm from the SciPy Python library (scipy.ndimage.label) to label individual classified vegetation clumps as connected components ( Hoshen and Kopelman, 1976 ). This labeling allowed us to count and reliably calculate the size of each detected clump. As a side note, KMeans (similarly to other clustering algorithms) is sensitive to the choice of k : higher choices often lead to the detection of more features, but at the risk of overfitting and loss of statistical power. To ensure that our choice of k = 2 was appropriate, we performed a silhouette analysis ( Rousseeuw, 1987 ) and calculated the percentage of unexplained variance for analyses using k ranging from 2 to 10 (not shown). For both the satellite data as well as the simulated data, we found that k = 2 was the most parsimonious choice. Download figure Open in new tab Fig S.4.1: Example of snapshot from our study area (left) and corresponding binarized version using the KMeans++ algorithm (right). Section 5. Fourier transform and Power Spectrum A Fourier transform shows one peak for every direction where there is periodicity in the corresponding image. The location of the peak represents how many times the motif is repeated within the picture (spatial frequency). For example, for the clump pattern ( Fig. S5.1a ), the Fourier transform ( Fig. S5.1b ) shows a peak at a distance from the origin that reflects the number of times that a clump is found in a given direction. Thus, tracing parallel lines that cross clumps in a given direction, the distance to the origin of each Fourier transform peak corresponds with the number of lines with such orientation within the picture; the orientation determines the angle at which the peak is located in the Fourier transform. The radial power spectrum integrates the Fourier transform across all possible angles and, in consequence, only depends on the distance to the origin (i.e. spatial frequency). Putting all together, any links connecting first neighbors in the original pattern (e.g. pink and brown lines in Fig. S5.1a ) lead to peaks of the Fourier transform at the shortest possible distance to the origin (pink and brown circles in Fig. S5.1b ) and contribute to the same peak in the radial power spectrum (brown solid box in Fig. S5.1c ). This peak is the so-called fundamental peak, which determines the main (spatial) frequency within the picture. Because the fundamental peak provides information about the typical distance to first neighbors, its presence and location are used to characterize regularity within images/patterns ( Deblauwe et al., 2011 ). Theoretically, this peak (and its harmonic peaks, located at twice, thrice, etc. the distance of the first peak) are typically the focus of research. The clump and gap patterns, however, also show an additional and typically ignored feature. When considering lines that connect secondary neighboring clumps (green and red lines in Fig. S5.1a ), because of the overlap of one row and the following one, the distance between lines is smaller, i.e. there are more lines than for first neighbors. These lines lead to farther peaks in the Fourier transform (red and green circles in Fig. S5.1b ), and thus produce a second fundamental peak in the power spectrum (dark red box in Figure S5.1c ) that is close to, but different from, the first harmonic of the fundamental peak (dashed-line box). Technically, this peak is a linear combination of modes (e.g. in Fig. S5.1b , the coordinates of the red circle can be obtained by adding the vectors of peaks from the first neighbors, specifically twice the vector of the purple circle plus the vector symmetric to the brown circle). Download figure Open in new tab Figure S5.1: Link between pattern and power spectrum for an idealized (i.e., perfectly hexagonal) clump pattern. Panel (a) is the original image with red, pink, brown, and green lines indicating directions/lines crossing the center of the clump, the number of lines corresponds to spatial frequency in panels (b) and (c). Panel (b) shows the corresponding 2D Fast Fourier Transform (FFT), with dark blue indicating lower values and light blue indicating higher magnitude values. The horizontal and vertical axes correspond to the spatial frequencies in the horizontal and vertical directions respectively, and the color of circles corresponds to those of the lines in panel (a). Panel (c) shows the radially averaged power spectrum calculated from the 2D FFT. The first fundamental peak, denoted by a brown solid box, represents the spatial frequency of the brown and pink lines in panel (a) (and circles in panel (b)), while its first harmonic peak is identified by the dashed brown box. The second fundamental peak, indicated by the red box, corresponds to the spatial frequency of the red and green lines in panel (a) (circles in panel (b)). To investigate the behavior of this second fundamental peak, we constructed perfectly hexagonal clump patterns and first varied the size of the clumps only (i.e., number and position of clumps remained constant, see Fig. S5.2 ). As expected, because the location of the clumps remains fixed, the location and height of the first fundamental peak in the power spectrum did not change. Our analyses showed that the location of the second fundamental peak also remained fixed at a given spatial frequency; however, its height changed non-monotonically: it decreased with clump size until clump size was ∼20 pixels, then started increasing again. At that size threshold, the first harmonic of the first fundamental peak started decreasing. These changes corresponded with changes in the intensity of the second line of peaks in the Fourier transform (see red boxes for the second fundamental peak and green boxes for the first harmonic of the first fundamental peak). Download figure Open in new tab Figure S5.2: Power spectra (left) and Fourier transform (center) for idealized clump patterns (right) with identical number and location, but different sizes (clump diameter changes from 22 to 18 pixels from top to bottom). The peaks in the power spectrum and associated peaks in the Fourier transform associated with the second fundamental peak are denoted by red boxes, and those associate with the first harmonic of the first fundamental peak by green boxes. To better understand these changes in the height of the second fundamental peak, we then reduced the number of clumps and monitored the power spectrum as a function of both clump size and number ( Figure S5.3 - S5.5 ). The power spectrum obtained with just five clumps showed now a well instead of a peak, with the location of the well increasing as clump size decreased and the depth of the well being slightly affected by size changes ( Fig. S5.3 ). We then removed four of the clumps to leave the central one, which produced a well at the same spatial frequency as that for the five clumps ( Figure S5.4 ), with a location that also increased as clump size decreased but a depth that barely changed. Finally, when we instead multiplied by nine the number of clumps (i.e. 45 clumps in total, Figure S5.5 ), the depth and position of the well was much more sensitive to clump size: as size decreased, the well moved toward higher spatial frequencies, leading to the decrease and eventual disappearance of the second fundamental peak for intermediate clump sizes, then its reappearance and decrease of the first harmonic of the first fundamental peak as clump size is decreased further. Taking all the above into consideration, we conclude that changing size and number in such patterns leads to changes in the location of the well and second fundamental peak and how much they overlap (and therefore how much they offset each other); the conclusion is also applicable (albeit for smaller clump sizes) to the first harmonic of the first fundamental peak. This explains why the heights of the second fundamental peak and the first harmonic change, and point out that motif size and number are important features to monitor when analyzing vegetation patterns such as clumps and gaps. Download figure Open in new tab Figure S5.3: Power spectra (left) and Fourier transform (center) for idealized clump patterns (right) with identical number and location, using only five clumps. Note that the second fundamental peak is not present, but instead a well between the first peak and its harmonic emerges (see green box). Download figure Open in new tab Figure S5.4: Same as in Fig. S5.3 but using a single clump in an otherwise empty system. Download figure Open in new tab Figure S5.5: Power spectra (left) and Fourier transform (center) for idealized clump patterns (right) with identical number and location, using with 45 clumps. The second fundamental peak is present in the power spectrum and marked with a red box, as well as the associated peaks in the Fourier transform. Note that, for the case with 20-pixel-diameter clumps (middle row), both the second fundamental peak and the well overlap, resulting in a shallower well in the power spectrum and lower intensity in the peaks of the Fourier transform (see red boxes). Power spectrum for images from the focal area We calculated the power spectrum associated with images from our study area, to explore whether they would allow us to identify a changing secondary peak and therefore the bistability region. Unfortunately, due to the heterogeneity and stochasticity associated with real patterns and environmental conditions, the second fundamental peak was not discernable (see Fig. S5.6 , which shows the radial-averaged power spectrum of clump patterns with rainfall levels ranging from 430.8 to 463.6 mm/year). Download figure Open in new tab Figure S5.6: Radially averaged power spectrum of the focal area’s clump patterns observed for rainfall ranging from 430.8 to 463.6 mm/year. Section 6. Results with a non-bistable model for vegetation patterns in semi-arid ecosystems To investigate which features of our results are specifically linked to bistability, we repeated the analysis with a model (Martinez-Garcia et al. 2013) that also replicates all three vegetation patterns observed in semi-arid ecosystems. Differently to the Klausmeier model, however, patterns emerge purely from vegetation growth and competitive interactions (i.e. it does not implement any facilitation). There is only one equation, which represents vegetation biomass, but the equation encodes interaction between plants through a non-local competition kernel that affects vegetation growth rate within a constant interaction radius. Rainfall or aridity are not directly included in the model, but rather indirectly represented through the maximum growth rate, which is the parameter whose change triggers desertification or pattern changes. The resulting desertification transition is smooth, i.e. there is a unique value of the average vegetation density expected for each maximum growth rate (see Fig.S6.1). Thus, there is no bistability when the system approaches the transition. Therefore, contrasting the results of this model with those obtained with the Klausmeier model can help understand the particularities associated with bistability. We used a generic parametrization for this model as we lacked information to tailor the parameter values to our study area. Thus, using the parametrization proposed in (Martinez-Garcia et al. 2013), which was inspired by vegetation in semi-arid ecosystems, we obtained the power spectra depicted in Fig. S6.2 . Similarly to our tailored Klausmeier model, there is a secondary peak next to the first harmonic, and both change heights as growth conditions change. Differently from the Klausmeier model, however, the location and height of the first peak did not change. All together indicates that clumps changed in size but not in location (and therefore not in number, see Section 5 above). Indeed, Figs. S6.3g-i show that both average and maximum clump size decreased with worsening growth conditions (lower maximum growth rate parameter, β ) whereas the number of clumps remained approximately constant. This behavior contrasts with that observed in the bistability region of the Klausmeier model, but agrees with the behavior observed outside of that bistability region. Therefore, these results suggest that bistability drives the decrease in the number of clumps, approximately constant in size, observed in the Klausmeier model and in the study area as aridity increases. Download figure Open in new tab Figure S6.1: Phase diagram for the model in (Martinez-Garcia et al. 2013). The desertification transition, which occurs at a maximum growth rate 0.5 is smooth and therefore there is no bistability; only one average vegetation density is possible for every value of the maximum growth rate, which makes this transition much more predictable than an abrupt one. Download figure Open in new tab Figure S6.2: Heatmap for the radial power spectrum for the model in (Martinez-Garcia et al. 2013). Colors and axes like in Fig. 5 (main text). Note that the frequency at which the first fundamental peak occurs does not change with growth conditions, differently from the Klausmeier model. Download figure Open in new tab Figure S6.3: Average clump size (left), maximum clump size (center) and number of clumps (right) as a function of growth conditions for the study area (top row), the Klausmeier model (middle), and the model in (Martinez-Garcia et al. 2013) (bottom). Panels d) and g) show not only the mean value (squares) but also full distribution of clumps for each replica and each rainfall level. The dashed line in the middle row separates the bistability region ( R 443 mm/yr). Section 7. Preliminary results with the Rietkerk model The Rietkerk model ( Rietkerk et al. 2002 ) is a three-equation model that, in addition to tracking vegetation density and surface water, also monitors the dynamics of soil water. Thus, its equations use terms that are more ecologically detailed than those in the Klausmeier model, although this increased realism comes with at the price expense of a much larger set of parameters. The Rietkerk model is a unique vegetation model because of its dynamical properties: on the one hand, like in the Klausmeier model, there is a range of rainfall values for which perturbing the desertic state triggers the emergence of clumps, evincing the bistability of vegetated and desert states within that range; on the other hand, when diffusivity is neglected and therefore space is not relevant, the homogeneous system does not show bistability (see Fig. S7.1 ). The latter property makes the Rietkerk model different from the Klausmeier model, for which bistability of vegetation and desert occurs both homogeneous and patterned versions of the model; the model thus shares this particular dynamical property with the model in (Martinez-Garcia et al. 2013), for which bistability does not occur in the homogeneous model (although, in the (Martinez-Garcia et al 2013) model, bistability does not occur even when patterns emerge). We ran preliminary analyses to test whether these dynamical differences between the Klausmeier and the Rietkerk models affected the behavior reported in the main text. The lack of information for the more complicated terms and parameters of the Rietkerk model prevented us from following the thorough bottom-up parameterization approach that we used with the Klausmeier model. Instead, we ran our tests utilizing parameter values that have been used in the past for arid regions in Africa ( Bonachela et al. 2015 ). As shown in Fig. S7.2 , the behavior of the power spectrum as rainfall changes is qualitatively similar to that of the Klausmeier model ( Fig. 5 ). First, the spatial frequency at which the first fundamental peak is located increases with rainfall. Moreover, the radial spectrum shows a second fundamental peak that differs from the first harmonic of the primary fundamental peak, and the height of this peak and the harmonic change as rainfall increases. In addition, these changes correlate with a change in clump number ( Fig. S7.3 ): as aridity increases (i.e. rainfall decreases), the number of clumps decreases. These similarities reinforce the idea of using the behavior of the first and second fundamental peaks, and that of the number of clumps, as indicator for bistability. The average and maximum clump size, on the other hand, both show a marked decline, which departs from the subtle changes observed for the Klausmeier model and seems closer to the behavior observed for the model with a smooth transition. We hypothesize that differences between the dynamical properties of the model in the presence and the absence of patterns affect the behavior of clump morphology when patterns are present, but we plan to investigate this technical aspect in future work. Note that this behavior is observed within the bistability regime which, for this model, extends up to rainfall levels R = 1.00 mm/d and includes all clump patterns. Outside the bistability regime, clumps converge to labyrinths, with the corresponding increase in clump size and decrease in number. Download figure Open in new tab Figure S7.1: Phase diagram for the Rietkerk model (left) and, for reference, the Klausmeier model (right), for both the full model (points) and the model without diffusivity (i.e. non-spatial model, for which the result of the interactions between vegetation and water is a homogeneous state, purple lines). For the Rietkerk model (left), in the absence of patterns (i.e. when space is not relevant) the system does not show bistability but rather a smooth transition from vegetated to desertic state; note how the purple line only shows one value of vegetation for each rainfall value, as opposed to the two possible values shown by the purple line of the Klausmeier model (right). When considering vegetation and water dispersal (points), the Rietkerk model behaves like a scale-dependent feedback model and patterns develop after the homogeneous state is perturbed, giving rise to an abrupt desertification transition and bistability for rainfall within [0.60,1.0] (shaded region). Differently from the Klausmeier model, in the Rietkerk model clumps are only observed within the bistability region. For the non-spatial models (purple lines), the continuous transition in the Rietkerk model is also referred to as a transcritical bifurcation, and the discontinuous transition in the Klausmeier model as a fold bifurcation (where the dashed line represents unstable values/states of the system). Download figure Open in new tab Figure S7.2: Heatmap for the radial power spectrum for the Rietkerk model. Colors and axes like in Fig. 5 . Note the similarities with that figure for the Klausmeier model, with frequency at which the first fundamental peak occurs changing with rainfall, and the existence of a secondary peak whose height decreases and ultimately becomes negligible as rainfall increases within the bistability region. Download figure Open in new tab Figure S7.3: Average clump size (left), maximum clump size (center) and number of clumps (right) as a function of rainfall for the Rietkerk model. As in Fig. S6.3 , panel a) shows both mean value (squares) and full distribution of clumps for each replica and each rainfall level. Note how, similarly to the Klausmeier model (and as expected from Fig. S7.2 ), the number of clumps decreases with lower rainfall; however, the decrease in clump size is marked, similarly to the behavior observed with the model in (Martinez-Garcia et al 2013). Acknowledgements We would like to thank R. Martínez-García for his comment on an early version of the manuscript and the code to run the model in Martínez-García et al. (2013) . This work was supported by NSF grant DMS-2052616 to J.A.B. Funder Information Declared U.S. National Science Foundation, https://ror.org/021nxhr62 , DMS-2052616 References ↵ Bastiaansen , R. , Doelman , A. , Eppinga , M. B and Rietkerk , M . ( 2020 ). The effect of climate change on the resilience of ecosystems with adaptive spatial pattern formation . Ecology Letters , 23 ( 3 ), 414 – 429 . OpenUrl PubMed ↵ Borgogno , F. , D’Odorico , P. , Laio , F. and Ridolfi , L . ( 2009 ). Mathematical models of vegetation pattern formation in ecohydrology , Rev. Geophys ., 47 , RG1005 . OpenUrl ↵ Dakos , V. , Kéfi , S. , Rietkerk , M. , van Nes , E. and Scheffer , M. ( 2011 ). Slowing down in spatially patterned ecosystems at the brink of collapse . The American Naturalist 177 : E153 – E166 . OpenUrl CrossRef PubMed Web of Science ↵ Deblauwe , V. , Barbier , N. , Couteron , P. , Lejeune , O. and Bogaert , J . ( 2008 ). The global biogeography of semi-arid periodic vegetation patterns . Glob Ecol Biogoegr 17 : 715 – 723 . OpenUrl ↵ Deblauwe , V. , Couteron , P. , Lejeune , O. , Bogaert , J. and Barbier , N . ( 2011 ). Environmental modulation of self-organized periodic vegetation patterns in Sudan . Ecography , 34 , pp. 990 – 1001 . OpenUrl CrossRef Web of Science ↵ ESRI . “World Imagery” . Dec 12, 2009 . Sources: Esri, DigitalGlobe, GeoEye, i-cubed, USDA FSA, USGS, AEX, Getmapping, Aerogrid, IGN, IGP, swisstopo, and the GIS User Community . Available online: https://www.arcgis.com/home/item.html?id=10df2279f9684e4a9f6a7f08febac2a9 (accessed on March 23, 2023). ↵ ESRI . “Slope map” [digital map]. “Terrain” . Dec 30, 2013 . Source: Airbus, USGS, NGA, NASA, CGIAR, NLS, OS, NMA, Geodatastyrelsen, GSA, GSI and the GIS User Community . Available online: https://www.arcgis.com/home/item.html?id=af25a795273440deb449b336543602be . (accessed on March 23, 2023) ↵ Gilad , E . ( 2004 ). Ecosystem engineers: From pattern formation to habitat creation , Phys. Rev. Lett ., 93 , 098105 . OpenUrl CrossRef PubMed Gilad , E. , von Hardenberg , J. , Provenzale , A. , Shachak , M. and Meron , E. ( 2007 ). A mathematical model of plants as Ecosystem Engineers . Journal of Theoretical Biology , 244 ( 4 ), 680 – 691 . OpenUrl CrossRef PubMed Web of Science ↵ Hagos , M.G. , ( 2001 ). The Influence of Tree Thinning and Subhabitat Differentiation on the Reproductive Dynamics of Acacia mellifera subsp . detinens. M.Sc. Agric-thesis , University of the Free State , Bloemfontein, South Africa . ↵ Harris , C.R. , Millman , K.J. , van der Walt , S.J. et al. ( 2020 ). Array programming with NumPy . Nature 585 , 357 – 362 . OpenUrl CrossRef PubMed ↵ HilleRisLambers , R. , Rietkerk , M. , van den Bosch , F. , Prins , HHT. and de Kroon , H. ( 2001 ). Vegetation pattern formation in semi-arid grazing systems . Ecology 82 : 50 – 61 . OpenUrl CrossRef Web of Science ↵ Leenaars , J. and Hengl , T . ( 2015 ) Africa SoilGrids - Texture . ISRIC . ↵ Leenaars , J. and Hengl , T . ( 2017 ) Africa SoilGrids - Root zone plant available water holding capacity aggregated at ERZD . ISRIC . ↵ Klausmeier , C.A . ( 1999 ). Regular and irregular patterns in semiarid vegetation . Science 284 , 1826 – 1828 . OpenUrl Abstract / FREE Full Text ↵ Kéfi , S. , Eppinga , M.B. , de Ruiter , P.C. and Rietkerk , M. ( 2010 ). Bistability and regular spatial patterns in arid ecosystems . Theoretical Ecology , 3 : 257 – 269 . OpenUrl ↵ Kéfi , S. , Rietkerk , M. , Alados , C. et al. Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems . Nature 449 , 213 – 217 ( 2007 ). OpenUrl CrossRef PubMed Web of Science ↵ Lejeune , O. , Tlidi , M. and Couteron P . ( 2002 ). Localized vegetation patches: a self-organized response to resource scarcity. Phys . Rev. E 66 : 010901R . OpenUrl ↵ Lejeune , O. , Tlidi , M. and Lefever , R . ( 2004 ). Vegetation spots and stripes: dissipative structures in arid landscapes . Int. J. Quantum Chem . 98 : 261 – 271 . OpenUrl CrossRef Web of Science ↵ Manning , R ., ( 1891 ). On the flow of water in open channels and pipes . Transactions of the Institution of Civil Engineers of Ireland . 20 : 161 – 207 . OpenUrl ↵ Martínez-García , R. , J. M. Calabrese , E. Hernández-García , and C. López ( 2013 ), Vegetation pattern formation in semiarid systems with-out facilitative mechanisms , Geophys. Res. Lett ., 40 , 6143 – 6147 . OpenUrl ↵ Martinez-Garcia , R. , Cabal , C. , Calabrese , J. M. , Hernández-García , E. , Tarnita , C. E. , López , C. and Bonachela , J. A . ( 2023 ). Integrating theory and experiments to link local mechanisms and ecosystem-level consequences of vegetation patterns in drylands . Chaos, Solitons & Fractals , 166 , 112881 . OpenUrl Meron , E . ( 2012 ). Pattern-formation approach to modelling spatially extended ecosystems . Ecol Model , 234 , pp. 70 – 82 . OpenUrl CrossRef Web of Science ↵ Meron , E. , Gilad , E. , von Hardenberg , J. , Shachak , M. and Zarmi , Y. ( 2004 ). Vegetation patterns along a rainfall gradient . Chaos Solitons Fractals 19 : 367 – 376 . OpenUrl CrossRef Web of Science ↵ NAU , ( 2010 ). The Effect of Bush Encroachment on groundwater resources in Namibia, a Desk Top Study . Windhoek, Namibia : Namibia Agricultural Union . unpublished. ↵ Okayasu , T. and Aizawa , Y . ( 2001 ). Systematic analysis of periodic vegetation patterns . Prog Theor Phys 106 : 705 – 720 . OpenUrl CrossRef ↵ Rietkerk , M. , Boerlijst , M.C. , van Langevelde , F. , HilleRisLambers , R. , van de Koppel , J. , Kumar , L. , Prins , H.H.T. ( 2002 ). Self-organization of vegetation in arid ecosystems . Am Nat 160 : 524 – 530 . OpenUrl CrossRef PubMed Web of Science ↵ Rietkerk , M. , Dekker , S.C. , de Ruiter , P.C. and van de Koppel , J. ( 2004 ). Self-Organized Patchiness and Catastrophic Shift in Ecosystems . Science , 305 , 1926 – 1929 . OpenUrl Abstract / FREE Full Text ↵ Rietkerk , M. and van de Koppel , J. ( 2008 ). Regular pattern formation in real ecosystems . Trends Ecol. Evol . 23 , 169 – 175 . OpenUrl CrossRef PubMed Web of Science ↵ Scheffer , M. , Bascompte , J. , Brock , W.A. , Brovkin , V. , Carpenter , S.R. , Dakos , V. , … Sugihara , G. ( 2009 ). Early-warning signals for critical transitions . Nature , 461 , pp. 53 – 59 , doi: 10.1038/nature08227 OpenUrl CrossRef PubMed Web of Science ↵ Sherratt , J.A ., ( 2016 ). When does colonisation of a semi-arid hillslope generate vegetation patterns? J. Math. Biol . 73 , 199 – 226 . OpenUrl PubMed ↵ Sherratt , J. A. and Lord , G. J . ( 2007 ). Nonlinear dynamics and pattern bifurcations in a model for vegetation stripes in semi-arid environments . Theor. Popul. Biol . 71 : 1 – 11 . OpenUrl CrossRef PubMed Web of Science ↵ Siteur , K. , Siero , E. , Eppinga , M.B. , Rademacher , J.D. , Doelman , A. and Rietkerk , M . ( 2014 ). Beyond Turing: The response of patterned ecosystems to environmental change . Ecol. Complex ., 20 , 81 – 96 . OpenUrl CrossRef Web of Science ↵ Snyder , K.A. and Tartowski , S.L . ( 2006 ). Multi-scale temporal variation in water availability: Implications for vegetation dynamics in arid and semi-arid ecosystems . Journal of Arid Environments , 65 : 219 – 234 . OpenUrl ↵ Thornthwaite , C . ( 1948 ). An Approach toward a Rational Classification of Climate . Geogr. Rev ., 38 ( 1 ), 55 – 94 . OpenUrl CrossRef Web of Science ↵ Ursino , N . ( 2005 ). The influence of soil properties on the formation of unstable vegetation patterns on hillsides of semiarid catchments , Adv. Water Resour ., 28 , 956 – 963 . OpenUrl ↵ Ursino , N . ( 2007 ). Modeling banded vegetation patterns in semiarid regions: Interdependence between biomass growth rate and relevant hydrological processes , Water Resour. Res ., 43 , W04412 . OpenUrl ↵ Vanderveken , L. , Martínez Montero , M. , and Crucifix , M . ( 2023 ). Existence and influence of mixed states in a model of vegetation patterns . Nonlinear Processes in Geophysics , 30 ( 5 ), 585 – 599 . OpenUrl ↵ Veldhuis , M.P. , Martinez-Garcia , R. , Deblauwe , V. and Dakos , V . ( 2022 ). Remotely-sensed slowing down in spatially patterned dryland ecosystems. Ecography , e 06139 . ↵ Von Hardenberg , J. , Meron , E. , Shachak , M. , Zarmi , Y. ( 2001 ). Diversity of vegetation patterns and desertification . Phys. Rev. Lett . 87 , 3 – 6 . OpenUrl References ↵ Ahmed , K.R. , Akter , S . ( 2017 ). Analysis of landcover change in Southwest Bengal Delta due to floods by NDVI, NDWI and K-means cluster with Landsat multi-spectral surface reflectance satellite data . Remote Sens Appl Soc Environ 8 : 168 – 181 . OpenUrl ↵ Alam , S. A. , Starr , M. and Clark , B. J. F . ( 2013 ). Tree biomass and soil organic carbon densities across the Sudanese woodland savannah: A regional carbon sequestration study . Journal of Arid Environments , 89 : 67 – 76 . OpenUrl ↵ Arcement , G. J. and Schneider , V. R. ( 1989 ). Guide for selecting Manning’s roughness coefficients for natural channels and flood plains . U.S. G.P.O . ↵ Arthur , D. and Vassilvitskii , S . ( 2007 ). K-Means++: The advantages of careful seeding . Proc 18th Annu ACM-SIAM Symp Discrete Algorithms 1027 – 1035 . ↵ Bein , E. , Habte , B. , Jaber , A. , Birnie , A. , and Tengnas , B . ( 1996 ). Useful trees and shrubs in Eretria . Technical handbook No. 12, Regional Soil Conservation unit, Nairobi . ↵ Bonachela , J.A. , Pringle , R.M. , Sheffer , E. , Coverdale , T.C. , Guyton , J.A. , Caylor , K.K. , Levin , S.A. , Tarnita , C.E . ( 2015 ). Termite mounds can increase the robustness of dryland ecosystems to climatic change . Science 345 : 651 . OpenUrl ↵ Celik , T . ( 2009 ). Unsupervised change detection in satellite images using principal component analysis and K-means clustering . IEEE Geosci Remote Sens Lett 6 : 772 – 776 . OpenUrl ↵ Donaldson , C. H. ( 1969 ). Bush encroachment and methods of control with special reference to the Blackthorn (Acacia mellifera) problem of the Molopo area . Armoedsvlakte research station, Vryburg . ↵ ESRI . “Slope map” [digital map]. “Terrain” . Dec 30, 2013 . Source: Airbus, USGS, NGA, NASA, CGIAR, NLS, OS, NMA, Geodatastyrelsen, GSA, GSI and the GIS User Community . Available online: https://www.arcgis.com/home/item.html?id=af25a795273440deb449b336543602be . (accessed on March 23, 2023) ↵ ESRI . “World Topographic Map” [digital map]. “Terrain” . June 13, 2013 . Sources: Esri, HERE, Garmin, Intermap, INCREMENT P, GEBCO, USGS, FAO, NPS, NRCAN, GeoBase, IGN, Kadaster NL, Ordnance Survey, Esri Japan, METI, Esri China (Hong Kong), © OpenStreetMap contributors, GIS User Community . Available online: https://www.arcgis.com/home/item.html?id=30e5fe3149c34df1ba922e6f5bbf808f . (accessed on March 23, 2023) ↵ Hoshen , J. , Kopelman , R . ( 1976 ). Percolation and cluster distribution. I. Cluster multiple labeling technique and critical concentration algorithm . Phys Rev B 14 : 3438 – 3445 . OpenUrl ↵ Leenaars , J. and Hengl , T . ( 2015 ) Africa SoilGrids - Texture . ISRIC . ↵ Leenaars , J. and Hengl , T . ( 2017 ) Africa SoilGrids - Root zone plant available water holding capacity aggregated at ERZD . ISRIC . ↵ Joubert , D.F. , Rust , A. , Smit , G. , and Hoffman , M ., ( 2017 ). Growth rates and mortality patterns of Acacia mellifera subsp. detinens in the semi-arid Highland Savanna, Namibia: Encroachment is not as rapid as previously believed . Namibian Journal of Environment , 1 , A-5 . OpenUrl ↵ Joubert , D.F. , Smit , G.N. and Hoffman , M.T . ( 2013 ). The influence of rainfall, competition, predation and germination and the establishment of an encroaching Acacia in an arid Namibian savanna . Journal of Arid Environments , 91 : 7 – 13 . OpenUrl ↵ Manning , R ., ( 1891 ). On the flow of water in open channels and pipes . Transactions of the Institution of Civil Engineers of Ireland . 20 : 161 – 207 . OpenUrl ↵ NAU , ( 2010 ). The Effect of Bush Encroachment on groundwater resources in Namibia, a Desk Top Study . Windhoek, Namibia : Namibia Agricultural Union . unpublished. ↵ Olaode , A. , Naghdy , G. , Todd , C . ( 2014 ). Unsupervised classification of images: A review . Int J Image Process 8 : 2014 – 325 . OpenUrl ↵ Rousseeuw , P.J . ( 1987 ). Silhouettes: A graphical aid to the interpretation and validation of cluster analysis . J Comput Appl Math 20 : 53 – 65 . OpenUrl CrossRef PubMed Web of Science ↵ Rietkerk , M. , Boerlijst , M.C. , van Langevelde , F. , HilleRisLambers , R. , van de Koppel , J. , Kumar , L. , Prins , H.H.T. ( 2002 ). Self-organization of vegetation in arid ecosystems . Am Nat 160 : 524 – 530 . OpenUrl CrossRef PubMed Web of Science ↵ Seed , G. and Malik , A. ( 1968 ). Some aspects of measured and estimated evaporation in the Sudan, Durham theses, Durham University . Available at Durham E-Theses Online: http://etheses.dur.ac.uk/10121/ ↵ Siteur , K. , Siero , E. , Eppinga , M.B. , Rademacher , J.D. , Doelman , A. and Rietkerk , M . ( 2014 ). Beyond Turing: The response of patterned ecosystems to environmental change . Ecol. Complex ., 20 , 81 – 96 . OpenUrl CrossRef Web of Science ↵ Thornthwaite , C . ( 1948 ). An Approach toward a Rational Classification of Climate . Geogr. Rev ., 38 ( 1 ), 55 – 94 . OpenUrl CrossRef Web of Science ↵ Ursino N ( 2005 ) The influence of soil properties on the formation of unstable vegetation patterns on hillsides of semiarid catchments . Adv Water Resour 28 : 956 – 963 OpenUrl ↵ Wan , L. , Li , Y. , Cen , H. , Zhu , J. , Yin , W. , Wu , W. , Zhu , H. , Sun , D. , Zhou , W. , He , Y . ( 2018 ). Combining UAV-based vegetation indices and image classification to estimate flower number in oilseed rape . Remote Sens 10 : 1484 . OpenUrl View the discussion thread. Back to top Previous Next Posted October 04, 2025. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Vegetation clump size and number as indicators for alternative stable states in semi-arid ecosystems Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Vegetation clump size and number as indicators for alternative stable states in semi-arid ecosystems Z. Chen , K. Halder , J. A. Bonachela bioRxiv 2025.10.03.680294; doi: https://doi.org/10.1101/2025.10.03.680294 Share This Article: Copy Citation Tools Vegetation clump size and number as indicators for alternative stable states in semi-arid ecosystems Z. Chen , K. Halder , J. A. Bonachela bioRxiv 2025.10.03.680294; doi: https://doi.org/10.1101/2025.10.03.680294 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Ecology Subject Areas All Articles Animal Behavior and Cognition (7635) Biochemistry (17697) Bioengineering (13894) Bioinformatics (41951) Biophysics (21455) Cancer Biology (18592) Cell Biology (25507) Clinical Trials (138) Developmental Biology (13380) Ecology (19903) Epidemiology (2067) Evolutionary Biology (24321) Genetics (15610) Genomics (22509) Immunology (17737) Microbiology (40398) Molecular Biology (17182) Neuroscience (88618) Paleontology (667) Pathology (2833) Pharmacology and Toxicology (4825) Physiology (7641) Plant Biology (15158) Scientific Communication and Education (2046) Synthetic Biology (4296) Systems Biology (9825) Zoology (2271)