Topological defects law for migrating banded vegetation patterns in arid climates

Self-organization and pattern formation are ubiquitous processes in nature. We study the properties of migrating banded vegetation patterns in arid landscapes, usually presenting dislocation topological defects. Vegetation patterns with dislocations are investigated in three different ecosystems. We show through remote sensing data analysis and theoretical modeling that the number of dislocations N(x) decreases in space according to the law N ∼ log(x/B)/x, where x is the coordinate in the opposite direction to the water flow and B is a suitable constant. A sloped topography explains the origin of banded vegetation patterns with permanent dislocations. Theoretically, we considered well-established approaches to describe vegetation patterns. All the models support the law. This contrasts with the common belief that the dynamics of dislocations are transient. In addition, regimes with a constant distribution of defects in space are predicted. We analyze the different regimes depending on the aridity level and water flow speed. The reported decay law of defects can warn of imminent ecosystem collapse.


INTRODUCTION
Self-organization phenomena leading to spatially periodic patterns are observed in complex or nonlinear systems (1)(2)(3)(4)(5).Vegetation population dynamics provide puzzling and notable examples of spatially periodic structures, generically called vegetation patterns, formed by large-scale botanical organizations controlled by a nonequilibrium symmetry-breaking instability (6)(7)(8)(9)(10)(11)(12)(13).The banded patterns, often called tiger bush (14), consist of dense vegetation bands alternating with sparsely covered or even bare soil, their wavelength ranges from decimeters to hundreds of meters.Banded vegetation patterns have probably been first reported by Macfadyen in the earlier fiftieth (15,16).The spontaneous symmetry-breaking instability causes their formation even when the topography is flat (6).The presence of the slope causes the migrating banded patterns (6,8,12).They grow by a few decimeters each year in the opposite direction of the water flow (8,12).Besides, a bibliography of empirical and scientific studies devoted to the origin of their formation and maintenance can be found in (6,12,(17)(18)(19)(20).
Most of the banded vegetation patterns observed in nature are disordered and present topological defects such as dislocations, as can be seen in Fig. 1.Dislocations in the banded vegetation patterns are indicated by red rings in the aerial photographs of Fig. 1.When two stripes join and transform into a single one, they form a defect called dislocation.Observations across large areas of numerous arid and semi-arid regions of Africa, Australia, America, and the Middle East show that topological defects are abundant.Banded vegetation is a well-documented issue that has been abundantly discussed and is by now fairly well understood.So far, however, the law governing the formation of such defects has neither been experimentally determined nor theoretically predicted.
Here, we establish a law governing the organization of dislocations.By analyzing satellite images taken from vast territories of the African and American continent, we show that the number of dislocations obeys the formula N ∼ log(x/B)/x, where x is the coordinate in the opposite direction of water flow and B is a suitable constant.Theoretically, we have considered three different ecological approaches describing the dynamics of topological defects.All these models quantitatively support this deterministic law.Furthermore, these ecological models predict an additional dynamical regime where the number of dislocations remains constant.In addition to the slope, which is the source of dislocation propagation, we show that boundary conditions play an essential role in their permanent creation; defect generation from boundaries is a documented phenomenon in nonlinear physics that appears in several situations, the most common being the dynamics of viscous flows (21).Therefore, with a source of dislocations through the boundaries, the dynamics of these topological defects can be permanent rather than transient.This fact strongly contrasts with previous work where dislocation formation is considered a transient dynamic due to their mutual annihilation interaction, leading at long times to a perfectly ordered banded pattern free of defects (8).The permanent dynamics of defects is the process of pairs of dislocations being created at the boundary with opposite topological charges, and then they move with the pattern migration velocity (toward x) at the same time they interact, approaching each other until annihilation; the process is repeated in time in an unpredictable way.This complex permanent dynamic leaves an imprint in the dislocation number as a function of the x direction.We demonstrate how a decaying number of dislocations in space may be used as an early indicator of an ecosystem's potential collapse under harsh environmental conditions.We conclude by showing how the measure of the dislocation distribution in space can be used as a noninvasive tool for diagnosing ecosystem health.The ecosystem transition to bare soil is a much-studied issue in which spatial vegetation models play a crucial role (22)(23)(24).Our theory complements the understanding of ecosystem adaptability and resilience until now, as we consider the role of sloped topography and boundary conditions in the dynamics.The predicted law is supported by field observations and can be crucial for identifying and comprehending the different spatiotemporal behaviors seen in complex systems other than ecological ones.

Remote sensing data analysis and the dislocation distribution decay law
To establish through field observation that the number of dislocations in the banded vegetation follows a logarithmic law, we perform an image analysis.Three regions of the world are considered: Chile, Sudan, and the United States.To do that, we use highresolution satellite images obtained from the Google Earth software (https://earth.google.com/web),together with the elevation database SRTM (Shuttle Radar Topography Mission) with 1-arc sec resolution (25).First, we select and create an adequate mask of the region where banded vegetation patterns are settled on sloppy landscape as shown in the satellite images of Fig. 2 (A to C).Second, we extract the mean orientation of the elevation gradient 〈θ〉 over the selected region as illustrated by Fig. 2 (D to F).We assume that the mean orientation of the elevation gradient is parallel to x.In the case of the banded vegetation pattern in hyper-arid landscapes of Chile, the x variable decreases with height, as water comes from the Eastto-West traveling fog (26,27).This means that the water bubbles move uphill, and therefore, the vegetation pattern migrates downhill.However, in arid landscapes of North America and Sudan, water is supplied by rainfall, and the x variable grows with height.
Once the x direction is defined, dislocation positions are marked.For the satellite images, because of the intrinsic fluctuations, the high anharmonicity, and the high variations in the wavelength in the banded vegetation, the dislocations could not be recognized with standard methods.To detect dislocations, we construct a skeleton of the banded vegetation pattern using the software for scientific image analysis Fiji (28) (see Materials and Methods section).This method allows us to identify the branch split points and the branch ends as points representing dislocations of the local pattern.The results are summarized in Fig. 2 (G to I).
Last, we select an area within the banded vegetation pattern in the plane (x, y), and we define the dislocation number N(x, y) over tiles of one wavelength side.Then, we average along the y direction.The obtained dislocation number N(x) is plotted as a function of x/λ where λ is the wavelength of the banded vegetation pattern.Note that N is the expected number of dislocations in a λ 2 surface tile centered on the (x, y) plane.The results are shown in Fig. 3.In the hyper-arid landscape of Chile and Sudan and the United States arid landscapes, the number of dislocations N(x) decreases with the x direction.From these results obtained from remote sensing observations, we can see that the spatial distribution  of defects is not uniform.Their distribution depends on the sloped direction along which the water flows.The fit of the observations is indicated by continuous orange curves in Fig. 3. Unexpectedly, the N(x) ∼ log(x/B)/x decay law fits well with the data obtained from Chile and Sudan.For the U.S. landscapes, the fit is excellent.
To understand the complex ecological phenomenon reported above and the role of the law dictating the number density of dislocations in space, mathematical modeling is indispensable.In the following subsections, we investigate the origin of the logarithm decay law through theoretical investigation and numerical simulations of three ecological models.

Theoretical modeling
To shed light on the observations of the previous section, we consider different standard approaches to explain biomass evolution.The dynamics of ecological systems are often described by either reaction-diffusion models that explicitly incorporate water transport or integrodifferential equations.The latter approach is grounded on nonlocal interactions associated with facilitative and competitive feedback and seed dispersion.Other models based on cellular automata have been first proposed ( 14) and also models based on environmental randomness (3,29).
where r = (x, y) and t are the spatial coordinates and time, respectively.m f and m c account for facilitation and competition plant-to-plant feedbacks.The nonlocal contributions read 2 fy;cy � are ellipsoidal coupling kernels with a shift in x with respect to the origin of magnitude x 0f,0c .The facilitative and the competitive ranges are l fx,cx and l fy,cy for the x and y direction, and the feedback strengths are measured by χ f,c .The Kernels ϕ f,c introduce an anisotropy and break the reflection symmetry x ↔ −x.The last term of the right-hand side of Eq. 2 models the seed dispersion with diffusive coefficient d.
In the weak gradient approximation, one can derive from model Eq. 1 a simpler partial differential equation (see Materials and Methods for details) of the form where α accounts for the translation parameter of the ellipsoidal kernel.The parameter η measures the decrease-to-growth rate ratio, called the aridity parameter.κ is the facilitation-to-competition strength difference, called the cooperativity parameter.γ is proportional to the difference of the squared competition-tofacilitation lengths and p plays the same role as d.
In addition to the integrodifferential and the weak-gradient models, we consider the water-biomass model describing the space-time evolution of the biomass (b) and water (w) density.This model reads (10) The slope effect is accounted for in the term α∂ x (w − vb), where α is the water speed, which flows opposite (in favor) to the x direction for α < 0 (α > 0).Because of the water absorption by plants, the biomass reduces the water advective transport mediated by the parameter v.The parameters γ and σ model the biomass production increase with water considering a saturable function, d models the seed dispersion, and μ accounts for mortality.The parameter p measures water input, ρ reduces the transpiration rate linearly with the biomass, and β models how plants affect water absorption by the soil.
Numerical simulations of the nonlocal model Eq. 1 with ellipsoidal translated kernels (l fx,cx ≠ l fy,cy ) display propagative banded patterns for small x 0f,0c values as shown in Fig. 4A.These results are obtained using Dirichlet boundary conditions with zero value in the flow direction edges (b = 0 for x = 0 and x = s, where s is the system size).Periodic boundary conditions are used in the y direction.Numerical simulations of all the models presented were conducted with a Runge-Kutta algorithm of fourth order for time integration and a finite difference scheme for space discretization.
As the translation parameter increases, the uniform banded patterns become unstable and the system generates permanent dislocations from the fixed edge x = 0, see Fig. 4B.Similarly, a permanent emission of defects can be sustained by environmental stochastic fluctuations (31).
The permanent dislocation dynamics are also obtained from the reduced model Eq. 2 (cf.Fig. 4, C and D) and the reaction-diffusion Eq. 3 (cf.Fig. 4, E and F).All models display a transition from a perfect traveling banded vegetation to a regime where dislocations are permanently emitted, as shown in Fig. 4.This transition occurs for α* < α, the system asymptotically tends to a regular banded pattern as x → ∞, cf.Fig. 5A(i), but with dislocations being created in the upstream boundary.The critical value α* is the threshold for the boundary layer instability, and below this value, the number of dislocations is zero.The α* parameter has no analytical expression and depends on the model considered.Hence, this parameter only is determined numerically.When dislocations are only created on the edge, numerical data follows where A, B, and x 0 are the fit parameters.The number of dislocations N(x) as a function of xis plotted in Fig. 5A(iii).This numerical result agrees with field observations using remote sensing image data analysis, as shown in the panels of Fig. 3. Table 1 summarizes the results of fitting law (4) to both observational and numerical data.
To understand analytically the origin of the logarithmic law, we perform a normal form analysis, which leads to the derivation of the well-known Ginzburg-Landau Eq. 9 (see Materials and Methods).Dislocations correspond to topological singularities in the phase of the Ginzburg-Landau equation (32)(33)(34).This reduction shows that in defect interaction, when the nonlinear phase correction β is small, the length l between the defects decays according to the law l 2 = t/log(t) (35)(36)(37).Then if the system is advected with speed α, one can interchange the role of time for space using the relation t = x/α.Hence, this characteristic length changes with distance as l 2 = x/αlog(x/α).Likewise, the average number of defects in a given area Π is N(x) = Π/l 2 = Παlog(x/α)/x.Again, the normal form analysis confirms the logarithmic law.
When, however, the advection parameter increases, i.e., large α, we identify a second transition where a permanent creation of dislocations occurs not only from the edge but also in the bulk, as shown in Fig. 5A(iv).This regime is well-known in nonlinear systems in general, and it is referred to as defects turbulence (32,38).In this regime, the averaged dislocation number is constant N(x) = c as a result of the continuous creation of defects in the bulk, see Fig. 5A(vi).This figure is obtained from numerical simulations of Eq. 2. The transition from nonturbulent to turbulent regime is also obtained from the Ginzburg-Landau Eq. 9, as shown in Fig. 5B.
The numerical analysis of ecological models indicates that by only measuring the number of dislocations in the pattern, one can infer if the semi-arid and arid ecosystems operate in the turbulent regime where N(x) = c or in the nonturbulent regime where N(x) obeys a logarithmic decay law.This law obtained from numerical simulations of the three models considered here is in good agreement with observations using remote sensing image analysis, as shown in the panels of Fig. 3. Therefore, the measure of the number of dislocations in the vegetation patterns and their spatial distribution can be used as a noninvasive tool for diagnosing the degree of complexity of arid landscapes and for identifying unexpected dynamical phenomena in ecological systems.

DISCUSSION
The transition between different regimes is investigated in terms of the speed of the water flows.We have shown that for a large speed, the ecosystem presents a turbulent behavior where the number of topological defects is constant.For a small value of the water flow speed, the number of defects decreases according to the logarithmic law.In what follows, we discuss the effect of the aridity level on dislocation formation.For this purpose, we fix the speed of the water flow and vary the aridity level.Figure 6 summarizes the different ecosystem operating regimes.For small aridity parameters, the system develops migrating banded pattern devoid of defects (cf.red curve).For a moderate level of aridity, the system exhibits a transition toward a turbulent regime where the number of dislocations is constant (see blue curve).However, for extreme aridity conditions, the system reaches a regime where the system undergoes self-organized dislocations with a logarithmic decay law (see yellow curve).Further increasing the aridity, the banded patterns exhibit a transition toward a state totally devoid of vegetation.Thus, for a given landscape with a homogeneous slope, the presence of a decaying number of dislocations can be an ecological indicator of imminent transition toward a bare state.This complements what is known about the catastrophic shift of ecosystems in flat topography, where different types of stationary patterns exist and where multistability of patterns with different wavelengths can be observed.The existence of many pattern branches permits the ecosystem to adapt to environmental changes, which allows a patterned ecosystem to survive past the tipping point compared to a homogeneous ecosystem (22).However, once a gentle slope is introduced into the system, the advective effects of water flow must be taken into account for the stability analysis of patterns and other solutions of the system (39).This changes the stable pattern branches compared to a flat territory case.In addition, complex and turbulent-like dynamics can emerge as a consequence of the slope.These complex dynamical regimes have their own relative stability compared to the different perfect migrating pattern regimes and homogeneous states.Multistability of the complex dynamical regimes and the perfect patterns can occur, as observed in Fig. 6, suggesting that in adaptation to change, ecosystems could transit to these complex regimes if in the presence of a slopped territory.Numerically, only the most stable branch of perfect patterns is accessible.Different stable branches of perfect patterns could be obtained from Eq. 9 analytically at the onset of pattern formation.
To summarize, we have investigated different transitions of migrating vegetation banded patterns: from zero defects, to constant, and to a decaying number of dislocations.We have shown analytically that the number of dislocations in space N(x) obeys a N(x) ∼ log(x/B)/x law.This formula is in good agreement with numerical simulations of the three ecological models and with remote sensing image data taken from three arid ecosystems of different continents.Furthermore, the dislocation law allows us to determine whether the self-organized response to the water scarcity of arid and semi-arid ecosystems favors uniform bands or ecological spatio-temporal complexity.
A usual approach to characterize the response of plants to changes in their environment is through studies of the plants themselves (local analysis).Characterizing dislocation distributions of migrating banded vegetation patterns (macroscopic analysis) opens a noninvasive diagnostic tool for determining the degree of aridity mediated by desertification and global warming processes.Likewise, a full characterization of the bifurcation diagram for models including reflection symmetry rupture and dislocations self-organization becomes relevant in designing conservation guidelines, preventing the further degradation of migrating patterned vegetation cover.
Last, the spatial distribution of defects is a consequence of their creation induced by the boundary condition and their annihilation through mutual interaction.The boundary induces inhomogeneities in the system.The inhomogeneous Ginzburg-Landau equation constitutes an ideal framework for investigating the dynamics of defects in banded vegetation patterns.It provides a unified and simple description containing the dynamics discussed.Thus, the analytical results can be easily extended to describe similar laws in other complex nonlinear spatially extended systems present in nature.

Detailed derivation of the weak gradient model with advection
We look for an approximation to Eq. 1 of the main text, in the form of a partial differential Eq. 2. To account for anisotropy, we consider that the interaction ranges associated with facilitation and the competition l cx,fx and l cy,fy are different.We seek corrections to the steady states close to μ = 1 and b = 0 that depend on time and space through the slow variables T = ϵt, X = ϵ 1/8 x, and Y = ϵ 1/8 y.We expand the parameters μ, χ f,c , l fx,fy , x 0c,0f , d, and the biomass b in terms of a small parameter ϵ (ϵ ≪ 1) that measures the distance from μ = 1 as follows Introducing these scalings and the above expansions in Eq. 1, we then obtain a sequence of linear problems for unknown functions.We analyze each problem and apply the solvability condition at each order.These conditions are automatically fulfilled at the orders ϵ 1/2 and ϵ.By applying the solvability condition at the higher order inhomogeneous problem (ϵ 3/2 ), we obtain the following partial differential equation for the biomass where the coefficients are χ cy χ 0 , and Λ ¼ 3l 2 fx l 2 cx .Model Eq. 5 has different homogeneous steady states which are b = 0 and b 0 = κ ± (κ 2 −2η) 1/2 .Note that the upper branch of b 0 is stable when l cx > l fx and γ y > 0. Otherwise, we need to consider higher ϵ orders in the equation.The condition γ y = 0 will be used throughout the work, as it does not change the qualitative behavior of the system.

Detailed derivation of the Ginzburg-Landau equation Derivation of amplitude equation in the bulk
The amplitude equation obtained using a normal form analysis constitutes an adequate tool for understanding pattern formation.For the boundary conditions considered, the system creates a boundary layer next to the upstream edge of the system.The effect of this boundary layer can be neglected when focusing on regions far from the edges.Let us consider first the linear problem for a perturbation of the homogeneous stable state u ≪ 1 as b = b 0 + u, where b 0 ¼ κ þ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi κ 2 À 2η p is the homogeneous cover.Introducing this ansatz in Eq. 2 of the main text for the field b yields the linear problem Linear stability analysis for finite wavenumber k perturbations leads to the growth rate of modes λ(k) . To obtain the amplitude equation for the critical mode, let us move slightly from the instability condition, using as the bifurcation parameter η, as η = η c + ϵ.Introducing the following expansion where A ≡ A(X, Y, T ) is the slowly varying envelope, with the scalings X = ϵ 1/2 x, Y = ϵ 1/2 y, T = ϵt, and the parameter Ω c ≡ −b c αk c .A [n]  accounts for the terms of order n in the amplitude A. At order (ϵ 1/2 ), we obtain At this order, the solvability condition is automatically satisfied.For the sake of simplicity, let us define d x , and d(k) ≡ iαk − γk 2 + k 4 .Then, performing expansions, up to order ϵ, limiting to the case of small group velocities v g = ∂ k Ω| k c ∼ O(ϵ 1/2 ), we get To solve the linear problem, the following inner product is introduced valid over the periodic functions in space and time of period 2π/k c and 2π/Ω c .The kernel of the operator (∂ t − L c ) † , defined as the solution of (∂ t − L c ) † ψ = 0, corresponds to ψ = e ±i(k c x+Ω c t) .Then, applying the solvability condition, we find where Last, at order ϵ 3/2 , the solvability condition yields

Boundary layer effect
To figure out the emission of dislocations from the boundary of the system in the regime of decaying number of dislocations, we need to consider the boundary layer effect arising from the Dirichlet boundary conditions.We use the method suggested in (40).For this purpose, we suppose that sufficiently near to the upstream edge, one can write b ¼ b 0 þ ϵMðXÞ where M(X ) is a function that helps to connect the population state b = b 0 with the nonpopulation state b = 0 at the boundary and satisfies M(X ) → 0 when X → ∞.The analytical solution close to the boundary is not known.Qualitatively, b(X) is a Monod function.On the basis of this nonuniform b, a modified amplitude equation is derived.Making straightforward calculations, one finds a similar amplitude equation compared to Eq. 7 but with inhomogeneous linear terms where the parameters depend on M(X ) as , we get Eq.9.

Fig. 2 .
Fig. 2. Remote sensing analysis: Determination of the x direction and defect recognition with remote sensing data.(A to C) show the vegetation patterns in Chile 20 β 29.5 0 S, 70°3.5 0 W, Sudan 11°9 0 N, 28°16.5 0 E, and the United States 31°2.5 0 N, 103°5.5 0 W, respectively.(D to F) exhibit the direction of the steepest variation in the altitude over the region of interest.(G to I) illustrate the pattern's skeletons, and insets show the patterns dislocations highlighted.

Fig. 3 .
Fig. 3. Dislocation number decay law obtained from remote sensing analysis in Chile, Sudan, and U.S. landscapes.Circles account for observed data, and the orange curves represent the fits.(A to C) correspond to a N(x) ∼ log(x/B)/x fit for patterns in Chile, Sudan, and the United States, respectively.Fit parameters in λ units are (A) x 0 = −2.9,B = 1.2, and A = 2.7; (B) x 0 = −2, B = 1.5, and A = 10.6;(C) x 0 = −2.4,B = 1.3, and A = 5.R 2 is the coefficient of determination of the fits.

Fig. 6 .
Fig. 6.Diagram of migrating banded vegetation pattern biomass as a function of the aridity parameter η.For Eq. 2, parameters κ = 0.3, p = 0.05, γ = 2, and α = 1.(A) illustrates the mean biomass 〈b〉 at the steady states of the model when changing the aridity.(B) corresponds to the branch of perfect patterns N(x) = 0. (C) shows the turbulent-like behavior where N(x) = c.(D) represents the branch of asymptotic patterns, where N(x) ∼ log(x/B)/x.

μ 1
ðXÞ ¼ ð2κ À 3b c þ γk 2 c À k 4 c ÞMðXÞ νðXÞ ¼ αk c MðXÞBoth terms are proportional to the slow inhomogeneity M(X ), so they asymptotically go to 0 as X → ∞.Hence, one recovers the homogeneous Ginzburg-Landau Eq. 7.With the change of parameters and variables as X ¼ ffi ffi ffi ffi ffi ffiD x p X, Y ¼ ffi ffi ffi ffi ffiffi D y p Y, A ¼ 1= ffi ffi ffi ffi ffi ffiffi À a p A, β = β/a,and α ¼ b c α= ffi ffi ffi ffi ffi ffi D x p

Table 1 . Summary of the best fit for the decaying spatial distribution of dislocations for mathematical models and remote sensing image analysis
. A, B, and x 0 are the fit parameters of Eq. 4, and R 2 is the coefficient of determination for the respective fits.