The Impact of Spatial and Temporal Dimensions of Disturbances on Ecosystem Stability

Ecosystems constantly face disturbances which vary in their spatial and temporal features, yet little is known on how these features affect ecosystem recovery and persistence, i.e., ecosystem stability. We address this issue by considering three ecosystem models with different local dynamics, and ask how their stability properties depend on the spatial and temporal properties of disturbances. We measure the spatial dimension of disturbances by their spatial extent while controlling for their overall strength, and their temporal dimension by the average frequency of random disturbance events. Our models show that the return to equilibrium following a disturbance depends strongly on the disturbance’s extent, due to rescue effects mediated by dispersal. We then reveal a direct relation between the temporal variability caused by repeated disturbances and the recovery from an isolated disturbance event. Although this could suggest a trivial dependency of ecosystem response on disturbance frequency, we find that this is true only up to a frequency threshold, which depends on both the disturbance spatial features and the ecosystem dynamics. Beyond this threshold the response changes qualitatively, displaying spatial clusters of disturbed regions, causing an increase in variability, and even a system-wide collapse for ecosystems with alternative stable states. Thus, spanning the spatial dimension of disturbances is a way to probe the underlying dynamics of an ecosystem. Furthermore, considering spatial and temporal dimensions of disturbances in conjunction is necessary to predict ecosystem responses with dramatic ecological consequences, such as regime shifts or population extinction.


the AE model
To understand the shape of the return time curve T (σ) for the AE model, as seen in bottom right panel of Figure 2, we consider the initial response of the system just after the disturbance has occurred. If we ignore the effect of diffusion, then the domain that was not disturbed does not react at all, while the disturbed domain either rebounds back to the vegetated state N = K, or falls further to the bare state N = 0. This depends simply on whether the current level of biomass N 0 (immediately after the disturbance) is higher or lower than α. If we note the system size as L, the spatial extent of the disturbance as σ, and the disturbance strength (percentage of biomass cut) as s, then the biomass value at the disturbed region just after the disturbance occurred is: We therefore have a critical value of σ when N 0 = α, which will decide whether the disturbed region would rebound without the effect of diffusion, as: As shown in Figure S1 with a vertical red line, this approximation gives a good indication to where the return time changes behaviour rapidly. For larger disturbances σ > σ c , each single point would recover from the disturbance on its own accord. Moreover, the contribution of spatial processes to this recovery would be negligible if the rate of isolated recovery r iso is much faster than the rate of rescue recovery r res , namely r iso r res . We can use simple dimensional considerations (Meron, 2015) to estimate what are the conditions for this situation.
First, we can estimate r iso by r, since it is the only parameter with units of time and no relation to space. This estimation will break down close to the unstable equilibrium (since the recovery rate will go to zero), i.e. when σ is close to σ c . To estimate r res we need to consider how fast a front between the disturbed and undisturbed domains moves. In general, the front's speed u between two domains in a model that can be written as N t = rf (N ) + d∇ 2 N would be proportional to both the rate of local dynamics r and the rate of diffusion d as u ∼ √ rd, due to dimensional considerations. The area of the disturbed region is σL, so that the overall time to recovery is the ratio between area and speed, τ ∼ σL/ √ rd. We can thus estimate the rescue recovery rate as r res ∼ √ rd/(σL), and compare this to r iso . This gives us the condition σ r/d/L, which is easily met in the parameters we choose since L is large enough. Therefore to approximate the rate of recovery for σ > σ c we can just use the return time when we ignore spatial effects (d = 0). Indeed, this estimation, shown by the right-more blue curve in Figure S1, works well, even for lower values of σ close to the critical value σ c .
For the opposite regime, σ < σ c we can easily assume that isolated recovery is negligible, since in the bistable AE model isolated recovery is not possible when locally N < α. To estimate the rescue recovery we again consider the speed of a front between the disturbed and undisturbed domains. Here, both of these can be assumed to be the two stable states, and under these conditions the front speed is a constant that depends only on the model parameters (Meron, 2015). In particular, the front speed u grows with the rate of local dynamics r and diffusion coefficient d as u ∼ √ rd, and for the AE model considered here, u also grows with α once α > K/2 (the front is stationary at α = K/2 due to the symmetry between the two stable states for these parameters). If we assume that once a disturbance occurs, the time until the disturbed region collapses to the alternative state is negligible, then the recovery time T in simply the time it takes the fronts to take over the disturbed region: This approximation is shown in the left-more blue curve in Figure S1. While there is a clear under-estimation of the return time in this case due to complex front dynamics, the overall trend and the value for very localized disturbances is well approximated.

Appendix B: Relation between return time and variability
In this appendix we derive a general relation between recovery dynamics following a single disturbance, and variability under repeated disturbances. In this manuscript we formalize disturbance events as a realization of a point process. We therefore begin by explicitly defining a point process that we use for our derivation as well as the assumptions we make for this derivation. We then detail the derivation itself, and finally we discuss the main results of this derivation.

Definitions
A point process (p.p.) Φ is a random sampling of points in R d such that the number of points sampled in a bounded set is always finite (Baccelli and Blaszczyszyn, 2010). We note Φ(A) as the restriction of a p.p. to a given set A and we note φ(A) = {x 1 , x 2 , ...} ⊂ A as any of its realizations. Let N (A) be the random variable representing the number of points in A sampled by Φ. The intensity measure Λ of Φ counts the expected number of point sampled from a given set, i.e. Λ(A) = E Φ N (A). Note that, if A and B are disjoint sets, N (A) and N (B) are independent random variables so that Cov Φ (N (A), N (B)) = 0. We call a Poisson process a p.p. such that, for any A, the mean number of points sampled is also equal to the variance so that E Φ N (A) = Var Φ N (A) = Λ(A). A homogeneous Poisson process satisfies Λ(A) = f |A| where |A| is the usual volume of A. In this case, this simply means that Φ(A) is a uniform sampling of points of A, with an expected number of points proportional to the volume of that set.
Our derivation, as given below, follows from the intuition that if the disturbances are not too frequent, so that they do not often interact in space, then we can use the superposition principle to find a relation between the response to one disturbance and the response to multiple disturbances. More concretely, our reasoning, as illustrated in Figure S2, requires the following assumptions: 1. We assume disturbance events to be realizations of a Poisson p.p. Φ in R of intensity measure Λ.
2. We assume no spatial interactions between disturbances, thus neglecting co-occurrence events, where multiple disturbances co-occur in a small region.

Derivation
Let us denote g(t) as the recovery function from a single disturbance, so that g(0) = sK, where s is the relative overall strength of the disturbance, and by our definition of return time T R we have: be a random realization of all past disturbance events. The displacement time series of overall biomass under a repeated disturbance Note that for Eq.(S4) to hold, assumption (2) is necessary only if the dynamical response to an individual disturbance is non-linear. If the response is linear, a superposition principle allows Eq.(S4) to be true (Arnoldi et al., 2016), even without assumption (2). In fact, for global disturbances (σ = 1), assumption (2) cannot hold, yet Eq.(S4) still holds because a global disturbance typically induces only small displacements from K, so that the recovery is essentially linear. By ergodicity of the point process, taking an average over long times for one realized time series is equivalent to taking a point average (say at time 0) over realizations, so that We are thus left to analyse random variables of the form T k ∈Φ([0,T ]) g(T k ), a realization of which is t k ∈φ([0,T ]) g(t k ) (note the change of variables −t k → t k , Figure S2: Illustration of how return time following a single disturbance translates into the variability under a regime of repeated disturbances. (a,b) Snapshot at time t * of the local and regional dynamics, following three different disturbances (D1, D2, D3) at times (t 1 , t 2 , t 3 ). The local dynamics (a) are shown using a potential (similar to Figure 2), with three balls and arrows corresponding to the local state at each disturbed region. A snapshot of the two-dimensional system shows the regional dynamics (b), as each disturbance occurred in a different location in space in a different time, and therefore has a different size and different biomass value in the disturbed region. (c) The trajectory of a system following a single disturbance is shown, where the current state at time t * of each disturbance (if it occurred alone) is noted by a red circle. (d) The overall biomass over time due to the three disturbances can be approximated well by adding the three trajectories following a single disturbance in appropriate times, as shown by the grey shading. The effect at time t * of each disturbance and all together is shown by the red circles. (e) A regime of repeated disturbances at different times results in a noisy time series (blue line) that can be measured by the variance (grey shade). This measure of variability V can be well approximated as the product of the average frequency of disturbance f and the second moment of the trajectory following a single disturbance A.
so as to restrict to positive times). We may approximate g(t) with arbitrary precision by a piece-wise constant functiong(t) = i g i 1 A i (t) where the A i 's are mutually disjoint intervals such that A i = [0, T ] and the functions 1 A i (t) are defined by 1 A i (t) = 1 when t ∈ A i and 0 otherwise. The numbers g i approximate the function g in A i . Importantly We thus see that Because the approximation of g(t) byg(t) is arbitrarily good this finally yields: In the main text we defined variability as the variance of h φ (t). If we repeat the above reasoning on h 2 φ instead of h φ we get that, with arbitrary accuracy, where the last two lines follow, respectively, from the fact that N (A i ) and N (A j ) are independent r.v. if i = j (no temporal correlations) and the fact that, for a Poisson process, Var (N (A i )) = EN (A i ) = Λ(A i ). In terms of the integral norms g n Λ,n = g n (t)dΛ(t), what we have proved is that, under assumption (1) and (2) h 2 = g 2 Λ,1 + g 2 Λ,2 (S10) so that variability becomes

Discussion
In the particular case of a homogeneous Poisson process of intensity f , so that Λ(A) = f |A| then we have that as illustrated in Figure S2. Note that homogeneity is not essential to our arguments. For instance, if t is expressed in years, then seasonality effects could be modelled by taking, dΛ(t) = f × (1 + α(2 cos 2 (πt) − 1))dt (S13) with 0 ≤ α ≤ 1. What is essential is the absence of temporal correlations between disturbance events, which assumes that a disturbance does not influence the occurrence of future disturbances. This is not always true in natural systems, (e.g. a drought can make a fire more likely), but is a reasonable assumption for most climatic events. Note that we assumed so far that all disturbances had the same spatial extent σ. This assumption can easily be relaxed by assuming σ to be a random variable, with p.d.f. p(σ). The contributions to the overall variability of disturbances of size σ ∈ [σ, σ + dσ] is dV (σ) = g σ 2 Λσ,2 p(σ)dσ; so that In this sense, the relationship Eq.(S14) is quite general. However, it requires assumption (2) to hold, especially for localized regimes of disturbances, for which the response will typically be highly non-linear. As we explain in the next section, departures from this prediction is the sign that disturbances are aggregating into long-lived disturbed regions, a process that increases variability and, in bistable systems, can lead to a regime shift.
Finally, it is interesting to note that had we assumed no temporal overlap between responses to disturbances (assuming a very low frequency of events, for instance) then we would have had where |φ([−T, 0])| is the number of disturbance events that have occurred in [−T, 0]. For a homogeneous p.p. and for n = 1 we recover Eq.(S8). However, for n = 2 we do not find Eq.(S10). This shows that assuming no temporal overlap between responses to compute V would induce an error equal to the squared mean displacement h 2 . An error that can easily be significant (but is of second order in f as the frequency of events goes to zero).

Appendix C: Aggregation of disturbances
As discussed in the Results section, the aggregation of disturbances is the main reason for the under-estimation of variability and for high collapse probability ( Figure 4). We describe here how this phenomenon affects variability and collapse probability.
We define a disturbance aggregation as a situation when multiple disturbances occur within the same time frame in different locations so that effectively a large contiguous region in the system is disturbed at a given time. This aggregation typically reduces the total biomass, and more significantly, leads to an overall slow recovery when compared to a simple addition of multiple non-interacting disturbances. As shown in Figure S3, occurrence of such aggregation changes the recovery trajectory considerably (top panel), which leads to higher variability (bottom panel). These events are more common for mid-sized disturbances, which is the reason that the under-estimation is most prominent for these values of σ. This is because an aggregation is more likely to occur when disturbances have both a longer time-span and a larger spatial extent. Larger values of σ mean that the size of each disturbed region is larger, and therefore it takes less disturbances to cover a certain region. On the other hand, more local disturbances (smaller σ) tend to lead to slower recovery (this can be quite different between models, depending on local dynamics), and hence the effect of each disturbance is retained for a longer time. Overall mid-sized disturbances are the middle-ground between these two conditions, and it is therefore there that disturbance aggregation occurs most frequently. For the specific case of a bistable system such as the AE model, mid-sized disturbances also have the slowest recovery, so the effect is even stronger.
We note that for bistable systems, a collapse to the alternative state is possible if enough disturbances push the system to the other state. This is essentially an extreme case of disturbance aggregation, where the disturbances cover the entire system. Therefore, when increasing the frequency of disturbances we first see a collapse occurring for mid-sized disturbances, and only for higher frequencies do we see it for smaller and larger extent of disturbances (see Figure S4).

Appendix D: Generalization of results
The results presented in this paper have focused on a simple model setup for the sake of clarity, but they can be easily generalized in various settings. For example, in eq. S13 in Appendix C we show how temporal seasonality can be taken into account and give the same overall results. We further demonstrate here that our results are not specific to a particular spatial structure, with two specific examples. We compare between one-dimensional (1D) and two-dimensional (2D) systems, looking at the results for return time, variability and collapse probability. We also consider a regime of repeated disturbances where the size of disturbances is not the same every time, but rather taken from a probability distribution.  Figure S3: Dynamics of single simulation with multiple disturbances, demonstrating the effect of disturbance aggregation. The dynamics of a single simulation of the AE model with disturbance frequency of f = 0.01 and mid-sized disturbances σ = 0.15 are shown in the three panels. This is compared with a "constructed time-series" where the effects of a single disturbance are added multiple times to form a time-series that behaves as a simulation where disturbances do not interact in space (so that their effect is always linearly additive). The top panel shows the overall biomass over time, for both the actual simulation (blue) and the constructed time-series (red). The middle panel shows a space-time simulation where the vertical axis shows the space of a one-dimensional simulation, and colour depicts the amount of biomass (darker green denotes more biomass). The bottom panel shows variability calculated over a running time-window of τ v = 4000, where for each point shown, the variability was calculated for part of the time-series with length τ v that is centred on that time point. The variability of the actual simulation is shown in blue, whereas the variability for the constructed time-series is shown in red. A large disturbance aggregation is seen starting t = 20000, which lowers the biomass and increases the variability considerably (blue), with no counterpart for the constructed time-series (red). Parameters used: r = 0.5, s = 0.1.
To demonstrate that the specific spatial structure of the ecosystem is not a fundamental part of our theory, we show here results for a two dimensional (2D) ecosystem, and compare them to results for a one dimensional ecosystem (1D) that we show in the main text. For the 2D case we consider an ecosystem with the size of 200 × 200 with periodic boundary conditions, and use disturbances with a circular shape. As can be seen in Figure S5, the qualitative properties for return time, variability and collapse probability, as described in the main text, are the same for the 1D and 2D ecosystems. Note that the specific frequencies that are compared are not the same, but rather chosen to show the similarity in trends between the 1D and 2D ecosystems.
The trends of variability and collapse probability as described in the main text hold under more general conditions. We consider here a more realistic scenario by allowing disturbance's size to randomly vary around a fixed average. As an example we set σ to be randomly chosen from a Gaussian distribution with a standard deviation of 0.05 around a given average ( Figure S6). Both the variability and collapse probability show a hump-shaped relationship with the spatial extent of the disturbance in the bistable AE model (top panels), while for models with a single equilibrium such as the SR model, no such behaviour is seen (bottom panels).  Figure 4). Left and middle columns show variability (for low and high frequency respectively), where the dashed line shows numerical estimations with grey shading noting error estimation. The analytical prediction (based on the response to a single disturbance) is shown in a solid black line, where deviation from this prediction implies some degree of interaction between disturbances. Right column shows collapse probability for the high frequency of disturbances. Disturbance extent is taken from a Gaussian distribution with a standard deviation of 0.01, where the mean is changed along the x-axis between 0.1 and 0.2. Comparing with Figure 4 we see that allowing random variations in disturbance extent leads to variability and collapse probability that are more of a smooth function of average disturbance extent while preserving the hump shape for the bistable AE model. Low and high frequencies used are f = 0.005, 0.02 for the AE model and f = 0.02, 0.05 for the SR model. Parameters used: r = 0.5, s = 0.1.