A dynamical approach to phytoplankton blooms

Algae in the ocean absorb carbon dioxide from the atmosphere and thus play an important role in the carbon cycle. An algal bloom occurs when there is a rapid increase in an algae population. We consider a reaction-advection-diffusion model for algal bloom density and present new proofs of existence and uniqueness results for the steady state solutions using techniques from dynamical systems. On the question of stability of the bloom profiles, we show that the only possible bifurcation would be due to an oscillatory instability.

1. Introduction. As algae play a significant role in the carbon cycle, it is crucial to understand the occurence and structure of blooms [1]. We will examine steady state solutions to the nonlocal reaction-diffusion-advection equation for phytoplankton density P (z, t) given by P t = P zz − CP z + A (F (r(z, t, P )) − B) P, P (z, t) ≥ 0, P z (0, t) − CP (0, t) = P z (L, t) − CP (L, t) = 0, r(z, t, P ) = e −z− z 0 P (y,t) dy , where z is depth, L ∈ (0, ∞), (z, t) ∈ [0, L] × [0, ∞), A, B, C are parameters, and F is a suitably smooth and monotone increasing function of r. This model is derived in [6] from the nutrient-phytoplankton-zooplankton (NPZ) model proposed in [13]. The surface of the water corresponds to z = 0 and the bottom of the layer is L > 0. The trivial state P = 0 is always a solution to (1)-(3), while stationary solutions to (1)-(3) satisfy the ODE 0 = P − CP + A (F (r(z, P )) − B) P, P (z) ≥ 0, r(z, P ) = e −z− z 0 P (y) dy .
The term r(z, P ) is a light source term which includes integration over the domain [0, z]; consequently, (4)-(6) is a nonlocal nonlinear system. This particular form of (6) includes the background turbidity of the water (e −z ) to decrease light with increasing depth. An alternative model is the "self-shading case" r(z, P ) = exp(− z 0 P (y) dy) in which light limitations are solely imposed by the population density.
Different models for phytoplankton blooms have been studied in depth by several authors since Riley et al. [17]. We focus here on models that assume that nutrients are in ample supply and therefore species growth is determined by light only. Shigesada and Okubo [18] proposed the self-shading case of (1) on an infinite domain (L = ∞). They use phase space techniques, but the assumption of self-shading allows them to reduce the problem to a two-dimensional phase plane. Our approach is close to theirs, but the turbidity introduces a z-dependent term that requires a three-dimensional phase space. Ishii and Takagi [11] followed [18] with a study of the existence and uniqueness of a positive steady state for both the self-shading case and in the case of turbid water, again on an infinite domain with the boundary condition lim z→∞ P (z, t) = 0. Using comparison methods, they found a bifurcation: either the trivial solution is the only steady-state solution and is globally stable, or it is unstable and there is a unique positive solution which is stable. As [11] restrict to the case of an infinite domain, they miss the existence of a critical depth.
More recently, Ebert et al [6] analyzed conditions leading to survival of algal bloom density profiles on a finite domain [0, L]. Their dimensional analysis results in the four nondimensional parameters A, B, C, L used in this paper, but with the assumption F (r(z, P )) = r(z, P ). They determined the existence of a critical depth L * and confirmed analytically the existence of a vertical turbulent diffusion found numerically by Huisman et. al. [9]. They include many numerical simulations and address the alternate case F (r(z, P )) = ar(z, P ) γ . Huisman et al. [8] explored parameter regimes and solved the boundary value problem numerically to show that under certain light conditions, the phytoplankton does develop a stationary density profile, which may attain a maximum density under the surface z = 0.
Using phase plane techniques, Kolonikov et al. [14] determined that the selfshading model has a unique positive steady state which is stable. Moreover, there is nontrivial solution at every depth; hence without turbidity there is no critical water depth L * . Du and Hsu [5] studied the stability of the advectionless model (C = 0) and established a critical death rate d * so that the density profile tends to a unique positive steady state if d < d * and tends to zero as t → ∞ if d ≥ d * . Hsu and Lou [10] studied (1)-(3) on a finite domain [0, L] with slightly different expressions for the parameters (the death rate, the sinking or buoyant coefficient, water column depth, and vertical turbulent diffusion rate). Using bifurcation and PDE techniques, they established a critical death rate d * : if the death rate is smaller than d * , then there is a unique positive steady state. Moreover, they found that the critical death rate is monotone decreasing with respect to the sinking/buoyant velocity coefficient. As in our approach, they use a general light source term which is monotonically increasing but which must satisfy an extra condition that we relax: F (r(z, P )) ≥ ar(z, P ) γ for I ∈ [0, I 0 ]. They found that there may or may not be a critical water depth and expressed this depth in terms of d * . Du and Mei answered the question of a critical depth in [4] and were able to relate existence of a critical depth to the advection coefficient C. Moreover, they established existence and uniqueness of the positive steady state for d ∈ (0, d * ).
Our model and parameters follow [6] with results similar to those in [4] and [10]. However, we approach the model with techniques from dynamical systems that work in the three-dimensional phase space setting. We set up the phase space for (4)- (6) to show that the existence and uniqueness of a positive steady state can be seen through the shape of a curve of solutions in phase space. This paper is organized as follows. In Section 2 we give the derivation of the model from the reduced Nutrient-Zooplankton-Phytoplankton (NPZ) model, following closely the approach in [6]. In 2.1, we make a generalization on the light source term and set up a system of first-order equations to illustrate the dynamics of the phytoplankton density equation. We describe how the geometry of the phase space determines both the existence of a steady-state phytoplankton bloom and the uniqueness of the surviving bloom's density profile. In Section 3, we establish conditions leading to the existence of a critical depth and prove that for subcritical depths, there must be a nontrivial stationary solution, while for supercritical depths the only solution is the trivial one. In Section 4, we prove that whenever a nontrivial stationary solution exists, it is unique. In Section 5, we address the nonlocal eigenvalue problem associated with (1)-(3) to prove that any instability in a nontrivial solution must result from a Hopf-type bifurcation.
2. Algal bloom model. The algal bloom model we consider is derived from a reduced NPZ model proposed by Klausmeier and Litchman [13]; see [12] for an overview by Kaper and Engler. The standard NPZ model includes a standard advection-diffusion PDE for each component N , P or Z; for phytoplankon, this equation is The first term is the overall change in the density of phytoplankton over time. The components u, v, w describe the velocity of the surrounding fluid in the horizontal (x and y) and vertical (z) directions, and w s is the vertical sinking speed. Hence the latter three terms on the left-hand describe horizontal and vertical advection. On the right-hand side, D h and D v are horizontal and vertical diffusion coefficients; the first two terms on the right-hand side describe turbulent diffusion. The final term is the change in phytoplankton due to the presence of light, nutrients, and zooplankton; this term is described by the NPZ model. This closed system of equations consists of five transfer functions relating the populations of the nutrients, phytoplankton, zooplankton, and the availability of light (f (I)): ∂N ∂t bio = −f (I)g(N )P + (1 − γ)h(P )Z + i(P )P + n(Z)Z.
When an algal bloom occurs, the population of phytoplankton has dramatically increased. To describe this phenomenon with the NPZ model, Klausmeier and Litchman [13] suppose the nutrients are in ample supply and the presence of zooplankton can be disregarded. Moreover, they reduce the dimensions of the model by assuming that the phytoplankton move passively and are uniformly distributed horizontally; thus P depends only on t and the depth z. Consequently, the terms in (7) depending on the horizontal components x and y vanish. Hence (7)-(10) reduce to the following advection-diffusion equation for P (z, t) Here the function g(N ) has been replaced by the constant a > 0 to model the ample supply of nutrients and i(P ) has been replaced with a constant loss rate . The light source is attenuated according to a function f (I) satisfying the Lambert-Beer law for light: where K bg is the turbidity and k is a coefficient determined by the particular type of phytoplankton. Light is maximal at z = 0 (the surface of the water) and monotone decreasing with increasing depth. Phytoplankton does not move across the boundaries at z = 0 or z = H > 0 (the unscaled depth of the water). Integrating over a thin layer results in the boundary conditions Ebert et. al. [6] rescaled and nondimensionalized (11)-(14) using with the expressions Under (15), the Lambert-Beer law for the light penetration through the layer is the nonlinear term Dropping the notation, the resulting advection-diffusion equation for P (z, t), with zero-flux boundary conditions depends on four dimensionless parameters A, B, C, and L. The parameter A is the ratio of the growth rate at the surface and the scales of absorption and diffusion, B is the ratio between the death and growth rates at the surface, C is a measure of the tendency of the algae to sink (C < 0 measures buoyancy, C > 0 measures sinking), and L is the dimensionless depth of the layer.
As A and B are both ratios, they must be positive. If the death rate is greater than the maximal growth rate at the surface, i.e. B > 1, there is no nontrivial solution [6]. (We rederive this result through the dynamics in Section 3, see the remarks following Lemma 3.4.) Thus we adopt the condition 0 < B < 1. Overall, the value range of each parameter is Lastly, for any choice of z > 0, the bounds of integration in (16) rely on knowledge of P (z, t) over the interval (0, z); therefore, (17) is a nonlocal, nonlinear equation.

Stationary solutions.
In this section, we set up (4)-(6) as a dynamical system to explore the geometry of steady state solutions to (1)-(3). We will make use of this system to give geometric proofs of existence and uniqueness of solutions to (4)- (6) in Sections 3 and 4. In addition, we include a generalization for the light term. Setting where it is understood from the integrand that r depends on the density profile P (z), we replace the term j(P ) in (17) with any F (r) satisfying Using (16) as in (17) is simply choosing F (r(z)) = r(z) or, after rescaling z and P , F (r(z)) = r(z) s , s > 0. Here F represents the phytoplankton growth rate as a function of the light source. A typical example where the growth saturates near the surface boundary is the Michaelis-Menten-type relation F (r) = ar 1 + cr (20) for appropriate choices of the growth rate a > 0 and c > 0. Another physical model is see [16] and [19]. With this light source term, (17) becomes with the natural requirement that P (z) ≥ 0 for all 0 ≤ z ≤ L, and with the mixed boundary conditions It is frequently convenient to express the boundaries of the layer in terms of r rather than z. Notice r (z) = (−1 − P (z)) r(z).
As r(0) = 1, and as P (z) ≥ 0 for solutions to (4)-(6), r(z) is a decreasing function of z that approaches 0 asymptotically as z → ∞. We can now view the boundaries 0 ≤ z ≤ L as moving from r = 1 to r = r(L) ∈ (0, 1). (As r depends on the phytoplankton density P (z) as well as z in (19), r(L) is different for different density profiles. If there is a nonphysical solution with P (z) < −1, then r(z) may increase; this mathematical possibility does not affect the results.) We introduce q = P to attain the following first-order system:  For any solution to the steady state problem (25)-(27) satisfying (23), we use a to refer to its initial condition at z = 0 so that (P (0; a), q(0; a), r(0; a)) = (a, Ca, 1). Using this notation, we define a curve formed by solutions to (25)-(27) at a chosen depth z ∈ [0, L] by Plotted in (P, q, r)-space, this curve tracks where each initial condition along the line q = CP is at depth z; essentially, it moves the line (a, Ca, 1) under the flow of (25)-(27) to depth z. For clarity, we also make use of a truncated version of γ(z) defined by Showing existence amounts to showing that the curve γ(L) intersects the plane {q = CP } away from the r-axis with P > 0. The point of intersection corresponds to a solution of (22) which satisfies the boundary conditions at the surface of the water (by virtue of being on γ(L)) and at the bottom of the layer (as it intersects the plane {q = CP }), see Figure 1. Uniqueness follows if there is no more than one such intersection.

Conditions for existence.
In this section, we use the phase space for (25)-(27) to establish that for a given light source term F (r), there is a critical value B * = F (1) for the parameter B ∈ (0, 1) that determines whether a nontrivial solution to (4)-(6) can exist. If B ≥ B * , then there can be no nontrivial solution.
If B < B * , then there is either always a nontrivial steady state solution to (1)-(3), or there exists a critical depth L * > 0 so that a nontrivial steady state only exists at subcritical depths. The existence of such a critical depth is shown in [6] for F (r) = r, we demonstrate its existence through an analysis of the dynamics of (25)-(27) for any F satisfying (H). To show that a nontrivial steady state solution exists at a certain depth L, it suffices to establish the validity of two criteria: 1. the solution curve γ(L) defined in (28) satisfies q(L, ) < CP (L, ) for small initial conditions > 0. In other words, the truncated solution curve γ(L, ) is under the plane {q = CP } (see Lemma 3.2). 2. As the initial condition a increases without bound, the curve γ(L) moves above the plane with P > 0. Hence there is some nontrivial initial condition α > 0 so that γ(L) intersects the plane at (P (L; α), q(L; α), r(L; α)). This is verified in Lemma 3.3. If these two statements hold, then the density profile P (z; α) is a solution to (4)- (6) for that particular depth L.
If B < F (1), we will show that the second statement is always true in Lemma 3.3. In the critical depth case, the first statement will only be true if L < L * , where L * may be infinite. In phase space, the existence of a critical depth means that for all L < L * , the curve γ(L) intersects the plane {q = CP } at some point (P (L; a 0 ), q(L; a 0 ), r(L; a 0 )) away from the r-axis (i.e. a 0 > 0). This point of intersection corresponds to a density profile P (z; a 0 ) which is a nontrivial steady state solution to (22)-(23). For all L ≥ L * , there is no such intersection, and thus no nontrivial solution to (22)-(23).
(If there is no critical depth L * , then γ(L) intersects the plane away from the r-axis at all depths L ∈ (0, ∞).) Remark 1. Existence of a nontrivial solution could theoretically occur if γ(L) is initially above the plane for small initial conditions and subsequently crosses below at some initial condition α > 0. This alternative configuration is proven impossible by Corollary 3 in Section 4 (the intersection can only occur in one direction).
To show that the solution curve initially moves under the plane, we use variational techniques to track tangent vectors to γ(z) as z increases. For any choice of z ∈ [0, L], the tangent vector field along γ(z) is computed by  27) is This set of equations describes the motion of the tangent vector field defined in (30) under the flow of (25)-(27). Along the trivial solution P = q = 0, which begins at r = 1 and descends to r = e −L < 1, the variational equations (31)-(33) reduce to Proof. We define a function W a (z) that measures the cross product of the constant vector (1, C) parallel to the plane {q = CP }, and the vector (δP, δq) tangent to γ(z) along any solution with initial condition a: The function W 0 (z) measures whether the trivial solution's tangent vector at any given depth z = − ln r, 0 < r ≤ 1, is under the plane (W 0 (z) < 0) or over the plane (W 0 (z) > 0). By the initial condition at r = 1, the curve γ(0) is the line (P, CP ); hence As A > 0, the sign of W 0 (z) is determined (B − F (r))δP . Initially, as 0 < B < F (1) by assumption. Hence for small positive depths z, W 0 (z) is a decreasing function. As a result the curve γ(z) must be initially under the plane {q = CP } for small enough initial conditions. Lemma 3.2 is stated for L sufficiently small as W 0 (z) is only a decreasing function when B − F (r(z)) < 0. The first criterion for existence is then satisfied for all L satisfying W 0 (z) < 0. For all L with W 0 (z) ≥ 0, existence will not occur, as we discuss in Lemma 3.4.
In particular, by Lemma 3.2 and Corollary 3 in Section 4, if B ≥ F (1), then γ(L) will be entirely above the plane {q = CP } at all depths. Thus we define a critical value B * = F (1) and obtain the following result. Henceforth we assume B < B * . The second requirement for existence, namely that for any L > 0 the curve γ(L) passes above the plane as the initial condition increases, is proven in the following lemma. Proof. We introduce a new term which compares the ratio of q/P to the plane {q = CP }. Let η(z; a) = q(z; a)/P (z; a) − C/2; this expression is useful as showing q > CP is equivalent to showing η > C/2. Initially η(0; a) = C/2, and where is differentiation with respect to z. Let y(z) = −C 2 /4 + A(F (r(z)) − B). By (H), y(z) < y(0) for all z > 0. Thus There are three cases to consider for the right-hand side of (40). We show that for each case, there is some K ∈ R and an interval [0, z 0 ] so that 1. Suppose y(0) > 0 and C ∈ R. Dividing both sides of (40) by η(z; a) 2 + y(0) and integrating the result from z = 0 to some choice of z =ẑ yields η(ẑ;a) Consequently, for all z ≥ 0 small, η(z; a) ≥ y(0) tan arctan C 2 y(0) − z y(0) .
By continuity, we choose z 0 > 0 small enough so that for some prescribed > 0. Consequently, in (39) the choice of a 0 and imply that η (z 0 ; a 0 ) is positive and bounded away from 0. As a result, there is az > z 0 such that η(z) > C/2. As we can make z 0 andz − z 0 as small as desired, we have shown that q(z; a 0 )/P (z; a 0 ) > C for small z if the initial condition a 0 is large enough. This proves the lemma.
Proof. As η(z; a 0 ) > C/2 with A(B − F (z)) > 0, then if η(z; a 0 ) = C/2 for any z >z, it follows that Thus the solution will not return below the plane at larger depths z. It is mathematically possible for (P (z; a 0 ), q(z; a 0 )) to solve (4)-(6) with P < 0; however, (43) forces (P (z; a 0 ), q(z; a 0 )) to cross from the first to the second quadrant in the (P, q)-plane, which contradicts P = q. Thus (P (z; a 0 ), q(z; a 0 )) is pinned above {q = CP } in the right-half of the (P, q)-plane.
Lemmas 3.2 and 3.3 demonstrate that at small depths, the solution curve emanating from the trivial state initially passes below the plane and then breaks above the plane for large initial conditions. Consequently, as long as W 0 (z) < 0, γ(z) behaves this way and a positive solution to (4)-(6) will exist.
There are two terms in (38) determining whether W 0 (z) is increasing or decreasing: B − F (r) and δP . By (H) and the assumption that B < F (1), at the surface B − F (1) < 0, B − F (r) increases as r decreases, B − F (0) = B > 0. Hence the term B − F (r) switches from negative to positive at a specific depth z = z F > 0.
If δP (z; 0) > 0 for all z, then W 0 (z) switches from a decreasing function to an increasing one at z F . The fact that W 0 increases is not enough to conclude that γ(z) is entirely over the plane for large enough z. To prove that this can happen, resulting in the existence of a critical depth L * , we establish Lemma 3.4 below.
Lemma 3.4. If there exists some value z F > 0 so that W 0 (z) > 0 for all z > z F , then either 1. the tangent vector along the trivial solution satisfies lim z→∞ δP = 0. In this case, any finite depth L has a nontrivial steady state solution. 2. there is a finite critical depth L * > 0 so that W 0 (L * ) = 0. Case 1 is equivalent to saying L * = ∞, and yields a nontrivial steady state solution at any depth. For Case 2 there is a critical depth L * such that W 0 (L) ≥ 0 for any L > L * , which implies the curve γ(L) does not go under the plane {q = CP } past a certain depth. The proof of Lemma 3.4 uses the dynamics of (δP, δq, δr) near the plane {r = 0}.
Proof. As z → ∞, the trivial steady state solution approaches the origin in (P, q, r)space along the r-axis, see Figure 1. The plane {r = 0} is invariant for the 6dimensional system formed by (P, q, r, δP, δq, δr) with a single fixed point at 0. The plane {r = 0} is invariant; moreover, by (30) δr must also be zero in {r = 0}.
For the final case, if lim z→∞ δP > 0, then in a neighborhood of the origin, the vector (δp, δq) tends to the straight-line solution parametrized by v + in the first quadrant. As λ + > C, this line lies above the line δq = CδP . In this case, the tangent vector to γ(z) at the trivial solution is above the plane in the limit z → ∞. This possibility is Case 2 of the lemma. Remark 2. Lemma 3.3 proves that the point where γ(z) intersects the plane q = CP has P coordinate tending to large values as L → 0. Hence arbitrarily dense phytoplankton profiles are sustained over correspondingly thin depths. (Conversely, thin phytoplankton profiles extend to deeper depths.) Remark 3. The parameter B is the ratio between the between the death rate and growth rate at the surface of the water. As indicated in [6], with F (r) = r there is no nontrivial solution if B > 1 (when the death rate is larger than the surface growth rate). The fact that there is no solution can now be viewed as a consequence of (38) together with Corollary 3 of Section 4: if B * = 1 and B ≥ 1, then W 0 (z) > 0 for all z. Hence the solution curve γ(z) is always above the plane {q = CP }.  To prove Theorem 4.1, we consider the same vector field (δP, δq, δr) defined by (30). As in Section 3, this vector field satisfies the equations of variations (31)-(33) associated with (25)-(27): A nontrivial solution on the domain [0, L] must be unique if γ(L) can only intersect the plane {q = CP } in one direction, see Figure 1. The direction that γ(L) intersects the plane is related to the behavior of an associated normal vector field. For this purpose, let ζ a (z) be the third component of the cross-product (P (z; a), q(z; a), r(z; a)) × (δP (z; a), δq(z; a), δr(z; a)).
The derivative of ζ a (z; a) along a solution with initial condition a is ζ a (z) = Cζ a (z) − AP 2 (z; a) dF dr δr(z; a).
As this is a linear differential equation with ζ a (0) = 0, we may write ζ a (z) analytically as ζ a (z) = −Ae Cz z 0 e −Cs P 2 (s; a) dF dr δr(s; a) ds.
The sign of δr depends on the behavior of δP in a manner similar to that of ζ a and δr. We first prove that for a given initial condition a > 0, δP (z; a) can only vanish after P (z; a) vanishes.
Lemma 4.2. For any fixed initial condition a > 0, if there is some z 0 > 0 so that z 0 is the shallowest depth at which δP (z; a) vanishes, then P (z; a) = 0 for some z ∈ (0, z 0 ).
However, we may explicitly compute ζ a (z 0 ) = P δq − q δP | z=z0 = P (z 0 ; a) δq(z 0 ; a) ≤ 0, which contradicts (51). Thus for a given positive initial condition, δP cannot vanish before P . This result on ζ a allows us to understand how the solution curve γ(L) intersects the plane {q = CP }. In particular, projecting into (p, q)-plane, when γ(L) intersects the plane at a solution to (4)-(6) the tangent vector (δP (z; a), δq(z; a)) must always point "up" to the region q > CP . See Figure 4 for two such intersections.
To conclude that any positive solution to (4)-(6) is unique, we must show that γ(L) cannot return below the plane. By Corollary 3, if ζ a (L) < 0 for some a > 0, then it is not a solution to (4)- (6), which means that P (z; a) must change signs before z = L. This situation is hypothetically possible as the restriction P (z; a) > 0 is based entirely on physical constraints. The following lemma eliminates such a configuration; throughout its proof we advise the reader to consult Figure 5. Proof. Suppose A 1 < A 2 are two positive initial conditions corresponding to solutions to (4)-(6) for some depth L > 0. By Corollary 3, both must intersect {q = CP } transversally in the same direction. Hence there is an α 0 ∈ (A 1 , A 2 ) so that δq(L; α 0 ) = C δP (L; α 0 ) with ζ α0 (L) < 0.
However, this contradicts Lemma 4.2 for P (z; α 4 ). Thus there can be no such α 0 . So if there is a solution with initial condition A 1 , it is necessarily unique.

5.
Stability. In this section we investigate the stability of steady state solutions to the two-point boundary value problem (17): for some chosen L > 0. The light source term is the standard r(z, t) = e −z− z 0 P (y,t) dy , and F satisfies (H). By Section 3, if B < B * , then there is either a unique positive stationary solution to (53)-(54) at every depth L > 0, or there is a critical L * > 0 so there is a unique positive stationary solution if and only if L < L * . We couple any steady state solution P 0 (z) to r 0 (z) := e −z− z 0 P0(y) dy . Consequently, (P 0 (z), r 0 (z)) satisfy P 0 (0) = CP 0 (0), P 0 (L) = CP 0 (L).
(57) We determine the stability of P 0 (z) in the usual way by applying (53)-(54) to the function Q(z, t) = P 0 (z) + e λt φ(z). Linearizing about (P 0 , r 0 ) results in the eigenvalue problem or equivalently with the boundary conditions [φ − Cφ] z=0,L = 0. Due to the integration over the domain [0, z], if P 0 is nontrivial then (58) is a nonlocal eigenvalue problem (NLEP). Along the trivial solution P 0 = 0, (58) reduces to a standard eigenvalue problem, and we characterize the stability of the trivial solution using the linear operator The spectrum of L 0 is determined by the Sturm-Liouville equation Therefore, the eigenvalues of L 0 form a real, discrete, descending sequence ν 0 > ν 1 > · · · , and the oscillations of each eigenfunction are well-understood, see [2].
Stability of the trivial solution over the domain [0, L] is assured if ν 0 > 0. We arrive at this conclusion by comparing the oscillation of the first eigenfunction φ 0 to the tangent vector component δP (z, 0) defined in (30). In Section 3, for the case δP (z, 0) > 0 for all z > 0, we defined L * geometrically as the depth at which the tangent vector (δP, δq) along the trivial solution is tangent to the plane {q = CP }. The following theorem provides an alternative way to define L * : as z passes through L * , a bifurcation in the stability of the trivial solution occurs.
Theorem 5.1. If there is no critical depth or if L < L * , the trivial solution is unstable. Otherwise, if L > L * , the trivial solution is stable.
Proof. Let P 0 denote the trivial solution and let φ be the eigenfunction with largest eigenvalue ν; by Sturm-Liouville theory, φ has no interior zeros on (0, L). Notice (31)-(33) implies that δP 0 satisfies (As it is not necessarily the case that δP 0 (L) = CδP 0 (L), the above does not allow us to conclude that δP 0 (z) is an eigenfunction of L 0 with eigenvalue 0.) We multiply (62) by φ and (61) by δP 0 , and take the difference of the two to obtain After simplifying the left-hand side, the above may written Integrating (64) from 0 to d, where d ≥ 0 is a chosen depth, and using the boundary conditions at z = 0 yields for the left-hand side. There are two possibilities: either δP (z; 0) vanishes at some smallest depth d, or δP (z; 0) > 0 for all z ∈ [0, ∞). In the first case, By hypothesis, δP 0 is positive on (0, d) while δP 0 (d) < 0 and φ(z) does not change signs. As a result, ν > 0 and P 0 is unstable.
Otherwise if δP (z, 0) > 0 for all z ∈ [0, ∞), then for any L we may say where W P0(0) (z) is defined in (37). As φ does not change signs and δP 0 > 0, the sign of the left-hand side is determined by W 0 (L): As ν > 0, ν = 0, and ν < 0, respectively, for each case, L * is the depth at which the stability of the trivial solution changes.
We consider the NLEP (59) by examining the linear operator It is not necessarily the case that (67) has a real or discrete set of eigenvalues. We can show, however, that given any nontrivial steady state solution P 0 , if an eigenvalue of L 1 is real, then it must be negative.
Proof. As in the proof of Theorem 5.1, it is still the case that δP 0 satisfies (62)-(63).
Repeating the initial steps of that proof, we obtain Although we can make no immediate assumptions about the oscillations of φ 0 (z), we suppose without loss of generality that φ 0 (0) > 0. If φ 0 (z) is positive for all z ∈ [0, L) (possibly zero at z = L), then setting d = L in (70) and integrating the right-hand side of (69) from z = 0 to z = L results in As the integrand on the left-hand side is positive on (0, L), ν ≤ 0.
This means that the only bifurcations that could occur from the trivial state would be associated with eigenvalues passing through the imaginary axis. Such an event, if it were to occur, would correspond to an oscillatory instability. 5.1. Summary. This work was inspired by the papers of Ebert et al. [6] and Huisman et al. [7]- [9]. The key question in this area is whether there is a critical depth for a phytoplankton bloom, and whether the associated depth profile is stable to various external influences. We follow many earlier investigators in looking at a simplified advection-reaction-diffusion model. Although this model has a non-local term, it can be readily dealt with using a simple augmentation of the system. We have introduced a new approach, based entirely on a geometric dynamical systems view of the equations and derived conditions for the existence and uniqueness of depth profiles. We have, moreover, showed that a consequence of our analysis is a restriction on any instabilities to be oscillatory.
The first natural extension is to consider either multi-species, as done in [15], or the inclusion of nutrients as in the work of Zagaris et al. [20]. In the latter case, an oscillatory instability (Hopf bifurcation) is found. A next step would be to extend our techniques to such more complex systems. The real challenging issue, however, is to see how these results hold up in the presence of more spatial effects.
A key assumption in almost all of the work to date is that the dynamics lies in a vertical column. This is, of course, unrealistic and a very interesting avenue of investigation would be to consider one, or even two, horizontal directions in the ambient fluid. Beyond this would lie the issues coming from considering an active fluid and explicitly coupling the phytoplankton dynamics with the fluid motion. This would almost certainly be only approachable by numerical means, but it would be hoped that the analytic results wherein a fluid whose motion turbulent mixing effects alone are encoded in the equations, such as done here, would form a useful guide as to what to expect.