Digital Object Identifier (DOI):
10.1007/s00285-002-0184-4
J. Math. Biol. 46, 333–346 (2003)
Mathematical Biology
Sandip Ghosal · Shreyas Mandre
A simple model illustrating the role of turbulence
on phytoplankton blooms
Received: 14 March 2002 / Revised version: 27 September 2002 /
Published online: 28 February 2003 – c© Springer-Verlag 2003
Abstract. The problem of the vertical distribution of phytoplankton is considered in the
presence of gravitational settling, turbulent mixing, population growth due to cell division
and a constant rate of loss due to predation and natural death. Nutrients are assumed to be
plentiful so that the production rate depends only on the light available for photosynthesis.
The non-linear saturation of plankton growth is modeled by allowing the attenuation rate
of light to be a linear function of the plankton density. The turbulent diffusivity is assumed
constant which corresponds to a mixed layer depth very much greater than the depth of light
penetration (euphotic depth). It is shown that an exact analytical solution of this non-linear
problem is possible for an idealized model in which the functional dependence of produc-
tion on light intensity is assumed to be a step function. Non-zero solutions are shown to
exist only if the parameters characterizing the system are above a certain critical curve in
a two dimensional parameter space. Numerical simulations using functional forms of the
production curve that resemble the measured photosynthetic response of plankton, show,
that the qualitative behavior of the system is similar to that of the idealized model presented.
Comparisons are made with other analytical approaches to the problem.
1. Introduction
The dynamics of plankton populations may be described [12] by the following par-
tial differential equation for the plankton density φ(z, t) as a function of the depth,
z and time, t :
∂φ
∂t
= Sφ − vp
∂φ
∂z
+
∂
∂z
(
kT
∂φ
∂z
)
. (1)
The water surface is at z = 0 and the z-axis is directed downwards. The “eddy
diffusivity coefficient” kT , represents the effect of turbulent mixing, and in general
could vary in space and time. The second term on the right hand side accounts for
gravitational settling at a speed vp relative to still water (vp > 0). The net growth
rate S is determined by various environmental factors. The boundary conditions
for φ are those of no flux at the water surface:
[
kT
dφ
dz
− vpφ
]
z=0
= 0, (2)
S. Ghosal, S. Mandre: Northwestern University, Department of Mechanical Engineering,
2145 Sheridan Road, Evanston, IL60208-3111 U.S.A.
e-mail: s-ghosal@northwestern.edu
Key words or phrases: Self-shading – Plankton – Nonlinear – Turbulence

334 S. Ghosal, S. Mandre
and either no flux at the lower boundary (if the water column is of finite depth, ‘h’)
[
kT
dφ
dz
− vpφ
]
z=h
= 0 (3)
or vanishing plankton density deep below the euphotic zone (if the water is con-
sidered infinitely deep)
φ(z → ∞) = 0. (4)
A reasonable model for the net growth rate is
S = P(I) − L (5)
where P(I) is the production, which in an eutrophic environment may be paramet-
rized by the local value of the light intensity, I , and L is a constant loss rate due to
grazing by zooplankton and other higher animals in the food chain and due to the
natural death of cells. The function P(I) has been measured for many species [9,
15]. It has the following qualitative behavior:
1. P(0) = 0,
2. P(I) increases monotonically with I until some saturation level characterized
by a value I = Ic,
3. P(I) is approximately constant for I > Ic, sometimes showing a slight decrease
with increasing values of I (photo-inhibition).
One of the earliest analytical models is due to Riley, Stommel & Bumpus
(henceforth RSB) [16], who assumed the form
S =
{
r if z < H ;
−L if z ≥ H (6)
where ‘H ’, the depth at which the net production is zero, was assumed known.
They looked for steady solutions of (1) with constant kT and the requirement that
φ(z) ≥ 0 for all z ≥ 0. When considered together with the boundary conditions
of zero flux at the surface and vanishing populations at infinite depth, it was found
that physically acceptable (non-zero) solutions exist only if (4rkT )/v2p > 1. Thus,
if the plankton reproduce too slowly, or sink too fast, or is subjected to too little
turbulence, steady populations cannot be sustained. Analytical solutions were also
shown to provide a reasonable fit to available oceanic data on the depth distribution
of plankton when the parameters were properly tuned. The formulation being lin-
ear, only the shape of the plankton concentration distribution could be determined
but not its amplitude.
A non-linear formulation capable of determining the amplitude was presented
by Shigesada & Okubo (henceforth SO) [17]. They allowed P(I) to be any positive
and non-decreasing function of the light intensity, I . The attenuation of light with
depth was described by
dI
dz
= −µI, (7)
where the total attenuation coefficient µ was assumed to be a linear function of
plankton density: µ = µ0 + µ1φ. The background attenuation coefficient, µ0, and