Integrodifference models for persistence in temporally varying river environments

To fully understand population persistence in river ecosystems, it is necessary to consider the effect of the water flow, which varies tremendously with seasonal fluctuations of water runoff and snow melt. In this paper, we study integrodifference models for growth and dispersal in the presence of advective flow with both periodic (alternating) and random kernel parameters. For the alternating kernel model, we obtain the principal eigenvalue of the linearization operator to determine population persistence and derive a boundary value problem to calculate it. For the random model, we establish two persistence metrics: a generalized spectral radius and the asymptotic growth rate, which are mathematically equivalent but can be understood differently, to determine population persistence or extinction. The theoretical framework and methods for calculations are provided, and the framework is applied to calculating persistence in highly variable river environments.


Background
Stream and river ecosystems are shaped by their physical environment of unidirectional water flow. Questions of population persistence in river ecosystems must necessarily consider the effect of the water flow on populations in space and time. To complicate matters, this water flow can vary tremendously with seasonal fluctuations of water runoff and snow melt.
How can populations persist in streams when they are being constantly washed downstream? This so-called drift paradox (Muller 1982) has engaged biologists and mathematicians in a series of modeling efforts using reaction-advection-diffusion equations to describe the population densities in space and time. The early paper of Speirs and Gurney (2001) uses a modification of Fisher's equation that includes advection to show the existence of a critical flow rate in the stream, below which the populations will persist, and above which the population will wash out, much as a chemostat population will persist or wash out in low flow and high flow conditions. The approach of Speirs and Gurney (2001) employs classical mathematical methods of population spreading speeds and critical domain size. The spreading speed for Fisher's equation, 2 √ rD where r is the intrinsic growth rate and D is the diffusion coefficient (Aronson and Weinberger 1975), yields the critical advection velocity v c , below which stream populations will persist, and above which they will wash out. This can be understood intuitively: when the advection velocity v exactly matches the spreading speed 2 √ rD, the population is washed downstream by water flow at the same speed it is moving upstream by the combined effects of growth (r ) and diffusion (D). Speirs and Gurney (2001) show the critical domain size L c exists for all advection velocities that lie below the critical value (0 < v < v c ) and that the critical domain size approaches infinity as v approaches v c . Biologically, this is interpreted as implying that stream populations will persist if the advection speed falls below a threshold value, and there is a sufficiently large stretch of stream available. This theory has been tested empirically by Walks (2007) who related the persistence of plankton in flowing water to stream advection velocities. A mathematical review of the ideas in Speirs and Gurney (2001) can be found in Lewis et al. (2009).
Extensions to the theory have focused on increasingly realistic models for the stream populations. These include stationary and mobile compartments to describe subpopulations on the benthos and in the stream , non-diffusive dispersal of stream populations that can include long-distance jumps , spatially varying stream environments (Lutscher et al. 2006), spatial interactions between competitors in the stream environment (Lutscher et al. 2007) and periodic fluctuations in environmental conditions Lewis 2011, 2012).
Despite these extensions to the theory, the models have been limited to the case where the stream environment is predictable. Although convenient from a modelling perspective, this is inaccurate. For example stream flows not only vary by an order of magnitude between spring and fall seasons (Abrahamsson and Hakanson 1998), they also vary unpredictably from year to year (Anderson et al. 2006). While some models exist for spreading populations in randomly fluctuating (Neubert et al. 2000) environments, none have investigated persistence and spread in environments such as streams, where unidirectional flow predominates.
In this paper we investigate persistence of populations in periodic and randomly fluctuating environments with predominantly unidirectional flow. Our mathematical model is based on a discrete-time and continuous-space dynamical system that takes the form of an integrodifference equation. In the next section we develop a modelling background for integrodifference models.

Integrodifference models
We consider organisms with separate growth and dispersal stages. Dispersal is assumed to be continuous in space and occurring over a fixed time interval (the dispersal event) while growth is independent of space, but depends on the local population density. Denoting n t (x) as the population density at stage t, the growth dynamics are modeled by f (n t ) = n t (x)g(n t (x)), ( where f is a nonnegative monotonically increasing function. The function g(n) is the per capita growth rate and we assume the maximum per capita growth rate is found as n approaches zero. The dispersal dynamics are modelled by the integral equation where the dispersal kernel K (x, y) models the probability density associated with an individual, that starts at y, settling at x during the dispersal event. We assume that K is a continuous nonnegative function with area one when integrated with respect to x over the real line for all fixed y. The combined model for growth and dispersal is then given by the the nonlinear integrodifference model Equations of the form (3) were first formulated to study gene flow and selection (Slatkin 1973), and were only later applied in ecological settings (Kot and Schaffer 1986). As written, Eq. (3) assumes an unstructured population that grows in the stream benthos, and disperses through the stream and settles back to the benthos each time step. Although the population is assumed to be unstructured, an extension of the model can be used to describe stage-structured populations where dispersal varies from stage to stage (Lutscher and Lewis 2004). We consider a habitat Ω = [x 0 , y 0 ] for some x 0 < y 0 . For such a bounded domain Ω the model (3) assumes that the organism can disperse across the boundary, but there is no source term from outside the boundary. This is the case if conditions outside Ω are unfavorable to growth and survival or if the organism cannot disperse back into the habitat Ω once it has left. This would be the case for a stretch of suitable habitat in a stream, surrounded by unsuitable habitat, where the organism cannot survive. A non-aquatic example is of a plant whose seeds are blown across the edge of a field into a parking lot or other unsuitable region.

The dispersal kernel
The dispersal kernel K can take a variety of forms. If an individual at y moves randomly for a fixed amount of time T and then settles, the dispersal kernel is a Gaussian with variance 2DT where D is the diffusion coefficient associated with the random movement. Alternatively, if the randomly moving individual settles at a constant rate β > 0 then, after a sufficiently long period of time, the kernel approaches a Laplace distribution where a = β D (Neubert et al. 1995). Figure 1a shows two sample Laplace kernels with a = 4 and a = 1.5 for the habitat Ω = [−1, 1]. Note that the kernel for the larger value a = 4 corresponds to a higher probability of the organism settling in Ω, which is consistent with the higher value of a arising from a larger settling rate or a smaller diffusion coefficient.
If, in addition, an organism experiences a unidirectional flow with velocity v (e.g., stream flow or wind) the kernel takes the form where the rate constants a i are defined in terms of the advection velocity v, settling rate β, and diffusion coefficient D by Lutscher et al. 2005). Figure 1b shows the same kernels in Fig. 1a, but now with an advective velocity v = 4. The derivation of (5) assumes a separation of time scales between settling and dispersal to employ a partial differential equation describing the time-dependent probability density z(t, x) of an individual that starts moving at point y at time t = 0 via advection (due to the flow), diffusion (as a first approximation to variability in the flow speed and direction), and settles at rate β: The density of settling points at x is then given by Integrating (6) over 0 < t < ∞ yields an ordinary differential equation for k defined on the real line −∞ < x < ∞ whose solution is K (x, y), as defined by (5). Here it is assumed that the dispersing individual does not modify its movement in response to the domain boundary. If such behaviour is included, it gives rise to boundary conditions for (8) which modifies the associated Green's function for K , Appendix E). Two quantities that are derived from the dispersal kernel and are relevant to our modelling considerations are the dispersal success function and the redistribution function. The dispersal success function s(y) indicates the probability that an individual starting at y successfully settles in the habitat Ω after the dispersal event: Since K (x, y) ≥ 0 and R K (x, y 0 ) dx = 1 for any y 0 we have 0 ≤ s(y) ≤ 1. The redistribution function r (x), on the other hand, corresponds to an area release experiment; if N individuals are released uniformly over the patch, then after one dispersal event the expected density of individuals is Nr(x) where By way of contrast with the dispersal success function s, the redistribution function r needs not be bounded above by 1 for all x. However, for symmetric kernels the two functions are identical. Figure 2 shows an example of the dispersal success and redistribution functions for an asymmetric Laplace kernel on Ω = [−1, 1]. The dispersal success function provides a means to approximate the principal eigenvalue of the linearization of (3), which itself can be used as a measure of population persistence (see Sect. 1.3.1). More specifically, let λ 1 (K ) denote the principal eigenvalue of the linearization of (3) at the equilibrium solution n * (x) = 0 and φ an associated positive eigenfunction, i.e., where R = f (0). It has been shown in Lutscher et al. (2005) that λ 1 (K ) is a strictly increasing function of L. If we assume Ω φ(y) dy = 1 and integrate (11) over the habitat we obtain the following relation between λ 1 (K ) and the dispersal success function s(y): Taking the approximation φ(x) ≈ 1 L , we can estimate λ 1 (K ) by which is known as the dispersal success approximation. To illustrate, Fig. 3 shows the dispersal success approximation (13) compared to the principal eigenvalue λ 1 (K ) as a function of stream length for a sample kernel.

Including temporal variation
In this paper we consider a generalization of model (3) that allows for the kind of temporal variation in growth and dispersal that is found in stream ecosystems. We focus on including the inter-annual variations in growth and dispersal that arise from temporally fluctuating environments. An organism in such an environment experiences different growth and dispersal dynamics, depending on the year. In this case the model The principal eigenvalue λ 1 (K ) and the dispersal success approximation λ a,1 for linearization of (3) with f (0) = 1.2 and asymmetric kernel K with D = 1, β = 1, and v = 0.1. Note that λ a,1 tends to underestimate λ 1 (K ) as the domain length increases where K t denotes the tth time step dispersal kernel and f t (n) the growth dynamics at time step t.

Mathematical setting
We briefly review the mathematical setting and known results for population persistence in the context of integrodifference equations. We restrict our attention to the case where K (x, y) = K (x − y) is a difference kernel, expressed in terms of the difference between the settling location x and the starting point y. Although this includes the case of symmetric distance kernels where K (x, y) = K (|x − y|) (such as (4)), we do not require symmetry since we are particularly interested in the case where K is an asymmetric advective kernel such as (5). We assume the population densities are given by elements of C(Ω), the Banach space of continuous real-valued functions defined on Ω. We first discuss the constant kernel case, then review some known results for temporally varying kernels.

Constant environments
We can rewrite Eq. (3) as where F : C(Ω) → C(Ω) is the nonlinear Hammerstein operator If f and K (x, y) are continuous, then it follows from the Arzelà-Ascoli theorem that F is a compact operator. Moreover, since K ≥ 0 and f ≥ 0, the operator F is a positive operator, mapping the cone of nonnegative functions C + (Ω) into itself. The linearization of (15) at the equilibrium n * (x) = 0 takes the form of a Fredholm equation of the first kind where L = F (0) is the Fréchet derivative of F at n * (x) = 0 and R = f (0) is the geometric growth rate for the population (Van Kirk and Lewis 1997). Since the Fréchet derivative of a compact operator is compact (Krasnoselskii 1964b), the operator L is a compact bounded linear operator. A further assumption regarding the positivity of L allows one to connect the population dynamics of (15) with the spectral properties of L . Namely, we say the operator L is strongly positive if for any continuous function n ≥ 0 there exists a power t = t (n) such that L t (n)(x) > 0 for all x ∈ Ω. Biologically, this condition implies that on a connected habitat repeated application of the kernel will allow an individual that starts at any point y ∈ Ω to eventually access all other points x in Ω (Lutscher and Lewis 2004, Condition A4). In this case, the Krein-Rutman theorem applies, and it follows that L has a principal eigenvalue λ 1 > 0 such that |λ| < λ 1 for all other eigenvalues and λ 1 is the only eigenvalue associated with a positive eigenfunction φ. Our assumptions on K imply that the zero solution of (15) is linearly stable when λ 1 < 1 and unstable when λ 1 > 1. Moreover, for λ 1 > 1 there exists a nontrivial equilibrium solution of (15) (Hardin et al. 1990, see also Van Kirk and Lewis 1997 for the L 2 (Ω) case).
Within this context, the principal eigenvalue λ 1 is equal to the spectral radius of the operator L which, by the Gelfand formula, can be expressed as and · ∞ is the sup norm on C(Ω) (Krasnoselskii 1964a,b). In particular, λ 1 ≤ L . Moreover, if we consider a positive eigenfunction φ 1 and integrate the associated equation where we used the fact that 0 ≤ s(y) ≤ 1. Thus, λ 1 ≤ R = f (0). If we assume there is dispersal loss from all points in the domain then s < 1 and we obtain a strict inequality We can interpret (20) biologically as indicating that dispersal loss will reduce the growth rate below the intrinsic growth rate for the non-spatial model.

Temporally Varying Environments
A series of papers by Hardin et al. (1988aHardin et al. ( ,b, 1990 consider the case where the Hammerstein operator is time-dependent so that where with K t having parameters that are chosen randomly from some set (defining the allowable range of environmental conditions). Here the linearized operator L t = F t (0) is time-dependent, so the eigenvalue analysis discussed in the previous section does not apply. However, they show that population persistence can still be understood via the limit of operator norms which acts as an effective spectral radius for the time-varying setting. Hardin et al. (1988a) derive conditions under which the limit (23) exists, and show how the quantity r determines long-term population persistence or extinction for solutions of (21). The limit r is similar to the dominant Lyapunov exponent (stochastic growth rate) of random matrix products, which has been applied to determine population persistence for a structured population in both correlated and uncorrelated random environments (e.g., Benaïm and Schreiber 2009) and coexistence for interacting structured populations living in a random environment (e.g., Roth and Schreiber 2014).
We consider (22) in the context of growth and dispersal in streams for the case of periodic and random dispersal parameters and show the Hardin et al. framework can be adapted to our setting. A range of different types of growth rates (and related elasticities) of populations have been introduced and used to study population dynamics in random environments (see e.g., Tuljapurkar 1990 andTuljapurkar et al. 2003). In this paper we consider population persistence via an asymptotic growth rate where n t (x) is defined by the system with nonzero initial condition n 0 (x) ≥ 0. We will prove that Λ = r , and hence, numerically, we can calculate Λ to determine population persistence or extinction.
We use this to consider several examples in the context of randomly fluctuating river populations.
Note that if A is the infinitesimal generator of a continuous semigroup For many generators A it is known that −∞ ≤ s(A) ≤ ω 0 < ∞ (see e.g., Hille and Phillips 1957) and the conditions for the equality of these two quantities have been studied (Greiner et al. 1981;Kato 1982;Thieme 2009). The quantities s(A) and ω 0 have also been used to derive persistence conditions for population models in temporally homogeneous or heterogeneous environments (see e.g., Thieme 2009). In this current work, r and Λ, as defined by (23) and (24), are analogous to the spectral bound and type for the infinitesimal generator A of a continuous semigroup T . Thus this work can be considered as a generalization of the idea of s(A) = ω 0 and using such quantities to determine population persistence.

Outline of the paper
In this paper, we study the integrodifference equation (15) for population persistence in temporally varying advective environments. In Sect. 2, we study the integro-difference equation with alternating kernels and growth rates in a periodically varying environment and obtain an explicit method to calculate the principal eigenvalue for the twostage process. We also give the approximation of the principal eigenvalue by virtue of the dispersal success function and the redistribution function. In Sect. 3, we study the model in a randomly varying advective environment, where both the growth rate and the dispersal kernel are random. This contrasts with the earlier work by Hardin et al. (1988a), where only the growth rate fluctuated randomly. We derive the persistence metric r , similar to (23) and obtain its equivalence to the asymptotic growth rate (24). We also provide exact formula for the asymptotic growth rate when kernels take an asymmetric advective form (5) with randomly chosen parameters. This allows us to explicitly calculate population growth rates in randomly fluctuating river environments. Our various methods for calculating persistence and growth metrics are illustrated using numerical examples for models describing randomly fluctuating rivers. A short discussion completes the paper in Sect. 4.

Alternating kernel model
We consider a deterministic case of time-varying kernels for the linearized model (17).
In particular, we consider the case of alternating kernels K 1 (x, y) and K 2 (x, y) with associated growth rates R 1 and R 2 . We can consider the two stages in succession using the linearized model (17) as: In this way, the two-stage model can be considered as a single model (for two stages) via where K is defined by (26) and R = R 1 R 2 . Let λ 1 be the principal eigenvalue associated with the two-stage model (27). The results in Sect. 1.3.1 imply the zero solution is unstable if λ 1 > 1 and it is stable if λ 1 < 1. Hence, the population persists if λ 1 > 1 and the population will be extinct if λ 1 < 1. It follows from estimate (20) that Since the principal eigenvalue λ 1 for the two-stage process models two years of population dynamics, the quantity √ λ 1 would be an effective single year growth estimate. In this sense, estimate (28) says the effective annual growth rate is bounded above by the geometric mean of the two growth rates.
We can gain further insight into the two-stage λ 1 by working directly with the eigenfunction equation. Suppose λ = 0 is an eigenvalue associated with an eigenfunction φ for (27). Then where K is the two-step kernel (26). Let ψ be defined by Then φ solves (29) if and only if Now suppose K 1 and K 2 are advective dispersal kernels of the form (5) for some β 1 , D 1 , v 1 and β 2 , D 2 , v 2 . Differentiating (30) and using (8) for K 1 we have Similarly, differentiating (31) and using (8) for K 2 we have It follows that the eigenfunction φ solves the fourth order equation where The boundary conditions can be determined as follows. Suppose Ω = [0, L]. First, differentiating (30) and (31) and using the definition of K we have and where a i, j denotes the constant a i in (5) for the kernel K j . Next, differentiating (33), using (37)-(38) and (33) we obtain two additional boundary conditions: The differential equation (34) together with the four boundary conditions (35)-(36) and (39)-(40) comprise a fourth-order boundary-value problem for the eigenpair (λ, φ) of (29) in the case of two-stage advective kernels. The general solution of equation (34) has the form where r i ∈ C are the roots of the associated characteristic polynomial. Applying the boundary conditions to φ yields a fourth-order linear system of the form Ac = 0 for the coefficients c = [c 1 c 2 c 3 c 4 ] T . This system admits a nontrivial solution, and hence φ in (41) is nontrivial, only if det A = 0, which, for a fixed domain length L, defines an implicit equation for λ. Solving det A = 0, we then obtain the principal eigenvalue λ 1 = λ 1,twostep of the two-stage operator. To illustrate, we derive the fourth order boundary value problem for (27) in the symmetric kernel case in Appendix A.1. We also note that this process can be used to determine the critical domain length by setting λ = 1 and determining conditions on L for which the fourth-order system admits a nontrivial solution (see e.g., Jin and Lewis 2011).
Example 1 To illustrate how the principal eigenvalue λ 1,twostep of (29) depends on flow velocities we consider the case of variable flow rates v 1 and v 2 , but with fixed As the flow rates v 1 and v 2 vary, while keepingv fixed, the value of λ 1,twostep varies within some interval. Figure 4 shows the possible values of λ 1,twostep as a function of the average of flow velocityv ∈ [0, 20]. Overall, λ 1,twostep decreases withv. When the averagev is sufficiently small (i.e.,v <v a ), then λ 1,twostep > 1 regardless of v 1 and v 2 , which corresponds to persistence of the population; whenv is sufficiently large (i.e.,v >v b ), then λ 1,twostep < 1 regardless of v 1 and v 2 , which corresponds to extinction. For moderate values ofv (i.e.,v a <v <v b ), different combinations of v 1 and v 2 may lead to λ 1,twostep > 1 or λ 1,twostep < 1, and hence, the population can persist or go extinct in different fluctuating flows even though the mean of flow velocity is constant.
For a fixed meanv, we can also consider λ 1,twostep as a function of the variation in flow |v 1 − v 2 |. Figure 5a illustrates this for the casev = 1.3. Note that λ 1,twostep is an increasing function of |v 1 − v 2 |. For this average flow velocity, the smaller the variation between v 1 and v 2 , the smaller the possibility that the population can persist in the river.
If we fix v 1 and vary only v 2 , then Fig. 5b shows that λ 1,twostep is a decreasing function of v 2 . This coincides with the fact that, when the flow velocity in one step is constant, then the larger the flow velocity in the second step, the harder it is for the population to persist in the river.
Example 2 To illustrate how one can study critical domain size questions in this setting we consider the case of fixing v 1 = 0.1 and determining the critical domain length as a function of v 2 (leaving the other parameters as in Example 1). We can study this by setting λ = 1 in (34) and determining conditions on L for which the fourthorder system admits a nontrivial solution. An example is shown in Fig. 6. As one might expect, as v 2 increases the critical domain length increases, with L approaching infinity The relation between the critical domain size and v 2 for the alternating kernel model with v 1 = 0.1, R 1 = 1.2, R 2 = 1.5, D 1 = D 2 = 1, and β 1 = β 2 = 1. As v 2 → 3, the critical domain size approaches infinity as v 2 tends to some value. Since the critical domain size represents the minimal length of the river such that population can persists, this observation implies that the higher the flow the more difficult it is for the population to persist in the river, consistent with earlier results in Lutscher et al. (2005) and Jin and Lewis (2011).
Finally, we note that, similar to the dispersal success approximation in (13), we can use redistribution and dispersal success to approximate λ 1,twostep as where r 1 is the redistribution function corresponding to K 1 and s 2 is the dispersal success function corresponding to K 2 . In Appendix A.2 we include an example comparing this approximation with λ 1,twostep for different domain lengths L.
The differential equation approach in this section for alternating kernels can be generalized to the case of a sequence of n kernels (Jacobsen and McAdam 2014). However, we will instead turn to the case of random kernels, and in particular, asymmetric Laplace kernels whose parameters are chosen from a given distribution.

Random kernel model
In this section, we consider the model where f t and K t denote "random" growth and dispersal kernels at step t. We apply the theory of Hardin et al. (1988aHardin et al. ( ,b, 1990 to the random difference kernel model (43) and show that the long-term population persistence or extinction can be determined by the generalized spectral radius r defined by (23). We also define another quantity Λ (as in (24)), which is mathematically equivalent to r , is computationally easier to work with, and most importantly has a biological meaning of the asymptotic growth rate of the population. We use this alternate framework to consider several examples for random kernels.

Persistence metrics
First, for notational clarity, we rewrite (43) as where {α t } t≥0 is a sequence of independent identically distributed random variables taking values in an index set A (representing the range of environmental conditions). We make the following assumptions on the kernels K α t and growth functions f α t : 2. There exists constants K > 0 and K such that is continuous with f α (u) = 0 for all u ≤ 0. 2. There exists m > 0, f > 0, and f > 0, such that for any α ∈ A, (a) f α (u) is an increasing function in u.
(d) f α is right differentiable at 0. For simplicity, we denote the right derivative as f α (0).
Under these assumptions it can be shown (see Appendix A.3) that the framework of Hardin et al. (1988a) can be applied to the nonlinear model (44), which yields the following result: For nonzero initial data n 0 ∈ C + (Ω), the population n t of (44) converges in distribution to a stationary distribution μ * , independent of n 0 , that is either concentrated at 0 ∈ C + (Ω) (extinction) or supported in C + (Ω)\{0} (persistence).
Moreover, let where α t ∈ A for all t ≥ 1 and L α t := F α t (0) is the linearization of F α t at n = 0.
The following results hold: (a) If r < 1, then the population will go extinct.
(b) If r > 1, then the population will persist.
The quantity (45) provides a means to study population persistence for our integrodifference stream model in the context of random dispersal and growth, within the framework of the hypotheses (C1)-(C3). We now show there is an alternate metric for (44), which can also be used to analyze persistence of the population, is numerically easier to work with, and has a clear biological interpretation.
Consider the linearization of (44) at the trivial solution where R α t = f α t (0). Let n t (x) be the solution of (46) for initial value n 0 ∈ C + (Ω)\{0}. Then the average growth rate of the population over the first t steps can be written as Define the limit of this average as In this sense, the limiting constant Λ represents the asymptotic growth rate of the population.
In Appendix A.4 we show the limit in (47) exists and is independent of the initial function n 0 ∈ C + (Ω)\{0}, so the definition in (47) can be simplified to for any n 0 ∈ C + (Ω)\{0}. Furthermore, as we state in the theorem below, the asymptotic growth rate Λ and the constant r are equal. The proof of this theorem, provided in Appendix A.4, also includes the proof of the existence of the limit in (47) and the equivalence of (47) and (48). (45) and (48), respectively. Then Λ = r.

Theorem 2 Let r and Λ be defined in
Remark 1 It follows from Theorems 1 and 2 that if Λ > 1, the population will be persistent and if Λ < 1, the population will go extinct. Thus population persistence or extinction can be studied by computing Λ for the iterates n t of the linear model (46) (recalling R α t and K α t change at each step).
We illustrate applications of Theorem 2 for several examples of (46), using Λ to determine population persistence or extinction. For simplicity, we use (48) to approximate Λ.
Example 3 (Random two kernel model) Consider (46) where the kernel K t is chosen at random from one of two asymmetric advective kernels K 1 and K 2 , with equal probability. For K i as in (5), we assume v 1 = 0.1, v 2 = 1, D 1 = D 2 = 1, and β 1 = β 2 = 1. Since we are effectively flipping a coin to determine the kernel K t we call this the "coin-flip kernel model" or CFK model. We assume R t = 1.2 if K t = K 1 and R t = 1.5, if K t = K 2 .
First, in Fig. 7 we illustrate sample rates of convergence for Λ for different initial conditions by plotting the tth-approximation of Λ defined by (48) for two different initial states n 0 = 1/20 and n 0 = (π/20) sin(π x/20), on a habitat Ω = [0, 20]. Notice Λ t converges to the same value of Λ ≈ 1.22 for each initial state and at roughly the same rate.
Next, we compare the principal eigenvalue for the alternating kernel model from Sect. 2 with the value of Λ for the random CFK model. Figure 8 shows a plot of the principal eigenvalue λ 1,twostep of the alternating kernel model (27) (using the same parameters from the CFK model) with Λ for the the CFK model (46). The principal eigenvalue of the alternating kernel model appears to match well with Λ for the random model.  Example 4 (Log-normal flow velocities) Our next example considers (46) with the flow rate for kernel K t chosen from a log-normal distribution (keeping the other parameters fixed). We consider the relation between the asymptotic growth rate Λ as a function of the variance in flow rate, while maintaining a fixed mean. First, to illustrate the log-normal distribution, Fig. 9 shows the probability density function for a log-normal distribution with a fixed mean for two different variances.  Figure 10 shows an example of how the asymptotic growth rate Λ depends on the variance of flow velocity, assuming a fixed mean. We see that Λ tends to increase as the variance increases. For a single step kernel with mean flow velocity v = 0.95, the associated principal eigenvalue λ 1 < 1, corresponding to extinction. However, if we increase the variance while keeping the mean fixed, the higher flow velocities are balanced by more frequent low flow velocities which provide favorable conditions for survival. The values of Λ are estimated by computing Λ t for t >> 1. We note that the numerical results are identical for Ω = [0, L] and Ω = [−L/2, L/2]. Figure 11 shows how Λ depends on the variance ν of the log-normal flow rate v for three different fixed means μ = 0.9, 0.95 and 1. Consistent with what one might expect, the higher the average flow rate the smaller Λ is, and hence, the harder it is for the population to persist. Again, the values of Λ are estimates based on Λ t for large t.

Explicit calculation for Λ t
In this section we compute an exact formula for the approximation Λ t of the asymptotic growth rate Λ. For notational simplicity, we write K α t as K t and R α t as R t . Beginning with a uniform initial density n 0 (x) = 1 on the habitat Ω, we have where s 1 (y) is the dispersal success function for the kernel K 1 . Similarly, where s 2 (y) is the dispersal success function for K 2 and r 1 (y) is the redistribution function for K 1 . Continuing, we have and, in general, for t > 3 we have We can compute this exactly when the kernels K i are advective kernels (5) with random parameters v i , β i and D i , provided no two sets of kernels parameters repeat themselves exactly (which is reasonable in the case of random parameter values). Roughly speaking, each population stage will be represented by a sum of exponentials, although the number of terms grows at each step. More precisely, for t ∈ N where a j := a 1 and b j := a 2 in definition (5) for kernel K j , and ρ r,t are certain computable coefficients that depend on the kernels up to step t (their precise form is inter-  55) on Ω we obtain an exact expression for Λ t : Although (56) provides an exact formula for Λ t for the random kernel case, it is not particularly stable for numerical calculations due to error accumulation in light of the many small divisors that appear in the coefficients ρ r,t (e.g., see (72)-(74) in Appendix A.5). This issue is compounded by the fact that our previous calculations of Λ t via numerical integration of (48) showed that the convergence to Λ is fairly slow (e.g., see Fig. 7). Figure 12 illustrates an application of (56) for a specific case of Example 4 with lognormal flow velocities.

Discussion
Even though classical ecological models assume environmental uniformity, the true natural environment shows a high degree of temporal variability. While the yearly specifics of the environmental variations rarely can be predicted, the general nature of the variability, as measured over many years, can be described statistically. One emerging challenge in mathematical biology has been to incorporate such measures of environmental variability into mathematical models for population persistence (see, for example, Benaïm and Schreiber 2009;Tuljapurkar 1990;Tuljapurkar et al. 2003;Schreiber 2010;Roth and Schreiber 2014). Although much recent mathematical attention has focused on this challenge, the pioneering work by Hardin et al. (1988a) actually provided mathematical tools to understand variability for integrodifference equations, as long as a quarter of a century ago.
While the mathematical foundation of our work rests on the seminal papers by Hardin et al. (1988a,b), Hardin et al. (1990), the methods we developed here have extend the approach significantly, and also transform the rather abstract results into concrete applications for the dynamics of river populations. Specifically, and perhaps most importantly, we have connected the Hardin et al. (1988a) metric r (45) to a more biologically reasonable equivalent metric, the asymptotic growth rate of the linearized operator, Λ (48). Indeed, we established the mathematical equivalence of r and Λ and interpreted this equivalence in terms of equivalent persistence metrics for the underlying stochastic, nonlinear dynamical system (44). Recall that r and Λ are defined as where L α t maps from C(Ω) to C(Ω), and We can rewrite Λ as which yields Λ = lim t→∞ L α t •· · ·•L α 1 1/t , but in this case each L α t is considered to map from L 1 (Ω) to L 1 (Ω). Therefore, we can interpret our result mathematically as stating that, when studying population persistence for our random model (44), it does not matter whether the function space is chosen as C(Ω) or L 1 (Ω). However, this result extends beyond the rather narrow mathematical interpretation given above; the asymptotic growth rate Λ has more biological significance and can be easier to calculate than r .
The connection between r and Λ has allowed us to infer persistence properties of the nonlinear stochastic dynamical system, describing population growth and dispersal in rivers, based on Λ. In particular, it means that our explicit calculations of the asymptotic growth rate for river systems with asymmetric exponential dispersal (56) can be rigorously connected to persistence in the associated nonlinear system.
In our analysis we also developed a connection between periodically fluctuating river system, with asymmetric exponential dispersal, and a differential operator describing growth of an associated eigenfunction. Numerical results show a close concordance of persistence thresholds for the alternating kernel model, where good and bad years alternate, and those for a related coin flip kernel model, where good and bad years are chosen randomly with equal probability (Fig. 8).
The class of models in this paper can be generally applied to river or stream populations, where unidirectional flow dominates. However, particular stream populations are likely to require more detailed and specific models. One advantage of a general model is the ability to draw general conclusions. What can be concluded, in general, from the models in this paper regarding the role of variability in persistence in streams and rivers? First, longer streams (Fig. 3) and lower flow rates (Figs. 4, 5b) increase the likelihood of persistence, and higher flow streams must be longer, providing more habitat, if populations are to persist (Fig. 6). These, by themselves are not new theoretical results, and have been understood theoretically since the work of Speirs and Gurney (2001). However, a closer look at Fig. 4 shows that the variability in the flow velocity, as given in the alternating kernel model, can determine persistence outcomes as much as the mean velocity. Specifically, increased variability gives an increased probability of persistence (Fig. 5a). Here the effects of flow rate variation do not simply average out, and the beneficial effect of a low-flow period more than compensates for the detrimental effect of a corresponding high-flow period. This relationship between flow variability and persistence holds over to the more complex case where the dispersal kernel is chosen from a family where the flow velocity is drawn from a continuous probability density function, such as a log-normal distribution (Figs. 9, 10, 11). We considered variations of the parameters in the positive space for the two-step alternating kernel model and the random model and made numerous simulations for the dependence of λ 1,twostep and Λ on the variance of the flow velocity v. In all our simulations, λ 1,twostep and Λ are increasing functions of the variance of v. While we are not able to theoretically prove this result for these two models, Figs. 4, 5, and 9, 10, 11 were typical numerical examples chosen to illustrate the calculations.
Although our model, with uncorrelated random environments, showed that increasing temporal variations can promote population persistence, this phenomenon may not hold in other settings. For example, it has been shown that for a given average population growth rate, temporal variations in the growth rate may increase the risk of extinction; see e.g., Lewontin and Cohen (1969), Turelli (1978), Lande (1993), Halley and Iwasa (1999). Positive temporal autocorrelations in environmental conditions can decrease or increase extinction risk depending on other features; see, e.g., Schwager et al. (2006), Heino et al. (2000), Ripa and Lundberg (1996). In particular, positive autocorrelations in temporal fluctuations can disrupt predator-prey coexistence (Roth and Schreiber 2014). In more general and realistic situations where there are environmental variations in space and time, the effect of interactions between temporal correlations, spatial heterogeneity and dispersal on population persistence becomes even more complex. For instance, metapopulations whose expected fitness in every patch is less than 1 can persist if there are positive temporal autocorrelations in relative fitness, sufficiently weak spatial correlations, and intermediate rates of dispersal between patches (Schreiber 2010). More recently, Roth and Schreiber (2014) develop a coexistence criterion for interacting structured populations in stochastic environments and show, among other applications, that autocorrelations in temporal fluctuations can interfere with coexistence in predator-prey models.
There is much further work that could be done. In this paper, we did not specifically address the critical domain size problem, other than illustrate how our method can be applied for an example with alternating kernels (Example 2). It is not our purpose here to study how the critical domain size is influenced by the variation of different factors, but this could be an interesting avenue for future work, especially for the random model, which would build upon the work for integrodifference equations in Kot and Schaffer (1986) for symmetric dispersal kernels and Hardin et al. (1988aHardin et al. ( ,b, 1990, Van Kirk and Lewis (1997Lewis ( , 1999, Latore et al. (1999) for more general dispersal kernels, including environmental heterogeneity both in space and in time.
By shifting the domain from [0, L] to [−L/2, L/2] we can use symmetry to reduce this to y(x) = c 1 cosh r x + c 3 cos kx.
Applying the boundary conditions (59) and (61) (which now hold at −L/2) we obtain the 2 × 2 system of equations for c 1 and c 3 : (a 1 r 2 + a 2 2 (a 2 − a 1 )) cosh This system admits a nontrivial solution only if the determinant of the coefficient matrix is zero, which yields an equation of the form which implicitly defines the principal eigenvalue λ 1,twostep . Furthermore, by setting λ = 1 in (64), we obtain an implicit equation for the critical domain length L = L(a 1 , a 2 , R 1 R 2 ).

A.2 Redistribution and dispersal success approximations for alternating kernel model
Suppose φ 1 is an eigenfunction associated with the alternating kernel model, i.e., φ 1 solves where K is the two-stage kernel (26). Integrating (65) over Ω yields where s(y) is the dispersal success function for K . If Ω = [0, L] and we use the approximation φ 1 ≈ 1/L we obtain Therefore, we obtain the approximation which approximates λ 1,twostep in terms of the dispersal success function (9) for K 2 and redistribution function (10) for K 1 .
To compare the estimates for the alternating kernel model we consider (27) on Ω = [0, L] with kernels K 1 and K 2 defined by (5). We denote the dispersal success approximations for the single stage case with either K 1 or K 2 as Figure 13 shows an example of the comparison of the approximations for the principal eigenvalues of the single stage and two stage operators with λ 1,twostep , the actual principal eigenvalue, obtained by solving the boundary value problem (34)-(40). The principal eigenvalues λ 1,twostep and (λ 1,twostep ) 1/2 both increase as L increases (as one would expect), λ a,1 underestimates λ 1,twostep but λ a,2 overestimates λ 1,twostep . The geometric mean (λ a,1 λ a,2 ) 1/2 of one year approximations appears to provide a good estimate for (λ 1,twostep ) 1/2 when L is not too large, but it overestimates (λ 1,twostep ) 1/2 when L is large. The geometric mean (λ a,1 λ a,2 ) 1/2 also tends to underestimate the two-step approximation (λ a,12 ) 1/2 , which itself overestimates the actual persistence measure (λ 1,twostep ) 1/2 . One should note that this conclusion is based on the results for the single model and the two-stage model with chosen parameters in the example. First we recall some preliminary facts from Hardin et al. (1988aHardin et al. ( , 1990. Consider the general model: where t = 0, 1, 2, . . ., ω ∈ Ω, Ω is a compact set in R n for some n ≥ 1, and x 0 ∈ C + (Ω). They consider the following assumptions: (C0) α t is a sequence of independent identically distributed random variables in some index set A.
(H1) For any α ∈ A, the continuity of F α in C + (Ω) follows from the continuity assumption on K α and f α . Moreover, (C1) and (C2) imply that F α (u) = 0 if and only if u = 0. (H2) By the monotonicity of f α , for any α ∈ A and u ≥ v in C + (Ω) we have Then d > 0 and for any α ∈ A, (H4) First we show that F α is compact for each α ∈ A. Let {u k } be a bounded sequence in C(Ω). Then for any α ∈ A we have which shows {F α (u k )} is uniformly bounded. Next, for any x 1 , x 2 ∈ Ω we have Since K α is continuous, it follows {F α (u k )} is equicontinuous. Thus, by the Arzelà-Ascoli theorem, {F α (u k )} has a convergent subsequence which implies F α is a compact map. Therefore, for any α ∈ A, where h = f K |Ω|.
(H10) (a) For any u, v ∈ C(Ω) with u ≤ v and x ∈ Ω, (b) If u ∞ ≤ 1 then u(x) ≤ 1 for all x, and it follows from monotonicity By a similar calculation as in (H5), it follows Thus for any x ∈ Ω, Since ξ u ∞ b ≤ 1, it follows from (H3) and (H7) that By (H1)-(H9) it follows from Theorem 4.2 in Hardin et al. (1988a) that for all n 0 = 0, n 0 ∈ C + (Ω) with probability one, n t converges in distribution to a stationary distribution μ * , independent of n 0 , such that either μ * ({0}) = 0 or μ * ({0}) = 1. By (H10)(c) we know L α ∞ ≤ h for all α ∈ A. Thus from Lemma 5.1 in Hardin et al. (1988a) it follows that lim t→∞ L α t • · · · • L α 1 1/t exists with probability one, and hence we can define a constant r as Finally, it follows from Theorem 5.3 in Hardin et al. (1988a) that if r < 1, then μ * ({0}) = 1 and n t → 0 with probability one and if r > 1, then μ * ({0}) = 0. In this sense, if r < 1 then, in the long run, the population density approaches zero with probability one (extinction), while if r > 1 the population density will not approach zero, and hence persists.
A.5 Explicit solution for random asymmetric exponential kernels Here we present the details for the dispersal success approximations (52) when the kernels K t are random advective kernels of the form (5), i.e., where the rate constants are defined by: where v is the advection velocity, β is the settling rate , and D is the diffusion coefficient for the kernel K t . The key step is the following lemma which can be shown by a direct integration: Lemma 1 Let Ω = (−L/2, L/2) and suppose K t (x − y) is an asymmetric advective kernel of the form (70).