A Fast-Slow Model of Banded Vegetation Pattern Formation in Drylands
Punit Gandhi, Sara Bonetti, Sarah Iams, Amilcare Porporato, Mary Silber
AA Fast-Slow Model of Banded Vegetation Pattern Formation in Drylands
Punit Gandhi a , Sara Bonetti b,c , Sarah Iams d , Amilcare Porporato f , Mary Silber a Department of Mathematics and Applied Mathematics, Virginia Commonwealth University, Richmond, VA 23284, USA b Department of Environmental Systems Science, ETH Z¨urich, 8092 Z¨urich, Switzerland c Bartlett School of Environment, Energy and Resources, University College London,WC1H 0NN London, UK d John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA e Department of Civil and Environmental Engineering and Princeton Environmental Studies, Princeton University,Princeton, NJ 08540, USA f Department of Statistics and Committee on Computational and Applied Mathematics, University of Chicago, Chicago, IL60637, USA
Abstract
From infiltration of water into the soil during rainstorms to seasonal plant growth and death, the eco-hydrological processes that are thought to be relevant to the formation of banded vegetation patterns indrylands occur across multiple timescales. We propose a new fast-slow switching model in order to cap-ture these processes on appropriate timescales within a conceptual modeling framework based on reaction-advection-diffusion equations. The fast system captures hydrological processes that occur on minute tohour timescales during and shortly after major rainstorms, assuming a fixed vegetation distribution. Theseinclude key feedbacks between vegetation biomass and downhill surface water transport, as well as betweenbiomass and infiltration rate. The slow system acts between rain events, on a timescale of days to months,and evolves vegetation and soil moisture. Modeling processes at the appropriate timescales allows parametervalues to be set by the actual processes they capture. This reduces the number of parameters that are chosenexpressly to fit pattern characteristics, or to artificially slow down fast processes by the orders of magnituderequired to align their timescales with the biomass dynamics. We explore the fast-slow switching modelthrough numerical simulation on a one-dimensional hillslope, and find agreement with certain observationsabout the pattern formation phenomenon, including band spacing and upslope colonization rates. We alsofind that the predicted soil moisture dynamics are consistent with time series data that has been collectedat a banded vegetation site. This fast-slow model framework introduces a tool for investigating the possibleimpact of changes to frequency and intensity of rain events in dryland ecosystems.
Keywords. pattern formation, vegetation bands, dryland ecohydrology, reaction-advection-diffusion equa-tions, fast-slow switching model.
1. Introduction
Dryland ecosystems are subject to infrequent and largely unpredictable rain inputs [1]. These rainevents act as pulses to the system that drive dynamics across multiple scales in time and space [2, 3]. Whenthere is insufficient rain to support uniform vegetation coverage, the vegetation may become patchy, and,in some water-controlled regions these patches exhibit striking regularity in spatial arrangement [4]. Thisis particularly the case for gently sloped terrain, where patterns often appear as bands of dense vegetationgrowth alternating with bare soil, each aligned transverse to the elevation grade. The vegetation bands canbe tens of meters wide with spacing on the order of a hundred meters, and form a regular striped patternthat often occupies tens of square kilometers on the landscape [5].
Email addresses: [email protected] (Punit Gandhi), [email protected] (Sara Bonetti), [email protected] (Sarah Iams), [email protected] (Amilcare Porporato), [email protected] (Mary Silber)
Preprint submitted to Elsevier May 8, 2020 a r X i v : . [ n li n . PS ] M a y a)(b) Figure 1: (a) Satellite image of banded vegetation patterns in the Haud region of Africa (8 ◦ ◦ Figure 1(a) shows a Sentinel-2A satellite image [6] of banded vegetation patterns from the Haud regionof the Horn of Africa taken in 2016. A time series of rainfall [7] at the location is shown in Fig. 1(b). Some ofthe earliest field studies of banded patterns were carried out within this region, and postulated a connectionbetween the formation of bands and the hydrology of the area [8, 9]. For example, early observations describesheet flow in the interband region during rain events, with the rain “arrested by the next vegetation arc downthe slope” [10]. The presence of vegetation promotes water infiltration into the soil both by interceptingrunoff and enhancing infiltration as a consequence of the increased biologically-induced soil macroporosities[11, 12, 13, 14]. Such feedbacks account for the significantly higher levels of soil moisture observed within thevegetation bands than in the bare soil regions after rain events [15]. Such interactions between vegetationand hydrological processes appear as nonlinear feedbacks in mathematical modeling studies, where theyact as mechanisms underlying the formation of regular vegetation patterns [16, 17, 18]. Vegetation patternmodels often also rely on a difference in timescales associated with fast water transport and slow biomassdispersal as a key ingredient for generating pattern-forming instabilities. This separation of scales has beenexploited in analysis of the simplest class of conceptual reaction-advection-diffusion models [19, 20, 21].Hydrology-vegetation feedbacks, where the presence of vegetation increases infiltration and reduces evap-oration, have been incorporated at an annually-averaged timescale into contemporary reaction-advection-diffusion models for vegetation patterns [22, 23]. We build on and extend conceptual models due to Rieterket al. [23] and Gilad et al. [22], by developing a fast-slow mathematical framework to capture these hydro-logical and ecological processes on their appropriate timescales. We also add an additional feedback to theconceptual modeling framework: the speed of the downhill surface water flow is reduced by surface rough-ness effects in the vegetated zones [24, 25]. This new fast-slow switching framework provides an alternativeboth to detailed mechanistic models and to conceptual models that do not resolve the rainstorm timescale.Mechanistic models encode many processes on the timescales on which they occur [26, 27] but often re-quire a large number of parameters, making it difficult to identify key underlying mechanisms. Conceptual2odels formulated on a long, annually averaged, timescale [22, 23, 28] require “effective” parameter valuesto capture faster processes [29]. For example, one way to encode water infiltration and biomass evolutionon the same timescale is to reduce the rate at which water infiltrates from the surface into the soil by afactor of 1000 or more [23]. This reduction of the infiltration rate then leads surface water to be presentyear-round at a height that would be observed only during rainstorms. In the framework we use here, wetake an intermediate path. Fast hydrological processes evolve on the timescale set by a rain event, with slowecohydrological processes evolving in between rain events with a timescale set by the vegetation dynamics.This results in a conceptual model, but one that is designed to operate on the timescales of key processesin the system.By capturing processes on the timescales on which they occur, we can use infiltration and water trans-port parameters consistent with the ecohydrology literature. The model proposed in this study is able toform patterns and capture certain pattern characteristics, with ecohydrologically consistent parameters, andwithout additional parameter fitting. The natural separation of the processes by their timescales admitsa significant simplification within each of the two systems, while still capturing the influence of both fastand slow processes on dynamics. Biomass-water feedbacks captured by the model, which act on the fastrain-event timescale, include increased infiltration of surface water into the soil in the bands, and a sloweroverland surface water flow in the vegetated zones.There have been other notable efforts to extend conceptual modeling frameworks to capture the effectsthat seasonality [30], pulse intermittency and intensity [31, 32, 33], and stochasticity [34, 35] of rainfall haveon vegetation patterns. However, many of these studies fail to appropriately adjust effective parameterswhen additional timescales are captured within the model. A few of the studies are similar in philosophyto the fast-slow framework we propose in that they handle the fast hydrological processes associated withrain events separately from the rest of the slower processes. The work of Siteur et al. [35] employs simplifiedinfiltration and soil moisture dynamics, compared to those in the model we are presenting. They treatrain events as delta-function pulses of surface water, instead of resolving surface water dynamics during rainevents. The benefit of their extreme simplification is the ability to make predictions about annually averagedquantities for the case that rain events are stochastic in time and intensity. Other work [32, 33] implementsa continuous-time model, based on [22], using a numerical scheme that also simplifies the computation ofthe fast hydrological processes during rain events. In particular, [32, 33] use a steady-state relationshipbetween surface water height and biomass while it is raining, instead of resolving the time dynamics of thesurface water.The paper is organized as follows. In Section 2 we first discuss some existing model predictions andground-based observations regarding the distribution of soil moisture relative to the vegetated and bareregions. There are significant qualitative differences in predictions between various classes of models, butavailable observations suggest that the soil moisture is dramatically increased within vegetation zones relativeto the bare soil during rain events and also depleted at a faster rate via transpiration during dry periods.We then present a reaction-advection-diffusion model for banded vegetation patterns that couples processesoperating on both fast and slow timescales in Section 3. We use this coupled model as motivation for thefast-slow switching model presented in Section 4. The precipitation in the coupled model is assumed tovary in time, parameters are chosen to be consistent with the hydrological literature when available, and weexplore predictions about pattern formation within the model based on linear stability analysis of spatially-uniform states. The time-dependence of precipitation has a substantial impact on pattern characteristicssuch as precipitation level at onset, band spacing and upslope colonization rate.We use the structure of the coupled model, along with chosen parameters, to identify a timescale sepa-ration in the system: the ratio of the timescale for surface water to infiltrate into the soil to the timescalefor biomass growth is small. This small ratio suggests a fast-slow switching model, which we describe inSection 4, and explore through numerical simulation in Section 5. The simulation results of the fast-slowmodel are consistent with linear predictions from the coupled model from Section 3 with the same temporalrain input profile, and qualitatively match the observed soil moisture dynamics from the literature that ispresented in Section 2. Finally, in Section 6, we discuss key predictions of the fast-slow model in the contextof other existing modeling frameworks, highlight advantages of the fast-slow model, and suggest potentialdirections for future research. Additional simulation results for the fast-slow switching model are presented3 igure 2: Adapted from Cornet et al. (1988) [12]. (Top) Rainfall data in 1985 at a site in the Chihuahuan Desert of Mexicoexhibiting banded vegetation patterns. (Bottom) Dashed line indicate measurements of soil water content in the vegetationband and dotted line indicates the inbetween bare soil region. in Appendix A of the online supplement and the linear stability calculations of the coupled model, whichare based on Floquet theory, are detailed in Appendix B of the online supplement.
2. Soil moisture: model predictions and data from literature
Surface water dynamics occur on a minute to hour timescale, while vegetation growth and death occurson a week to month timescale. Because soil moisture is replenished by fast infiltration of surface waterand depleted by slow evapotranspiration, its dynamics occur on multiple timescales. Models of vegetationpatterns vary widely in both how they capture these multiscale dynamics, and in their predictions about thesoil moisture distribution in the vegetation bands relative to the bare soil regions. In this section, we discussdifferences in predictions among different classes of models ranging from conceptual, annually averagedones to very detailed mechanistic ones. We start the discussion by presenting ground measurements of soilmoisture at a banded vegetation site in Mexico [12].During rainfall events biomass feedback on infiltration acts to increase soil moisture where plants are byallowing surface water to enter the soil at a higher rate relative to the bare soil regions [15, 11]. However,soil moisture may be lost faster where biomass is present because of transpiration. While time-resolved soilmoisture data at sites that exhibit banded vegetation is limited, the work of Cornet et al. [12] provides atime series of soil water content for a vegetation band and the bare soil region just uphill at a site in theChihauahuan Desert for 1985. In data from [12], after rain events, the soil water content within the bandincreases dramatically relative to that of the bare soil, as shown in Fig. 2. It is only after extended periodswithout rain that the soil water content within the vegetation band becomes comparable, or even slightlybelow, the soil moisture of the bare ground. Even with these fluctuations, it appears from Fig. 2 thatthe soil moisture is, on average, more concentrated in the vegetation bands than in the bare soil regions.Qualitatively similar features occur for a site with banded vegetation patterns in Niger [36]. Data on gappedpatterns collected by Barbier et al. [37] also indicates increased soil moisture levels in vegetated regions,although the shorter-timescale dynamics are qualitatively different. More comprehensive field studies woulddetermine if similar observations occur at additional sites beyond these observational studies, and overlonger time-periods. Many predictions from models are inconsistent with these observations of soil moisturedistribution.Conceptual reaction-advection-diffusion models can make qualitative predictions about the system overcentury timescales or longer. Such long-time simulation results ensure that the asymptotic behavior of the4odel is captured, and not just transient behavior that may be sensitive to the choice of initial conditionson biomass distribution. While these conceptual models typically encode a separation of scales into theadvection and diffusion terms to capture slow biomass dispersal and fast water transport, the reaction termsassociated with local ecohydrological processes are typically formulated on the slow biomass timescale.Because these models do not resolve individual rain events, the resulting predictions are interpreted asthose of the annually averaged system, and therefore miss important details about soil moisture on shortertimescales.The simplest of these conceptual models lump processes associated with surface and subsurface hydrol-ogy into a single “water” field, leading to a pair of reaction-advection-diffusion equations that capturesinteractions between biomass and water. Such “two-field” models, exemplified by Klausmeier [28], predictthe highest concentration of water in the regions of bare soil between the vegetation bands. So-called “three-field” reaction-advection-diffusion models separately track the dynamics of surface and soil water, in additionto the biomass field. These models predict the soil moisture, on average, will peak within the biomass fortypical parameter choices [22, 23]. However, it has also been shown, in the context of the three-field modelby Gilad et al. [22], that the soil moisture spatial profile can be either in or out of phase with the spatialprofile of biomass [38], depending on which mechanism for pattern formation dominates. While such modelsare capable of qualitatively capturing annually averaged dynamics observed in Fig. 2, they fail to provideuseful predictions on shorter timescales.Detailed mechanistic ecohydrological models for general water-limited environments [26], which captureprocesses on the fast time scales of individual rain events, have also been used to study the banded vegetationpatterns. Such models predict that soil moisture can switch between being more highly concentrated underthe vegetation bands to being more highly concentrated in the bare soil region for extended periods oftime (a year or more). In these models, decadal time-averaged soil moisture is more highly concentrated inthe bare soil region [27]. The annually averaged predictions of these models, like those of the Klausmeiermodel, are in contrast with the measurements shown in Fig. 2. Moreover, the complexity of such modelslimits numerical simulations to timescales of a few decades or less, and prohibits detailed exploration of theextensive parameter space.Any prediction that soil moisture is, on annual average, more concentrated in the bare soil region, isinconsistent with the limited time-series data that is available for vegetation patterns [12, 36, 37]. Becauseof competition between increased infiltration and increased transpiration at vegetation bands relative to thebare soil regions, modeling choices about biomass growth based on plant physiology may play an importantrole in predictions about soil moisture distribution [39]. While we focus on the modeling of hydrologicalprocesses at appropriate timescales in this work, we expect that a critical analysis of model details of biomassdynamics, especially combined with more field measurements, could provide additional insights.
3. A three-field coupled timescale model
Before introducing the fast-slow switching model, the main focus of this study, we first present a reaction-advection-diffusion model in Section 3.1 with a seasonally-varying precipitation input that captures processesacross the relevant fast and slow timescales. This ‘three-field coupled timescale model’ provides motivationfor the switching model, and we discuss how parameters are chosen in Section 3.2. In Section 3.3, weexplore the linear stability of the spatially uniform states within the seasonal model. Using a dimensionlessversion of the coupled model (Section 3.4) we demonstrate that the separation of timescales leads to a smallparameter. The fast and slow systems that comprise the switching model are introduced in Section 4 byconsidering the limit in which this small parameter approaches zero. While the coupled model motivatesthe switching model, we treat the switching model as distinct from the coupled model instead of consideringit as a means to approximate solutions of the coupled model.
We first consider a model in which we input a time-dependent precipitation P ( T ) uniformly on a one-dimensional spatial domain. The model evolves a biomass field B ( X, T ) ( kg/m ), a soil moisture field5 ( X, T ) ∈ [0 , H ( X, T ) ( cm ). In restricting the description to these threedynamical variables, the model has its roots in the three-field conceptual reaction-advection-diffusion modelsof Rietkerk et al. [23] and Gilad et al. [22]. However, those models typically evolve the system on thetimescale of biomass and input the rain at a constant rate given by the mean annual value. Although thebasic structure of those models is a starting point for ours, we make modifications to capture details of thehydrological processes on the timescales at which they occur.We will refer to our version of the three-field model as a “three-field coupled timescale model” or simply“coupled model” in this paper. It takes the following form: ∂H∂T = P ( T ) − I ( H, s, B ) + K V ∂∂X (cid:18) √ ζ H δ N B (cid:19) (1a) φZ r ∂s∂T = I ( H, s, B ) − Ls − Γ Bs (1b) ∂B∂T = C (cid:16) − BK B (cid:17) Γ Bs − M B + D B ∂ B∂X (1c)where the infiltration rate of water from the surface into the soil is given by I ( H, s, B ) = K I (cid:16) B + f QB + Q (cid:17)(cid:16) HH + A (cid:17) (1 − s ) β . (2)This empirical infiltration model captures feedback from biomass, surface water and soil saturation in ahighly simplified way as compared to more physical-based approaches that rely on modeling the verticaldistribution of soil moisture [40]. Such simplified approaches have proven successful in other contexts [41].In the next section, we describe each of the terms in the coupled model (1) to highlight differences fromprevious models and to justify our parameter choices which are summarized in Table 1. A summary ofsome of the most significant differences between the coupled model (1) and, for example, the Gilad et al.model [22, 42] are: Precipitation is time-dependent, there is no soil diffusion in the coupled model, there isno root augmentation feedback in the coupled model, there is biomass feedback on surface water transportin the coupled model and, in the coupled model, the infiltration slows as soil moisture saturates and becomesindependent of surface water height for large enough surface water height. We take the Manning formula for gravity-driven open channel flow [43, 44] as abasis for our surface water transport model. According to this formula, the average velocity of the fluid in awide channel is empirically determined to be
V ∝ n − ζ / H / , where n characterizes the surface roughness, ζ is the local elevation grade and H is the height of the water. By modifying the surface roughness termbased on literature values for vegetated surfaces and bare soil, we include the influence of biomass on runoffinterception leading to a slow-down in the surface water speed within vegetated regions. In particular, wetake the water flow speed to be given by V ( H, B ) = K V √ ζ H δ − (1 + N B )where the term 1 +
N B allows us to describe the change in surface roughness from bare to vegetated soil.The value of the exponent that is consistent with the Manning formula corresponds to δ = 5 /
3. With δ = 5 / ζ = 0 .
005 (i.e. an 0 .
5% grade) we estimate theremaining parameters using Manning roughness coefficients n [43, 44]. Typical values of n for bare soil( n g = 0 . s/m / ) and dense vegetation ( n v = 0 . s/m / ) lead to K V = 2 × ( m/day ) / ( cm / ) and N = 20 m /kg . Estimations using a typical surface water height around H = 1 cm result in the samevalues of K V = 2 × m/day (note different units) and N = 20 m /kg when using δ = 1. These valuesand 1 cm surface water when δ = 5 / V ≈ cm/s on bare soil andjust over a factor of 10 slow down through vegetation with biomass value B = 0 . kg/m . While the value6arameter units default value(s) description K I cm/day
500 infiltration rate coefficient f – 0.1 bare/vegetated infiltration contrast Q kg/m A cm H -independent for H (cid:29) Aβ – 4 infiltration exponent (soil moisture) ζ – 0.005 elevation grade δ – 5 / ∗ surface water transport exponent K V m/day/cm / (or m/day ) ∗ × surface water transport coefficient N m /kg
20 surface roughness coefficient D B m /day .
01 biomass diffusion φZ r cm ∗∗ soil water capacity L cm/day cm/day ) / ( kg/m ) 0.67 transpiration coefficient K B kg/m C ( kg/m ) /cm M day − Table 1: Summary of the coupled model (1) default parameter values, as described in Section 3.2. ∗ in numerical simulations,we often use δ = 1, which changes the units of K V , although not the default numerical value we use for it. ∗∗ This value isbased on typical soil porosity φ = 0 .
45 and a root depth Z r = 60 cm .
7f the exponent that is consistent with the Manning formula is δ = 5 /
3, we take δ = 1 for most of theresults in Section 5 for computational convenience. This makes the water flow speed independent of H andthe advective transport term in Eq. (1a) linear in the surface water height H . (We present some numericalcalculations and simulations with δ = 5 / δ = 5 / δ = 1.) Soil Moisture Dynamics.
We use a “bucket model” that tracks a depth-averaged soil moisture within theplant rooting zone according to the soil moisture balance equation (1b), where s is the relative soil moisture[45]. Denoting the rooting depth by Z r and the soil porosity by φ , the quantity sφZ r represents the volumeof water contained in the root zone per unit ground area, and has units of column height. In the modelhere, we assumed typical values of Z r = 60 cm and φ = 0 .
45 [45]. Water enters the soil from the surfacevia the infiltration model (2) and is lost through evapotranspiration as ( L + Γ B ) s . The linear evaporationrate coefficient L = 0 . cm/day is set so that the evaporation rate is approximately 0 . cm/day when soilmoisture is at s = 0 .
5, which is approximately the field capacity [45]. We neglect the shading effects onevaporation by making the evaporation rate independent of biomass B . We set the transpiration parameterto be Γ = 0 .
67 ( cm/day ) / ( kg/m ) in order that the transpiration rate with biomass density B = 1 kg/m reaches approximately 0 . cm/day when soil moisture is s = 0 . K I in Eq. (2) represents an effective infiltration rate, which was here assumed equalto K I = 500 cm/day in order to ensure typical values of the infiltration rates, e.g. ∼ cm/day , in thenumerical simulation. We note, however, that the value of K I depends on the soil type. Similar to previousmodels, we assume a factor f = 0 . Q = 0 . kg/m . This biomass feedback models the presence of a soil crust that reduces infiltrationin bare soil and root systems that enhance infiltration in vegetated areas. We additionally assume thatthe infiltration rate becomes independent of surface water height for H significantly greater than A , where A = 1 cm (many descriptions of infiltration have this H -independence, see e.g. [25]). . Finally we include afactor of (1 − s ) β with β = 4 to account for the reduced infiltration values as the moisture content reachessaturation. While soil saturation effects are typically negligible for the case of constant precipitation, we areinterested in considering large, concentrated rain events where the soil moisture can increase dramaticallyin a very short time.We neglect subsurface water flow, under the assumption that surface water flow will be the dominantmode of water transport. In other models, subsurface transport has been included as linear diffusion butwe do not notice significant differences in the limited simulations where we have included it (see AppendixA of the online supplement for more details). We also neglect leakage of water from the root zone to deeperareas for similar reasons: it did not have a large impact on the limited simulation results where we haveincluded it. However, a more comprehensive exploration of the role of these two subsurface soil processes,and how they are modeled, could provide a more detailed understanding of when they can and cannot beneglected. Biomass Dynamics.
We leave the biomass dynamics largely unchanged in form from a simplified versionof the model by Gilad et al. [42]. The plant growth rate is taken to be proportional to transpiration rateΓ Bs , except with a logistic term that can be thought of as capturing a decrease in water use efficiency asbiomass increases. Specifically, the growth term is C (1 − B/K B )Γ Bs , where the proportionality constant C = 0 . kg/m ) /cm is set by the typical water-use efficiency for plants. Increased biomass begins tosignificantly decrease growth rates when B ∼ K B , where we take K B = 4 kg/m (however, biomass levelsstay well below K B in all simulations presented here). One notable difference from the biomass growth termof the Gilad et al. [42] model is that we do not include a root augmentation feedback associated with thelateral spreading of roots. Mortality is modeled by linear loss with coefficient M = 0 . day − , which wasestimated by Mauchamp et al. [47] using typical carbon maintenance costs.We model seed dispersal by linear biomass diffusion with rate D B = 0 . m /day . This parameter isoften chosen in models of this type without clear justification; D B ranges from ∼ − m /day in Giladet al. [42] to ∼ − m /day in Klausmeier [28] to 0 . m /day in Rietkerk et al. [23]. (It may ultimately8ave been chosen to help the model match pattern characteristics; see discussion in Gandhi et al. [48].) Ourlinear stability results, presented in Fig. 3, show that this parameter, for our model, controls a cut off forshort-wavelength linear instabilities of the uniform state within the seasonal model. While this parameteris unconstrained and diffusion is likely a poor model of seed dispersal, we make no adjustments to D B inour simulations, leaving it fixed at a value of 0 . m /day . We note that this model does not attempt tocapture drought resistant behavior that plants in dryland ecosystems are known to exhibit. Nor does itinclude transport of organic material or seeds via water during rain storms. If we were to set the precipitation rate to a constant, P ( T ) = P , then we could explicitly determinetwo possible types of spatially uniform steady state solutions of (1). One is the bare soil state, defined by B bs = 0, and with s bs = P L , I bs = K I f (1 − s bs ) β , H bs = A P I bs − P , (3)which exists for all P . This trivial solution becomes unstable at a transcritical bifurcation point P c ≡ M L/C
Γ, and for P > P c , there also exists a uniform vegetation solution with nonzero biomass B uv = K B CP − M L/ Γ CP + M K B , s uv = P L − Γ B uv , (4) H uv = A P I uv − P , I uv = K I B uv + f QB uv + Q (1 − s uv ) β . We find that these two uniform solutions, and the associated transcritical bifurcation at P = P c , carryover to the case that the precipitation comes into the system during n r equally-spaced rain events per year.We consider rain events with constant precipitation rate P ( cm/day ), each of duration T r . Specifically, ifwe denote the Heaviside unit step function by H , and we let P ( T ) = PH ( T r − T ) , T ∈ [0 , /n r ] , (5)then a periodic seasonal rainfall is obtained by repeating this pattern, i.e. by letting P ( T + 365 /n r ) = P ( T ).The total annual precipitation is n r T r P , and it plays the role of a mean annual precipitation rate M AP inour investigations. For the default parameters given in Table 1, the transcritical bifurcation P c occurs whenthe total annual precipitation is ∼ mm/yr , and we find that this value does not change when we usethe pulsed, periodic rain input (5) rather than constant precipitation P .The stability of the spatially-uniform vegetation state to heterogeneous perturbations is determined bycomputing Floquet multipliers associated with the linear variational equations, given perturbations propor-tional to e ikX (See, e.g., the book by Meiss [49]). Specifically, we let H k ( T ) = H ( T ) + H ,k ( T ) e ikX ,s k ( T ) = s ( T ) + s ,k ( T ) e ikX ,B k ( T ) = B ( T ) + B ,k ( T ) e ikX , where ( H ( T ) , s ( T ) , B ( T )) is a spatially-uniform, temporally periodic solution of (1) for the periodicrain input (5). In practice, we compute ( H ( T ) , s ( T ) , B ( T )) by integrating the ODEs associated withspatially uniform states of (1) forward in time until the system converges sufficiently close to this stableperiodic orbit. For the linear stability calculation, we linearize in the perturbations ( H ,k , s ,k , B ,k ), and(numerically) compute the Floquet multipliers as eigenvalues of the associated Monodromy matrix M .This calculation is carried out on a two-dimensional grid of values for mean annual precipitation ( M AP )associated with Eq. (1) and the wavenumber k of the perturbation. The stability boundary in the ( M AP, k )parameter plane is then determined by interpolating where the leading Floquet multiplier crosses the unitcircle; we refer to the instability region in the (
M AP, k )-plane as the “Turing-Hopf bubble”. (See AppendixB of the online supplement for more details of the linear stability computations.)9 AP ( mm/yr )
MAP, k )-parameter plane for default parameters of Table 1 in the case thatprecipitation consists of a six-hour long rain event every six months. Computations were done with transport exponents δ = 1(larger red bubble) and δ = 5 / δ = 1 case). The large blue bubble has constant precipitation P , the red has two equally-spacedrain events of 6 hour duration, and the small yellow bubble is associated with one rainy season of 12 hour duration. Notethat the scale has changed for this panel; the red bubble is the same as in figures (a), (c) and (d). (c) Impact of varying thetransport parameters D B , associated with biomass diffusion, and K V , associated with overland surface water flow in the casethat δ = 1 in Eq. (1). The bubble associated with default parameters of Table 1 is in red, and unless otherwise indicated D B = 0 . m /day and K V = 2 × m/day . The short-wavelength (high k ) cut-off vanishes if there is no biomass diffusion(blue region). The solid black line indicates the impact on the Turing-Hopf bubble of decreasing K V by an order of magnitude;this suppresses the longest-wave (small k ) instabilities. (d) Comparing the Turing-Hopf bubble for default parameters in red( δ = 1) with those obtained when infiltration feedback is removed (cyan, f = 1), or biomass feedback on transport is removed(yellow, N = 0 m /kg ). δ = 1; the precipitation function Eq. (5) had n r = 2 and T r = 0 . days . This means that there were two big rain events per year, separated by sixmonths (two rainy seasons), and each lasted for six hours. Here we summarize the findings from our linearstability calculations. • Figure 3(a) shows that increasing the exponent from δ = 1 to δ = 5 / δ = 1 astaking δ = 5 / • Figure 3(b) shows the impact on the instability region of changing the duration or frequency of therain events. With constant rain input P the instability region is significantly greater as shown bythe large blue Turing-Hopf bubble that extends all the way to M AP ≈ mm/yr . Alternatively, ifwe keep the intensity of rain the same but decrease the frequency from two rain events per year tojust one rain storm, the instability region shrinks quite significantly to the small yellow Turing-Hopfbubble. This suggests that the MAP level at onset of instability of the uniform state decreases as rainevents become less frequent. In fact, if we decrease the frequency of large rain events to one every 2years, then there is no pattern-forming instability of the uniform state for these parameter values. • Figure 3(c) shows how the transport parameters K V and D B in Eq. (1) impact the shape of the Turing-Hopf bubble. While there is an instability of the uniform state even when we neglect biomass diffusion( D B = 0 m /day ), there is no short-wavelength (high k ) cut-off as indicated by the asymptoticallyvertical boundary to the blue region at M AP ≈ mm/yr . As D B increases to D B = 0 . m /day we find that the short-wavelength modes are stabilized. The long-wavelength boundary is apparentlyset by the overland transport parameter K V , as suggested by the change of shape of the Turing-Hopfbubble when K V is decreased by an order of magnitude to 2 × m/day . • Figure 3(d) shows how the Turing-Hopf bubble changes when one of the two biomass-hydrology feed-backs is shut off. The N = 0 m /kg ( f = 0 .
1) Turing-Hopf bubble is quite similar in shape, althoughsmaller, to that for N = 20 m /kg ( f = 0 . only the transport feedback (i.e. when f = 1) is more circular, with an expected wavelength ofpattern that is smaller than that when the infiltration feedback is turned on.The Turing-Hopf bubbles in Fig. 3 indicate regions of linear instability of the uniform vegetation state. Theyare not equivalent to the regions of stability of nonlinear patterns which are typically referred to as “Busseballoons” [50]. Even still, the numerical simulation results presented in Section 5 suggest that these linearcalculations provide some insight into characteristics of fully nonlinear states. (cid:15) We now introduce a rescaling of the three-field coupled timescale model Eq. (1) to put it into dimension-less form. The scaling for time is chosen based on the infiltration timescale T I = A/K I which corresponds tothe characteristic time for a surface water column of height A to infiltrate into the soil under optimal condi-tions: large reservoir of surface water, dense vegetation, and dry soil. T I ≈ min for the parameter valuesof Table 1. In contrast the characteristic time associated with maximal biomass growth within this model is T G = 1 /C Γ, which corresponds to ∼ days for our default parameters. The spatial scale is chosen basedon the characteristic advection distance of surface water before infiltrating into the soil X A = K V A δ √ ζ/K I .This distance is ∼ m for our parameters. Biomass is scaled by K B = 4 kg/m and surface water heightby A = 1 cm . For this choice of scalings, namely, B = K B b, H = Ah, T = AK I t, X = K V A δ √ ζK I x p P/K I η N K B δ – f – 0.1 q Q/K B β – 4 α A/φZ r σ L/ ( φZ r C Γ) 0.11 γ K B / ( φZ r C ) 1.48 (cid:15) AC Γ /K I . × − µ M/ ( C Γ) 0.15 δ b D B K I / (cid:0) C Γ K V √ ζ (cid:1) . × − Table 2: Definitions and parameter values of the dimensionless form of the coupled model given by Eq. (6). the coupled model (1) becomes ∂h∂t = p ( t ) − ι ( b, s, h ) + ∂∂x (cid:18) ν ( b, h ) h (cid:19) (6a) ∂s∂t = α ι ( b, s, h ) − (cid:15) ( σs + γbs ) (6b) ∂b∂t = (cid:15) (cid:18) sb (1 − b ) − µb + δ b ∂ b∂x (cid:19) (6c)where ι ( b, s, h ) = (cid:18) b + qfb + q (cid:19) (cid:18) hh + 1 (cid:19) (1 − s ) β , ν ( b, h ) = h δ − ηb , and the dimensionless parameters are defined in Table 2. The parameter (cid:15) = T I /T G represents the ratio ofthe infiltration timescale to the growth timescale and is of order 10 − . We exploit the smallness of this ratioin the following section to motivate the switching model that is the focus of the rest of the paper.
4. A fast-slow switching model
Rainfall initiates fast hydrological processes associated with overland waterflow and, with it, infiltrationfrom the surface into the soil. These processes occur on timescales of minutes to hours. Evapotranspirationand plant growth occur on much longer timescales, days to months, and we assume can be neglected duringrain events. This separation of scales is evident in the coupled model presented in Section 3 from the smallparameter (cid:15) in Eq. (6) which represents the ratio of the water infiltration timescale T I to the biomass growthtimescale T G . In this section, we present a switching model that captures dynamics of fast processes initiatedby rainfall and switches to the dynamics of slower processes in the absence of surface water. Figure 4 showsa schematic diagram of such a model, which we describe in detail in Section 4.1. We then, in Section 4.2,discuss numerical methods employed to explore the model along with some potential issues associated withour choice of periodic boundary conditions. 12 ast System : h ( x, t ), s ( x, t )Equation (7): precipitation,surface water transport, infil-tration from surface into soilsoil moisture from Fast , after h is negli-gible, initializes Slow Slow System : s ( x, τ ), b ( x, τ )Equation (9): evapotran-spiration; biomass growth,mortality and dispersalSoil moisture andbiomass from Slow ,once rain starts, ini-tialize
Fast . h ( x, t f ) ∼ O ( (cid:15) ) s f ( x ) p ( τ s ) > s s ( x ), b s ( x ) Figure 4: The fast system is initialized using the current biomass profile and soil moisture distribution. It evolves the surfacewater and soil moisture on the fast time until surface water is no longer present. The final soil moisture distribution andunchanged biomass profile are then used to initialize the slow system. It evolves the soil moisture and biomass until the startof the next rain event, and the cycle repeats. (cid:15) → (cid:15) (cid:28) (cid:15) = 0, corresponding to the assumption thatthe dynamics occurring on the biomass growth timescale are completely negligible. The (cid:15) = 0 assumptionleads to a fixed biomass profile b ( x ) since ∂b/∂t = 0 from Eq. (6c). The remaining two equations ofsystem (6) become: ∂h∂t = p ( t ) − (cid:18) b + qfb + q (cid:19) (cid:18) hh + 1 (cid:19) (1 − s ) β + ∂∂x (cid:18) h δ ηb (cid:19) (7a) ∂s∂t = α (cid:18) b + qfb + q (cid:19) (cid:18) hh + 1 (cid:19) (1 − s ) β . (7b)We refer to Eq. (7) as the “fast system,” which models the dynamics during and shortly after rain events.(For convenience, we substituted in the expressions for infiltration rate ι ( b, s, h ) and transport speed ν ( b, h )that are defined below Eq. (6).)We start evolution of the fast system when a rain event begins, and continue until the surface water h becomes negligible, i.e. below a threshold of order (cid:15) . We expect these fast dynamics to occur over an hourstimescale, after which we obtain an updated soil moisture profile s fast system −−−−−−−→ s + θ ( b, s, p ) (8)that has been replenished by the rain event. The profile of the soil moisture replenishment term θ isdetermined by the fast system (7), and depends on the biomass b ( x ) and soil moisture s ( x ) profiles prior tothe rain event along with the rain event p ( t ) itself.The fast dynamics of Eq. (6) are terminated once h ∼ O ( (cid:15) ), after which the soil moisture and biomassevolution proceeds on a slow timescale given by τ = (cid:15)t under the assumption that the surface water heightremains fixed at h = 0. Rescaling time in Eqs. (6b) and (6c) then leads to ∂s∂τ = − ( σs + γsb ) (9a) ∂b∂τ = sb (1 − b ) − µb + δ b ∂ b∂x . (9b)13 igure 5: Spacetime plot of surface water height in units of cm shown with the blue indicating 2 cm . The green indicateswhere vegetation is above 0.01 kg/m . The surface water, initialized at 2 cm along the top 15 m of the domain falls to below0.2 mm by the trailing edge of the vegetation band at X ≈ m . Parameter values used in Eqs. (7) and (9) are given byTable 2 with δ = 1. We refer to Eq. (9) as the “slow system” and use it to evolve biomass and soil moisture in the absenceof surface water. The slow system (9) is initialized with the biomass profile b prior to the rain event andthe replenished soil moisture s + θ , and evolves until the next rain event occurs. In results and discussionsections that follow, we numerically explore this fast-slow model that has been motivated by the coupledmodel (1). We discretize in space with finite differences using first-order upwinding for advection and second-ordercentered differences for diffusion. For time integration we employ Matlab’s ode15s, an implicit variable-ordersolver designed for stiff equations [51]. At the start of each rain event, the fast system is initialized with theoutput biomass (which will remain fixed), soil moisture of the slow system and no initial surface water. Thefast system is run for twice as long as the rain event, with a constant rain input for the first half of this time.The surface water height typically has fallen to below 10 − cm by this point, and the simulation proceedsto the slow system initialized with the updated soil moisture and previously output biomass. However, ifthe surface water is not below a threshold value of H thresh = 0 . cm , the fast system is run again with norain input for fixed intervals of time until the threshold is reached and the simulation can move on to theslow system.The simulations are carried out on a periodic domain that typically corresponds to a few wavelengthsof the banded pattern. The underlying assumption for this choice of boundary conditions is that we areconsidering a small section in the middle of a very long hillslope. One potential issue that may arise withperiodic boundary conditions is that, during the fast system, surface water may repeatedly reach the bottomof the domain and get re-injected at the top of the domain. We don’t expect surface water to travel morethan a few vegetation bands before infiltrating into the soil, and find that our simulations are consistentwith this expectation. For example, a simulation of the fast system (7) with parameters from Table 1 showsthat 2 cm of surface water at the top 15 m of a hillslope is decreased to approximately 1% of its originalamount after traveling through an 80 m wide band of vegetation with peak biomass value of 0 . kg/m .The water travels with a speed of approximately 4 cm/s through the vegetation band and 13 cm/s over thebare soil. This simulation, which is shown in Fig. 5, is initialized with the state reached just prior to a rainevent at 200 years from the simulation shown in Fig. 8 of Section 5.2.
5. Numerical results for fast-slow model
In this section, we explore pattern formation within the fast-slow switching model, Eqs. (7) and (9),through numerical simulation. Inspired by the biannual rainy seasons shown in Fig. 1(b), we take as a base-case scenario two major rain events each year, spaced six months apart, with a mean annual precipitationof 160 mm/year . In this case, each rain event lasts for six hours, and deposits half of the mean annual14
100 200 300 400 5000120 100 200 300 400 500050100 0 100 200 300 400 50000.5 0 100 200 300 400 50000.5 0 100 200 300 400 5000120 100 200 300 400 500050100 0 100 200 300 400 50000.5 0 100 200 300 400 50000.5 (a) (b)
Figure 6: Spatial profiles of patterns with initial perturbation of wavenumber (a) k = 2 π/ m − and (b) k = 2 π/ m − .For each wavenumber, surface water height H , Infiltration rate I , soil moisture s , and biomass density B at t = 3000 years with mean annual precipitation of 160 mm/year are shown. The solid line is the annually averaged profile while the dottedlines show the pointwise minimum and maximum values over the course of the year. Parameter values used in Eqs. (7) and (9)are given by Table 2, and the rain is input uniformly over two evenly-spaced six-hour rainstorms per year. precipitation at a constant rate during that short time-frame. Unless otherwise noted, we use parametervalues from Table 2 with δ = 1, corresponding to dimensional quantities given in Table 1, throughout thissection. We also report all quantities in dimensional units, with the conversions described in Section 3.4.(See also parameters in Table 1.)We first study characteristics of stable periodic solutions obtained using sinusoidal perturbations of theuniform state as initial conditions in Section 5.1. These solutions provide a basis for interpreting the resultsof simulations initialized with random initial perturbations, which we discuss in Section 5.2. We provide adetailed picture of transients that persist on ecologically relevant timescales, e.g. decades to centuries, forthis case and describe the observed model behavior on the timescale of millennia. We then, in Section 5.4,study the influence that the amount of rainfall has on pattern characteristics. Finally, in Section 5.5, weexplore the role that biomass feedback plays in infiltration and surface water transport in the pattern formingprocess. Additional supporting simulation results are reported in Appendix A of the online supplement. We consider a spatial domain of length L = 500 m and take the positive x -direction to be uphill. Weinitialize the fast-slow system with a 1% sinusoidal perturbation of wavenumber k j = 2 πj/L , j = 1 . . .
30, tothe uniform solution with constant precipitation given by Eq. (4). For 1 ≤ j ≤
19, the simulation convergesto a periodic traveling wave solution with wavenumber k j , and for j >
19 the simulation converges toa periodic pattern with different wavenumber (typically k j with 4 ≤ j ≤ ∼ m ( j = 19) all the way to the size of the domain 500 m ( j = 1).Figure 6 shows examples of spatial profiles of two different wavenumber solutions for the fixed parameterset considered here. The solid line in each plot represents the annual average while the dotted lines indicatethe pointwise maximum and minimum. The average biomass, the surface water and soil moisture all peakat the leading edge of the vegetation bands. While the surface water and infiltration rates are near zero onaverage, the maximum values are at around 2 cm and 100 cm/day , respectively, during rain events. The soilmoisture within the vegetation bands also varies significantly in time, going from s = 0 . (a)(b)(c) Figure 7: Characteristics of solutions from simulation as a function of pattern wavenumber. (a) Biomass averaged over timeand space, and peak value of annually-averaged biomass are indicated by solid and dotted lines, respectively. (b) Fraction ofwavelength covered by vegetation band. (c) Uphill migration frequency of biomass pulse in wavelengths of the pattern percentury. Parameter values used in Eqs. (7) and (9) are given by Table 2, and the rain is input uniformly over two evenly-spacedsix-hour rainstorms per year with mean annual precipitation of 160 mm/year . Figure 7 indicates that many pattern characteristics (e.g. average biomass, fraction of wavelength coveredby vegetation) show remarkably little variation with respect to wavenumber. Panel (a) indicates a monotonicdecrease in the maximum biomass value as a function of wavenumber. The fraction of the wavelength withvegetation cover, shown in Panel (b), varies little with k , on the order of 5% with a maximum at around n = 5. To compute the fractional cover, we take the edges of the vegetation band to be the points of steepestincrease/decrease in biomass. Using a threshold value of Q = 0 . kg/m to define the edges leads to nearlyidentical fractional cover values.Unlike the average biomass and fractional cover, the average migration speed of the bands does varysignificantly with wavenumber: it decreases monotonically from about 200 cm/year for the largest bandspacing of 500 m ( k = 2 π/ m − ) to about 10 cm/year for a small spacing of 28 m ( k = 36 π/ m − ),and has a value of 65 cm/year for patterns with a wavelength of 100 m . As a possible reference point,observational studies of banded patterns in the Horn of Africa record typical wavelengths of 40 − m andmigration speeds measured in 10s of centimeters per year, on slopes with typical elevation grades 0 . − . . − . m . Deblauwe et al. [5], drawingon pattern observations around the globe, report slightly smaller mean migration speed values in the rangeof 0 . − . constant rain input tothe coupled model (1) are way off, both in comparison to the reported observations and to our fast-slowsimulations that incorporate short rain events. For instance, the onset of Turing patterns for a constantrain input have wavelength of 26 m , and travel rapidly at a speed of ∼ m/yr (i.e. ∼
13 wavelengths percentury). Moreover, the associated Turing-Hopf bubble of Fig. 3(b) extends all the way to a
M AP ≈ mm/yr , giving a rather high lower bound to where patterns would occur. Using linear predictions from thecase of constant rain input, and trying to adjust the parameters to fit those observed pattern characteristics,16
100 200 300 400 5000120 100 200 300 400 500050100 0 100 200 300 400 50000.5 0 100 200 300 400 50000.20.4 (b)(a)
Figure 8: (a) Spatial profile of surface water height H , infiltration rate I , soil moisture s , and biomass density B at t = 200 years . The solid line is the annually averaged profile while the dotted lines show the pointwise minimum and maximum valuesover the course of the year. (b) Spacetime plot of annually averaged biomass B in units of kg/m over the course of the 200 year simulation. Yellow indicates low biomass while green indicates high biomass, and the color scale range is 0 < B < kg/m .All subsequent biomass spacetime plots use the same color scale described here. Parameter values used in Eqs. (7) and (9) aregiven by Table 2 with δ = 1, and the rain is input uniformly over two two evenly-spaced six-hour rainstorms per year withmean annual precipitation of 160 mm/year . would lead to different parameter choices. While the periodic states summarized in Figs. 6 and 7 persist for at least 5000 years and are stable tosmall perturbations, long-time simulations with random initial perturbations eventually approach the long-wavelength patterns. We first provide a detailed account of long-lived transients that persist on ecologicallyrelevant timescales of centuries before discussing this asymptotic behavior on the scale of millennia.Figures 8 to 10 show details of a 200-year simulation on a 500 m domain with mean annual precipitationof 160 mm/year . These are the same parameters we used for the periodic patterns in the previous section.However, we modify the initial condition to be a 1% random perturbation of the uniformly vegetated statethat’s obtained for constant precipitation, given by Eq. (4). These values, B = 0 . kg/m and s = 0 .
100 200 300 400 50000.10.2 (a) (b)(c) (d) (cid:54) r a i n H ( cm ) I ( cm/day ) Figure 9: Fast system during last year of simulation shown in Fig. 8. (a) Fixed biomass profile during simulation of fast system.(b) Initial (final) soil water content ( φZ r s ) in units of cm shown with dotted (solid) line. (c) Spacetime plot of surface waterheight H in units of cm during and shortly after a 6-hour storm with 80 mm of rainfall in the final year of a 200-year simulation.(d) Spacetime plot of instantaneous infiltration rate in units of cm/day during the same time period. The simulation evolves the fast system, Eq. (7), during and shortly after each rain event. Specifically,constant precipitation, at a rate of about 13 mm/hour , is input during the first six hours and none duringthe next six hours, during which time the surface water height drops to below 10 − cm . The biomass profile,which remains fixed during this time, is shown in Fig. 9(a). Figure 9(b) shows the initial and final soil watercontent Φ Z r s in units of cm . A total of 8 cm of water is input uniformly across the domain by the rainstorm and the vegetation bands collect approximately 10.6 cm on average, while the bare soil regions collectapproximately 1.7 cm on average. The vegetated regions therefore absorb, on average, approximately 6.2times the amount of rainfall because of runoff from the bare soil regions. Figure 9(c,d) show spacetimediagrams for surface water height H and instantaneous infiltration rate I over the 12 hour period that thefast system is run. During the rain, both the surface water height and the infiltration rate are enhanced bythe presence of biomass. The surface water is increased in the vegetation bands because increased surfaceroughness slows water transport. This increased height, along with a positive feedback of biomass directlyon infiltration rate, lead to the heightened accumulation of soil water content in the vegetation bands thatis seen in Fig. 9(b).At the end of the 12-hour period of the fast system, we then initialize the slow system with the updatedsoil moisture from the fast system and the same biomass profile. Figure 10 tracks soil moisture s andbiomass B at locations in the vegetation band (blue and green) and in bare soil (orange) over time for theslow system during the final year of the simulation shown in Fig. 8. The soil moisture in the vegetationband initially increases significantly, but eventually falls to nearly the same level as the bare soil just beforethe next rain event. The biomass peak is delayed by approximately two months after the rain event.While the behavior in Figs. 8 to 10, initialized with random perturbation, persists for a millennium, 3500-year simulations (see Fig. A2) reveal it to be transient. This is in contrast to the cases where a periodicperturbation of specific wavenumber is introduced, as described in Section 5.1. The periodic states in Figs. 6and 7 persist as traveling wave solutions (on annually averaged timescales) for at least 5000 years.The long transients persist on ecologically relevant timescales and have wavelengths consistent withpredictions from linear stability analysis of the coupled model, Eq. (1), with identical rain input (i.e. six-hour constant rain intensity storms every six months). Those calculations indicate an onset of the pattern-forming instability occurs for a total annual rain of ∼ mm (input in two large rain events), and are18
50 100 150 200 250 300 35000.5 (a)(b)
Figure 10: Time series of (a) soil moisture and (b) biomass during the last year of the simulation shown in Figs. 8 and 9.The blue (green) lines indicate soil moisture (biomass) at the peak of the vegetation band ( x ≈ m ), and the orange linesindicate the center of the uphill bare soil region ( x ≈ m ). The cyan vertical lines indicate the time of the rain events. associated with a wavelength of ∼ m . For mean annual precipitation of 160 mm/year , the “Turing-Hopf bubble” has a width that indicates the uniform state is unstable to perturbations with wavenumbersbetween k = 0 . m − and k = 0 . m − , corresponding to wavelengths greater than ∼ m and allthe way up to a kilometer scale; the most unstable mode for these conditions occurs for k ≈ . m − (awavelength of ∼ m ), where the Floquet multiplier reaches a peak value of 1 .
15. Simulations of thefast-slow system on a 500 m domain, with an initial periodic perturbation of fixed wavenumber within thisTuring instability range, appear to be stable, at least for wavelengths greater than ∼ m , i.e. they do notundergo band-merging events on the millennial simulation timescale. (See results presented in Fig. 7.) The theory for open channel flow suggests a nonlinear dependence of surface water height on transport,e.g. δ = 5 / δ = 1 out of numerical convenience for the majority of simulationspresented here in Section 5. The predictions of the linear theory for the three-field coupled model presentedin Fig. 3(a) show relatively modest differences between the δ = 1 and δ = 5 / δ = 5 /
3, compared to δ = 1. To illustrate the minor quantitative differences in periodicsolutions for the two cases, we take the solution with k = 2 π/ m − for the case of δ = 1, shown inFig. 6(b), as an initial condition for a simulation with δ = 5 /
3. A comparison of this initial condition tothe final state that the simulation converges to is shown in Fig. 11. The resulting periodic state exhibitsslightly less biomass for δ = 5 /
3: the bands are slightly narrower with slightly lower peak values. Despitethese minor quantitative differences, the overall qualitative character of the solution profiles is the same.(Another example, presented as Fig. A3 in Appendix A of the online supplement, starts with a randomperturbation of the uniform state, and shows a much slower development of patterns for δ = 5 /
3, indicatinganother aspect of the numerical speed-up we achieve for δ = 1.) In this section we explore the impact on the patterns of ramping down the mean annual precipitationand then ramping it back up, with details on our protocol for this provided below. We find that the range ofexistence for patterns is considerably greater than suggested by the linear theory of Section 3.3; rather thanpatterns only existing in the range of 109 − mm/yr , they exist as highly nonlinear states in the range of ∼ − mm/yr . We also find that on gradually decreasing the precipitation there are a number of discretesteps that reduce the number of bands on the 500 m domain until only one band persists. However, as wethen increase the precipitation, these jumps in the number of bands do not reverse themselves. Instead thesingle band persists but with ever increasing width until it eventually fills the domain as a spatially uniform19
100 200 300 400 500010 100 200 300 400 500050 0 100 200 300 400 50000.5 0 100 200 300 400 50000.20.4
Figure 11: Nonlinear dependence of transport on surface water height ( δ = 5 / δ = 5 / δ = 1. Colored lines show, from top tobottom: spatial profile of surface water height H , infiltration rate I , soil moisture s , and biomass density B at t = 200 years .The solid line is the annually averaged profile while the dotted lines show the pointwise minimum and maximum values over thecourse of the year. The profiles of the initial condition, a solution when δ = 1, are shown in light gray for reference. Parametervalues used in Eqs. (7) and (9) are given by Table 2, and the rain is input uniformly over two evenly-spaced six-hour rainstormsper year with mean annual precipitation of 160 mm/year . vegetation state. We find that the average biomass decreases approximately linearly with precipitationlevel, but that the peak biomass values in the bands are relatively constant. In particular, it is the fractioncoverage per wavelength that contributes most significantly to the change in average biomass on the domain,and not the peak biomass level in the bands.Results associated with decreasing mean annual precipitation are summarized in Fig. 12. In order todetermine where the uniform state loses stability, we begin with the stable uniform vegetation state at240 mm/year and decrease the precipitation amount in steps of 2 mm/year . At each step, we add 1%perturbations and allow the simulation to run for 300 years (or shorter if it converges to a uniform statemore quickly). We find that the uniform state of the fast-slow system remains stable until 174 mm/year .Perturbing the uniform state at 174 mm/yr with 1% sinusioidal perturbations of wavenumber k j = 2 πj/L with j = 1 . . .
19, we find that the uniform state is unstable to the 3 ≤ j ≤ m . This is quantitatively consistent with the predictions of lineartheory for the three-field coupled model presented in Section 3. Interestingly, the average biomass on thedomain for these highly nonlinear patterned states is approximately equal to the maximum biomass obtainedfor the uniform vegetation state.We continue decreasing the mean annual precipitation from 174 mm/yr , in 2 mm/year steps, startingwith the j = 5 wavenumber pattern, corresponding to a wavelength of 100 m . This state persists until theprecipitation reaches 52 mm/year , well below the transcritical bifurcation that creates the uniform state at109 mm/year . Moreover, the patterned state sustains more biomass on average than the uniform state, forprecipitation values where the uniform state does exist. While the patterns maximum biomass value doesnot change very much over this range, the fraction of the domain covered by vegetation decreases linearlywith precipitation, i.e. there is a linear decrease in the width of the vegetation bands with M AP . At 52 mm/year , a pattern with 250 m wavelength emerges from the simulation and transitions to a single bandon the 500 m domain at 40 mm/year . This so-called ‘oasis state’ finally collapses to the bare soil state at34 mm/year .As shown in Figure 13, we also scan up in precipitation, starting with the oasis state at 35 mm/year and moving in increments of 2 mm/year as before. This state with a single vegetation band on the domainpersists all the way up to 201 mm/year before loosing stability to the uniform vegetation state. As theprecipitation increases, the width of the vegetation band increases. We note that at around 179 mm/year ,20
100 200 300 400 50000.20.4
20 40 60 80 100 120 140 160 180 200 220 24000.10.20.30.4
20 40 60 80 100 120 140 160 180 200 220 24000.51
Figure 12: Decreasing Precipitation. The green solid (dotted) line indicates average (minimum and maximum) values ofbiomass for uniformly vegetated state. At 174 mm/year , this spatially uniform state is unstable to a periodic perturbationwith wavenumber k = 2 π/ m − . The red left arrow (dot) indicate the maximum (average) value of the pattern thatemerges. At 52 mm/year this pattern with wavelength 100 m is unstable and a pattern with wavelength 250 m emerges,and at 40 mm/year a pattern with a single wavelength on the domain is selected. At 34 mm/year this single-wavelegnthpattern collapses to the bare soil state. Average and min/max values of the biomass profiles are shown for the three differentwavelengths. When it exists, the uniform vegetation state is included in gray for the given precipitation value. A plot of thefraction of a wavelength covered by biomass shows linearly decreasing coverage with precipitation. Parameter values used inEqs. (7) and (9) are given by Table 2, and the rain is input uniformly over two evenly-spaced six-hour rainstorms per year.
20 40 60 80 100 120 140 160 180 200 220 24000.10.20.30.4
Figure 13: The green lines and red left-arrows, from Fig. 12, indicate the uniform vegetation state and the patterns obtained bydecreasing precipitation with steps of 2 mm/year . The blue right-arrows indicate the maximum of the biomass for increasingprecipitation. At each step of 2 mm/year , the state remains a single biomass pulse on the domain but with increased width.Parameter values used in Eqs. (7) and (9) are given by Table 2, and the rain is input uniformly over two evenly-spaced six-hourrainstorms.
100 200 300 400 5000120 100 200 300 400 5000500 100 200 300 400 50000.50 100 200 300 400 50000.20.4 (a) transport feedback offinfiltration feedback on (b) both feedbacks on (c) transport feedback oninfiltration feedback offCorresponding soil moisture time series s ( T ) for the bare and vegetated soil: Figure 14: Influence of biomass feedback on surface water transport and infiltration rate with mean annual precipitation of140 mm/year . For each of the following choices of ( f, N ), we show spatial profiles during the last year of the simulation abovethe spacetime plot of the biomass from a 1000-year simulation. Below those, we present the time series of soil moisture duringthe last year at the location associated with peak biomass in a vegetation band (blue) and in the center of a bare soil region(orange). Standard parameter choices from Table 1 are used in each case except that: (a) transport feedback is turned off, N = 0 m /kg and f = 0 . N = 20 m /kg and f = 0 . N = 20 m /kg and f = 1. All other parameter values used in Eqs. (7) and (9) are given by Table 2 with δ = 1, and therain is input uniformly over two evenly-spaced six-hour rainstorms per year. there is a qualitative change in the shape of biomass profile, in terms of its concavity. Above 179 mm/year ,the biomass values at the trailing edge of the vegetation band correspond to those of the uniform vegetationstate for the given precipitation level. If we again decrease precipitation, starting now with a one-pulsesolution, we find an apparent bistability of two distinct single-pulse states for precipitation between approx-imately 175 and 177 mm/year , which is in a region where the uniform state is also (barely) stable. SeeFig. A4 in Appendix A of the online supplement for more details. The fast component of our model (7) incorporates two important feedbacks between the biomass andthe water resource re-distribution, and in this section we explore the relative impact of those feedbackson the pattern characteristics. The presence of biomass is known to enhance infiltration of surface waterinto the soil within the dryland ecosystems that exhibit banded vegetation patterns. This positive feedbackplays a central role in the pattern formation process within many conceptual models [22, 28, 46]. Themodel we consider here captures an additional feedback of biomass on the hydrology, namely an increase in22urface roughness that can slow down water flow where biomass is present [24, 25]. Our standard choice ofparameters includes both kinds of biomass feedback but, as Fig. 14 illustrates, patterns are possible with either feedback alone. These results are consistent with linear predictions from the three-field coupled modelshown in Fig. 3(d).The simulations for Fig. 14 are initialized with a random perturbation of the uniform state, and carriedout with mean annual precipitation of 140 mm/year . The parameters characterizing the biomass feedbackon surface water transport, N in Eq. (1a), and infiltration rate, f in the infiltration function I in Eqs. (1a)and (1b), are varied with all other parameter values taking the default parameters from Table 1. Simulationresults with the standard parameters, N = 20 m /kg and f = 0 .
1, are shown in Fig. 14(b). These results areto be contrasted with those of columns (a) and (c) for which either the transport or the infiltration feedbackis turned off by setting N = 0 or f = 1, respectively. The top row shows plots of the spatial profiles ofsurface water, infiltration rate, soil moisture and biomass during the last year of the simulation. The middlerow shows spacetime diagrams of the annually-averaged biomass. The bottom plots show time series of thesoil moisture at the peak biomass location within the vegetation band (blue) and in the center of a baresoil region (orange) during the last year of the simulation. The soil moisture plot provides a qualitativecomparison for field measurements, as shown in Fig. 2.We make the following observations about Fig. 14: • While the width of the bare soil region is somewhat comparable between all three cases, the widthsof the vegetation bands are significantly greater when there is no feedback in the transport, i.e. N =0 m /kg . This leads to far fewer bands on the domain in this case, compared to the cases withtransport feedback. • In absence of the infiltration feedback i.e. f = 1, the bands are more regular in spacing and travelapproximately twice as fast across the domain and with uniform speed. In contrast, when there is notransport feedback, i.e. N = 0, then we see a lag between the leading edge colonization and trailingedge death. • The asymmetry in the profiles of the vegetation bands, between leading uphill edge and trailingdownhill edge, is much more pronounced in the presence of the transport feedback. In fact, theprofiles are almost triangular when the infiltration feedback is turned off, and almost rectangularwhen there is no transport feedback. This difference in profiles is also quite pronounced in the soilmoisture profiles. • The transport feedback is required to concentrate the surface water in the bands. Morever, the peaksof the surface water are lower without an infiltration feedback, and greatest when both feedbacks arepresent. • The time series of the soil moisture in the lower plots shows then when there is no infiltration feedbackthere is a decreased contrast in soil moisture between the bare soil and vegetation band immediatelyfollowing the rain event. • When there is no infiltration feedback the soil moisture at the peak of the biomass actually falls belowthe soil moisture in the bare soil zone. Specifically, a switching point at around 100 days after therain event is seen in the time series of the soil moisture, where s of the bare soil becomes higher than s within the vegetation band. This switching behavior is also observed in some drought phases of thefield observations illustrated by Fig. 2.
6. Discussion
In this paper we have motivated and then investigated a fast-slow mathematical modeling framework,Eqs. (7) and (9), for vegetation pattern formation. This fast-slow framework exploits the inherent separationof timescales between fast hydrological processes and slow biomass dynamics. Our fast-slow model switchesbetween a fast system, Eq. (7), that resolves aspects of the surface water and soil moisture dynamics during23ain events at the fast timescale of minutes to hours, and a slow system, Eq. (9), that captures the dynamicsof biomass and soil moisture interactions on the slow timescale of weeks to months when no surface wateris present. By modeling key processes on the timescales at which they occur, we are able to employ aconceptual-level model with parameter values that are consistent with observations of the processes beingmodeled. Soil moisture dynamics occur on the fast timescale, through infiltration, as well as the slowtimescale, through evapotranspiration. A comprehensive set of appropriately resolved time series data onsoil moisture for banded vegetation patterns would be useful to probe this modeling framework on itsmultiple timescales.We do not choose parameters to fit predicted pattern characteristics to observations of banded vegetationpatterns, as has been done for models that were developed for the slow timescale of the biomass dynamicsalone (See, e.g. [37]). Even still, the fast-slow switching model is able to accurately capture certain details ofthe phenomenon with hydrologically-informed parameter values. Numerical simulation of the model revealsthat the spacing of the vegetation bands, and their up-slope colonization speeds are of the right order ofmagnitude. Another focus of our investigation relates to the distribution of the soil moisture relative to thebiomass. Our expectation, based on limited ground-based studies of this, is that the soil moisture should betrapped for much of the year where the biomass is located. Our simulation results are consistent with this.However, we make a number of simplifying modeling assumptions, including to neglect subsurface watertransport and the influence of lateral spreading of roots. More careful thought must also be given to thebiomass model in future work to ensure that their key processes are being captured.Many model studies suggest that bare-soil areas should increase in size to compensate for a decreasein precipitation level [54]. The fast-slow switching model further predicts that whether this compensationoccurs by adjusting the band spacing or the band width depends on whether the rainfall is increasing ordecreasing. Our most extensive model simulation results focus on a situation where there are two heavy rainevents of equal strength per year, which drive the pattern formation process. From those simulations, thereis strong evidence for multi-stability, history-dependence, and hysteresis. In some of the longest simulationswe gradually decrease and then increase the annual precipitation levels. As precipitation levels decrease, agradual decrease in band width is observed. This narrowing is accompanied by occasional losses of bandsleading to increases in band spacing. For low precipitation levels, the vegetation bands are narrow, withlarge swaths of bare soil in between them. As the precipitation level increases, the widths of the bandsincrease without any change to the wavelength of the pattern. These results are consistent with predictionssuggested to Hemming by his detailed observations of banded patterns in Northern Somalia that he reportedin 1965 [10]. In particular, he predicted that increasing rainfall might widen individual arcs while decreasingrainfall might reduce the number of arcs.The fast-slow modeling framework is well-suited for exploring the impact that the temporal rainfallpattern has on spatial vegetation pattern formation in these ecosystems, and this will be a major focus offuture efforts. Preliminary simulation results, which are described further in the online supplement, indicatethat the temporal distribution of rain events influences whether patterns will occur at a given level of meanannual precipitation, but more work is needed to understand this in detail. With fewer rain events (at thesame annual precipitation level), patterns occur within a smaller
M AP range. For one twelve-hour rainevent each year, at our standard parameters (see Table 2), patterns occur at a range of 89 to 147 mm/year .With the same rain amount spread between two six hour events each year, the range grows to 34 to 201 mm/year (Fig. 13). A previous study by Guttal et al. [30] on the effect of seasonal rain input withinthe three-field model by Rietkerk et al. [23] found the opposite trend. However they did not attempt toincorporate timescales associated with individual rain events and instead focused on capturing the impactof drought on the biomass dynamics.It may be tempting to interpret the presence of patterns as a sign of stress in the ecosystem, with theirexistence at higher
M AP levels indicating reduced resilience. However, we caution against this. In the fast-slow model presented here, the patterned states have a higher total biomass than the corresponding uniformvegetation state for the same precipitation level. In our ramped precipitation simulations summarizedby Fig. 13, we find that the pattern state only disappears (around 201 mm/yr ) when it’s no longer “anadvantage” to the system over the uniform vegetation case.While the fast-slow switching model is motivated by a small parameter in a three-field coupled timescale24odel, Eq. (6), we do not consider the fast-slow model as an approximation of the coupled model in thesense of singular perturbation theory. However we do show that linear predictions of the coupled modelare consistent with simulations of the fast-slow model when the same temporal rain profile is used. Themathematical relationship between these two models will be explored at a more formal level in future work.Our goal is to be able to start from a more detailed model that captures fast hydrological processes andformally make the reduction to a two-field model that operates on the slow timescale. While this has beendone in the high infiltration limit where there is no overland water flow [55], capturing the influence ofsurface water transport on soil moisture through an effective transport could help make a direct link tomodels like those of Klausmeier. Moreover, it may provide insight into how well rainfall is approximatedby the inclusion of an additive constant to the soil moisture rate equation, as is commonly done, or suggestalternate approaches that may be more appropriate.In future work, we are particularly interested in considering the influence of stochastic rainfall withinthe fast-slow switching framework. One approach for doing so, conceptually similar to the work of Siteuret al. [35], is to develop the interpretation that the fast system provides a kick to the slow evolution ofsoil moisture and biomass. The kicks from the fast system, under certain simplifying approximations, canbe solved in closed form via the method of characteristics. This may lead to progress towards analyticinsights as well as an increased speed-up of numerical simulations. Alternatively, numerical speed-up couldbe achieved by machine learning approaches [56], or by more parsimonious modeling of the fast system.The fast-slow framework introduced in this paper is a conceptual model that captures the naturaltimescales of important processes. In particular, it resolves the fast hydrological processes that contributeimportant positive feedbacks between the biomass and the water resource. It produces physically reasonablepredictions in the absence of parameter tuning and provides a testbed for exploring the impact of differentprecipitation regimes on patterns. At a practical level, this framework’s advantage is two-fold: (1) capturingboth the fast and slow timescales in the model allows for the parameters to be chosen based on informationabout the particular processes being modeled, and (2) separating fast and slow processes in a switchingmodel provides a computational speed-up over a coupled timescale model. We also believe that there is aconceptual advantage in that the fast-slow switching framework captures the qualitatively different behav-iors of the system when surface water is and is not present in a very transparent way.
Acknowledgements.
The work by PG was supported in part by the National Science Foundation grantDMS-1440386 to the Mathematical Biosciences Institute. The research of MS was supported by NSF-DMS-1517416. PG and MS also benefited from participating in the wonderful birthday conference,
Advances inPattern Formation: New Questions Motivated by Applications , honoring Ehud Meron’s pioneering work inthe field of pattern formation.
References [1] I. Noy-Meir. Desert ecosystems: environment and producers.
Annu. Rev. Ecol. Syst. , 4(1):25–51, 1973.[2] S. Schwinning, O. E. Sala, M. E. Loik, and J. R. Ehleringer. Thresholds, memory, and seasonality: understanding pulsedynamics in arid/semi-arid ecosystems.
Oecologia , 141(2):191–193, 2004.[3] S. L. Collins, J. Belnap, N. B. Grimm, J. A. Rudgers, C. N. Dahm, P. D’odorico, M. Litvak, D. O. Natvig, D. C. Peters,W. T. Pockman, R. L. Sinsabaugh, and B. O. Wolf. A multiscale, hierarchical model of pulse dynamics in arid-landecosystems.
Annu. Rev. Ecol. Evol. Syst. , 45:397–419, 2014.[4] V. Deblauwe, N. Barbier, P. Couteron, O. Lejeune, and J. Bogaert. The global biogeography of semi-arid periodicvegetation patterns.
Global Ecol. and Biogeogr. , 17(6):715–723, 2008.[5] V. Deblauwe, P. Couteron, J. Bogaert, and N. Barbier. Determinants and dynamics of banded vegetation pattern migrationin arid climates.
Ecol. Monogr. , 82(1):3–21, 2012.[6] M. Drusch, U. Del Bello, S. Carlier, O. Colin, V. Fernandez, F. Gascon, B. Hoersch, C. Isola, P. Laberinti, P. Martimort,A. Meygret, F. Spoto, O. Sy, F. Marchese, and P. Bargellini. Sentinel-2: ESA’s optical high-resolution mission for GMESoperational services.
Remote Sensing of Environment , 120(Supplement C):25–36, 2012. The Sentinel Missions - NewOpportunities for Science.[7] R. Gelaro, W. McCarty, M. J. Su´arez, R. Todling, A. Molod, L. Takacs, C. A. Randles, A. Darmenov, M. G. Bosilovich,R. Reichle, et al. The modern-era retrospective analysis for research and applications, version 2 (MERRA-2).
J. Climate ,30(14):5419–5454, 2017.
8] W. A. Macfadyen. Vegetation patterns in the semi-desert plains of British Somaliland.
Geograph. J. , 116(4/6):199–211,1950.[9] J. E. G. W. Greenwood. The development of vegetation patterns in Somaliland Protectorate.
Geogr. J. , pages 465–473,1957.[10] C. F. Hemming. Vegetation arcs in Somaliland.
J. Ecol. , pages 57–67, 1965.[11] S. B. Boaler and C. A. H. Hodge. Observations on vegetation arcs in the northern region, Somali Republic.
J. Ecol. ,52(3):511–544, 1964.[12] A. Cornet, J.-P. Delhoume, and C. Montana. Dynamics of striped vegetation patterns and water balance in the ChihuahuanDesert.
Diversity and pattern in plant communities , pages 221–231, 1988.[13] D. L. Dunkerley. Infiltration rates and soil moisture in a groved mulga community near Alice Springs, arid central Australia:evidence for complex internal rainwater redistribution in a runoff–runon landscape.
J. Arid Environ. , 51(2):199–219, 2002.[14] A. Casenave and C. Valentin. A runoff capability classification system based on surface features criteria in semi-arid areasof West Africa.
J. Hydrol. , 130(1-4):231–249, 1992.[15] G. A. Worrall. The Butana grass patterns.
J. Soil Sci. , 10(1):34–53, 1959.[16] E. Meron.
Nonlinear physics of ecosystems . CRC Press, 2015.[17] E. Meron. Pattern formation–a missing link in the study of ecosystem response to environmental changes.
Math. Biosci. ,271:1–18, 2016.[18] E. Meron. From patterns to function in living systems: Dryland ecosystems as a case study.
Annu. Rev. Condens. MatterPhys. , 9:79–103, 2018.[19] Robbin Bastiaansen, Paul Carter, and Arjen Doelman. Stable planar vegetation stripe patterns on sloped terrain indryland ecosystems.
Nonlinearity , 32(8):2759, 2019.[20] R. Bastiaansen and A. Doelman. The dynamics of disappearing pulses in a singularly perturbed reaction–diffusion systemwith parameters that vary in time and space.
Physica D , 388:45–72, 2019.[21] R. Bastiaansen, A. Doelman, M. B. Eppinga, and M. Rietkerk. The effect of climate change on the resilience of ecosystemswith adaptive spatial pattern formation.
Ecol. Lett. , 2020.[22] E. Gilad, J. von Hardenberg, A. Provenzale, M. Shachak, and E. Meron. Ecosystem engineers: from pattern formation tohabitat creation.
Phys. Rev. Lett. , 93(9):098105, 2004.[23] M. Rietkerk, M. C. Boerlijst, F. van Langevelde, R. HilleRisLambers, J. van de Koppel, L. Kumar, H. H. T. Prins, andA. M. de Roos. Self-organization of vegetation in arid ecosystems.
Amer. Nat. , 160(4):524–530, 2002.[24] E. Istanbulluoglu and R. L. Bras. Vegetation-modulated landscape evolution: Effects of vegetation on landscape processes,drainage density, and topography.
J. Geophys. Res. , 110(F2), 2005.[25] S. E. Thompson, C. J. Harman, P. Heine, and G. G. Katul. Vegetation-infiltration relationships across climatic and soiltype gradients.
J. Geophys. Res. Biogeosci. , 115(G2), 2011.[26] S. Fatichi, V. Y. Ivanov, and E. Caporali. A mechanistic ecohydrological model to investigate complex interactions in coldand warm water-controlled environments: 1. theoretical framework and plot-scale analysis.
J. Adv. Model. Earth Sys. ,4(2), 2012.[27] A. Paschalis, G. G. Katul, S. Fatichi, G. Manoli, and P. Molnar. Matching ecohydrological processes and scales of bandedvegetation patterns in semiarid catchments.
Water Resour. Res. , 52(3):2259–2278, 2016.[28] C. A. Klausmeier. Regular and irregular patterns in semiarid vegetation.
Science , 284(5421):1826–1828, 1999.[29] N. Ursino. The influence of soil properties on the formation of unstable vegetation patterns on hillsides of semiaridcatchments.
Adv. Water Resour. , 28(9):956–963, 2005.[30] V. Guttal and C. Jayaprakash. Self-organization and productivity in semi-arid ecosystems: Implications of seasonality inrainfall.
J. Theor. Biol. , 248(3):490–500, 2007.[31] N. Ursino and S. Contarini. Stability of banded vegetation patterns under seasonal rainfall and limited soil moisturestorage capacity.
Adv. Water Resour. , 29(10):1556–1564, 2006.[32] A. Y. Kletter, J. Von Hardenberg, E. Meron, and A. Provenzale. Patterned vegetation and rainfall intermittency.
J.Theor. Biol. , 256(4):574–583, 2009.[33] M. Baudena, J. von Hardenberg, and A. Provenzale. Vegetation patterns and soil-atmosphere water fluxes in drylands.
Adv. Water Resour. , 53:131–138, 2013.[34] P. D’Odorico, F. Laio, and L. Ridolfi. Vegetation patterns induced by random climate fluctuations.
Geophys. Res. Lett. ,33(19), 2006.[35] K. Siteur, M. B. Eppinga, D. Karssenberg, M. Baudena, M. F. P. Bierkens, and M. Rietkerk. How will increases in rainfallintensity affect semiarid ecosystems?
Water Resour. Res. , 50(7):5980–6001, 2014.[36] J. Seghieri, S. Galle, J.-L. Rajot, and M. Ehrmann. Relationships between soil moisture and growth of herbaceous plantsin a natural vegetation mosaic in Niger.
J. Arid Environ. , 36(1):87–102, 1997.[37] N. Barbier, P. Couteron, R. Lefever, V. Deblauwe, and O. Lejeune. Spatial decoupling of facilitation and competition atthe origin of gapped vegetation patterns.
Ecology , 89(6):1521–1531, 2008.[38] S. Kinast, Y. R. Zelnik, G. Bel, and E. Meron. Interplay between turing mechanisms can increase pattern diversity.
Phys.Rev. Lett. , 112(7):078701, 2014.[39] N Ursino. Modeling banded vegetation patterns in semiarid regions: Interdependence between biomass growth rate andrelevant hydrological processes.
Water Resour. Res. , 43(4), 2007.[40] H. Vereecken, L. Weiherm¨uller, S. Assouline, J. ˇSim˘unek, A. Verhoef, M. Herbst, N. Archer, B. Mohanty, C. Montzka,J. Vanderborght, et al. Infiltration from the pedon to global grid scales: An overview and outlook for land surface modeling.
Vadose Zone J. , 18(1), 2019.[41] J. R. Rigby and A. Porporato. Simplified stochastic soil-moisture models: a look at infiltration.
Hydrol. Earth Syst. Sci. , J. Theor. Biol. , 244(4):680–691, 2007.[43] V. T. Chow. Open-channel hydraulics. In
Open-channel hydraulics . McGraw-Hill, 1959.[44] Sara Bonetti, Gabriele Manoli, Costantino Manes, Amilcare Porporato, and Gabriel G Katul. Mannings formula andstricklers scaling explained by a co-spectral budget model.
Journal of Fluid Mechanics , 812:1189–1212, 2017.[45] I. Rodr´ıguez-Iturbe and A. Porporato.
Ecohydrology of water-controlled ecosystems: soil moisture and plant dynamics .Cambridge University Press, 2007.[46] R. HilleRisLambers, M. Rietkerk, F. van den Bosch, H. H. T. Prins, and H. de Kroon. Vegetation pattern formation insemi-arid grazing systems.
Ecology , 82(1):50–61, 2001.[47] A. Mauchamp, S. Rambal, and J. Lepart. Simulating the dynamics of a vegetation mosaic: a spatialized functional model.
Ecol. Model. , 71(1-3):107–130, 1994.[48] P. Gandhi, S. Iams, S. Bonetti, and M. Silber. Vegetation pattern formation in drylands. arXiv preprint arXiv:1905.05115 ,2019.[49] J.D. Meiss.
Differentiable Dynamical Systesm . SIAM, 2017.[50] FH Busse. Nonlinear properties of thermal convection.
Rep.Prog. Phys. , 41(12):1929, 1978.[51] L. F. Shampine and M. W. Reichelt. The matlab ode suite.
SIAM J, Sci. Comput. , 18(1):1–22, 1997.[52] K. Gowda, S. Iams, and M. Silber. Signatures of human impact on self-organized vegetation in the horn of africa.
Sci.Rep. , 8(1):3622, 2018.[53] R. Bastiaansen, O. Ja¨ıbi, V. Deblauwe, M. B. Eppinga, K. Siteur, E. Siero, S. Mermoz, A. Bouvet, A. Doelman, andM. Rietkerk. Multistability of model and real dryland ecosystems through spatial self-organization.
Proc. Natl. Acad.Sci. , 115(44):11256–11261, 2018.[54] E. Meron. Vegetation pattern formation: The mechanisms behind the forms.
Physics Today , 72(11):30–37, 2019.[55] Y. R. Zelnik, E. Meron, and G. Bel. Gradual regime shifts in fairy circles.
Proc. Nat. Acad. Sci. , 112(40):12327–12331,2015.[56] O. Crompton, A. Sytsma, and S. Thompson. Emulation of the saint venant equations enables rapid and accurate predictionsof infiltration and overland flow velocity on spatially heterogeneous surfaces.
Water Resour. Res. , 55(8):7108–7129, 2019. Fast-Slow Model of Banded Vegetation Pattern Formation in Drylands:Supplemental Information
Punit Gandhi , Sara Bonetti , Sarah Iams , Amilcare Porporato , Mary Silber This supplement is divided into two parts: Appendix A provides additional simulation results in supportof the main text, and Appendix B provides details of the linear stability calculations presented in the maintext. In what follows, references labels that begin with “A” and “B” correspond to equations and figurescontained within this supplement. All other labels correspond to equations, figures, tables, etc. in the maintext.
Appendix A. Additional simulation results
We present additional simulation results to show more details of certain interesting aspects of the model,predictions and to support comments made in Sections 5 and 6 of the main text.
Appendix A.1. Subsurface water transport
We neglect subsurface water transport throughout the main text. In this section we show results fromlimited simulations indicating that including linear diffusion of soil moisture does not qualitatively changethe results. We fix the value of the diffusion coefficient associated with subsurface transport, which is notvery well constrained, and further parameter exploration is required.In order to model subsurface transport, we consider the addition of a linear diffusion term of the form δ s ∂ s∂x in the soil moisture equation (9a) of the slow system in the fast-slow switching model, Eqs. (7) and (9).The value of the coefficient δ s , like the biomass diffusion constant δ b in Eq. (9b), is largely unconstrained.We assume δ s = 10 δ b based on previous modeling efforts that take the ratio of these coefficients to bebetween 1 and 100 [Rietkerk et al. (2002) Amer. Nat. 160(4):524530; Gilad et al. (2004) Phys. Rev. Lett.93(9):098105].We carry out the simulation shown in Figs. 12 and 13 in which precipitation is slowly decreased and thenincreased with soil moisture diffusion, and the results are shown in Fig. A1. For reference, the dotted linesin this figure show the maximum and average biomass values of the patterned states from the simulationwith no soil moisture diffusion. There is no qualitative difference, and the simulations are only very slightlyqualitatively different. Appendix A.2. Random initial conditions and asymptotic behavior
Figure A2 illustrates the asymptotic behavior for two different random initial conditions on a one-kilometer spatial domain. The simulation shown in Fig. A2(a) is initialized with a 1% random noise on topof the predicted uniform vegetation state for constant precipitation of P = 160 mm/year , given by Eq. (4).The simulation shown in Fig. A2(b) is initialized with a 1% random noise on bare soil. Department of Mathematics and Applied Mathematics, Virginia Commonwealth University, Richmond, VA 23284, USA Department of Environmental Systems Science, ETH Z¨urich, 8092 Z¨urich, Switzerland John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA Department of Civil and Environmental Engineering and Princeton Environmental Studies, Princeton University, Princeton,NJ 08540, USA Department of Statistics and Committee on Computational and Applied Mathematics, University of Chicago, Chicago, IL60637, USA Figure A1: Including Soil Moisture diffusion. The solid (dotted) green lines indicate the average (maximum and minimum)biomass values of the uniform vegetation state. The red left arrows indicate the maximum biomass values of two-band patternsobtained by decreasing precipitation with steps of 2 mm/year . The blue right-arrows indicate the maximum biomass valuesfor single-band patterns for increasing precipitation. In both cases the dots indicate the (time and space) averaged biomassvalues, and the vertical lines indicate the range of biomass values within the pattern. The corresponding results from Fig. 13are indicated by the blue and red dotted lines. (a) (b)
Figure A2: spacetime diagram of biomass from simulations initialized with 1% random perturbation of (a) uniform vegetationstate (b) bare soil state. For scale, the spacetime simulation window of Fig. 8(b) is indicated by a dotted box.
At approximately year 1500 the pattern formed in Fig. A2(a) begins undergoing a series of band-mergingevents. At around 2500 years, the bands begin to exhibit a periodic breathing dynamic in which the widthof the bands oscillate as the trailing edge periodically moves downhill for a short time. The simulation inFig. A2(b) is initialized with the same perturbation of the bare soil state, given by Eq. (3). In this casethe band merging events are initiated sooner, but the transient still persists for around two millennia beforeeventually reaching a pattern with a single band on the domain.
Appendix A.3. Nonlinear dependence of transport on surface water height
For a majority of the results presented in the main text, we take the exponent δ = 1 in the surface watertransport speed term ν ( b, h ) = h δ − / (1 + ηb ) of Eq. (7a). While δ = 5 / k = 2 π/ m − for δ = 1and δ = 5 /
3. Figure A3 shows a comparison between δ = 1 and δ = 5 / δ = 5 /
3, thecontrast in surface water height between vegetation bands and bare soil regions is less pronounced at peakvalues, which occur during rain storms. The δ = 5 / t = 500. Despite these quantitative differences, the results are qualitatively similar.2
100 200 300 400 5000120 100 200 300 400 500050100 0 100 200 300 400 50000.5 0 100 200 300 400 50000.20.4 (a) (b)(c) (d)
Figure A3: Nonlinearity in surface water advection, Left (a,c): δ = 5 / δ = 1. The upper plots (a,b) showspatial profile of surface water height H , Infiltration rate I , soil moisture s , and biomass density B at t = 500 years . The solidline is the annually averaged profile while the dotted lines show the pointwise minimum and maximum values over the courseof the year. The lower plots (c,d) show spacetime diagrams of annually averaged biomass B in units of kg/m over the courseof the 500 year simulation.
100 200 300 400 50000.5
165 170 175 180 18500.10.20.30.40.5
Figure A4: Bistability of one-band patterns. The solid (dotted) green lines indicate the average (maximum and minimum)biomass values of the uniform vegetation state. The red left arrows indicate the maximum biomass values of one-band patternsobtained by decreasing precipitation with steps of 2 mm/year . The blue right-arrows indicate the maximum biomass valuesfor single-band patterns for increasing precipitation. In both cases the dots indicate the (time and space) averaged biomassvalues, and the vertical lines indicate the range of biomass values within the pattern. The solution profiles are obtained at175 mm/year during the precipitation scan up (boxed in blue) and the scan down (boxed in red).
100 200 300 400 5000240 100 200 300 400 500050100 0 100 200 300 400 50000.5 0 100 200 300 400 50000.5 0 100 200 300 400 5000240 100 200 300 400 500050100 0 100 200 300 400 50000.5 0 100 200 300 400 50000.5 (a) (b)(c) (d)
Figure A5: Rain Intensity, Left (a,c) 3-hour rain events and Right (b,d): 12-hour rain events. (a,b) Spatial profile of surfacewater height H , Infiltration rate I , soil moisture s , and biomass density B at t = 200 years . The solid line is the annuallyaveraged profile while the dotted lines show the pointwise minimum and maximum values over the course of the year. (c,d)Spacetime plot of annually averaged biomass B in units of kg/m over the course of the 200- year simulation. Appendix A.4. Dependence on mean annual precipitation, rain intensity and timing
In Fig. 13 of the main text, the single-band solutions undergo a jump in terms of maximum biomassvalue at around 177 mm/year as precipitation is increased. Scanning down in mean annual precipitationwith a single-band solution reveals a hysteresis loop associated with this jump. Figure A4 shows biomassprofiles of the two distinct single-band solutions obtained at 175 mm/year . The solution found by decreasingprecipitation has a lower peak biomass value, a wider vegetation band, and has a bump in biomass at thetrailing edge.The simulation results presented in the Section 5 use a fixed rain input scheme of two six-hour eventsspaced 6 months apart, inspired by rainfall data shown in Fig. 1. In order to explore the influence of theintensity of the rain events, we consider a fixed mean annual precipitation of 160 mm/year and vary thelength of the biannual rain events. Figure A5 shows spacetime diagrams of biomass and spatial profiles forthe case of 3 hour and 12 hour rainstorms.Linear stability analysis for the coupled model, Eq. (1), shown in Fig. 3(b) indicates a drastic shrinkingof the precipitation range over which the uniform state is unstable to patterns when a single 12-hour rainevent per year is used instead of two equally spaced 6-hour rain events per year. Figure A6 shows thatsimulations of the fast-slow switching model, are consistent with these predictions. With one rain event peryear, patterns are found from 89 to 147 mm/year with precipitation scans. This is a significantly narrowerrange than the 34 to 201 mm/year that is found with two rain events per year.In order to explore the influence of moving towards seasonal rain input, we consider a case in whichtwelve rain events of 30 minutes each are evenly distributed over a rainy seasons that occurs every 6 months.5
Figure A6: One rain event per year. The solid (dotted) green lines indicate the average (maximum and minimum) biomassvalues of the uniform vegetation state. The red left arrows indicate the maximum biomass values of two- and one-band patternsobtained by decreasing precipitation with steps of 2 mm/year . The blue right-arrows indicate the maximum biomass valuesfor single-band patterns for increasing precipitation. In both cases the dots indicate the (time and space) averaged biomassvalues, and the vertical lines indicate the range of biomass values within the pattern.
Initially each rainy season spans the entire 6-month period, but the length is gradually decreased in stepsof 12 hours every decade until the 12-event season occurs over a period of 12 days. The first row of Fig. A7shows the initial and final patterns of instantaneous rainfall rates, from right to left. We note that, for the12 day rainy season, there are still 12 distinct rain events of 30 minutes even though they are not resolvedin the figure. For the simulation, we take the parameters given in Table 2 and a MAP of 160 mm/year .The second row shows the average biomass profiles (solid) along with the pointwise minimum and maximum(dotted) values. As the rainfall becomes more seasonal, the vegetation bands become wider with lower peakbiomass value and larger fluctuations. The last two rows of Fig. A7 show the peak and averaged biomass,and the fraction of the domain covered by biomass as a function of the length of the rainy seasons. Whilethe more seasonal rain pattern leads to a larger fraction of the domain covered by biomass, the less seasonalrain leads to more biomass on the domain. The peak biomass and average biomass values with a 12-dayrainy season are approximately 0.36 and 0.79 times the corresponding values with a 6-month rainy season.
Appendix B. Linear stability of the uniformly vegetated state in the three-field coupledtimescale model
Here we summarize the equations used to compute the linear stability of the spatially uniform vegetationsolution ( h ( t ) , s ( t ) , b ( t )) of the non-dimensionalized version of the three-field coupled timescale model (6),in the case of temporally periodic precipitation p ( t ) = p ( t + τ ), where τ is the (minimal) period associatedwith the rain events. The uniform solution satisfies the following system of ordinary differential equations:˙ h = p ( t ) − ι ( h , s , b )˙ s = α ι ( b , s , h ) − (cid:15) ( σ + γb ) s (B.1)˙ b = (cid:15) ( s (1 − b ) − µ ) b , where, recall, ι ( h, s, b ) ≡ (cid:16) b + qfb + q (cid:17)(cid:16) hh + 1 (cid:17) (1 − s ) β . (B.2)The linear variational equations, associated with a perturbation of the form e ikx are obtained by insertingthe ansatz h ( t, x ) = h ( t ) + h ,k ( t ) e ikx ,s ( t, x ) = s ( t ) + s ,k ( t ) e ikx ,b ( t, x ) = b ( t ) + b ,k ( t ) e ikx ,
100 200 30002040
20 40 60 80 100 120 140 160 18000.20.40.60.811.2
20 40 60 80 100 120 140 160 18000.20.40.6