From nonequilibrium initial conditions to steady dryland vegetation patterns: How trajectories matter

The multiscale nature of ecohydrological processes and feedbacks implies that vegetation patterns arising in water‐limited systems are directly linked to water redistribution processes occurring at much shorter timescales than vegetation growth. This in turn suggests that the initially available water in the system can play a role in determining the trajectory of the system, together with the well‐known role of the rainfall gradient. This work explores the role of initial hydrological conditions on vegetation dynamics and vegetation patterns. To do so, the HilleRisLambers–Rietkerk model was solved with different rainfall amounts and a large range of initial hydrological conditions spanning from near‐equilibrium to far‐from‐equilibrium conditions. The resulting vegetation patterns and ecohydrological signatures were quantitatively studied. The results show that not only do initial hydrological conditions play a role in the ecohydrological dynamics but also they can play a dominating one even resulting in divergent vegetation patterns that exhibit convergent mean‐field properties, including a new set of hybrid patterns. Our results highlight the relevance of assessing both global ecological and hydrological signatures and quantitatively assessing patterns to describe and understand system dynamics and in particular to determine if the systems are transient or steady. Furthermore, our analysis shows that the trajectories the system follows during its transient stages cannot be neglected to understand complex dependencies of the long‐term steady state to environmental factors and drivers.

There is a growing awareness that the long-term asymptotically steady states of ecological systems are often not representative of the short-term dynamics (Hastings, 2004;2010;Larsen et al., 2016;Arnoldi, Bideault, Loreau, & Haegeman, 2018) and relevant classical concepts (stability and resilience) require a better understanding of transient dynamics (Arnoldi, Loreau, & Haegeman, 2016;Siteur, Eppinga, Doelman, Siero, & Rietkerk, 2016), as well as assessing when the system can be considered steady (Hastings, 2004;Hamann, Schmickl, & Crailsheim, 2012). Furthermore, the role of initial conditions on the transient dynamics may not be captured by the asymptotically steady regime (Arnoldi et al., 2018). One example of this is that how near or far the initial conditions of a system are from equilibrium has implications on the rates and timescales of the dynamics and timescales, a topic yet to be thoroughly addressed (Hastings, 2016;Larsen et al., 2016). Additionally, dryland ecosystems frequently experiment far-from-equilibrium conditions thus staying in a quasipermanent transient condition and exhibit nonlinear responses to boundary conditions and forcing (Puigdefábregas, 1998). That is, the systems of interest here are prone to be strongly affected by initial conditions, which may also lead to long, complex transient trajectories.
The stability of these ecohydrosystems to perturbations and their regions of bistability (Zeng, Shen, Zeng, & Dickinson, 2004;Borgogno et al., 2009) has been a focus of research. This has lead to simulating the system evolution from arbitrary initial conditions into steady states (most often tacitly) defined when mean-field signatures-a quantity characterizing a complex behaviour (Addor et al., 2018)-such as biomass, converge to steady conditions (e.g., Baudena, von Hardenberg, & Provenzale, 2013;Kletter, von Hardenberg, Meron, & Provenzale, 2009;Thompson, Katul, & McMahon, 2008). It must be acknowledged that alternative analytical and mathematical approaches, not requiring an ad hoc definition of a steady state or equilibrium, have also been used (Lejeune, Tlidi, & Lefever, 2004;van der Stelt, Doelman, Hek, & Rademacher, 2012;Zelnik et al., 2013). Many studies have focused on patterns at steady state and their response or dependence to environmental conditions or system properties (e.g., Rietkerk et al., 2004;, and others have been concerned also with monitoring domain-wide integrated quantities such as total biomass and vegetation cover (e.g., Kletter et al., 2009;Pueyo, Kefi, Diaz-Sierra, Alados, & Rietkerk, 2010). The formulation of these models allows to achieve steady states only in the presence of steady rainfall forcing (Kletter et al., 2009). Natural time-varying rainfall is unlikely to force the system into steady states (Guttal & Jayaprakash, 2007;Kletter et al., 2009). Nonetheless, using idealized set-ups with steady rainfall, many authors have analysed the transient dynamics resulting from these inherently transient partial differential equation (PDE) models and drawn insights on different issues, such as the rates of transitions (Yizhaq et al., 2014), transitions of patterns on the basis of rainfall (Roitberg & Shoshany, 2017;Zelnik & Meron, 2018), response to gradual environmental change (Dagbovie & Sherratt, 2014;, succession dynamics (Pueyo et al., 2010), the pseudosteady migration of vegetation bands (Rietkerk et al., 2002;Sherratt & Lord, 2007), colonization of a slope (Sherratt, 2016), or the relevance of the transient history of patterns (Sherratt, 2013). This has further been studied using the concepts of pattern selection and wavelength selection (Dagbovie & Sherratt, 2014;Siero et al., 2015;Zelnik & Tzuk, 2017) used to describe the process under which a given pattern (with a particular characteristic wavelength) is established. Pattern (wavelength) selection has been shown to be history dependent (Dagbovie & Sherratt, 2014;Zelnik & Tzuk, 2017).
From the model perspective, its time evolving nature requires the definition of both the initial vegetation distribution and the initial surface and subsurface water distributions, henceforth initial hydrological conditions (IHCs). Random vegetation distributions have often been selected as initial conditions throughout the literature. Some researchers have reported little qualitative sensitivity to different random distributions (Guttal & Jayaprakash, 2007;Jenerette et al., 2011;Thompson & Katul, 2009). In contrast, Meron, Gilad, von Hardenberg, Shachak, and Zarmi (2004) initialized simulations with random perturbations over unstable uniform vegetated states, showing that different patterns can arise while using the von Hardenberg, Meron, Shachak, and Zarmi (2001) model. Zelnik et al. (2013) showed that different initial vegetation disturbances can lead to different patterns under the same environmental conditions. Similarly, in a two-species context, Baudena and Rietkerk (2012) showed that different two-species patterns can be obtained from different initial species distributions.
The role of IHC on asymptotic steady-state patterns has received less attention. IHCs have been mostly taken as the so-called plantless-equilibrium conditions (e.g., Rietkerk et al., 2002), that is, setting the IHC as if the domain is homogeneously bare. This has been proven useful when exploring stability properties of such systems , although the assumption of near-equilibrium IHC can be viewed as untested and has been increasingly challenged (Vetter, 2005;Mueller et al., 2014). A few notable studies are the exception. Meron et al. (2004) obtained similar patterns from random perturbations over the uniformly vegetated state, or from isolated vegetation spots, while experiencing different transient states including ring-like patterns. Sherratt and Lord (2007) showed that properties of vegetation bands can vary in response to small random variations of initial biomass and initial water. Both studies, relying on randomized small perturbations, concluded that the spatial patterns can depend on initial conditions, including IHC, which is most relevant for this work.
In contrast, Guttal and Jayaprakash (2007) used both well-developed ecohydrological patterns and zero initial water depths and found that the final patterns did not depend on initial conditions. The apparent contradiction between the results of Meron et al. (2004) and Sherratt and Lord (2007) compared with those of Guttal and Jayaprakash (2007) prompts the systematic study of the possible effects of IHC, and the dependency on how specific IHC could generate specific patterns or trajectories remains a gap that motivates this work.
Mathematically, initial conditions are the required, known, and prescribed initial states. In natural systems, initial conditions may be impossible to identify or retrieve (Sherratt, 2015), although systems that have experienced large natural or anthropogenic disturbances-e.g., overgrazing, mining, agriculture, fires, and volcanic eruptions-may allow to identify characteristic initial conditions prior to a (natural or assisted) restoration process. Initial states in such cases may strongly constrain or determine the long-term state of complex natural systems (Perron & Fagherazzi, 2012;Maurer & Gerke, 2016).
IHCs can also be viewed as the result of intermittent rainfall or a single (pulsed) event (Schwinning & Sala, 2004). In the context of semiarid vegetation patterns, it has been seldom asked what is the nature of the pulsed events that may be responsible for initializing systems (Eddy, Humphreys, Hart, Mitchell, & Fanning, 1999) or the transients triggered after individual events . We argue that IHC deviating from ecohydrological equilibrium can be idealized as singular hydrological disturbances, and studying the transient response to it may lead to insights of the transient dynamics and equilibrium states, without requiring complex and stochastic rainfall signals.
The transient trajectories from the initial conditions to the converged steady state have received less attention due to a wide range of reasons, starting with the inherent difficulties of observing such transient processes in nature and difficulties in model parametrization to allow meaningful comparisons among others (Karssenberg, Bierkens, & Rietkerk, 2017). Modelling efforts to address transient dynamics do exist, and the importance of transients has been highlighted, relating to many issues. Observed patterns have been reported not to coincide with simulated steady-state patterns (Meron, 2016). Ring and spiral patterns have been identified as transient (Bonanomi et al., 2014;Fernandez-Oto, Escaff, & Cisternas, 2019;Sheffer, Yizhaq, Shachak, & Meron, 2011;Tlidi et al., 2018;Yizhaq et al., 2019), likely triggered by excitable behaviour (Fernandez-Oto et al., 2019;Tlidi et al., 2018).
Discrete rainfall events have been observed to produce transients . The temporal variability of rainfall has been pointed as partially responsible for the rates at which the systems react to environmental change (Zelnik et al., 2013), and the transient behaviour of the systems interact with the rates of variation of external drivers . Ecological succession has been explained through modelling ecohydrological interactions (Pueyo et al., 2010). Importantly, patterns have been shown to be history dependent (Sherratt, 2013;, and steady patterns can be used to infer the historical origin of the system (Sherratt, 2015). Additionally, the so-called hybrid patterns (Meron, 2012;Zelnik, Meron, & Bel, 2015)-that is, stable patterns characterized by anomalies in pattern periodicity-have been identified arising from different system trajectories and also strongly affect the response trajectory of the system to perturbations. Mueller et al. (2014) include the dependency on previous states as a central element of the complexity of these ecohydrosystems.
In light of this, we hypothesize that initial water availability can have a long-lasting effect and contribute to determine the ecohydrological steady state (typically defined as a steady biomass signal). We propose that the extent to which such initial water availability can play a significant role on the long-term steady state can be explained by the evolution trajectory of the system. Therefore, the following research questions are proposed: 1. Which role does initial water availability play on determining the trajectory of the system and the long-term ecohydrological steady state?
2. Which observable variables and behaviours of the system, that is, global signatures (biomass, cover, and hydrological balance) and spatial patterns and quantitative indicators of their properties, can be affected by different IHCs?
We address these questions, performing a numerical sensitivity study along a rainfall gradient and an initial surface water gradient.

Mathematical ecohydrological model
The spatially distributed ecohydrological model is based on two hydraulic submodels (surface and subsurface flow) and a biomass submodel, as proposed by HilleRisLambers, Rietkerk, van den Bosch, Prins, and de Kroon (2001) and Rietkerk et al. (2002). The system is comprised by three reaction-diffusion PDEs and has been widely used in literature for similar studies (e.g., Roitberg & Shoshany, 2017;Sherratt, 2016;Thompson et al., 2008;Yizhaq et al., 2014).
where b is biomass density (g m −2 ), h is surface water depth (mm), W is subsurface water depth (mm), x are the two-dimensional spatial coordinates (m), t is time (day), r is rainfall intensity (mm day −1 ). The infiltration rate i (mm day −1 ) is a function of biomass in order to model infiltration enhancement due to vegetation.
The rate at which vegetation uptakes water from the subsurface  is All model parameters are defined in Table 1, along with their selected values for this study. It is noteworthy that Systems (1)- (3) are a set of three diffusion-reaction equations: all water and biomass propagate horizontally by linear diffusion processes. Biomass growth is modelled as linearly dependent on uptake, and biomass death is a first-order decay term. Infiltration couples both water systems directly, and indirectly to the biomass field, whereas uptake directly couples the subsurface water field to biomass.
Systems (1)-(3) are integrated in space using a first-order finite volume approach and in time with a first-order, explicit, and forward Euler method. Because all three PDEs in Systems (1)-(3) are linear diffusion equations, this numerical approximation is conditionally stable and must satisfy the time-step size restriction: The solver was implemented in C, with OpenMP-based parallelization, and run in parallel on 20 central processing units i7-6950X at 3.00 GHz.
Here, the simulation duration is set significantly larger than most of the reported simulations to ensure that the final states are steady or nearly steady.
Such IHC is in fact the steady bare soil solution of the system, that is, the water depth h and subsurface water depth W fields are set equal to their steady states solutions when biomass is 0. This can be shown to be IHCs are a function of rain and satisfy that the infiltration-rain ratio I∕R = (rA) −1 ∫ A i(x, y) dA is equal to 1. The reader should note that, despite that this condition is termed plantless, it does not mean that the initial vegetation condition is set to b(x, t = 0) = 0. In fact, initial biomass conditions are set in a similar fashion as previous studies (e.g., Pueyo et al., 2010) as randomly located biomass peaks of biomass where  is a set of randomly determined areas (computational cells), is the vegetated area and A is the total domain area. Initial biomass density was set to b o = 90 g m −2 , and initial vegetation cover was set to V = 0.01. All simulations use the same initial biomass distribution. The reasoning for favouring varying surface IHC over subsurface IHC is that the system's evolution is expected to be less sensitive to the latter. There are two reasons for this: the infiltration feedback has a higher interaction with surface flow, and the first-order decay percolation/evaporation sink term allows for a net loss of water directly from the subsurface, which strongly pushes the system back into equilibrium (refer to the supporting information for some illustrative results on this). For every rainfall scenario, twelve different surface IHC are tested (Table 2). Boldface values in the table correspond to plantless equilibrium IHC. Column colours in Table 2 describe the initial surface water depths, and row colours describe precipitation. These colours are systematically used throughout the paper. Nonzero surface IHCs must be interpreted within the context of a closed system (periodic boundaries), from which water cannot escape. It is therefore possible to conceptually interpret the initial water depths as the accumulated rainfall of events occurring before the (arbitrary) initial time. Although this study is mainly concerned with surface IHC, illustrative examples of the effects of different initial subsurface water and initial biomass b o are included as a supporting information, showing that the sensitivities to either subsurface and biomass initial conditions are much lower than for surface initial conditions. In consequence, henceforth, only surface IHC variations are explored.

Postprocessing and indicators
The results are assessed in of terms mean-field (i.e., integrated in the whole domain) ecological signatures (total biomass [B], and vegetation cover, and mean density) and mean-field hydrological signatures-for example, infiltration/rain ratio I∕R and surface water fraction-together with the spatial patterns themselves and quantitative statistical descriptors of the spatial patterns (number of patches, patch area, area-perimeter ratio, patch circularity, and patch separation). These metrics do not require any base assumption on the shapes or periodicity of patterns. We purposely avoid using metrics that assume periodicity (e.g., wavelength and frequency), which are often used to describe periodic patterns (e.g., Deblauwe, Couteron, Lejeune,

Description Parameters Value Units
Water-biomass conversion rate c b 10 g m −2 /mm Maximum uptake per unit biomass g b 0.05 mm /(g day m 2 ) Half-saturation water-uptake constant k b 5 m m Death rate coefficient  porting information for details). The sensitivity of the final states to rainfall and initial water availability is studied, as well as the temporal evolution of selected representative cases.

RESULTS AND DISCUSSION
3.1 Transient results: The 0.75-mm day −1 case Figure 1 shows, for the 0.75-mm d −1 rainfall (as an example), the evolution of the key domain-wise integrated (mean-field) ecohydrological signatures. Colours in Figure 1 follow Table 2 and represent each IHC.
Biomass evolution (Figure 1a) for all IHCs asymptotically converge to the same final biomass, except for the cases that result in zero biomass (bare states). Total biomass is nearly steady from the first year onward and is clearly steady by the end of simulations with vegetated final states. Biomass evolution is clearly different early on (t ≲ 3 years), with trajectories describing overshooting and undershooting relative to the final biomass, after which the trajectories quickly converge. Biomass dynamics can be argued to be (too) fast and comparable with the rates of water dynamics-that is, large changes in the first year. It also seems implausible that natural systems might achieve a steady state in such short timescales. This is nonetheless ubiquitous in the VSO modelling literature and has been recognized as a modelling limitation (Guttal & Jayaprakash, 2007;Thompson et al., 2008), arguably related to the simplicity of both the model and the selected parametrization (following common parameters in the literature).
The steady forcing and periodic boundaries, all cases necessarily converge to the same I∕R = 1 global equilibrium steady state, despite being highly variable at the beginning. The reference case of h 0 = 18.75 mm (in black) is the least variable one, showing that the system is initially-and remains always-close to global hydrological equilibrium.
Complementarily, more water is accumulated on the surface in bare cases (Figure 1c), which shows that there is more surface water (5%) in the bare steady cases than in vegetated cases (0.02%). Figure 1 also shows that the infiltration ratio, and the surface water ratio vary at significantly different rates across IHC, and that the initial period of IHC dominated water availability (in which the systems is not forced exclusively by rain) lasts roughly one tenth of the time that it is required for the system to restabilise.
This suggests that IHC can have a long period of sustained influence, until eventually the global signatures of the system become independent of it. Importantly, it must be noted that the time at which steady conditions can be said to exist, depends on the observed signature. Figure 1 shows that signatures describe steady behaviour at the end    Figure 2, suggesting that the excitable behaviour leads to a hybrid pattern that is stable under r = 0.75 mm day −1 (as bands/spots are stable for such rainfall). However, for a lower rainfall r = 0.60 mm day −1 (Figure S3), rings only appear for h 0 =60 mm, and the final pattern is one of strongly periodic circular spots. This suggests that, under a lower rainfall, the ring-generating excitable behaviour then leads into an unstable hybrid pattern that breaks back to a stable classical spotted pattern. That is, excitable behaviour will only have an impact on the long-term asymptotically steady state if the hybrid states are stable.
In order to better assess the pattern trajectories, Figure 3 shows Altogether, from the evolution of the system under different IHC as described by Figures 1, 2, and 3 lead to some insights.  (see supporting information). Although, without explicitly computing for longer times, it cannot be empirically ruled out that an extremely long evolution of the pattern will lead to spotted patterns; the point holds in that long term, stable hybrid patterns are a consequence of certain IHCs. The stable hybrid patterns are likely to be related to homoclinic snaking that results in hybrid patterns (Zelnik et al., 2013), contrary to ring patterns that are inherently transient and not a form of homoclinic snaking (Sheffer et al., 2011). 4. Ring patterns can emerge as part of far-from-equilibrium trajectories triggered by nonequilibrium IHC, extending upon previously identified possible causes such as parametrization  and rainfall (Sheffer et al., 2011). Furthermore, we show transient ring patterns emerging in the Rietkerk model, without a root-augmentation feedback present in previously reported ring structures (e.g., Sheffer et al., 2011;Meron, 2016) or other nonhydraulic nonlocal effects (Sheffer et al., 2011;Tlidi et al., 2018;Yizhaq et al., 2019). Nonequilibrium IHC may be viewed as a single pulsed event and therefore these results support that rainfall variability favours ring formation, through a surface water redistribution mechanism (Sheffer et al., 2011;Yizhaq et al., 2019). These observations support and extend the current modelling-based understanding of ring patterns, as (to the authors' knowledge) the transient nature of rings has still not been clearly shown empirically in field observations.

Final state indicators and signatures for all cases
In this section, the final state results for all rainfalls are summarized and assessed. Figure 4 shows For the 0.85-mm day −1 case, it exhibits hybrid patterns of spirals and spots towards the extremes of the IHC gradient. For 1.0 mm day −1 , the patterns also differ, but the pattern structure remains much more similar. The symmetry along the IHC gradient also manifests in that bare steady states can be obtained for both very dry low and very wet IHC, that is, the relation between IHC and long-term productivity is nonmonotone. This is a novel observation in VSO, which is the result of the different system trajectories spawning from far-from-equilibrium IHC. The patterns vary along the rainfall gradient consistently with the literature along the equilibrium cases (spots to labyrinths with increasing rainfall). Nonetheless, for nonequilibrium cases (e.g., 1/16), because of the hybrid patterns, the well-established variation following rainfall does not hold. Additionally, note that ring structures are not present in the final patterns in Figure 4, although they are present in the transient states (e.g., Figures 2, S3, and S4).  (1/16) and the 300 mm (16) cases is similar. These two hybrid patterns also behave similarly in terms of area and ratio.

SUMMARY AND CONCLUSIONS
The results of the study allow the formulation of the following conclusions regarding the posed research questions, as well as additional insights and implications:  (Sherratt, 2013(Sherratt, , 2016 (Hastings, 2010;Konings, Dekker, Rietkerk, & Katul, 2011;Zelnik et al., 2013;Karssenberg et al., 2017). If state transitions are to be properly understood, and their rates and characteristic times well represented, this study suggests that both a better assessment and treatment of IHC is necessary and that the transient dynamics still require significant study. Theoretical reasoning and observations in broader ecohydrological settings, which also highlight the need for better transient understanding and assessing nonequilibrium dynamics (Larsen et al., 2016;Ratajczak et al., 2017), are concurrent with this conclusion.
2. Spatial features of the patterns can be dependent on initial hydrological conditions. However mean-field signatures such as total biomass yield, vegetation cover, and hydrological partitioning still converge either to bare states or to the same vegetated steady state for different IHCs for a given rainfall. This implies that mean-field estimations of biomass and cover on the basis of total rainfall are history independent and that they are robust relative to the IHC. However, total rainfall alone is not a robust predictor of the spatial distribution of such biomass.
3. The study also raises interesting questions on how to define or assess if an ecohydrological system (in the semiarid VSO context) has reached a steady state. Mean-field (global) signatures can achieve steady states quickly in contrast to spatial patterns and their indicators, which can exhibit transient behaviour long after the mean-field indicators suggest a steady state. This implies that mean-field signatures are insufficient to assess steadiness, and spatial patterns and indicators are a necessary complement, complementing recent insights arising in broader studies (Bastiaansen et al., 2018). How relevant this is for particular natural ecohydrosystems is difficult to assess given that-as eloquently stated by Puigdefábregas (1998) (Dagbovie & Sherratt, 2014;Zelnik et al., 2013), amplitudes and frequencies of Fourier decompositions (Deblauwe et al., 2011;Yizhaq & Bel, 2016), and the ring index (Sheffer et al., 2011). A complete study of how these metrics behave is out of the scope of this work but may provide further insights.
The main limitation of our study is that the conclusions are model based and cannot be blindly extrapolated into natural reality. The main reason is that the resulting temporal scales and rates are too short and too fast compared with natural systems for several reasons. First, the interpretation of what the IHC mean in natural systems is difficult and which IHC actually can be expected is not well defined. Second, some sources of complexity have not been considered (e.g., rainfall variability and heterogeneity). Third, the accepted parametrization has only accounted for steady state, qualitative, pattern-matching validation (Bastiaansen et al., 2018;Bokulich, 2014), therefore not guaranteeing correctly reproducing rates. It may be that alternative parametrizations, accounting for different physiological rates of vegetation, may affect these results (see supporting information for a minimal example).
It remains open to explore if the trajectories observed in this work can be generalized to other parametrizations or under which conditions these behaviours are likely to occur or not. Nonetheless, this study sheds some light on these issues. The interpretation of the systems' response to an arbitrary IHC can be understood as the transient response of the ecohydrosystem to a discrete hydrometeorological event. Therefore, the magnitude of rainfall events, which may trigger colonization processes in natural systems, may determine the resulting pattern to the same extent-or even more-than annual rainfall.
Semiarid regions do not exhibit constant, steady precipitation but rather a highly pulsed and intermittent rainfall regime (Modarres & de Paulo Rodrigues da Silva, 2007;Dunkerley, 2018;Pilgrim, Chapman, & Doran, 1988). This work points out that each of those pulsed events may act as a new initial condition, which triggers a new transient process whose final steady state may be determined by the magnitude of the event together with yearly rainfall, as suggested by Ratajczak et al. (2017) and Schwinning and Sala (2004). In light of this, the challenge of accurately modelling water redistribution in response to such pulsed events, and analysing the resulting transient signals, is all the more relevant. This may prove essential to reproducing rates and timescales of natural responses to hydrometeorological events, which are not as dramatic as the modelling results suggest.
An interesting question relating to the resilience of systems with patterns rising from nonequilibrium and far-from-equilibrium initial conditions is hinted by this study. The new set of hybrid patterns, qualitatively and quantitatively different from those spawning from near-equilibrium trajectories, allows to ask if their stability (in terms of resilience and resistance) to environmental perturbations (e.g., drought and climate change) or disturbances (e.g., grazing and fires) might be different from patterns generated by near-equilibrium initial conditions for which stability and resilience have been extensively studied along rainfall gradients. The literature shows that resilience is sensitive to the spatial structure of vegetation (Rietkerk et al., 2004;von Hardenberg et al., 2001;Baudena et al., 2013;Realpe-Gomez et al., 2013); particular patterns, for example, spots, have also been identified to favour facilitation over other patterns that favour resilience, for example, bands , and nonperiodic patterns induced by heterogeneity have been noted to be more resilient (Yizhaq & Bel, 2016). The diversity of patterns reported warrants to explore their stability, as it may behave differently from periodic patterns, and also suggest that not only the spatial structure at any given time may affect resilience but the instantaneous trajectory and the ongoing rate of change at the time of disturbance may also play a role. Alternatively, by not assuming that VSO ecohydrosystems are steady, it becomes natural to ask if particular trajectories are more or less prone to be affected by environmental changes.