High-integrity human intervention in ecosystems: Tracking self-organization modes

Humans play major roles in shaping and transforming the ecology of Earth. Unlike natural drivers of ecosystem change, which are erratic and unpredictable, human intervention in ecosystems generally involves planning and management, but often results in detrimental outcomes. Using model studies and aerial-image analysis, we argue that the design of a successful human intervention form calls for the identification of the self-organization modes that drive ecosystem change, and for studying their dynamics. We demonstrate this approach with two examples: grazing management in drought-prone ecosystems, and rehabilitation of degraded vegetation by water harvesting. We show that grazing can increase the resilience to droughts, rather than imposing an additional stress, if managed in a spatially non-uniform manner, and that fragmental restoration along contour bunds is more resilient than the common practice of continuous restoration in vegetation stripes. We conclude by discussing the need for additional studies of self-organization modes and their dynamics.


Illustration of concepts and additional figures
Instabilities for vegetation rehabilitation model  (Fig Aa) destabilizes the bare-soil state as the precipitation rate P is increased beyond a threshold value (see Fig 3a). As the spectral density shown in the rightmost panel indicates, the SO mode associated with this instability is a stripe Fourier mode of the form A(t) cos(kx + φ) with a wavenumber k = k f even though the natural wavenumber that the system tends to form, k 0 , is generally different from k f . The growth of this mode directs the system towards a resonant pattern of parallel stripes consisting of a vegetation band at each modulation stripe.
The second instability (Fig Ab) destabilizes the resonant stripe pattern as P is decreased below another threshold value (see Fig 3a). The SO modes associated with this instability are a symmetric pair of oblique Fourier modes, a(t) cos(k x x±k y y +φ ± ), where k x = k f /2 and k y is such that the total wavenumber k is equal to the natural wavenumber k 0 , i.e. satisfying k 2 = k 2 x + k 2 y = k 2 0 , as the rightmost panel in Fig Ab  shows [1]. The growth of these modes directs the system towards a resonant rhombic pattern (frame t = 125 in Fig Ab).  (2) and (3) in SI beyond the instability of bare soil to a stripe pattern (a) and the instability of a stripe pattern to a rhombic pattern (b). The rightmost panels show spectral densities of the asymptotic stripe and rhombic patterns in the Fourier plane (k x , k y ). They reveal the SO modes (dark dots) that grow beyond the instabilities: (a) a resonant stripe mode (k x = k f , k y = 0), (b) resonant oblique modes (k x = k f /2, ±k y ), lying on a circle of radius k 0 , the natural wavenumber in the absence of ground modulations. Note that each mode (k x , k y ) is accompanied by a conjugate mode (−k x , −k y ), necessary to guarantee the real-valuedness of the state variables.

Glossary
Most of the terms appearing below are mathematical concepts. In describing them we favored intuitive clarity over mathematical rigor, given the broad and mostly nonmathematically oriented readership.
Basin of attraction. The set of points in phase space, representing initial conditions that converge in time to the same stable state. Bifurcation diagram. A diagram that shows the existence ranges and stability properties of possible systems states. The vertical axis represents a state variable (e.g. biomass) while the horizontal axis represents a control parameter. Solid (dashed) lines represent stable (unstable) states.
Growing eigenmode. The direction in phase space along which a system changes following an instability. More technically, a growing eigenmode is a vector that describes the relationships between state variables that result in the growth of small perturbations as an instability threshold is traversed. Fourier modes. Functions that describe sinusoidal periodicities in time or space, and are characterized by their frequency or wavenumber values, respectively. Instabilities of stationary uniform states. Several instability types can be distinguished according to the eigenmode (self-organization mode) that grows beyond the instability threshold: • Uniform stationary instability. An instability that is driven by the monotonic growth of a uniform eigenmode. Such an instability generally results in a different stationary uniform state.
• Uniform oscillatory instability. An instability that is driven by the oscillatory growth of a uniform eigenmode. Such an instability leads to uniform timeperiodic oscillations (also called 'Hopf instability').
• Nonuniform stationary instability. An instability that is driven by the monotonic growth of a spatially periodic eigenmode. Such an instability leads to a stationary periodic pattern (also called 'Turing instability').
• Nonuniform oscillatory instability. An instability that is driven by the oscillatory growth of a spatially periodic eigenmode. Such an instability leads to a travelingwave pattern. (also called 'Turing-Hopf instability').
Instability and bifurcation. A threshold phenomenon in which a combination of feedback processes act in concert to amplify small perturbations about a system state along a particular eigenmode and lead to a state change. A state that lost its stability may still exist as an unstable state. Instabilities are often referred to as bifurcations since the qualitative changes they induce often involve the emergence of new additional states. Linear stability analysis. A method of pattern-formation theory by which instabilities of system states can be identified by analyzing the dynamics of infinitesimally small perturbations. The analysis provides information about the instability threshold and the feedback processes that drive the instability (through the eigenmode).
Numerical continuation. A method of computing approximate solutions of differential equations by iterative methods and following them as a control parameter is varied. The method is particularly useful for calculating bifurcation diagrams that include unstable solution branches. Pattern formation theory. The mathematical theory of spatially extended nonlinear dynamical systems, which accounts, among other things, for the spontaneous emergence of spatially non-uniform states in homogeneous systems through symmetrybreaking instabilities. The reader is referred to Ref. [2] for an introduction to patternformation theory with applications to ecology. Phase space (also state space). The space spanned by the state variables of a dynamical system ('phase plane' in the case of two state variables).
Phase-space trajectory. The trajectory emanating from a point in a phase space that describes how the state of a dynamical system changes with time. Spatial resonance. Self-adjustment of spatial periodicity (wavelength) of a patterned state to an imposed external periodicity. When the external periodicity, λ f , is close to the system's inherent periodicity, λ 0 , the adjustment is 1:1, i.e. the actual system periodicity, λ, is equal to λ f . When λ f ≈ λ 0 /2, the adjustment is 2:1, i.e. λ = 2λ f . Resonant pattern. A periodic pattern with a wavelength that matches exactly the wavelength of an externally imposed periodicity. Saddle point. An unstable system state that has both stable and unstable manifolds, that is, directions of convergence and departure in phase space. Saddle-node (fold) bifurcation. A collision and disappearance of two states of a dynamical system, as a control parameter traverses a threshold value. Spectral density. The distribution of power among the Fourier modes of an oscillating state or of a spatial pattern. The power of a Fourier mode is measured by the square of the absolute value of its amplitude. Periodic patterns are characterized by spectral densities with distinct peaks at wave-vectors that describe the pattern's periodicity.
Stable manifold of a system state. The set of points (initial conditions) in phase space that converge at long times to the system state. Stable state. A system state that recovers from any sufficiently small perturbation or disturbance. The state can be stationary or oscillatory, spatially uniform or periodic. Stationary state. A state of a dynamical system for which all state variables are independent of time. A stationary state can be spatially uniform, spatially periodic or disordered. The latter two are also called periodic or disordered stationary patterns. Unstable manifold of a system state. The set of points (initial conditions) in phase space that converge to the system state when time is run backwards. The eigenmode that grows at an instability point provides an approximation to the unstable manifold in the vicinity of the system state. Unstable state. A system state that evolves towards a different state when subjected to small perturbations of a certain form, e.g. perturbations that contain a particular wavenumber.
Wavenumber. Spatial frequency of a periodic pattern. The wavenumber k of a periodic pattern in a one-dimensional (1d) system is related to the pattern's periodicity or wavelength λ through the relation k = 2π/λ. In 2d systems the wavenumber is the magnitude of the wavevector k = (k x , k y ), i.e. k = k 2 x + k 2 y , where k x = 2π/λ x , k y = 2π/λ y and λ x , λ y represent the periodicities of the pattern in the x, y directions. 3 Analysis of empirical afforestation project

Precipitation data
In order to demonstrate the precipitation regime for the afforestation project in the Northern Nevgev, Israel, we collected precipitation information from Israel's meteorological service (https://ims.data.gov.il). We chose data from the nearby stations in the city of Be'er Sheva (a few kilometers south of the afforestation area, see next subsection) in the time-frame of 1980-2019. Three stations had usable information, in the following years: A) Be'er-Sheva, South (station #251900) during 1980-1995; B) Be'er-Sheva, Negev Institute (station #251690) during 1980-2009; C) Be'er-Sheva (station #251691) during 2004-2019. For each station, we collected monthly precipitation data, and summed over the winter seasons, where we included the first 7 months in a year with the last 5 months of the previous year in a single year (precipitation is negligible between June and August). For instance, the precipitation between, and including, August 2018 to July 2019 was summed and noted as the precipitation for the year 2019. The annual precipitation of each station and the averaged is given in table 3.1. The averaged precipitation data is also visualized in Fig D,

Analysis of remote sensing data
We give here further technical details on the analysis and categorization of the aerial images. We also show in Figs F1-F12 results of the analysis for the 100 regions studied. All regions are located within the area shown in Fig E. The process of choosing the regions, as well as their normalization (i.e. so as each region shows 4 vegetation stripes along the same direction) is described in the main text (Materials and Methods). The data of these regions is given in: https://doi.org/10.5281/zenodo.3902397.
The spectral analysis of the images (initially as RGB true color) was preformed as follows. Each image was converted into a intensity map (grayscale), to facilitate the FFT analysis. This was done using several image processing functions of Matlab. The red channel (of the RGB image) was taken, and the function adapthisteq was used to normalized intensities. The functions imerode and then imreconstruct were used (both with a disk kernel of 5x5), to smooth out the result. This image was then converted to spectral density map using FFT, with the central pixel set to 0 (i.e. ignoring average intensities).
In this way, we now have six representations of each region: a color (RGB) image, a grayscale image, and a spectral density map (i.e. after applying FFT), each for the year 2010 and the year 2019. These six representations are shown, for each of the 100 regions, in Figs F1-F12. For each region, based on the six representation, a manual analysis of the transition between 2010 and 2019 was preformed. The analysis categorized the transition that occurred in each region into one of 5 categories (shown in different colors in Fig E): • Category I: 16 regions show a clear oblique modes (hexagonal pattern in the FFT plane), implying a transition to a rhombic vegetation pattern (Figs F1-F2).
• Category II: 48 regions show a substantial change, with some marks of oblique modes (Figs F3-F8).
• Category III: 9 regions only show a minimal change in vegetation cover (Fig F9).
• Category IV: 9 regions appear to show a collapse of the vegetation, with more than half the trees disappearing (Fig F10).
• Category V: 18 regions show a change that is not well detected by the FFT analysis (Fig F11), or are difficult to categorize (Fig F12).
On the top of each region analysis (in the figures below), basic information is given. A score is noted for each region, where a higher score indicates that the 2010 region had a more regular stripe pattern (calculation of region score is defined in the Materials and Methods section of the main text). Besides the score, a region number is given (numbered from highest to lowest score), x&y coordinates (in pixel terms of the image shown in Fig E), its size L (length of the square region), and the angle θ of the region, relative to the image in