Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://chaos.phys.msu.ru/ryabov/pdfs/Ryabov_Blasius_MMNP2008.pdf Äàòà èçìåíåíèÿ: Thu Oct 30 15:23:30 2008 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 20:12:39 2012 Êîäèðîâêà:

Math. Model. Nat. Phenom. Vol. 3, No. 3, 2008, pp. 42-86

Population growth and persistence in a heterogeneous environment: the role of diffusion and advection
A. B. Ryabov
a a1

and B. Blasius

a

Institute for Chemistry and Biology of Marine Environment (ICBM), Oldenburg University, 26111 Oldenburg, Germany

Abstract. The spatio-temporal dynamics of a population present one of the most fascinating aspects and challenges for ecological modelling. In this article we review some simple mathematical models, based on one dimensional reaction-diffusion-advection equations, for the growth of a population on a heterogeneous habitat. Considering a number of models of increasing complexity we investigate the often contrary roles of advection and diffusion for the persistence of the population. When it is possible we demonstrate basic mathematical techniques and give the critical conditions providing the survival of a population in simple systems and in more complex resource-consumer models which describe the dynamics of phytoplankton in a water column. Key words: diffusion, advection, survival, population, phytoplankton. AMS subject classification: 35K57, 76R10, 76R50, 92D25, 92D40

1. Introduction
For long times field biologists and naturalists have been fascinated by the richness and beauty of complex patterns that can be observed in spatially extended populations. However, the same observations also constitute a challenge for theoreticians who aim to explain this complexity by means of mathematical models. At first glance one might be tempted to argue that the spatial diversity of natural populations mainly originates from some underlying abiotic heterogeneity of the environment. If growth conditions vary between different locations then this spatial variation should be reflected in the density distribution of natural populations. Thus, it is reasonable to assume that a large part of the observed richness in the patterning of biotic landscape can be
1

Corresponding author. E-mail: a.ryabov@icbm.de

42


A. B. Ryabov and B. Blasius

The role of diffusion and advection

attributed to the spatial heterogeneity of growth conditions. On the other hand in nearly all systems populations at different locations are coupled and are able to interact with each other, be it by random, diffusive mixing or by directed, advective flows or currents, which are able to transport individual organisms from one region to another. Note that the term "advection" here is used in a very general sense, including different mechanisms of transport, such as drift caused by a current of water, sinking and floating-up in the gravity field, chemotaxis, etc. The impact of such spatial interactions on the fate of a local population can be very diverse, depending on whether the flow brings-in new immigrants into the habitat or if it takes them away. Furthermore, diffusion and advection can have opposing influences, where, for example, one process may be beneficial for the population, while another process has adverse effects. Usually, the role of an advective flux is evident and mainly depends on its strength and direction. In contrast, the role of diffusivity can be of two kinds. On the one hand, diffusivity accelerates the spread of a population in a habitat and is necessary to provide the population's persistence in a flux. On the other hand, intensive mixing may transport too many organisms into unfavourable zones, resulting in the extinction of the population. Moreover, in consumer-resource models usually also the resource fluxes are driven by diffusivity and advection, a fact which leads to even new patterns and dynamical behaviour. As a consequence, it is quite possible to find stable populations in locations where growth conditions alone would not permit persistence. On the other hand, seemingly well-being habitats may not be able to sustain a stable population, if they suffer population outflows. From all these effects, populations are spatially structured in an intricate interplay between local growth and its geographical variations on the one hand and spatial transport by diffusion and advection on the other hand. This article is devoted to present an introductory review of these topics and to introduce the reader some of the most commonly used models for the population dynamics in a non-uniform turbulent environment. Even though we have in particular the specific case of phytoplankton dynamics in mind, the main concepts considered here hold for any population subjected to mixing and advection. When it is possible, to reproduce the whole picture we perform derivations of the main mathematical relations. Otherwise we aim at least to demonstrate some basic ideas and refer to detailed discussions, providing information about the main techniques that are useful in this field of research. There exist many excellent reviews about spatial population dynamics (see e.g. [32, 47, 65, 72]). Nevertheless to our knowledge, the interplay of advection and diffusion for a heterogeneous population, as elaborated in this text, has never been described. So, while the text will not provide much new insights for the specialist, we hope that many readers will find this text a valuable introduction into the basic mechanism and critical conditions providing the survival of a population in a heterogeneous environment. To focus on the main ideas, we had to confine the review to topics related with population survival. However, we would like to briefly mention other important aspects of spatial population dynamics that had to be omitted here. First, we had to restrict to single species populations, while one of the most important and still debatable challenges concerns the competition of spatially structured populations [104, 25, 98] and the high diversity of phytoplankton species (see e.g. [88] and the references therein). Second, for the sake of simplicity we included only Fickian diffusion models, whereas in the ocean diffusivity depends on the scale of phenomena [72]. Third, we had 43


A. B. Ryabov and B. Blasius

The role of diffusion and advection

to omit the interesting aspects concerning conditional persistence in situations where the growth is influenced by an Allee effect [100]. Fourth, the interaction of turbulent currents with other processes can produce complex spatial structures in phytoplankton distributions [57], whereas locally oscillatory behaviour may lead to spatio-temporal chaotic waves [79, 80]. Finally, we did not even touch such important problems as model validation and the simulation of multicomponent natural ecosystems [60]. The review is structured as follows. In chapter two, we shortly discuss some non-spatial models, which will be used later to build-up spatial explicit models. In chapter three we consider critical patch models, starting from the well known KiSS model, and proceed with more sophisticate models, which draw out the influence of boundary conditions, finiteness of the favourable patch, heterogeneity of the mixing, etc. Furthermore, we consider the Fisher-Kolmogorov equation and its extensions to advection, non-uniform mixing and heterogeneous environments. In the last chapter, we consider 1D models of the vertical phytoplankton distribution, starting from those that include only light limitation and then considering models including many limiting factors. However, even in the latter models we focus on the basic theoretical aspects and general conclusions.

2 . N o n - s p a t ia l m o d e ls
Logistic growth As a warming-up, in this chapter we shortly describe some simple population models in a non-spatial context. These models will mainly be used as a base for developing spatially extended models, but they may also serve as a benchmark for exploring the properties of their spatial analogues. Let P (t) denote the density of a population of interest at time t. In the simplest way, the dynamics of P can be modelled in terms of an ordinary differential equation of first order dP = µ (P )P , dt (2.1)

where the growth rate µ(P ) depends on the population density P . The function µ(P ) takes into account for density dependent regulatory factors which usually are associated with resource depletion or conditions of over-crowding. Consequently, (neglecting the possibility of an Allee effect) we may assume that µ(P ) is monotonically decreasing and eventually becomes negative for large densities. The simplest form to put this into a model is the logistic growth dP = µ0 P dt 1- P K , (2.2)

where K is the carrying capacity of the system and describes the maximal population density that can be sustained by the system, i.e. µ(K ) = 0. Equation (2.2) can be thought to arise from a Taylor expansion of equation (2.1) for small densities P K , where µ0 µ(0) is the growth rate at small density and K -µ0 ( µ/ P )-1 is the carrying capacity. Resource-limited growth As a secondary model class, we investigate resource limited population growth, where the dynamics of the resource is explicitly included into the model equations. 44


A. B. Ryabov and B. Blasius
1.6 1.4 1.6 1.4 1.2 1 0.2 P(t) 0.4 0 0.2 P(t)

The role of diffusion and advection
1.6 1.4 1.2 1 0.4 0 0.2 P(t) 0.4

A

B

C

N(t)

1.2 1 0

Figure model of t he poi nt . c = 5,

1: Typical phase trajectories of the non-spatial models. (A) the dynamics of the conservative (2.3) is reduced to a straight line; (B) the non-conservative model (2.7) leads to the extinction population; (C) the trajectory of system (2.10) spirals in phase space to an isolated fixed The initial values P(0) = 0.1, N(0) = 1.5 are marked by filled circles, other parameters are HN = 5, m = 1, = 0.005, = 0.5, and Ni = 15.

Possibly the most simple way to model an isolated consumer-resource system (i.e., in the absence of any fluxes out of the system) is given in the following form dP dt dN dt = µ(N )P - mP , (2.3) = -µ(N )P + mP .

Here N (t) denotes the resource concentration (e.g., a limiting nutrient such as nitrate or phosphate) at time t and the population abundance P (t) is measured in units of the organism's resource content. As usual, in contrast to equation (2.1), here we explicitly separated the growth into a term proportional to the resource uptake, µ(N )P , and a death process with mortality m. Note that in this review the symbol µ has a somewhat different interpretation in the context of a resource-implicit (2.1) or resource-explicit model (2.3). We hope that the reader is not confused by this notation. Usually, µ(N ) is assumed to vanish for small resource concentrations (µ(0) = 0) and to be saturated for large N . A convenient parametrisation is given by the Monod form µ (N ) = c N HN + N (2.4)

with the half-saturation constant HN . Equation (2.3) represents a closed system, where recycling is assumed to be 100% efficient. Thus the biomass that is taken out of the system due to death processes is remineralised into the nutrient pool, i.e. no material is lost in the conversion from consumer to nutrients. As a consequence, this system can be reduced to a one-dimensional model, despite the fact that it contains two equations. The reason is that the total resource concentration, either free or bound to the population, is a first integral of motion: S = N + P = const. After plugging this into the equation for 45


A. B. Ryabov and B. Blasius

The role of diffusion and advection

P we obtain a single equation of the form (2.1) for the time dynamics of the population density, dP = P [µ(S - P ) - m] . dt (2.5)

If the saturation of the uptake rate is neglected, µ(N ) = cN , we retrieve the logistic growth model (2.2) with m µ0 = cS - m and K=S- . (2.6) c The phase trajectory of system (2.3) in (P , N ) coordinate is a straight line from the initial conditions to an equilibrium state (Fig. 1A). We note though, that the reduction to a single equation works only for a closed system. In general, recycling efficiency will never be fully effective. Thus, one may assume that only a fraction of the dead biomass is remineralised, leading to a non-conservative model dP dt dN dt = µ(N )P - mP ,

(2.7)

= -µ(N )P + mP .

The two approaches introduced in this section will be used in chapter 4 as building blocks to model the local reaction of a spatial extended population. However we should note that both (2.3) and (2.7) have some shortcomings, which ultimately are related to the fact that the models do not contain in- and out-flows. For example, the model (2.7) does not even admit isolated stationary states in the (N , P )-phase plane and for any initial condition leads to the extinction of the population, P (t) 0 (Fig. 1B). While this is no problem in the spatial extended counterparts of the models, which necessarily have exchanges with surroundings specified by the boundary conditions, in local models these exchanges should be introduced as additional terms. Chemostat models To show that these problem disappear as soon as external fluxes are incorporated into the models, we shortly introduce the most simple and frequently used chemostat model. Consider a well-stirred reactor that contains phytoplankton cells with concentration P (t) and a limiting nutrient with concentration N (t). The chemostat is supplied with the nutrient at an input concentration Ni from an external nutritive medium. The outflow contains both medium and phytoplankton cells. Inflow and outflow are characterised by the dilution rate . Then the chemostat equations take the following form dP dt dN dt = µ (N )P - (m + )P , = (Ni - N ) - µ(N )P + mP .

(2.8)

In this model, the total amount of free and bound nutrients S = N + P is not conserved. Thus, the chemostat sustains a stable feasible equilibrium and the system spirals in phase space to an isolated 46


A. B. Ryabov and B. Blasius

The role of diffusion and advection

fixed point (Fig. 1C), see e.g., [38, 10, 39]. Furthermore, for perfect nutrient recycling ( = 1) one can easily show that S (t) follows the simple dynamics dS = (Ni - S ) . dt (2.9)

Therefore, after some transient time, determined by the exchange rate , the system settles to a state where S = Ni . This means that the total amount of nutrients in the system, either bound to the population or free, asymptotically goes to the external nutrient concentration, and one asymptotically retains a one-dimensional equation, similar to (2.5) dP = P [µ(Ni - P ) - (m + )] . dt (2.10)

3. Critical patch models for an extended population
In this section we review some simple one-dimensional spatio-temporal models which demonstrate the intricate interplay between mixing, advection, boundary conditions and other factors which are critical for the persistence of a population. Consider a population growing on a one-dimensional (1D) axis, which may represent either a horizontal or vertical spatial extension. Let P (x, t) denote the population density at time t and position x. We assume that the change of density occurs as a result of the local death-birth processes and of the spatial movement of organisms due to diffusion and advection. The dynamics of such a population can be written in terms of a reaction-diffusionadvection equation [72, 62] P (x, t) = reproduction - advection + mixing t = µ(x, P )P - v P P + D , x x x (3.1)

where µ represents the growth rate (compare to equation (2.1)), v is the advection velocity and D is the diffusivity, which can depend on x and other variables. The advective term v P describes x the drift of the population with the constant velocity v . There are many physical or biological factors which can give rise to advective processes. For instance, in a vertical system the organisms (e.g. phytoplankton cells) may sink (or rise) because they are more heavy (or light) than their surrounding medium. Another example arises for a population in a vertical/horizontal stream (e.g., in a horizontal flow of a river or in an ocean current). Advection and diffusion can be written in terms of a derivative x J (P , x) so that the whole flux of biomass is given as J = vP - D P . x (3.2)

Additionally, model (3.1) has to be complemented by appropriate boundary conditions which specify the environmental influence on the spatially extended system. One can distinguish between 47


A. B. Ryabov and B. Blasius

The role of diffusion and advection

three important types of boundary conditions. The first type (or Dirichlet) boundary condition merely specifies the value of the solution at the boundary. For example, the condition P (0) = 0, P (L) = 0 (3.3)

assumes that the population density vanishes at the habitat edges, which corresponds to an absolutely hostile environment outside the segment [0, L]. The second type (or Neumann) boundary condition specifies the flux across the boundary. For instance, the condition J (x)|
x=0,L

=

vP - D

P x

=0
x=0,L

(3.4)

specifies impenetrable boundaries, in other words there is no in- and outflow of biomass across the habitat edges. The boundary condition of the third (Robin) type is a linear combination of the Dirichlet and Neumann conditions. Consider, for instance, a penetrable barrier, where the flux of organisms across the barrier is proportional to the difference of the population density outside, Pout , and inside, P (x), yielding J (x)|
x=0,L

=

vP - D

P x

x=0,L

= h(P (x) - P

out

)|

x=0,L

.

(3.5)

If the permeability of the barrier h = 0 this condition is equivalent to the zero flux boundary condition, and as h it approaches the Dirichlet boundary conditions P (L) = Pout .

3.1. Persistence on a finite patch: critical patch size and mixing
For simplicity, we begin our investigation with models which implement only diffusion but no advection. We concentrate on the population dynamics on a finite favourable patch and show that in this case diffusion plays a negative role for the populations' persistence. In the absence of advection, v = 0, equation (3.1) reduces to P (x, t) 2P = µ(x, P )P + D 2 (3.6) t x complemented by appropriate boundary conditions. Perhaps the most important biological question that one may ask in this model concerns the conditions under which the population will be able to persist or die out. This question can be trivially answered for an infinite homogeneous environment, where the growth rate is independent of the spatial coordinate, µ(x, P ) = µ(P ). In this case (and assuming that there is no multistability, i.e. the equation µ(P ) = 0 has a unique positive solution P ) the fate of the population entirely depends on the sign of the linearisation of the growth rate around P = 0. The population can invade if and only if the linearised growth rate µ0 = µ(P = 0) > 0. Furthermore, it is easy to verify that the final distribution must be uniform, P (x) = P . In reality the growth conditions will differ between various locations so that µ(x, P ) will depend on x. We refer to such a situation as a "heterogeneous environment". In a heterogeneous system, obviously also the linearised growth rate will be a function of the spatial position, µ0 = µ0 (x), 48


A. B. Ryabov and B. Blasius

The role of diffusion and advection

which makes the problem of the population persistence more complicated. Even if locally the population faces bad environmental conditions, µ0 (x) < 0, due to the population inflow from other locations it may still be able to persist. On the other hand, a favourable patch with µ0 (x) > 0 is no longer save, because diffusion out off the patch can cause additional losses which may even lead to local extinction. KiSS model One simple and elegant model to study this problem has been independently introduced by Kierstead and Slobodkin [44] and Skellam [96]. The main idea is to separate the landscape into a favourable area (which will be denoted as the species' habitat) adjoining some hostile environment from both sides. For its simplicity, Akiro Okubo assigned to the model the name KiSS, which on the one hand includes the authors' initials, and from the other hand can be deciphered as "Keep it Simple Stupid" [97]. The KiSS model suggests a population growing on a finite patch of size L, surrounded by an absolutely hostile environment with infinite mortality (Fig. 2). Thus, per definition the population density outside of the favourable patch equals zero, and it is sufficient to consider the dynamics inside the segment 0 x L, upon the condition that the solution vanishes at the boundaries. To simplify the situation even further, the growth term in equation (3.6) is linearised for small density and is set constant µ0 within the habitat, which yields the following equations P t = µ0 P + D 2P , x2 0xL (3.7) (3.8)

P (0) = P (L) = 0 .

This model could represent, for instance, the vegetation patterns of coastal plants growing on some island in the ocean, where diffusion of plants takes place via seed dispersal. Seeds which are transported out of the island into the water cannot survive, which is equivalent to infinite mortality outside of the favourable patch. A detailed analysis of the KiSS model can be found in [96, 44, 72]. Here we just present one simple solution. Using the method of separation of variables we assume the existence of a solution of the form P (x, t) = X (x) T (t) . (3.9) Substituting this expression into (3.7) we obtain the eigenvalue problem
Tn (t) = n Tn (t) D Xn (x) + (µ0 - n )Xn (x) = 0

(3.10)

where the eigenfunctions Xn (x) should satisfy the boundary conditions (3.8), i.e. Xn (0) = Xn (L) = 0. The solution to the last equation of system (3.10) is well known Xn (x) = An sin µ 0 - n x + Bn cos D 49 µ 0 - n x. D


A. B. Ryabov and B. Blasius

The role of diffusion and advection

Figure 2: Influence of diffusion on the persistence of an isolated population in the KiSS model (3.7). (A) D = 3.2 < D0 , the population growth is unbounded; (B) D = 3.3 > D0 , the population goes extinct. Parameters are µ0 = 2, L = 4, K = 400, thus D0 = 3.24. The horizontal thick black line shows the location of the favourable patch. It satisfies the boundary conditions only if Bn = 0 and the argument of the sine is divisible by n when x = L. From the last condition we obtain the eigenvalues 2 n2 D n = µ 0 - , L2 n = 1, 2 , . . . (3.11)

Solving the first equation of system (3.10), we find that Tn (t) = C en t . The problem is linear, therefore using the superposition principle, we can write P (x, t) as the sum


P (x, t) =
n =1

An sin

n x n t e, L

where the coefficients An depend on the initial conditions. Thus, with the separation ansatz (3.9), we were able to present the dynamics of the full model as a sum of simple "modes". The dynamics of mode n is determined by its eigenvalue n . A mode grows if n > 0, decays if n < 0, and can also oscillate if the imaginary part of n is not vanishing. Moreover, as t the full solution approaches to the mode that corresponds to the largest (principal) eigenvalue , because this mode has the fastest time exponent T (t) e t . From (3.11) we find that the largest eigenvalue = 1 . Therefore, x 1 t - - e. P (x, t)-- A1 sin t L The population grows, and so is able to persist, if 1 is positive. By contrast, for 1 < 0 the population dies out exponentially on the whole habitat. Thus using expression (3.11) for n = 1, we obtain the fundamental relation between the diffusivity, the growth rate and the critical (minimal) size of the favourable patch, which provides the survival of the population L L0 = 50 D . µ0 (3.12)


A. B. Ryabov and B. Blasius

The role of diffusion and advection

Okubo [71, 72] demonstrated that this equation can be obtained by means of simple dimensional analysis. Suppose that L0 = f (D , µ0). Dimensionally, [D ] = m2 /s, [µ0 ] = 1/s, and [L] = m. Thereby, the simplest combination of parameters yields L0 = c D , µ0

where c is a non-dimensional constant which equals in the 1D model. The same expression holds in 2D systems [96, 44], and it can be shown [32, 9] that the constant c depends on the principal eigenvalue, thereby it represents the geometry of the model. Equation (3.12) can also be written in the form D D0 = L2 µ0 , 2 (3.13)

yielding the critical, maximum mixing D0 under which the population can survive. Fig. 2A and Fig. 2B demonstrate the population dynamics for different diffusivities. The population density grows exponentially, if the diffusivity is smaller than D0 (Fig. 2A) and decays exponentially otherwise (Fig. 2B). These results highlight the negative aspects of diffusivity for a population and reveal an important ecological insight, namely that of a critical patch size. A finite population under the influence of random mixing must be larger than a minimal extension L0 in order to sustain a stable population. This critical size L0 simply scales as the square root of the strength of mixing divided by the growth rate, (3.12). Despite the model's simplicity the fundamental results (3.12) and (3.13) have an important ecological message that prevails in more realistic settings. If a favourable patch adjoins unfavourable areas, then the higher the diffusivity, the higher will be the loss rate across the boundaries, thus the larger should be the internal area and the growth rate on the habitat to compensate for these losses. Moreover, one large patch is better for the persistence of a species than two smaller patches of the same total size, since two patches would have four ends, that would double the loss rate. Diamond and May [15], McMurtrie [58], Cantrell and Cosner [8] applied this concept of a critical patch size to the design of national parks and natural reserves of optimal size and form. Logistic growth As the growth of biomass in the KiSS model is not limited, one natural extension of the KiSS model is to consider a logistic growth function P (x, t) = µ0 P t 1- P K +D 2P , x2 for 0 x L (3.14)

and the same (hostile) boundary conditions (3.8). For this model only approximate solutions are available (see e.g. [52], [2], [61], [72]). The stationary solution can be expressed in terms of elliptic functions and was investigated by Skellam [96], Levandowsky and White [54], and Ludwig et al. [56]. However, applying the invasibility criteria, we can conclude that this system possesses the same critical values as the KiSS model. Indeed, as the population density approaches zero, 51


A. B. Ryabov and B. Blasius

The role of diffusion and advection

Figure 3: Dynamics of the population with logistic growth, Eq. (3.14), under different values of diffusion. (A) Conditions are nearly critical, D = 3.2 D0 . (B) Conditions are sufficiently far from the survival-extinction transition D = 0.1 D0 . Other parameters are µ0 = 2, L = 4, K = 400, thus D0 = 3.24. The horizontal thick black line shows the location of the favourable patch. equation (3.14) goes over to equation (3.7). This means that if a population can invade in the KiSS model, it will also invade in the model with logistic growth, and vice versa. Figs. 3 show examples of the population dynamics. In both figures, the resulting population density has a characteristic shape, with large densities in the central part of the patch and a decay of population numbers, the closer one comes to the hostile border. Note, that if the conditions are close to the critical values (3.12) or (3.13), the growth of biomass becomes limited by the diffusive transport, and the maximum of density reached by the population can be much smaller than the carrying capacity. This is illustrated in Fig. 3A where max(P (x)) 6 even though K = 400. However, if the conditions are sufficiently far from the critical region, the carrying capacity is almost reached in the middle of the patch (Fig. 3B). To obtain an intuitive understanding into the influence of the density dependence in equation (3.14) consider a KiSS model with an effective growth rate of µ = µ0 (1 - P (x)/K ) x , which equals the average growth rate of the logistically growing population. The new effect now is that the effective growth rate µ decays with an increase of the population density P (x). Therefore, also the maximal diffusivity D0 , admitting a survival of the population (3.13), will be a decaying function of P (x). This means, however, that the biomass can grow only until the critical diffusivity is reached D0 (P ) = D . Thereby, if initially (when P is negligible) D is close to the critical diffusivity of the KiSS model, a small increase of density is enough to decrease D0 and to reach the balance between production and loss. Finite mortality and other extensions A second unrealistic feature of the KiSS model is the assumption of infinite mortality outside of the favourable habitat. This assumption does not hold in many important cases and limits the applicability of the KiSS model. For example, it hardly could be applied for phytoplankton simulations. Ludwig et al. [56] investigated an extension of the KiSS model, where this assumption of an absolutely hostile environment was relaxed. In this 52


A. B. Ryabov and B. Blasius

The role of diffusion and advection

study the authors assumed a positive growth rate µ0 inside and a finite mortality m outside of the habitat µ0 , for 0 < x < L µ(x) = (3.15) -m, for x 0 or x L . As the mortality is finite, the organisms can survive outside of the habitat and diffuse back. This effect reduces the losses at the edges and the population can survive on a smaller favourable patch compared to the KiSS model. In the model (3.15) the population can persist if L2 D arctan µ0 m . µ0

(3.16)

As the mortality approaches infinity this value approaches the critical patch size L0 of the KiSS model. Many other examples of critical patch models can be found in the books by Okubo and Levin [72] and Murray [62]. We just briefly mention some extensions. Okubo [68, 70] considered a model for growth and diffusion under an attractive force toward the centre of a patch. He found that the modified critical size equals Lc = L0 f (v 2 /4D), where f (0) = 1 and the function f (v 2 /4D) monotonically decreases with v toward zero. Gurney and Nisbet [27] considered a model in which the growth rate parabolically depends on the distance from the habitat centre. This approach is more realistic since it includes a gradual transition from favourable to unfavourable areas. Wroblewski et al. [116], Wroblewski and O'Brien [114] and Platt and Denman [81] included the effect of grazing and obtained an expression similar to (3.12) for the critical patch size, in which, however, µ0 is replaced by µ0 - g , where g characterises the grazing rate. Influence of boundary conditions and spatial arrangement It is interesting to extend the analysis for more complicated spatial geometries. Seno [93] and Cantrell and Cosner [7] investigated the influence of a spatial sequence of favourable and unfavourable habitats and boundary conditions. Seno considered a population of organisms migrating between n patches of different quality. Cantrell and Cosner [7] used a mean field representation of this problem. They compared the total size of a population living on a finite habitat of size L, surrounded by either completely hostile environments or by impenetrable boundaries. Using the logistic model (3.14), they assumed that the habitat possesses patches of different quality, that is, the growth rate changes between positive or negative values, µ(x) = ±1, provided that the favourable and unfavourable patches have the same total area (Fig. 4). In this model the best spatial configuration of favourable and unfavourable patches depends on the boundary conditions. If the exterior region is completely hostile then the location of a favourable patch in the middle of the habitat (Fig. 4A) provides the best conditions for the population, as two buffer zones separate the favourable region from the completely hostile regions, which decreases the population losses. Therefore, the total biomass will be higher than that on the habitat shown in Fig. 4B, which in turn is more preferable than the habitat shown in Fig. 4C, where two favourable patches adjoin the hostile surroundings. 53


A. B. Ryabov and B. Blasius
Most favorable 3 P(x) 2 µ (x) 1 0 -1 0 0.5 x Most favorable 1 Intermediate 3 2 1 0 -1 0 0.5 x Less favorable

The role of diffusion and advection
Less favorable 3 2 1 0 -1 1 0 0.5 x Less favorable 1

A

B

C

3

D

3 2 1 0 -1 0.5 x 1 0

E

3 2 1 0 -1 0.5 x 1 0

F

µ (x) 1
0 -1 0

P(x) 2

0.5 x

1

Figure 4: Comparison of several spatial arrangements of favourable (µ(x) = 1) and unfavourable patches (µ(x) = -1), assuming either lethal boundaries (top) or impenetrable boundaries (bottom) at the borders of the habitat (x = 0 and x = 1) according to [7]. Plotted is the final population density P (x) in model (3.6) assuming logistic growth, µ(x, P ) = µ(x)(1 - P /4) (solid line). The patch quality, µ(x), is shown as dashed line. Each case is classified either as "most favourable", "intermediate" or "less favourable", reflecting the total amount of biomass. If however, the boundaries act as a barrier, the maximum biomass is achieved if a favourable patch directly adjoins one of the impenetrable boundaries (Fig. 4D). Even though the derivation of this result requires rather sophisticated calculations, it is evident that in this configuration there is a single favourable patch which has only one border with an unfavourable environment, and so all possible losses are minimised. By contrast, any splitting of the favourable patch leads to a worsening of the habitat quality (Fig. 4E and 4F). It is clear, that the choice of the best arrangement of habitat quality is very important for the design of national parks, natural reserves, etc. Further, the same mechanism also captures some aspects of relevance for the vertical distribution and competition of phytoplankton species. In a water column the surface acts as an impenetrable barrier. Thus, in an incompletely mixed water column, a more light limited or buoyant species obtains a competitive advantage to another species, whose favourable patch is located in subsurface layers [34, 89]. Arbitrary spatial dependence of the growth rate The behaviour of model (3.6) in the case where the growth rate µ(x) is an arbitrary function of the coordinates was investigated by Cantrell and Cosner [9]. Applying separation of variables to the linearised problem it is possible obtain a

54


A. B. Ryabov and B. Blasius

The role of diffusion and advection

system of equations, similar to (3.10)
Tn (t) = n Tn (t)

(3.17) D X (x) + µ(x)Xn (x) = n Xn (x) , where X (x) should again satisfy certain boundary conditions. Now, however, µ(x) = const and a solution to this equation is known only for some partial forms of the function µ(x). If the principal eigenvalue of these equations is positive, the stationary solution P (x) = 0 is not stable, implying that the population can invade and establish on the habitat. It can be shown that this is always the case if the average growth rate is positive, µ(x)dx > 0. However, this condition is not necessary and does not hold even in the simple models considered before. General conditions of uniqueness and the existence of positive eigenvalues were obtained by Hess and Kato [29], Senn and Hess [92], Cantrell and Cosner [9]. The analysis of 2D and 3D models raises even more questions. In higher dimensions the persistence of a population may depend on the geometric form of the favourable patches [15], the form of the edges separating the patches and finally it may depend on the behaviour of individuals, moving across or along the edges [19].
n

3.2. Persistence on an infinite habitat
Travelling fronts While in the previous section we have highlighted some negative aspects of diffusion for the fate of a population, in this section we show that depending on the circumstances diffusion may as well support population growth. Maybe the most drastic example is the possibility to generate the spread or geographic expansion of a population into a new area. Consider, for simplicity, an infinite homogeneous habitat which provides a positive growth rate µ0 everywhere. If, as it is commonly assumed (i.e., there is no Allee effect), the stationary state P = 0 is unstable, the appearance of organisms in one spot will lead to their expansion over the whole habitat. Assuming logistic growth the spatio-temporal dynamics of the population can be represented by the following equation P (x, t) 2P P +D 2 . = µ0 P 1 - (3.18) t K x which is known as the Fisher-Kolmogorov equation after Fisher [20], who considered the logistic dynamics of advantageous genes and Kolmogorov et al. [48], who investigated a general form of this problem (see also [62, 51, 108] ). Fisher [20] and Kolmogorov et al. [48] have shown that if initially some part of the habitat is not occupied, P (x) = 0, then the population will propagate into this part with the constant velocity vf = 2 µ0 D . (3.19)

This solution is illustrated in Fig. 5A which shows a population propagating from left to right through a 1D habitat. As will be shown below, the magnitude of the propagation velocity is a 55


A. B. Ryabov and B. Blasius

The role of diffusion and advection

Figure 5: Front-propagation process in the Fisher-Kolmogorov model (3.18). (A) Propagation of a travelling front, if initially P (x, 0) = 0, for x > 0. (B) Propagation of a "pseudo-wave" in the limit of zero diffusivity, if initially P (x, 0) = 1/(x + 1)5 for x > 0. In both plots P (x, 0) = 1 for x 0. crucial factor not only for an invasion process but also for the survival of a population in the presence of drift, sinking or other advective processes [14, 99, 101]. There is a variety of ways to derive relation (3.19). The more common and rigorous approach suggests to assume a travelling solution of the form P (x - v t) and then to prove that this solution is stable only if v = vf (see e.g. [51]). Below we will provide another heuristic derivation confirming the validity of this expression. Note that in large aquatic basins the horizontal turbulent mixing increases with the scale of phenomena [73, 67, 66, 74]. Petrovskii [77, 78] showed that this should result in the increase of the front propagation velocity. Moreover, this velocity should grow with the size of the area occupied by a population. Pseudo waves We should stress a condition which is sometimes missed. The velocity vf is the minimal possible propagation velocity, which is realised if the population invades into an empty area or if, at least, in this area P (x, 0) decays more rapidly than a Gaussian distribution (see Fig. 5A). Another pattern of so-called "pseudo-waves" may occur for special initial conditions if from the start a small amount of biomass is distributed over the whole space. This is illustrated in Fig. 5B where the initial distribution of biomass is algebraically decaying for x > 0 and so is visually indistinguishable from that in Fig. 5A. As shown this initial configuration yields a much faster wave-like spread of the population. Strictly positive initial conditions lead to simultaneous logistic growth toward the carrying capacity at every position and, even in the limit of zero diffusivity, a wave-like pattern arises because the capacity is reached at slightly different times in different points and not due to the transport of biomass by diffusion. Pseudo-waves appear if the time scale of diffusive transport, D , is slower than the time difference dem of demographic processes in the neighbouring points D x2 x P < dem - , 2D µ (P )P x (3.20)

where dem can be found expanding the equation P (x, t) = P (x + x, t + dem ) into a series. It is easy to check, that for a Gaussian distribution of P (x, 0) both sides of (3.20) are proportional. 56


A. B. Ryabov and B. Blasius

The role of diffusion and advection

v P(x, t)

v

f

x
Figure 6: Schematic presentation illustrating the spread of a population into an oppositely directed flow of velocity v . Initially the 1D habitat is populated on the left, while the right is not occupied by the species. In a system without advection the population would propagate into the right with the front velocity vf , given by equation (3.19). This spread is hindered by the advective flow, which is aimed into the opposite direction. As a result the population front propagates with the reduced velocity vf = vf - v . The population can persist only if vf v , otherwise it is washed out. In contrast, for algebraic or exponential distribution of P (x, 0) the right hand side will be larger if x exceeds a threshold value. Consequently, the diffusive transport will be always slower than the time scale dem of the demographic process, giving rise to a "pseudo-wave". Advection Consider now an extension of the Fisher-Kolmogorov equation for a population which is additionally subjected to an advective flow with constant velocity v [63, 64, 14, 99, 101, 4, 50] P (x, t) = µ0 P t 1- P K -v P 2P +D 2 . x x (3.21)

To investigate the role of advection let us again suppose initial conditions, such that only the left part of the habitat is populated, while the remaining part is not occupied by the species (see Fig. 6). In the absence of advection in this system, the population would propagate into the right with the velocity vf . Now assume that this spread occurs in an advective flow, which is going into the opposite direction with the velocity v . This can be analysed best by considering a frame of reference moving parallel to the flow with the same velocity v , so that in this frame the flow velocity is zero. In the moving reference frame the population dynamics obey equation (3.18) and the front propagates with the velocity vf given by (3.19). Therefore, in the fixed reference frame, the propagation velocity is reduced by the velocity v of the advective flow, so that the front propagates with the velocity vf = vf - v . Obviously, the sign of the reduced propagation velocity determines the persistence of the species. A negative propagation velocity leads to a population wash-out, whereas a positive propagation velocity results in the invasion of the empty habitat. Thus, on an infinite habitat in an 57


A. B. Ryabov and B. Blasius

The role of diffusion and advection

advective flow the population can survive only if v vf = 2 µ0 D , or, written in a different form DD = (3.22)

v2 , (3.23) min 4µ 0 which means that in a flow a population can persist only if the diffusivity exceeds the threshold value Dmin . Equations (3.22) and (3.23) relate the critical values of diffusivity D and flow velocity v and illustrate the constructive interplay that can arise in the presence of both processes, advection and diffusion. The region in the (v , D )-coordinate plane, where a population can outgrow an advective flow, is visualized in Fig. 7 (area under the solid line). Note that if the habitat is unlimited then for any flow velocity v there is a diffusivity which provides the survival of the population. This transition plays an important role in many ecological situations and constitutes, for example, a necessary condition for the persistence of a population in a river [99], as well as for the persistence of sinking phytoplankton species in a vertical water column [87, 95, 33]. Derivation of the front propagation velocity These results about the spread of a population in an advective flow can be elegantly used to derive the population spread, Eq. (3.19), in a system without advection. For these aims we consider the general model

P 2P P (x, t) = µ(x)P - v +D 2 . (3.24) t x x Here, the local growth term has been linearised for small densities, however we allow for an arbitrary spatial dependence of the growth rate, µ = µ(x). This model can be simplified with the following transformation vx ~ P (x, t) = exp P (x, t) . (3.25) 2D The aim of this `trick' is to eliminate the advective term and as a result we obtain the following equation of a system without advection, however with a modified growth term ~ ~ P (x, t) 2P v2 ~ P +D 2 . = µ(x) - (3.26) t 4D x Note that the transformation (3.25) should be applied also to any boundary conditions of the original system and may alter them. However, if we can neglect the influence of boundaries (e.g. if the habitat is infinite, but still heterogeneous) advection obtains a very simple ecological interpretation as an additional mortality of strength v 2 /4D . This means that the presence of advection effectively reduces the growth rate µ(x), and this effect increases as D 0. Now we can "derive" the velocity vf of the front propagation in the Fisher-Kolmogorov equation. On a homogeneous infinite habitat (µ(x) = µ0 ) the population can survive only if the growth rate is positive. Thus, using (3.26) the velocity should be constrained However, as we showed before, the maximal possible advection velocity, v , is equal to the front propagation velocity, vf , which finally leads us to formula (3.19). 58 v2 µ0 D .


A. B. Ryabov and B. Blasius

The role of diffusion and advection

3.3. Finite habitats in an advective flow: the "drift paradox"
While the Fisher-Kolmogorov equation (3.18) assumes an infinite homogeneous habitat, such conditions are a strong idealisation for most natural populations. To describe some more realistic situations, in the following we investigate the role of advection for a population in a heterogeneous environment which is additionally constrained by certain boundary conditions. Separation of variables Consider again the general model (3.24) of the previous section, which now is supposed to be complemented by some boundary conditions. Again we use the change of variables (3.25) to eliminate the advection term. The linear form of equation (3.26) allows ~ a separation of variables P (x, t) = X (x) T (t) and by comparison to Eq. (3.17) we obtain the eigenvalue problem for the time-independent eigenfunctions D X (x) + µ(x)X (x) = v + v2 4D X (x) , (3.27)

while the equation for T (t) remains unchanged. Here, v denotes the eigenvalues for the problem in the presence of advection. Furthermore, X (x) should satisfy certain boundary conditions, which ~ depend on the boundary conditions obtained for P (x, t). If we introduce = v + v 2 /4D this equation will take the form of the second equation of system (3.17) for the problem without flux. Thus we can conclude that the presence of advection for a population under boundary conditions simply reduces the eigenvalues [24, 64, 14] v = - v2 . 4D

With the same arguments as in Section 3.1., to provide the persistence of a population, the largest eigenvalue v must be positive. Therefore, a population is only able to persist in a flow if the proper model without flow has the principal eigenvalue v2 . 4D (3.28)

KiSS model with advection As a simplest example, consider the KiSS model (3.7) in the presence of an advective flow. Recall, that this model describes a population on a finite favourable patch provided that the population density vanishes at the boundaries. These boundary conditions (P (0, t) = P (L, t) = 0) are not changed by the coordinate transform (3.25), so that the functions ~ P (x, t) and X (x) also must vanish at the boundaries. Therefore, we will obtain the same expression (3.11) for the eigenvalue spectrum with the dominant eigenvalue = µ0 - 2 D /L2 . Taking into account (3.28), we can easily derive the conditions for the persistence of a population on a finite favourable patch in an advective flow 2D v vf = 2 D = 2 D µ0 - 2 L 59 . (3.29)


A. B. Ryabov and B. Blasius

The role of diffusion and advection

This expression describes a semi-circle in the (D , v ) parameter plane (see dashed line in Fig. 7). Note that as L approaches infinity or D approaches zero, the maximal velocity (3.29) goes over to the condition (3.22). This means that in the limit of small diffusivity (D µL2 / 2 ) the conditions for survival on a finite patch are almost the same as those on a homogeneous infinite habitat. Thus, for small D the two critical curves in Fig. 7 almost coincide and show a similar behaviour. For larger diffusivity, however, these curves diverge because a finite habitat provides a slower propagation velocity (3.29) than an infinite habitat (3.22). Finally, for large values of D the increase of the losses across the habitat edges results in the upper diffusivity limit Dmax on a finite habitat. Solving this inequality for D , we obtain a limiting interval Dmin (v ) D Dmax (v ), with D
min/max

=

D0 2

1±

1-

2v2 L2 µ2 0

.

(3.30)

Here, D0 is the critical (maximal) diffusivity (3.13) in the KiSS model. The interval [Dmin , Dmax ] specifies the limits of mixing intensity which prevent the population wash-out (D > Dmin ), but still enable the persistence of the population on a finite patch (D < Dmax ). These values are real onl y i f Lµ v vmax = . (3.31) Note that the critical velocity vmax is a threshold when the characteristic time scale of growth gr = 1/µ becomes slower than the time scale of advection v = L/v . If v > vmax the population cannot persist, because on the one hand the large advection requires a strong mixing intensity to provide the expansion of organisms upstream, but on the other hand such mixing increases the transport of organisms into unfavourable areas and the population becomes extinct. Critical patch size In the following we investigate the influence of advection on the critical patch size in the KiSS model (3.7) and in the model by Ludwig et al. (3.15). In the first model the population vanishes at the habitat edges, whereas in the second the population density should vanish when x ±. Therefore, in both models the transformation (3.25) does not alter the boundary conditions and the presence of an advective flow simply reduces the growth rate µ. This leads to the effective growth rate v2 µv = µ0 - 0 4D in the KiSS model and to 2 µ-v , for 0 < x < L 0 4D µv (x) = 2 - m+ v , for x 0 or x L , 4D in the model by Ludwig et al. 60


A. B. Ryabov and B. Blasius

The role of diffusion and advection

Figure 7: Parameter regions which permit the persistence of a population in an advective flow. Persistence of the population is possible below the critical curves in the (D , v )-parameter plane. Two different scenarios are compared. The solid line shows the persistence regime for an infinite uniform habitat (cross and diag al hatching). In this case the limiting velocity scales as the on square root of diffusivity v D , equation (3.22). In comparison, the dashed line shows the result for the persistence of a population on a finite favourable patch with lethal boundaries (cross hatching). Here, the persistence regime yields a semi-circle in the parameter plane, equation (3.30). Assuming intermediate growth conditions, e.g. similar to the model by Ludwig et al. (3.15), one should expect the critical curve to be located somewhere in-between the solid and the dashed line. Substituting the modified growth rate µv into equation (3.26), we find similar expressions for the critical patch size as in Eqs. (3.12) and (3.16), where µ should be replaced by µv . Thus, for the KiSS model in an advective flow the critical patch size equals Lv = 0 D = µv 0 D = µ 0 - v 2 /4D 2 D
2 vf - v 2

,

(3.32)

and for the model by Ludwig et al. we obtain Lv = 4D
2 vf

arctan
2

where vf = 2 D µ0 . Note that both Lv and Lv approach infinity as v vf . In particular, 0 this can happen if D Dmin = v 2 /4µ0 which is of importance in marine biology as climate models predict that the ongoing global warming may result in a higher stratification of the ocean water [6, 91], increasing thereby the requirement on the critical (vertical) patch size for sinking phytoplankton species. 61

-v

4mD + v 2 , 2 vf - v 2

(3.33)


A. B. Ryabov and B. Blasius

The role of diffusion and advection

Persistence in a river Speirs and Gurney [99] investigated the conditions for the population survival in a river of length L and with a flow velocity v . This problem is known as the "drift paradox" because any advection will ensure that the average location of a population will move downstream, so at first glance it seems counter-intuitive that a population can persist in a river [28]. In their study Speirs and Gurney used the linear model (3.24) and assumed an impenetrable boundary upstream and a totally hostile environment downstream. In this model the critical size Lv of the favourable patch depends on the ratio between the flux velocity v and the front propagation velocity vf L=
v

2D
2 vf - v 2

arctan

2 vf - v 2 . v2

Note that Lv again increases with an increase of the advection velocity v and approaches infinity as v vf . Furthermore, because the favourable patch in this model has only one boundary with the hostile environment, we obtain the limit Lv L0 /2 as v 0, where L0 is the critical patch size (3.12) of the KiSS model. Extending this model, Pachepsky et al. [75] derived conditions for the persistence and spread of a population of organisms living and reproducing on the sediment and occasionally entering the water flow where they can drift and disperse. Locally elevated growth rate Dahmen et al. [14] considered another extension of the model by Ludwig et al. (3.15), assuming an advective flow and periodic boundary conditions. Furthermore they suggested a simple experimental set-up in which a light limited colony of bacteria grows on a ring. Almost the whole ring is shaded and unfavourable for the population. Only a small area is illuminated through a window of length L, which moves around the ring with a constant velocity v . Up to a change of reference frame, this set-up is equivalent to an advective flow and, changing the speed v , one can easily regulate the "advection" rate. Depending on the parameters and the velocity of the light supply, the population can become extinct, can be localised (the maximum of density tracks the location of the favourable patch), or delocalised [64] (due to the periodic boundary conditions the population can survive even if the advection velocity v exceeds the propagation velocity vf provided that the average growth rate is positive). Fig. 8 shows an example for such a localisation. Solving the eigenvalue problem Dahmen et al. exploited a similarity of equation (3.27) to the well studied "square well potential" problem of quantum mechanics [53]. They showed that depending on the dimensionless parameter x = (2/L) D /(µ0 + m), which characterises the growth rate in relation to the diffusivity and the habitat length, the critical flow velocity can be

62


A. B. Ryabov and B. Blasius

The role of diffusion and advection

Figure 8: Locally elevated growth rate provides the persistence of a population in an advective flow with velocity v , Dahmen et al. [14]. A favourable patch (µ(x) = µ0 > 0) is indicated by the horizontal thick line. expressed as vc = 2 D D µ0 - 2 L
2

,

if x 1 (3.34)
2

vc = 2D

(µ0 + m)L 2D

-

m D

,

if x 1

where all parameters have the same meaning as in (3.15). The first solution also corresponds to the limit of high mortality (m ) and coincides with equation (3.29), which was derived for the KiSS model with advection. The second solution describes a strongly mixed system. If the advection velocity is higher, the population becomes extinct. The value of the critical patch size and the critical diffusivity can be easily expressed from equation (3.34). Examining bacterial growth, Lin et al. [55] confirmed the results of Dahmen et al. [14] experimentally and by means of numerical simulation. Joo and Lebowitz [43] carried out computer simulations in a stochastic spatially discrete population model and obtained similar results, confirming the robustness of the model. Locally elevated diffusivity Consider a population growing on an infinite favourable patch in an advective flow. If the diffusion is too low this patch will not provide a proper propagation velocity (3.19) and the population will become extinct. Straube and Pikovsky [101] noted that locally increased diffusivity can drastically change the situation. A sufficiently large and well mixed patch can stop the drift of biomass, stabilising the population dynamics downstream (Fig. 9). Straube and Pikovsky considered the Fisher-Kolmogorov equation with advection (3.21), assuming that a patch of length L has an elevated level of diffusivity D1 , which exceeds the diffusivity D in the rest 63


A. B. Ryabov and B. Blasius

The role of diffusion and advection

Figure 9: Locally increased diffusivity stops the population wash-out in an advective flow with velocity v , according to Straube and Pikovsky [101]. The horizontal thick black line indicates the location of the patch with increased diffusivity D1 . of the habitat. Note that if the diffusivity depends on the coordinate, D = D (x), equation (3.21) takes the form P P P P (x, t) -v D (x) . = µ0 P 1 - + t K x x x Straube and Pikovsky [101] showed that the population survives, if the size of the intensively mixed patch is larger than 2 v 2 - vf 2 D1 Lcr = arctan . 4 µ 0 D1 - v 2 4 µ 0 D1 - v 2 As the diffusivity D1 approaches infinity the critical patch size approaches lim Lcr =
2 v 2 - vf

D1

2µ

.

0

Thus, even for the infinite mixing intensity the critical patch should be of finite size. Techniques from quantum mechanics Birch et al. [4] considered the Fisher-Kolmogorov equation with variable growth rate and advection on a 2D plane. In this study they made use of the structural similarity between equation (3.27) and the time-independent Schrodinger equation. Based on ¨ this similarity, they demonstrated a few examples where the application of perturbation theory, the method of Wentzel, Kramers, and Brillouin (WKB), and other quantum mechanical techniques are beneficial for the analysis of equations similar to (3.27). In particular, Birch et al. were able to determine a trade-off between the critical diffusivity and growth rate, which provide the persistence of the population in this system. We believe that the application of such well-known techniques from quantum mechanics has not been fully exhausted yet. Such methods could be a promising direction for the further development of the theory of extended populations. 64


A. B. Ryabov and B. Blasius

The role of diffusion and advection

Summary At the end of this section we want to summarize the main results. On an infinite homogeneous patch a stationary state is achieved when the growth rate is equal to the mortality. On a finite patch, which adjoins some unfavourable environment, the growth rate should exceed mortality to compensate for losses across the patch edges. The internal patch production grows with the patch size and with the growth rate, whereas the losses at the edges increase with the diffusivity. The spread of organisms over a habitat is defined by the propagation velocity of the front of the population density. This velocity scales as a square root of the growth rate and diffusivity. A population can survive in a flow only if the propagation velocity is higher than the velocity of the flow. Thus, an increase of the growth rate improves the conditions for the survival on a finite favourable patch and in a flow. The increase of the patch size increases the internal patch production, whereas the losses across the patch edges show only a weak dependence. The critical patch size is achieved, when the internal production is balanced by the external losses. The size of a large habitat has a small influence of the propagation velocity. However, an unfavourable environment around a small favourable patch truncates the propagating front, this can essentially reduce the propagation velocity and worsen the conditions in a flow. Finally, an increase of diffusivity on a finite habitat increases the losses into the unfavourable environment, which can lead to the population extinction. However, a minimal level of diffusivity is necessary to prevent the population wash-out in a flow. These two opposing processes result in a diffusivity window, which provides the population survival on a finite habitat. This window exists only if the flow velocity is less than some critical value.

4. Vertical phytoplankton distribution
In the previous chapter, using simple one-dimensional models, we discussed the influence of mixing and advection on the population survival in a heterogeneous environment. To model a nonuniform environment, we simply used an explicit form for the growth rate, µ(x), however we did not discuss the origin of this heterogeneity. In this chapter we extend these results by considering more complex models which describe resource limited population growth. These models are based on simple physical and biological laws and consistently describe the dynamics of a population and its limiting resources. Thus, they demonstrate natural mechanisms which lead to the appearance of favourable patches. The intention of this chapter is twofold: on the one hand we aim to illustrate how the main findings of the previous chapter apply to the context of consumer-resource models. On the other hand we show new effects arising due to new properties, which are not present in the simple models. The main object, which we will use for illustration, is the dynamics of a vertical phytoplankton distribution. This is important as phytoplankton are the primary producers in almost all aquatic food webs with a major influence on nearly all freshwater and marine ecosystems. The two main factors limiting the production of phytoplankton are the availability of nutrients and light. To understand how these resources affect the phytoplankton biomass, consider their distributions in a water column. In general, the light intensity reduces with depth and, in nutrient rich regions of the

65


A. B. Ryabov and B. Blasius

The role of diffusion and advection

ocean, the well illuminated surface layer constitutes a favourable area for photosynthetic phytoplankton species. By contrast, the nutrient concentration can behave just in the opposing way. The sedimentation of dead biomass (detritus), with the successive remineralisation in the deep layers or in the sediment [16] causes an increase of the nutrient concentration with depth [117]. Thus, while light limitation may lead to the formation of a surface phytoplankton maximum, a lower nutrient concentration favours a phytoplankton build-up in deeper layers. This tension between light and nutrient limitation from two opposite sides frequently causes optimal growth conditions in subsurface layers (e.g. [1, 13, 31]). This fact often leads to the appearance of maxima of chlorophyll or biomass distributions at approximately 30-100 m depth. So called deep chlorophyll maxima (DCM) [1, 111, 13, 110, 31] and deep biomass maxima (DBM) [49, 5] are ubiquitous phenomena and can be observed in many oligotrophic regions in the ocean, marine systems, and deep lakes. Another important component of a stratified water column is an upper mixed layer (UML). A UML commonly occurs in oceans and lakes due to mechanical perturbation of the surface waters (e.g. due to wind, waves, and storms). This layer is separated from the deep layers by a thermocline [12], which is defined as a relatively thin layer below a UML characterised by an strong change in temperature with depth. Mixing in a UML is much stronger than in the layers below it. As a result, the distributions of nutrients, temperature, salinity, etc. are nearly uniform in a UML and have gradients below it. The depths of a UML can usually vary from 10 m to 100 m, see e.g. [110]. Compared to the models of the previous chapter, the system behaviour in a water column can be further complicated due to a feedback loop between the biomass and resource distributions. The growing biomass shades light, consumes nutrients and is remineralised, which ultimately changes the total resource distribution. This, in turn, can lead to a new biomass distribution, which will generate a new resource profile and so on. As will be shown below, these complicated, selforganised dynamics can lead to new phenomena and diverse behaviour. For example, if the mixing is small, the final solution becomes non-stationary and oscillates [36], whereas in the presence of an upper mixed layer the system may exhibit bistability and the solution may be sensitive to the initial conditions [118, 89]. Equation of growth To formulate a mathematical framework for this chapter, let us consider a vertical water column of depth ZB . Let P (z , t) denote the density of phytoplankton at time t and depth z . Note that z = 0 denotes the sea level surface and the z -axis is directed downward. For the sake of simplicity, assume that phytoplankton growth is limited only by the availability of light and a nutrient (the model can easily be extended to take into account multiple nutrient limitation [26]). In our approximation the dynamics of a phytoplankton population obey a reaction-diffusionadvection equation, similar to equation (3.1) considered in the previous chapter (see [85, 95, 46, 36] among others) P (z , t) P P = µ(N , I )P - mP - v +D , (4.1) t z z z where m is the mortality (compare to equation (2.3)), v is the phytoplankton sinking velocity, and D is the diffusivity, which in general can depend on z . Furthermore, the growth rate µ(N , I ) depends on the local values of light intensity I (z , t) and nutrient concentration N (z , t) at each vertical position. If both nutrients are essential, µ(N , I ) can 66


A. B. Ryabov and B. Blasius

The role of diffusion and advection

be represented in the form of Liebig's law of minima µ(N , I ) = µ0 min [fI (I ), fN (N )] . It can be also written in the multiplicative form µ(N , I ) = µ0 fI (I ) fN (N ) . (4.3) (4.2)

Here µ0 is the maximum growth rate and fI (I ) and fN (N ) describe the limitation by light and the nutrient. The specific form of these functions depends on many factors, for instance, strong light may photoinhibit photosynthesis and reduce the growth rate for large values of I [83, 22, 86]. However, usually it is suggested that f (x) 1 as x , that is, the maximum growth rate is achieved when all resources are unlimited. For phytoplankton modelling, the most frequently used form is the Monod (or Michaelis-Menten) kinetics [107] fI (I ) = I , HI + I fN (N ) = N , HN + N (4.4)

where HI and HN are the half-saturation constants for nutrient-limited and light-limited growth, respectively. However, this non-linear form often admits only numerical investigation. Analytical solutions are commonly possible only for a linear or algebraic form of f (x) [95, 17]. Boundary conditions phytoplankton By default, we assume that the surface and bottom are impenetrable for v P (z , t) - D To mo and below As s u m i n g bot t om of P z =0.
z =0,Z
B

(4.5)

del a stratified water column, one can either separately solve the equations in a UML it, supposing infinite mixing within the UML and a small diffusivity DD in deep layers. continuity of the flux across the thermocline, we obtain the boundary condition at the a UML (see e.g. [35]) v P (z )|
z =Z T -0

=

v P (z ) - D

D

P z

,
z =Z T +0

where ZT is the depth of the thermocline. On the other hand, to simulate the water column in a single framework, one can assume a gradual transition from a UML to the deep layers [89] D ( z ) = DD + DU - D D , 1 + e(z -ZT )/w (4.6)

where DU and DD are the diffusivities within and below a UML, respectively, and the parameter w characterizes the width of the thermocline. So far, we did not specify any equations for the distribution of nutrients N (z , t) and light I (z , t). In the next two sections we will review several models which couple the light and nutrient dynamics with phytoplankton growth (4.1). First we will consider models in which phytoplankton growth is limited only by the light availability, whereas in a second step, we will review more complex models, incorporating both light and nutrient limitation. 67


A. B. Ryabov and B. Blasius

The role of diffusion and advection

4.1. Light limitation
In this section we will consider theoretical models in which the light gradient is a key factor. Such models give an adequate description for eutrophic aquatic environments (as observed in many regions), where the nutrients are in ample supply and light becomes a crucial factor which determines the distribution and the dynamics of phytoplankton [59, 90, 11]. So in the following we assume that the nutrient dependence of the growth rate is saturated, fN (N ) 1, and we can neglect the limitation of growth by nutrients. The spatial profile of light intensity in a water column is described by Lambert-Beer's law (see e.g. [45]) which states that the gradient of light intensity at depth z is proportional to the light intensity at this depth dI = -I . (4.7) dz The coefficient includes both the absorption of light by water and the attenuation by the phytoplankton cells = Kbg + k P (z , t) , where Kbg is the background turbidity and k is the per capita attenuation coefficient of the algae cells. Integrating (4.7) from surface to depth z , we obtain
z

I (z ) = Iin exp -Kbg z -

k P (t, z )dz
0



,

(4.8)

where Iin is the light intensity at surface. Equations (4.1) and (4.8), being coupled by means of the growth rate (4.2) or (4.3), yield an integro-differential system of equations. It is not straightforward to obtain rigorous or analytical results for such a system and even a numerical solution encounters certain difficulties [35, 103]. Nevertheless, without solving any equations, it is clear that the light intensity in the water column is reduced with increasing depth. Thus the light limitation forms a favourable area close to the surface, and the dynamics of the phytoplankton population should be related with the results of the previous chapter, obtained for heterogeneous environments in the presence of advection and diffusion. Critical values for phytoplankton growth Depending on the depth of a water column or on the diffusivity, a light limited phytoplankton population can survive or become extinct. Huisman et al. [33, 35] combined the conditions for survival into a single conception of the critical conditions for phytoplankton blooming in a closed water column (Fig. 10A) and within an UML (Fig. 10B). The main difference between these models is that in a closed system, the sinking of cells is stopped at the bottom, whereas in an UML the biomass can sink across the thermocline to the deep aphotic layers. To determine these conditions, consider first an unstratified water column with a constant diffusivity and assume impenetrable boundary conditions (4.5) for the phytoplankton biomass. Using the system of equations (4.1) and (4.8) in the limit of zero background turbidity, Kbg = 0, and for 68


A. B. Ryabov and B. Blasius

The role of diffusion and advection

Figure 10: Critical conditions for phytoplankton blooming: (A) in a closed water column and (B) in an upper mixed layer. Both figures are reprinted from Journal of Sea Research, 48, Huisman and Sommeijer, pp. 83-96 [35], Fig. 4 and 6, with permission of the author. a general monotonic growth rate µ(I ), Shigesada and Okubo [95] showed that a sinking phytoplankton species can establish a population only if DD
min

=

v

2

4(µ(Iin ) - m)

.

(4.9)

This expression coincides with condition (3.23) from the previous chapter and implies that the minimal diffusivity should provide a front propagation velocity which is larger than the sinking velocity. The same condition was derived earlier by Riley [87] and other authors from the FisherKolmogorov equation. It is interesting to note [95] that if Kbg = 0 and a non-trivial solution exists, then the total biomass in this model does not depend on the sinking velocity. The sinking just shifts the bulk of biomass downward, preserving, however, the total amount of biomass in the water column (Fig. 11A). Ishii and Takagi [40] relaxed the condition Kbg = 0 and proved some existence, stability and uniqueness results for this system. Assuming an algebraic form of the growth rate, µ(I ) I , Ebert et al. [17] have found some approximation for the minimal diffusivity, Dmin , and for other critical parameters. If a water column is sufficiently deep and Kbg = 0 then the net production rate is positive only above the compensation depth, Zc , which is defined as the depth z at which the local production rate is zero in the absence of biomass (Fig. 11B). From the light attenuation curve (4.8) we find Zc = ln Iin - ln Ic , Kbg

where Ic is defined as the compensation light intensity at which the growth rate is equal to mortality, µ(I ) = m, thereby the compensation depth is species specific. 69


A. B. Ryabov and B. Blasius

The role of diffusion and advection

A

B

m µ(I)

I(z)

favourable

Depth

Zc Zcr unfavourable

Figure 11: (A) Vertical profiles of phytoplankton for different values of velocity = v /µ0 D [95]. While the shape of the profile is changed, th in the system remains unchanged. The figure is reprinted from Journal 12, Shigesada and Okubo, pp. 311-326 [95], Fig. 4, with permission of representation of the compensation depth Zc and the critical depth Zcr .

the dimensionless sinking e total amount of biomass of Mathematical Biology, the author. (B) Schematic

If a water column is shallower than the compensation depth and we assume that the bottom is impenetrable for the biomass, then the population can survive even if mixing is less than a minimal diffusivity (4.9) because the cell settling will be stopped at the bottom (Fig. 10A). In general, the compensation depth divides a water column into a favourable and an unfavourable regime (Fig. 11B). In a well mixed water column the losses in the deep layers can lead to the population extinction. However, they can be compensated by the production in the euphotic zone, if the unfavourable region is relatively small. Considering a simple mathematical model of a well mixed water column Sverdrup [102] defined a critical depth as the depth of a water column at which the total growth is equal to the total loss of biomass. Similar to the compensation depth, the critical depth can be reinterpreted in terms of the critical light intensity [37] Zc r = ln Iin - ln Iout , Kbg

where Iout is the light intensity at the bottom of a sufficiently shallow (ZB < Zcr ) closed water column after the light limited population of phytoplankton have reached an equilibrium state. The critical depth depends on many parameters, it increases with the incident light intensity and with the phytoplankton growth rate, and it decreases with the mortality rate [37]. In a well mixed water column, an excess of the critical depth over the compensation depth determines the maximal possible losses in dark layers, which can be still compensated by the production in the euphotic zone. However, similar to the KiSS model or to the model by Ludwig et al. (see Sec. 3.1.) these losses diminish with a decrease of mixing. Extending this research, Huisman et al. [37] showed that if the depth of a water column or a thermocline exceeds the critical depth, the population survival still is possible if turbulent mixing is less than a maximal diffusivity. 70


A. B. Ryabov and B. Blasius

The role of diffusion and advection

This critical condition is similar to the existence of the maximal diffusivity (3.13) and (3.30) in the KiSS model, however here it describes the behaviour of a more realistic system. Therefore similar to the KiSS model in advection (3.30), a sinking population can survive in a water column of any depth if mixing remains between a minimal and a maximal value (Fig. 10A). We now turn to a stratified water column with a UML. We assume that the mixing in the deep layers is less than a critical value (4.9), so that the population survival will depend on the characteristics of a UML. Condie and Bormans [12] showed that if a UML is shallower than Z
T ,min



v , µ(Iin ) - m

a population cannot survive (compare with (3.31)). In other words, for the survival in a UML, the demographic time scale should be faster than the characteristic time of advection. However usually, ZT ,min is sufficiently small and this criterion is satisfied. In this case the population can persist if the strength of mixing remains within the turbulent window [Dmin , Dmax ] (see Fig. 10B). Furthermore, if the diffusivity exceeds a maximal value the population survives if the depth of a thermocline is smaller than a maximal depth, which is defined as the maximal depth of a well mixed upper layer at which losses and production are equal. This depth is slightly smaller than the critical depth in a closed water column, owing to additional losses of biomass across the thermocline. The fact that a deep upper layer can prevent phytoplankton blooming was noted experimentally in 1935 by Gran and Braarud [23], who investigated the conditions of phytoplankton blooming in the upper mixed layer. They reported that until there exists a deep UML, phytoplankton production cannot exceed the destruction by respiration and phytoplankton blooming is not possible. The concept of the maximal diffusivity is also consistent with field experiments, see e.g. [106, 18, 76].

4.2. Light and Nutrient limitation
In the last section we will discuss models which take into account both light and nutrient limitation of phytoplankton growth. These models are more difficult to analyse and often admit only numerical investigation. However, they are more realistic and provide some understanding of the processes occurring in deep waters of many regions where surface layers are nutrient depleted [85, 113, 82, 115, 112, 109, 36]. Furthermore, in the tension of two opposing resource gradients the location and the size of a production layer becomes a function of the phytoplankton abundance and the initial conditions, that can lead to new patterns and new dynamical behaviour. A coupled system of reaction-diffusion equations describing nutrient-phytoplankton cycling was probably first investigated by Okubo [69, 71]. Radach and Maier-Reimer [85] suggested a mathematical model of phytoplankton growth which included light-nutrient-phytoplankton dynamics. This model was extended by Jamart et al. [41] who considered limitation by two nutrients, grazing and the variability of the parameters with depth and time. This approach (see also [113, 82, 115, 112]) gave rise to a growing set of ecological models, which include cycling of many chemicals [117], coupling with meteorological data [42], interplay of different phytoplankton groups, and 3D simulations [60, 21]. Here, however, we will focus on the theoretical aspects and consider only simple conceptual models. 71


A. B. Ryabov and B. Blasius

The role of diffusion and advection

Conservative models The nutrient dynamics include uptake by phytoplankton, remineralisation of dead biomass back into a nutrient pool and diffusion. Assuming absolutely effective recycling we obtain P 2P P (z , t) = µ(N , I )P - mP - v +D 2 , t z z (4.10) N (z , t) N = -µ(N , I )P + mP + D 2 , t z where the biomass is measured in terms of its nutrient content (compare to the non-spatial version (2.3) of this model). We do not include advection in the second equation, as nutrients, which are dissolved in water, are only slightly influenced by the gravity force. Nevertheless, this term should appear, if advection is caused by a vertical or horizontal stream. Furthermore, we assume that the nutrient cannot diffuse across the surface and a large nutrient pool in the sediment or in deep ocean layers sustains a constant concentration, NB , the bottom of the water column N (0, t) =0, N (ZB , t) = NB . (4.11) z Fig. 12 shows typical final distribution of phytoplankton and nutrient given by model (4.10), supplemented by equation (4.8) for light. Hodges and Rudnick [30] pointed out that, independent of the functional form of the growth rate and of the light distribution (assuming that light decreases with depth), this model can reproduce a deep stationary phytoplankton maximum only if v > 0. In other words, the presence of opposing resource gradients is not sufficient to reproduce a deep phytoplankton maximum. To prove this, let us define the total concentration of the nutrient as S = P + N . Consider an equilibrium state, when the left-hand-side of (4.10) equals zero. By adding both equations (4.10) we obtain D 2S P -v =0. 2 z z
2

Assuming v = 0 and integrating this equation over z we find S P N = const = + =0, z z z owing to the boundary condition (4.5) and (4.11) at the surface. Thus S = N + P = const and a deep phytoplankton maximum should be accompanied by a deep minimum of nutrient. However, if the light intensity reduces with depth, this profile is unstable because there is no factor limiting phytoplankton growth in the upper layer. Thus this system should exhibit a surface maximum (Fig. 12A). However, similar to the model without nutrient limitation (Fig. 11A), the phytoplankton sinking shifts the maximum of biomass downwards (Fig. 12B). Extending this model, Hodges and Rudnick [30] included a detrital pool, T (z , t), as the third

72


A. B. Ryabov and B. Blasius

The role of diffusion and advection

A

B

Figure 12: (Colour online) Typical distributions of phytoplankton (green), nutrient (blue), and light (red) in the conservative model (4.10) without sinking (A) and with sinking (B), according to Hodges and Rudnick [30]. compartment P (z , t) = t T (z , t) = t µ(N , I )P - mP -vP mP - r T -vT P 2P +D 2 , z z T 2T +D 2 , z z (4.12)

2N N (z , t) = -µ(N , I )P + r T +D 2 , t z where vP and vT are the sinking velocities of phytoplankton and detritus respectively. Note that usually detritus sinks much faster than phytoplankton [94, 84]. In this model the cycle of chemicals includes three stages: the transfer of biomass to detritus with mortality m, the remineralisation of detritus back into nutrients with remineralisation rate r , and finally the consumption of nutrient by biomass. While this model can exhibit deep maxima if v = 0, the change of phytoplankton concentration is very small and cannot represent real data. An apparent maximum can be observed only if one assumes sinking of detritus or phytoplankton. Hodges and Rudnick extended this statement to any number of compartments, which however do not include depth dependent parameters. Thus, sinking is a major component of this system. The sedimentation of organic matter removes the nutrient fixed in phytoplankton cells from the upper layer, which leads to the formation of deep phytoplankton maxima. Beckmann and Hense [3] performed numerical simulations and analytical evaluations of model (4.12), assuming that detritus sinks relatively fast, whereas the phytoplankton sinking is negligible. Fig. 13 reproduces a typical distribution of physical characteristics in this model. Furthermore, Beckmann and Hense suggested to extend the concept of compensation depths. Instead of the static definition in the absence of biomass they suggested to use two dynamical depths at which the in situ production rate of phytoplankton is zero, owing to the light or nutrient limitation. If phytoplankton sinking velocity is zero then in equilibrium (see (4.12)) these values can be expressed from the 73


A. B. Ryabov and B. Blasius

The role of diffusion and advection

Figure 13: (Colour online) Typical distributions in the three compartments model (4.12). (A) Vertical profiles of phytoplankton (P ), detritus (D ), nutrient (N ), and light (I ). (B) The upper and lower limits of the production layer (dashed lines) are the species compensation depths. The figures are reprinted from Progress in Oceanography, 75, Beckmann and Hense, pp. 771-796 [3], Fig. 2, with permission of the author. phytoplankton distribution 2P z2 where Z
(N ) c

=0,
Z
(N ) c

2P z2

=0,
Z
(I ) c

and Z

(I ) c

are the compensation depths due to nutrient and light limitation (Fig. 13B).

Non-conservative models Model (4.12) contains three reaction-diffusion equations and an equation for the light distribution. This makes further analysis difficult. However, since detritus sinks relatively fast [94, 84], we can simplify the model assuming that a part of the dead biomass is instantly remineralised in situ, whereas the rest sinks until it reaches the bottom and sustains a constant nutrient concentration at the bottom. Thus we obtain the following non-conservative system of equations P 2P P (z , t) = µ(N , I )P - mP - v +D 2 , t z z (4.13) N (z , t) N = -µ(N , I )P + mP + D 2 . t z 74
2


A. B. Ryabov and B. Blasius

The role of diffusion and advection

This model (compare to equation (2.7)) can reproduce deep phytoplankton maxima even if the sinking velocity is zero, owing to the fact that a part (1 - ) of the fixed nutrient is implicitly transferred from the upper layer to the bottom. Even though this model is non-conservative and has apparent disadvantages, it or similar models were successively applied to reproduce field data [41, 109, 36, 118]. However, we are not aware of any comparison of the two model classes (4.10) and (4.13), which might be interesting. In the case of zero sinking velocity, Klausmeier and Litchman [46] performed analytical calculations for model (4.13). Assuming that the phytoplankton distribution can be approximated by a Dirac -function and further that an infinitely small production layer should be located to balance the light and the nutrient limitation, Klausmeier and Litchman found an equation for the position Z of a deep maximum, which for boundary conditions (4.5) and (4.11) reads as µ0 D (NB - Nc ) ln (Iin /Ic ) Kbg - Z= , k k m(1 - )(ZB - Z ) where Nc and Ic are the critical values of light and nutrient intensity for which the growth rate is equal to the mortality rate. Oscillations and chaos Huisman et al. [36] pointed out that system (4.13) exhibits oscillations of biomass if the mixing is reduced below a critical value. Fig. 14 shows the behaviour of biomass and of nutrient in two typical cases. In the first case (Fig. 14a) the mixing intensity is high enough to provide a stable distribution of biomass. If however the level of diffusivity is reduced, then only oscillatory, or even chaotic patterns, can appear (see Fig. 14b and 14c). As noted by Huisman et al. these oscillations are caused by the difference in the time scales of the rapid transport of phytoplankton, consuming the nutrients, and the slow upward transport of nutrients. Furthermore, as shown in Section 3.3., for the survival of a population in an advective flux the diffusivity should exceed a minimal level (3.30), which increases with the reduction of the habitat and of the growth rate. In the absence of biomass the nutrient can be nearly uniformly distributed over the water column, thereby the growth rate becomes only light limited and the production layer extends from the surface to the compensation depth, which is usually sufficiently large. Thus, without biomass, the level of mixing might be sufficient to induce population growth. However, the consumption of nutrients and self-shading of light reduce both the growth rate and the width of the production layer. That, in turn, increases the value of the minimal diffusivity (3.30) and finally the sinking may lead to the population wash-out if the diffusivity in the water column becomes insufficient. Koszalka et al. [50] noted that the periodic oscillations of phytoplankton biomass will be most probably disguised by currents and horizontal inhomogeneity in a real ecosystem. Upper mixed layer Hodges and Rudnik [30] and Beckmann and Hense [3] showed that if selfshading of light can be neglected in equation (4.8), then an upper mixed layer does not lead to any qualitative changes in the system dynamics. However, Yoshiyama and Nakajima [118] pointed out that a UML can lead to bistability of phytoplankton profiles. Ryabov et al. [89] generalised this result by taking into account the competition of two species and relaxing other assumptions. They considered the model (4.13), assuming a gradual change 75


A. B. Ryabov and B. Blasius

The role of diffusion and advection

Figure 14: (Colour on line) Evolution with time. (a) A stable DCM (D = 0.5 and (c) large-amplitude oscillations in figures are reprinted from Nature, 439 of the author.

of the phytoplankton density and the nutrient concentration cm2 /s), (b) small oscillations in the DCM (D = 0.2 cm2 /s), the DCM, with double periodicity (D = 0.12 cm2 /s). The , Huisman et al., pp. 322-325 [36], Fig. 2, with permission

of diffusivity (4.6) from a UML to the deep layers, and showed that under certain parameters, depending on the initial conditions the production layer can be steadily located either within a UML or below it. Fig. 15 provides a rough insight into the system dynamics. In the absence of an upper mixed layer the difference in the locations of biomass and of the production layer drives the bulk of biomass towards the production layer (Fig. 15A). The shift of biomass can lead to the redistribution of resources, which in turn can change the location of the production layer. This process repeats until the system reaches an equilibrium configuration (Fig. 15B), when the centre of biomass coincides with the centre of production. Now consider a system with an UML. In a certain range of parameter the UML does not affect distributions with a deep maximum of biomass (Fig. 15C). However, the initial growth of biomass within the UML begets another stable solution with a maximum of biomass located within the UML. The biomass is almost uniformly distributed within the UML and its location is uncoupled from the location of the production layer (Fig. 15D). As a result, a gradual shift of the bulk of biomass into deep layers is no longer possible and the transition to a deep biomass maximum can only take place if the light intensity below the UML is sufficiently large to provide positive net growth in deep layers ­ otherwise the phytoplankton

76


A. B. Ryabov and B. Blasius

The role of diffusion and advection

Figure 15: Typical vertical phytoplankton profiles P (z ) in a system without a UML (top) and with a UML (bottom), assuming a gradual change of diffusivity (4.6) in model (4.13). Without a UML, a non-stable phytoplankton distribution (A) evolves to a single stable solution (B). Under the same conditions in the system with a UML, we observe two stable distributions: with a maximum in the deep layers (C) and with a maximum in the UML (D). The dot-dashed and dashed lines show the limitations of growth (4.4) by light and by nutrient, respectively, vertical dotted line shows the level of mortality. Black and grey arrows show the centers of biomass and net production, respectively. remains trapped in the UML. Thus the production layer can occupy different parts of the water column, depending on the current system state and on initial conditions.

5 . D is c u s s io n
Concluding this review we would like to make a few notes. First, let us compare the behaviour of the critical patch models and those based on the consumer-resource dynamics. The latter models can be divided into two large groups. In the first group we would include those systems in which the location of a favourable patch is constrained by some environmental conditions. For instance, the limitation of phytoplankton growth by light leads to the formation of the favourable patch in the upper level of a water column. The dynamics of this group and of the critical patch models demonstrate many general traits and many effects can by predicted and evaluated on the basis of the minimal models. The second group consists of the models in which the location of the favourable patch is determined by the dynamical interplay of different factors. For example, we can consider the growth of phytoplankton biomass driven by two opposing resource gradients. In this group, the location of the favourable patch is not predefined. Moreover, the system dynamics becomes

77


A. B. Ryabov and B. Blasius

The role of diffusion and advection

very sensitive to the implementation of the consumer-resource cycling. This complexity leads to the arising of new patterns and new dynamical behaviour, which can hardly be reproduced in the framework of the critical patch models. The second remark concerns the advantages and disadvantages of partial differential equations (PDEs) for modelling ecological systems. PDEs provide very convenient and powerful tools for the investigation of population dynamics. First, in the same framework, we can consider such different and complex phenomena as, for instance, the vertical distribution of sinking phytoplankton cells or the survival of a population drifting in a flow. Second, analytical solutions in many cases provide important predictions and understanding of the main effects, which can appear in more realistic systems. Third, one can perform an exhaustive numerical simulation of a model, determining all possible bifurcation points. Finally, seemingly the pool of methods developed for the analysis of partial differential equations is not played out yet and this approach can still gain a lot of useful techniques from quantum mechanics and statistical physics. However, we would like to mention as well some restrictions of this approach. Intrinsically it is always suggested that this approach is suitable for systems containing many organisms, so that the relative fluctuations of density become negligible and all function are continuous. However this statement does not hold if we consider the survival-extinction transition. As the system approaches its critical state, the population density declines and the fluctuations of density (demographic stochasticity) start to play a crucial role [105]. Thus, in reality, the extinction of a population might occur under conditions which still allow for the population survival in a deterministic PDE framework. Therefore, the development of a theory including stochastic effects is necessary for the correct representation of the transient behaviour.

Acknowledgements
We are grateful to Aike Beckmann, Jef Huisman, and two anonymous referees for their helpful comments on the manuscript.

References
[1] M. R. Abbott, K. L. Denman, T. M. Powell, P. J. Richerson, R. C. Richards, and C. R. Goldman. Mixing and the dynamics of the deep chlorophyll maximum in lake tahoe. Limnol. Oceanogr., 29 (1984), 862­878. [2] R. Barakat. A note on the transient stage of the random dispersal of logistic populations. B. Math. Biol., 21 (1959), 141­151. [3] A. Beckmann and I. Hense. Beneath the surface: Characteristics of oceanic ecosystems under weak mixing conditions - a theoretical investigation. Prog. Oceanogr., 75 (2007), 771­796.

78


A. B. Ryabov and B. Blasius

The role of diffusion and advection

[4] D. A. Birch, Y.-K. Tsang, and W. R. Young. Bounding biomass in the fisher equation. Phys. Rev. E, 75 (2007), 066304­14. [5] J. K. B. Bishop. Transmissometer measurement of POC. Deep-Sea Res. Pt. I, 46 (1999), 353­369. [6] L. Bopp, P. Monfray, O. Aumont, J. L. Dufresne, H. L. Treut, G. Madec, L. Terray, and J. C. Orr. Potential impact of climate change on marine export production. Global Biogeochem. Cycles, 15 (2001), 81­100. [7] R. S. Cantrell and C. Cosner. The effects of spatial heterogeneity in population dynamics. J. Math. Biol., 29 (1991), 315­338. [8] R. S. Cantrell and C. Cosner. On the effects of spatial heterogeneity on the persistence of interacting species. J. Math. Biol., 37 (1998), 103­145. [9] R. S. Cantrell and C. Cosner. Spatial heterogeneity and critical patch size: Area effects via diffusion in closed environments. J. Theor. Biol., 209 (2001), 161­171. [10] S. Clodong and B. Blasius. Chaos in a periodically forced chemostat with algal mortality. P. Roy. Soc. Lond. B Bio., 271 (2004), 1617­1624. [11] F. Colijn and G. C. Cadee. Is phytoplankton growth in the wadden sea light or nitrogen ´ limited? J. Sea Res., 49 (2003), 83­93. [12] S. A. Condie and M. Bormans. The influence of density stratification on particle settling, dispersion and population growth. J. Theor. Biol., 187 (1997), 65­75. [13] J. J. Cullen. The deep chlorophyll maximum: comparing vertical profiles of chlorophyll a. J. Fish. Aquat. Sci, 39 (1982), l­803. [14] K. A. Dahmen, D. Nelson, and N. M. Shnerb. Life and death near a windy oasis. J. Math. Biol., 41 (2000), 1­23. [15] J. M. Diamond and R. M. May. Island biogeography and the design of natural reserves. In R. M. May, editor, Theoretical Ecology: Principles and Applications, pages 163­186. Saunders, Philadelphia, USA, 1976. [16] W. Ebenhoh, C. Kohlmeier, and P. J. Radford. The benthic biological submodel in the ¨ european regional seas ecosystem model. Netherlands J. Sea Res., 33 (1995), 423­452. [17] U. Ebert, M. Arrayas, N. Temme, B. Sommeijer, and J. Huisman. Critical conditions for ´ phytoplankton blooms. B. Math. Biol., 63 (2001), 1095­1124. [18] H. C. Ellertsen. Spring blooms and stratification. Nature, 363 (1993), 24. [19] W. E. Fagan, R. S. Cantrell, and C. Cosner. How habitat edges change species interactions. Am. Nat., 153 (1999), 165­182. 79


A. B. Ryabov and B. Blasius

The role of diffusion and advection

[20] R. A. Fisher. The wave of advance of advantageous genes. Ann. Eugenics, 7 (1937), 355. [21] M. J. Follows, S. Dutkiewicz, S. Grant, and S. W. Chisholm. Emergent biogeography of microbial communities in a model ocean. Science, 315 (2007), 1843­1846. [22] R. J. Geider, H. L. MacIntyre, and T. M. Kana. Dynamic model of phytoplankton growth and acclimation: responses of the balanced growth rate and the chlorophyll a: carbon ratio to light, nutrient-limitation and temperature. Mar. Ecol-Prog. Ser, 148 (1997), 187­200. [23] H. H. Gran and T. Braarud. Quantitative study of the phytoplankton in the bay of fundy and the gulf of maine (including observations on hydrography, chemistry and turbidity). J. Fish. Res. Board. Can., 1 (1935), 279­467. [24] P. Grassberger and I. Procaccia. Diffusion and drift in a medium with randomly distributed traps. Phys. Rev. A, 26 (1982), 3686­3688. [25] J. P. Grover. Resource competition in a variable environment: Phytoplankton growing according to monod's model. Am. Nat., 136 (1990), 771­789. [26] J. P. Grover. Resource Competition. Springer, 1997. [27] W. Gurney and R. Nisbet. The regulation of inhomogeneous populations. J. Theor. Biol., 52 (1975), 441­457. [28] A. E. Hershey, J. Pastor, B. J. Peterson, and G. W. Kling. Stable isotopes resolve the drift paradox for baetis mayflies in an arctic river. Ecology, 74 (1993), 2315­2325. [29] P. Hess. Periodic-Parabolic Boundary Value Problems and Positivity. Longman Sc & Tech, 1991. [30] B. A. Hodges and D. L. Rudnick. Simple models of steady deep maxima in chlorophyll and biomass. Deep Sea Res. Part I, 51 (2004), 999­1015. [31] O. Holm-Hansen and C. D. Hewes. Deep chlorophyll-a maxima (DCMs) in Antarctic waters - I. Relationships between DCMs and the physical, chemical, and optical conditions in the upper water column. Polar Biol., 27 (2004), 699­710. [32] E. E. Holmes, M. A. Lewis, J. E. Banks, and R. R. Veit. Partial differential equations in ecology: Spatial interactions and population dynamics. Ecology, 75 (1994), 17­29. [33] J. Huisman, M. Arrayas, U. Ebert, and B. Sommeijer. How do sinking phytoplankton species ´ manage to persist? Am. Nat., 159 (2002), 245­254. [34] J. Huisman, J. Sharples, J. M. Stroom, P. M. Visser, W. E. A. Kardinaal, J. M. H. Verspagen, and B. Sommeijer. Changes in turbulent mixing shift competition for light between phytoplankton species. Ecology, 85 (2004), 2960­2970.

80


A. B. Ryabov and B. Blasius

The role of diffusion and advection

[35] J. Huisman and B. Sommeijer. Population dynamics of sinking phytoplankton in lightlimited environments: simulation techniques and critical parameters. J. Sea Res., 48 (2002), 83­96. [36] J. Huisman, N. N. P. Thi, D. M. Karl, and B. Sommeijer. Reduced mixing generates oscillations and chaos in the oceanic deep chlorophyll maximum. Nature, 439 (2006), 322­325. [37] J. Huisman, P. van Oostveen, and F. J. Weissing. Critical depth and critical turbulence: Two different mechanisms for the development of phytoplankton blooms. Limnol. Oceanogr., 44 (1999), 1781­1787. [38] A. Huppert, B. Blasius, and L. Stone. A model of phytoplankton blooms. Am. Nat., 159 (2002), 156­171. [39] A. Huppert, B. Blasius, and L. Stone. What minimal models can tell: A reply to van Nes and Scheffer. Am. Nat., 163 (2004), 927­929. [40] H. Ishii and I. Takagi. Global stability of stationary solutions to a nonlinear diffusion equation in phytoplankton dynamics. J. Math. Biol., 16 (1982), 1­24. [41] B. M. Jamart, D. F. Winter, K. Banse, G. C. Anderson, and R. K. Lam. A theoretical study of phytoplankton growth and nutrient distribution in the pacific ocean off the northwestern U.S. coast. Deep Sea Res., 24 (1977), 753­773. [42] K. D. Johnk, J. Huisman, J. Sharples, B. Sommeijer, P. M. Visser, and J. M. Stroom. Summer ¨ heatwaves promote blooms of harmful cyanobacteria. Glob. Change Biol., 14 (2008), 495­ 512. [43] J. Joo and J. L. Lebowitz. Population dynamics in spatially heterogeneous systems with drift: The generalized contact process. Phys. Rev. E, 72 (2005), 036112­7. [44] H. Kierstead and L. B. Slobodkin. The size of water masses containing plankton bloom. Mar. Res., 12 (1953), 141­147. [45] J. T. O. Kirk. Light and Photosynthesis in Aquatic Ecosystems. Cambridge University Press, 1994. [46] C. A. Klausmeier and E. Litchman. Algal games: The vertical distribution of phytoplankton in poorly mixed water columns. Limnol. Oceanogr., 46 (2001), 1998­2007. [47] C. A. Klausmeier and D. Tilman. Spatial models of competition. In U. Sommer and B. Worm, editors, Competition and coexistence, volume 161, pages 43­78. Springer, 2002. [48] A. N. Kolmogorov, I. G. Petrovskii, and N. S. Piskunov. Study of the diffusion equation with growth of the quantity of matter and its application to a biology problem. Moscow Univ. Math. Bull. (Engl. Transl.), 1 (1937), 1­25. 81


A. B. Ryabov and B. Blasius

The role of diffusion and advection

[49] K. Kononen, M. Huttunen, S. Hallfors, P. Gentien, M. Lunven, T. Huttula, J. Laanemets, M. Lilover, J. Pavelson, and A. Stips. Development of a deep chlorophyll maximum of Heterocapsa triquetra Ehrenb. at the entrance to the Gulf of Finland. Limnol. Oceanogr., 48 (2003), 594­607. [50] I. Koszalka, A. Bracco, C. Pasquero, and A. Provenzale. Plankton cycles disguised by turbulent advection. Theor. Popul. Biol., 72 (2007), 1­6. [51] M. Kot. Elements of Mathematical Ecology. Cambridge University Press, 2001. [52] D. H. Landahl. A note on population growth under random dispersal. B. Math. Biol., 21 (1959), 153­159. [53] L. D. Landau and E. M. Lifshitz. Quantum Mechanics (Non-Relativistic Theory), volume 3. Butterworth-Heinemann, 1981. [54] M. Levandowsky and B. S. White. Randomness, time scales, and the evolution of biological communities. Evolutionary Biology, 10 (1977), 69­169. [55] A. L. Lin, B. A. Mann, G. Torres-Oviedo, B. Lincoln, J. Kas, and H. L. Swinney. Localization and extinction of bacterial populations under inhomogeneous growth conditions. Biophys. J., 87 (2004), 75­80. [56] D. Ludwig, D. G. Aronson, and H. F. Weinberger. Spatial patterning of the spruce budworm. J. Math. Biol., 8 (1979), 217­258. [57] A. P. Martin. Phytoplankton patchiness: the role of lateral stirring and mixing. Prog. Oceanogr., 57 (2003), 125­174. [58] R. McMurtrie. Persistence and stability of single-species and prey-predator systems in spatially heterogeneous environments. Math. Biosci, 39 (1978), 11­51. [59] B. G. Mitchell, E. A. Brody, O. Holm-Hansen, C. McClain, and J. Bishop. Light limitation of phytoplankton biomass and macronutrient utilization in the southern ocean. Limnol. Oceanogr., 36 (1991), 1662­1677. [60] A. Moll and G. Radach. Review of three-dimensional ecol. model. related to the north sea shelf system part 1: models and their results. Prog. Oceanogr., 57 (2003), 175­217. [61] E. W. Montroll. Lectures on nonlinear rate equations, especially those with quadratic nonlinearities. In A. O. Barut and W. E. Brittin, editors, Lectures in Theoretical Physics, volume X­A, pages 531­573. Interscience Publishers, 1968. [62] J. D. Murray. Mathematical Biology. Springer, 2003. [63] J. D. Murray and R. P. Sperb. Minimum domains for spatial patterns in a class of reaction diffusion equations. J. Math. Biol., 18 (1983), 169­184. 82


A. B. Ryabov and B. Blasius

The role of diffusion and advection

[64] D. R. Nelson and N. M. Shnerb. Non-hermitian localization and population biology. Phys. Rev. E, 58 (1998), 1383­1403. [65] C. Neuhauser. Mathematical challenges in spatial ecology. Not. Am. Math. Soc., 48 (2001), 1304­1314. [66] J. Nihoul, editor. Modelling of Marine Systems Elsevier Oceanography Series, volume 10. Elsevier, Amsterdam, 1975. [67] A. Okubo. Oceanic diffusion diagrams. Deep Sea Res., 18 (1971), 789­802. [68] A. Okubo. A note on small organism diffusion around an attractive center: A mathematical model. J. Oceanogr., 28 (1972), 1­7. [69] A. Okubo. Diffusion-induced instability in model ecosystems. Chesapeake Bay Institute Tech. Rept, 86 (1974), 17. [70] A. Okubo. Advection-diffusion in the presence of surface convergence. Oceanic fronts in coastal processes, pages 23­28, 1978. [71] A. Okubo. Horizontal dispersion and critical scales for phytoplankton. In J. Steele, editor, Spatial pattern in plankton communities, pages 21­42. Springer, 1978. [72] A. Okubo and S. A. Levin. Diffusion and Ecological Problems. Springer, 2001. [73] R. V. Ozmidov. Horizontal Turbulence and Turbulent Exchange in an Ocean. Nauka, Moscow, 1968. [74] R. V. Ozmidov. Phytoplankton patches in the ocean under various modes of ocean turbulence. Oceanology, 38 (1998), 3­11. [75] E. Pachepsky, F. Lutscher, R. Nisbet, and M. Lewis. Persistence, spread and the drift paradox. Theor. Popul. Biol., 67 (2005), 61­73. [76] F. Peeters, D. Straile, A. Lorke, and D. Ollinger. Turbulent mixing and phytoplankton spring bloom development in a deep lake. Limnol. Oceanogr., 52 (2007), 286­298. [77] S. V. Petrovskii. On the diffusion of a plankton patch in a turbulent ocean. Oceanology, 39 (1999), 737­742. [78] S. V. Petrovskii. On the plankton front waves accelerated by marine turbulence. J. Marine. Syst., 21 (1999), 179­188. [79] S. V. Petrovskii and H. Malchow. A minimal model of pattern formation in a prey-predator system. Math. Comput. Model., 29 (1999), 49­63. [80] S. V. Petrovskii and H. Malchow. Wave of chaos: New mechanism of pattern formation in spatio-temporal population dynamics. Theor. Popul. Biol., 59 (2001), 157­174. 83


A. B. Ryabov and B. Blasius

The role of diffusion and advection

[81] T. Platt and K. L. Denman. A general equation for the mesoscale distribution of phytoplankton in the sea. In Mem. Soc. Roy. Sci. de Liege, 6e serie, pages 31­42, 1975. [82] T. Platt, K. L. Denman, and A. D. Jassby. Modeling the productivity of phytoplankton. The Sea. Marine modeling, 6 (1977), 857­890. [83] T. Platt, C. L. Gallegos, and W. G. Harrison. Photoinhibition of photosynthesis in natural assemblages of marine phytoplankton. J. Mar. Res, 38 (1980), 687­701. [84] H. Ploug, M. H. Iversen, and G. Fischer. Ballast, sinking velocity, and apparent diffusivity within marine snow and zooplankton fecal pellets: Implications for substrate turnover by attached bacteria. Limnol. Oceanogr., 53 (2008), 1878­1886. [85] G. Radach and E. Maier-Reimer. The vertical structure of phytoplankton growth dynamics; a mathematical model. In Mem. Soc. Roy. Sci. de Liege, 6e serie, pages 113­146, 1975. [86] C. S. Reynolds and A. E. Irish. Modelling phytoplankton dynamics in lakes and reservoirs: the problem of in-situ growth rates. Hydrobiologia, 349 (1997), 5­17. [87] G. A. Riley, H. Stommel, and D. F. Bumpus. Quantitative ecology of the plankton of the western north atlantic. B. Bingham Oceanogr. C., 12 (1949), 1­169. [88] S. Roy and J. Chattopadhyay. Towards a resolution of `the paradox of the plankton': A brief overview of the proposed mechanisms. Ecol. Complex., 4 (2007), 26­33. [89] A. B. Ryabov, L. Rudolf, and B. Blasius. The influence of an upper mixed layer on the vertical distribution and composition of phytoplankton biomass. In preparation. [90] E. Sakshaug, D. Slagstad, and O. Holm-Hansen. Factors controlling the development of phytoplankton blooms in the antarctic ocean: a mathematical model. Mar. Chem., 35 (1991), 259­271. [91] J. L. Sarmiento, R. Slater, R. Barber, L. Bopp, S. C. Doney, A. C. Hirst, J. Kleypas, R. Matear, U. Mikolajewicz, and P. Monfray. Response of ocean ecosystems to climate warming. Global Biogeochem. Cycles, 18, 2004. [92] S. Senn and P. Hess. On positive solutions of a linear elliptic eigenvalue problem with neumann boundary conditions. Mathematische Annalen, 258 (1982), 459­470. [93] H. Seno. Effect of a singular patch on population persistence in a multi-patch system. Ecol. model., 43 (1988), 271­286. [94] A. L. Shanks and J. D. Trent. Marine snow: sinking rates and potential role in vertical flux. Deep-Sea Res., 27 (1980), 137­143. [95] N. Shigesada and A. Okubo. Analysis of the self-shading effect on algal vertical distribution in natural waters. J. Math. Biol., 12 (1981), 311­326. 84


A. B. Ryabov and B. Blasius

The role of diffusion and advection

[96] J. G. Skellam. Random dispersal in theoretical populations. Biometrika, 38 (1951), 196­ 218. [97] L. Slobodkin. Akira okubo and the theory of blooms. Oceanography, 12 (1999), 9­14. [98] U. Sommer and B. Worm, editors. Competition and Coexistence. Springer, 2002. [99] D. C. Speirs and W. S. C. Gurney. Population persistence in rivers and estuaries. Ecology, 82 (2001), 1219­1237. [100] P. A. Stephens and W. J. Sutherland. Consequences of the allee effect for behaviour, ecology and conservation. Trends Ecol. Evol., 14 (1999), 401­405. [101] A. V. Straube and A. Pikovsky. Mixing-induced global modes in open active flow. Phys. Rev. Lett., 99 (2007), 184503­4. [102] H. U. Sverdrup. On conditions for the vernal blooming of phytoplankton. ICES J. Mar. Sci., 18 (1953), 287­295. [103] N. N. P. Thi, J. Huisman, and B. P. Sommeijer. Simulation of three-dimensional phytoplankton dynamics: competition in light-limited environments,. J. Comput. Appl. Math., 174 (2005), 57­77. [104] D. Tilman. Resource Competition and Community Structure. Princeton University Press, 1982. [105] D. Tilman and P. M. Kareiva, editors. Spatial Ecology: The Role of Space in Population Dynamics and Interspecific Interactions. Princeton University Press, 1997. [106] D. Townsend, M. D. Keller, M. E. Sieracki, and S. G. Ackleson. Spring phytoplankton blooms in the absence of vertical water column stratification. Nature, 360 (1992), 59­62. [107] D. Turpin. Physiological mechanisms in phytoplankton resource competition. In C. Sandgren, editor, Growth and reproductive strategies of freshwater phytoplankton. pages 316­ 368. Cambridge University Press, 1991. [108] W. van Saarloos. Front propagation into unstable states. Phys. Rep., 386 (2003), 29­222. [109] R. A. Varela, A. Cruzado, and J. Tintore. A simulation analysis of various biological and physical factors influencing the deep-chlorophyll maximum structure in oligotrophic areas. J. Mar. Syst., 5 (1994), 143­157. [110] E. L. Venrick. Phytoplankton seasonality in the central north pacific - the endless summer reconsidered. Limnol. Oceanogr., 38 (1993), 1135­1149. [111] E. L. Venrick, J. McGowan, and A. Mantyla. Deep maxima of photosynthetic chlorophyll in pacific ocean. Fishery Bulletin, 71 (1973), 41­52. 85


A. B. Ryabov and B. Blasius

The role of diffusion and advection

[112] R. A. Walters. A time- and depth-dependent model for physical, chemical and biological cycles in temperate lakes. Ecol. model., 8 (1980), 79­96. [113] D. F. Winter, K. Banse, and G. C. Anderson. The dynamics of phytoplankton blooms in puget sound a fjord in the northwestern united states. Mar. Biol., 29 (1975), 139­176. [114] J. S. Wroblewski and J. O'Brien. A spatial model of phytoplankton patchiness. Mar. Biol., 35 (1976), 161­175. [115] J. S. Wroblewski. A model of phytoplankton plume formation during variable oregon upwelling. J. Mar. Res, 35 (1977), 357­394. [116] J. S. Wroblewski, J. O'Brien, and T. Platt. On the physical and biological scales of phytoplankton patchiness in the ocean. Mem. Soc. Roy. Sci. de Liege, 6e serie, pages 43­57, 1975. [117] E. Yakushev, F. Pollehne, G. Jost, I. Kuznetsov, B. Schneider, and L. Umlauf. Analysis of the water column oxic/anoxic interface in the black and baltic seas with a numerical model. Mar. Chem., 107 (2007), 388­410. [118] K. Yoshiyama and H. Nakajima. Catastrophic transition in vertical distributions of phytoplankton: Alternative equilibria in a water column. J. Theor. Biol., 216 (2002), 397­408.

86