Tail‐dependent spatial synchrony arises from nonlinear driver–response relationships

Abstract Spatial synchrony may be tail‐dependent, that is, stronger when populations are abundant than scarce, or vice‐versa. Here, ‘tail‐dependent’ follows from distributions having a lower tail consisting of relatively low values and an upper tail of relatively high values. We present a general theory of how the distribution and correlation structure of an environmental driver translates into tail‐dependent spatial synchrony through a non‐linear response, and examine empirical evidence for theoretical predictions in giant kelp along the California coastline. In sheltered areas, kelp declines synchronously (lower‐tail dependence) when waves are relatively intense, because waves below a certain height do little damage to kelp. Conversely, in exposed areas, kelp is synchronised primarily by periods of calmness that cause shared recovery (upper‐tail dependence). We find evidence for geographies of tail dependence in synchrony, which helps structure regional population resilience: areas where population declines are asynchronous may be more resilient to disturbance because remnant populations facilitate reestablishment.


S1 General theory
To substantiate our theoretical claim that, under fairly general conditions, appropriately nonlinear, synchronous environmental influences on population growth rates can produce tail dependent spatial synchrony (asymmetric tail associations) between population time series in different locations, we first formulate a general model. Suppose N i (t + 1) = N i (t)λ(N i (t)) exp(e s (t) + e l,i (t)) (S1) for i = 1, 2, where N i (t) is population density at location i and time t, λ is a densitydependent growth rate, e s (t) is a spatially synchronous environmental effect on growth rates, and e l,i (t) represents local-noise effects on growth rates. We assume the e l,i (t) are normally distributed with mean e l and standard deviation σ l , and that they are independent across both space and time. We use either and overbar or E(·) to denote expected value of an expression. We assume the e s (t) are independent and identically distributed (iid) across time, and are independent of the e l,i (t).
We let e s (t) = f (δ(t)), where δ(t) is normal with mean δ and standard deviation σ s , and where f (x) = c 1 + exp(−( x a + b)) is a sigmoid function. The random variable δ(t) represents a spatially synchronous environmental variable and f represents its influence on population growth rates. The function f is shaped similarly to the logistic function, but it has been "stretched" both horizontally (controlled by the parameter a) and vertically (controlled by the parameter c), and also translated both horizontally (controlled by b) and vertically (controlled by d).
We assume λ(N ) is a monotonically decreasing continuous function of N , and that lim N →∞ λ(N ) = 0. Defining e = e s + e l , we furthermore assume λ(0) > 1/ exp(e). We can then uniquely define N * = λ −1 (1/ exp(e)), and note that N * = N * λ(N * ) exp(e s + e l ), i.e., the deterministic one-patch model obtained by setting e s (t) = e s and e l,i (t) = e l has an equilibrium at N * . We assume this equilibrium is stable, i.e., After some calculations (see next paragraph), this reduces to the equivalent assumption that This completes the presentation of the model and associated assumptions.
To see that (S4) is equivalent to (S3), we proceed as follows. First, by computing deriva-tives, (S3) is equivalent to But we know from above that N * = N * λ(N * ) exp(e), so 1 = λ(N * ) exp(e), so (S5) becomes This is equivalent to But N * > 0, dλ dN < 0, and exp(e) > 0, so the first inequality in (S7) is automatically satisfied, and therefore (S3) is equivalent to But because λ(N * ) = 1/ exp(e) (see above), this is equivalent to (S4). The right side of (S4) is the instantaneous proportional change in λ(N ) per unit change in N , at N = N * .
Model parameters are e l , σ l , δ, σ s , a, b, c, d, but some of these parameters turn out to be redundant. We now carry out simplifications and parameter reductions, starting with the observation that the parameters δ and σ s are redundant with the parameters a and b.
Multiple combinations of these parameters can lead to the same distribution of f (δ), and the distribution of f (δ) is the only way these parameters influence model dynamics. It is straightforward to see that it suffices to set a = 1 and b = 0 -all distributions of f (δ) that can be obtained by utilizing full flexibility in the values of the parameters δ, σ s , a, and b can also be obtained with a = 1 and b = 0.
• The synchronous environmental driver, e s (t) is of the form e s (t) = c 1+e −δ(t) + d. Here the δ(t) are independent of the e l,i (t) and are iid through time, and are normally distributed with mean δ and standard deviation σ s . And d is selected so that e s = 0.
• The density-dependent growth rate, λ(N ), is a monotonically decreasing, continuous function of N such that lim N →∞ λ(N ) = 0, and such that λ(0) > 1. Furthermore, the model equilibrium is stable, which amounts to assuming dλ The free parameters of the model are then σ l , c, δ, and σ s . The functional form of λ is also flexible, subject to the listed assumptions/constraints; though we will find, below, by linearizing, that for stochasticity of modest intensity, only certain features of λ make a difference.
We next linearize the model. Letting g(N, e s , e l ) = N λ(N ) exp(e s + e l ) and Here, Therefore, (S17) becomes where We know −1 < A < 1, by model assumptions, so this linearized model is a stable ARMAlike process of autoregressive order 1. Any solution x i (t) to (S25) is paired with a solution x i (t)/B to the modified model and the two solutions will have the same synchrony and tail association properties, so it suffices to study (S29). For similar reasons, we can assume, without loss of generality, that c = 1. For "weak noise" (i.e., sufficiently small standard deviations of e s (t) and e l,i (t)), the dynamics of (S29) are essentially the same, for our purposes, as the dynamics of the original model. So we henceforth consider (S29), which has only the parameters −1 < A < 1, σ l > 0, δ, and σ s > 0. We next perform a simulation study to see if outputs of (S29) show asymmetric tail association, and how that depends on model parameters. and σ l = 0.1, 0.2, 0.3, we plotted this difference against δ and σ s (Fig. S1).
Results showed clearly that when the distribution of δ substantially overlapped the left "shoulder" of the logistic sigmoid, resulting populations were predominantly upper-tail associated; whereas when the distribution of δ substantially overlapped the right "shoulder"  Figure S1: Results of the simulation study in section S1. See text for details.
of the logistic sigmoid, populations were predominantly lower-tail associated. These patterns were mitigated for strong density dependence (i.e., large |A|), especially for strongly over-compensatory dynamics (i.e., strongly negative A), but were detectable even in those cases. u + v = 2b u , which intersect the unit square. The partial Spearman correlation associated with the band circumscribed by these lines is

S2 Example of the theory
where sample variances and means are computed over t = 1, . . . , T but the summation is only computed over t such that u(t) + v(t) > 2b l and u(t) + v(t) < 2b u . These are the t corresponding to the points in the band described above.

S4 Details of spatial linear regression and model assumptions
To test for statistical effects on tail dependence in giant kelp spatial synchrony of tail dependence in the relationship between kelp biomass and, respectively, wave calmness, the NPGO, and seawater nitrate concentration, we used multiple regression fit using generalized least squares (Pinheiro & Bates, 2000). The generalized least squares framework has similar assumptions to ordinary least squares regression, including linearity of relationships between predictor and response, and normality of residuals, but allows for various functional forms of correlation in the residuals. Using the 'gls' function in the 'nlme' R package, we fit multiple regression models with the error term having spatial correlation that decays exponentially with distance. The model is: Here, y is the tail dependence strength of giant kelp synchrony averaged over all sites within some threshold distance of the central site (we consider a threshold of 25 km in the main text; see tables S1-S4 for additional thresholds from 10 to 200 km). The β n are regression parameters including the intercept. The x n are the predictors tail dependence strength in, respectively, wave calmness, the NPGO, and seawater nitrate concentration; these are also averaged over the same distance thresholds as tail dependence strength of giant kelp synchrony. is the error term, the spatial covariance of which is specified to decay exponentially Figure S2: Examples of tail association in empirical data analogous to theoretical illustrations in Figure 2b- Table S1: Spatial multiple linear regression relating tail dependence between giant kelp biomass and driver variables to tail dependence in giant kelp spatial synchrony. The distance threshold for averaging was 10 km.
with distance d; a is the nugget and b is the range.
We tested whether model assumptions were met by examining normal QQ plots of model residuals and asking whether spatial autocorrelation in model residuals approximately follows a negative exponential distribution. Spatial autocorrelation was examined using spline correlograms (Bjørnstad & Falck, 2001) using the 'ncf' package in R. The assumption of normality of residuals was well-met by the statistical models (see, e.g., Figure S3a), but spatial autocorrelation in model residuals was not particularly well-described by a negative exponential function (e.g., Figure S3b). However, we proceeded using this model for reasons enumerated below. The observed relationship between spatial autocorrelation and distance is non-monotonic, but existing approaches including other correlation structures available in 'nlme' and notable alternatives like spatial lag models assume monotonic distancedecay in autocorrelation. Hence, there is no clear optimal solution, but by accounting-even imperfectly-for spatial autocorrelation, our models are more conservative than ordinary least squares regression.
We used the same procedures to evaluate the relationship between tail dependence in kelp spatial synchrony and site mean wave calmness. Similarly, model residuals were consistent with the assumption of normality of residuals, but distance-decay of spatial autocorrelation was non-monotonic but we proceeded with using a negative exponential correlation for the error term for the same reasons as above.   Table S5: Spatial multiple linear regression relating tail dependence between giant kelp biomass and driver variables to tail dependence in giant kelp spatial synchrony. The distance threshold for averaging was 200 km. Figure S3: Regression diagnostics for a model of tail dependence in giant kelp spatial synchrony as a function of tail dependence in the relationship between giant kelp abundance and, respectively, wave calmness, the NGPO, and seawater nitrate concentration. The distance threshold for averaging was 25 km. a) normal QQ plot. b) spline correlogram showing spatial autocorrelation as a function of distance. Figure S4: Regression diagnostics for a model of tail dependence in giant kelp spatial synchrony as a function of site mean wave calmness. The distance threshold for averaging was 25 km. a) normal QQ plot. b) spline correlogram showing spatial autocorrelation as a function of distance.