Summer heatwaves promote blooms of harmful
cyanobacteria
KLAUS D. JO
¨
HNK
*
,1, 2
,JEFHUISMAN
*
,2
, JONATHAN SHARPLESw, BEN SOMMEIJERz,
PETRA M. VISSER
*
andJASPER M. STROOM§
*
Aquatic Microbiology, Institute for Biodiversity and Ecosystem Dynamics, University of Amsterdam, Nieuwe Achtergracht 127,
1018 WS Amsterdam, The Netherlands, wProudman Oceanographic Laboratory, University of Liverpool, 6 Brownlow Street,
Liverpool, L3 5DA, UK, zCenter for Mathematics and Computer Science, Postbus 94079, 1090 GB Amsterdam, The Netherlands,
§Water Board Rijnland, PO Box 156, 2300 AD Leiden, The Netherlands
Abstract
Dense surface blooms of toxic cyanobacteria in eutrophic lakes may lead to mass
mortalities of fish and birds, and provide a serious health threat for cattle, pets, and
humans. It has been argued that global warming may increase the incidence of harmful
algal blooms. Here, we report on a lake experiment where intermittent artificial mixing
failed to control blooms of the harmful cyanobacterium Microcystis during the summer
of 2003, one of the hottest summers ever recorded in Europe. To understand this failure,
we develop a coupled biological–physical model investigating how competition for light
between buoyant cyanobacteria, diatoms, and green algae in eutrophic lakes is affected
by the meteorological conditions of this extreme summer heatwave. The model consists
of a phytoplankton competition model coupled to a one-dimensional hydrodynamic
model, driven by meteorological data. The model predicts that high temperatures favour
cyanobacteria directly, through increased growth rates. Moreover, high temperatures also
increase the stability of the water column, thereby reducing vertical turbulent mixing,
which shifts the competitive balance in favour of buoyant cyanobacteria. Through these
direct and indirect temperature effects, in combination with reduced wind speed and
reduced cloudiness, summer heatwaves boost the development of harmful cyanobacter-
ial blooms. These findings warn that climate change is likely to yield an increased threat
of harmful cyanobacteria in eutrophic freshwater ecosystems.
Keywords: buoyant cyanobacteria, climate change, competition model, harmful algal blooms, heat-
wave, mixing, sinking phytoplankton, turbulence
Received 6 January 2006; revised version received 8 June 2007 and accepted 20 August 2007
Introduction
The summer of 2003 was probably the hottest summer
in Europe of the past 500 years (Levinson & Waple,
2004; Luterbacher et al., 2004). Temperature extremes
were locally 5 1C higher than on average (Beniston,
2004). Climate models indicate that the summer heat-
wave of 2003 might offer a glimpse of summers to be
expected in Europe in the later part of this century, as a
result of global warming (Beniston, 2004; Scha
¨
r et al.,
2004; Stott et al., 2004). During this hot summer, we ran a
large-scale lake experiment to control surface blooms of
the harmful cyanobacterium Microcystis.
Microcystis is a cosmopolitan cyanobacterium of eu-
trophic freshwaters (Reynolds et al., 1981; Zohary et al.,
1996; Chen et al., 2003; Verspagen et al., 2006). Cells of
Microcystis and several other cyanobacteria contain gas
vesicles (Walsby, 1994), providing them with buoyancy.
These buoyant cyanobacteria float upwards during
weak vertical mixing, and can form dense blooms at
the water surface (Zohary & Robarts, 1990; Visser et al.,
1996a). Moreover, Microcystis can produce microcystins,
a family of toxins damaging the liver of birds and
1
Present address: Leibniz-Institute of Freshwater Ecology and
Inland Fisheries, Alte Fischerhu
¨
tte 2, 16775 Neuglobsow, Germany.
2
KDJ and JH contributed equally to this work.
Correspondence: Jef Huisman, Aquatic Microbiology, Institute for
Biodiversity and Ecosystem Dynamics, University of Amsterdam,
Nieuwe Achtergracht 127, 1018 WS Amsterdam, the Netherlands,
tel. 131 20 5257085, fax 131 20 5257064,
e-mail: jef.huisman@science.uva.nl
Global Change Biology (2008) 14, 495–512, doi: 10.1111/j.1365-2486.2007.01510.x
r 2008 The Authors
Journal compilation r 2008 Blackwell Publishing Ltd 495
mammals, including humans. Hence, Microcystis
blooms provide a serious threat for water quality
(Chorus & Bartram, 1999; Codd et al., 1999; Carmichael,
2001; Huisman et al., 2005). Previous studies have
indicated that warm summers with little vertical mixing
provide ideal conditions for surface blooms of harmful
cyanobacteria (Reynolds, 1997; Ibelings et al., 2003;
Robson & Hamilton, 2003; Mooij et al., 2005). However,
the lack of large-scale lake experiments and detailed
mechanistic models has thus far limited rigorous pre-
dictions on the impact of weather conditions on cyano-
bacterial blooms.
Changes in vertical mixing may shift the competitive
balance between buoyant cyanobacteria and sinking
phytoplankton species (Walsby et al., 1997; Huisman
et al., 2004). In this paper, we report on Lake Nieuwe
Meer, a hypertrophic lake in the Netherlands, where
artificial mixing successfully controls Microcystis
blooms since 1993 (Visser et al., 1996a; Huisman et al.,
2004). The entire lake is vertically mixed throughout
summer by air bubbling. Artificial mixing led to major
changes in phytoplankton species composition. Surface
blooms of Microcystis dominating the lake before 1993
were replaced by a mixture of green algae and diatoms.
Artificial mixing of Lake Nieuwe Meer is accompanied
by an extensive monitoring program. Phytoplankton
and nutrients are sampled frequently. Temperature
loggers are moored in the lake to obtain online tem-
perature records. A microstructure profiler is employed
to monitor vertical mixing (Sharples et al., 2001; Huis-
man et al., 2004). We complemented the monitoring
program by a model study of competition between
buoyant cyanobacteria and sinking phytoplankton
species. We coupled this phytoplankton competition
model (Huisman et al., 1999a, 2004, 2006; Klausmeier &
Litchman, 2001) to a one-dimensional hydrodynamic
model (Rodi, 1993; Mohammadi & Pironneau, 1994;
Jo
¨
hnk & Umlauf, 2001; Hutter & Jo
¨
hnk, 2004). The
coupled biological–physical model is driven by meteor-
ological fields such as irradiance, wind speed, and air
temperature.
During the summer of 2003, we aimed to test a new
mixing regime in Lake Nieuwe Meer. We implemented
a scheme in which artificial mixing was switched on
and off with a 1- or 2-week periodicity, to reduce the
energy costs of artificial mixing without inducing sur-
face blooms of harmful cyanobacteria. However, when
we launched this research program, we were not aware
that the summer of 2003 would become the hottest
summer ever recorded in Europe. In August, at the
peak of the summer heatwave, there were almost in-
stant outbursts of Microcystis as soon as artificial mixing
was switched off. In this paper, we report on the failure
of intermittent mixing to suppress Microcystis blooms in
Lake Nieuwe Meer. We take advantage of our coupled
biological–physical model to unravel which mechan-
isms promoted the development of surface blooms of
these harmful cyanobacteria during the summer heat-
wave of 2003.
Lake experiment
Lake Nieuwe Meer is a recreational lake in the city of
Amsterdam, with a surface area of 1.3 km
2
, a mean depth
of 18 m, and a maximum depth of 30 m. The lake is
connected to the canals of Amsterdam at one side and
to a canal leading through the agricultural areas of the
Haarlemmermeer Polder on the other side. The lake has
very high nutrient concentrations (Table 1). Mean sum-
mer concentrations of total nitrogen in the epilimnion
were up to 3800 mgNL
1
in the late 1980s, and gradually
declined to 2400 mgNL
1
in 2006. Summer concentra-
tions of total phosphorus amounted to 450 mgPL
1
in the
late 1980s, and gradually declined to 260 mgPL
1
in 2006.
Dissolved inorganic nitrogen and soluble reactive phos-
phorus were never depleted to limiting values (Table 1;
compare with Sas, 1989; Wetzel, 2001). Summer concen-
trations of silica usually ranged from 700 to
3000 mgSiL
1
, with a lower concentration of 100 mgSiL
1
in August 1984. Accordingly, Lake Nieuwe Meer can be
classified as a hypertrophic lake, where nutrient limita-
tion of phytoplankton growth is negligible.
Dense surface blooms of the buoyant cyanobacterium
Microcystis aeruginosa dominated the lake during sum-
mer for many years. In 1993, Water Board Rijnland
installed a system of seven perforated air tubes just
above the lake sediment to start artificial mixing of the
lake by means of air bubbling. Artificial mixing led to
major changes in phytoplankton species composition
(Visser et al., 1996a; Huisman et al., 2004). Microcystis
was replaced by a mixture of diatoms and large green
Table 1 Summer nutrient concentrations (in mgL
1
) in Lake
Nieuwe Meer, measured monthly at 0.5 m depth in May–
August
Nutrient Mean SD Minimum Maximum
Total nitrogen 3595 920 1800 6450
Dissolved inorganic
nitrogen
2255 762 380 4100
Total phosphorus 406 130 200 740
Soluble reactive
phosphorus
343 109 130 610
Silicate 1532 1268 100 4000
Data show the mean, SD, minimum concentration and max-
imum concentration measured over the years 1980–2006.
496 K. D. JO
¨
HNK et al.
r 2008 The Authors
Journal compilation r 2008 Blackwell Publishing Ltd, Global Change Biology, 14, 495–512
algae. The diatoms that became dominant by artificial
mixing of Lake Nieuwe Meer consisted mainly of
Cyclotella and Stephanodiscus species. The green algae
that flourished with artificial mixing were dominated
by Scenedesmus species, with smaller amounts of Mono-
raphidium, Kirchneriella, and Dictyosphaerium. Also,
small flagellates (mainly Cryptomonas spp.) have be-
come more abundant after artificial mixing was in-
stalled in the lake.
Mixing experiment
Artificial mixing of an entire lake requires quite some
energy, and is, therefore, rather expensive. To reduce the
costs, we therefore proposed an intermittent mixing
regime, in which artificial mixing was alternately
switched on and off for periods of 14 days during the
summer of 2003. Artificial mixing of the lake was al-
ways off in winter. Artificial mixing was switched on at
the end of April 2003. The intermittent mixing regime
started on June 17, when artificial mixing was switched
off for 2 weeks. This period was followed by only 1
week of artificial mixing (July 2–9), due to a failure in
the air bubbling system. After this, the mixing was
switched off for a period of 2 weeks (July 10–24), and
subsequently switched on again for 2 weeks (July 25–
August 7). The next mixing-off period, in the second
week of August, was in phase with the peak of the
extreme heatwave across Europe. Within a few days,
Microcystis rapidly increased. Although Microcystis con-
centrations were still far below those experienced dur-
ing years without artificial mixing, the Water Board
decided to switch mixing on permanently from August
14 onwards, thus suppressing the unexpected Micro-
cystis bloom. The measurements ended on August 26.
Temperature structure
The temperature structure of the lake was monitored
using a permanently moored vertical array of six tem-
perature loggers (8-bit Minilog; VEMCO, Shad Bay,
Nova Scotia, Canada), placed equidistantly from the
water surface to 20 m deep with a spacing of 4 m. The
vertical array was mounted to a buoy near the middle of
the lake. The temperature loggers recorded the tem-
perature every 5 min over a period of 3 months, includ-
ing the full experimental period.
Turbulence structure
The vertical turbulence structure can be derived from
detailed profiles of the temperature microstructure (Im-
berger & Ivey, 1991; Kocsis et al., 1999; Sharples et al.,
2001). We employed a Self-Contained Autonomous
MicroProfiler (SCAMP; Precision Measurement Engi-
neering Inc., Carlsbad, CA, USA) to measure high-
resolution profiles of temperature, pressure (depth),
irradiance, chlorophyll fluorescence and conductivity
(Stevens et al., 1999). The free-falling SCAMP was
calibrated to sink downwards at a constant velocity
of 10 cm s
1
, collecting data at a frequency of
100 Hz. This yields a vertical resolution of 1 mm.
SCAMP profiles were run at three sampling stations
in the lake, on 16 sampling days spread over the
experimental season. At each sampling station, we took
10 consecutive SCAMP profiles within 80 min to
capture the variable nature of turbulent mixing. Turbu-
lent dissipation rates and vertical turbulent diffusivities
were calculated from the high-resolution temperature
profiles according to Sharples et al. (2001).
Phytoplankton development
Water samples were taken from 1 m and 8 m depth at
three sampling stations in the lake, on 17 days spread
over the experimental season. The samples were fixed
with Lugol’s Iodine. Phytoplankton cells were counted
microscopically using a Sedgewick-Rafter counting
chamber and identified to the species level. Population
abundances of the phytoplankton species, including
colonial species like Microcystis, were expressed as
number of cells per unit volume.
Model study
Complementary to the lake experiment, we developed a
model to predict changes in temperature structure,
turbulence structure, and phytoplankton population
dynamics. The model consists of a one-dimensional
hydrodynamic model coupled to a one-dimensional
phytoplankton competition model. This coupled biolo-
gical–physical model is forced by meteorological con-
ditions, thus enabling model studies of different climate
scenarios.
Hydrodynamic model
Simulations of temperature and turbulent diffusivity
are based on a one-dimensional ke turbulence model
(Rodi, 1993; Mohammadi & Pironneau, 1994; Jo
¨
hnk &
Umlauf, 2001; Hutter & Jo
¨
hnk, 2004). Because Lake
Nieuwe Meer is a relatively small lake, the Coriolis
effect is insignificant and pressure gradients generated
by internal seiching can be neglected. The turbulence
model is based on a system of five partial differential
equations, describing the dynamics of momentum, heat,
turbulent kinetic energy, and turbulent dissipation rate,
respectively. The depth of the lake is indicated by the
HEATWAVES PROMOTE HARMFUL CYANOBACTERIA 497
r 2008 The Authors
Journal compilation r 2008 Blackwell Publishing Ltd, Global Change Biology, 14, 495–512
variable z, where z runs from 0 at the surface to a
maximum depth z
m
at the bottom of the lake. The
bathymetry of the lake is implicitly taken into account
via the area–depth relation, A(z).
Horizontal momentum related to the horizontal ve-
locity vector, u, is driven by wind stress at the surface,
and vertically distributed by diffusion. The change in
horizontal momentum can be described by a system of
two partial differential equations (for momentum in the
x and y direction, respectively)
@u
@t
¼
@
@z
ðD
m
þ D
z
Þ
@u
@z
þ c
b
ujuj
1
A
@A
@z
; ð1Þ
where D
m
is the molecular diffusivity of momentum for
water and D
z
is the vertical turbulent diffusivity which
is space- and time-dependent. The second term on the
right-hand side describes the loss of momentum at each
depth induced by the lake’s boundaries. This boundary
stress term is formulated as a sliding law with a drag
coefficient c
b
, chosen according to the lake’s settings.
The two boundary conditions for the momentum equa-
tion are the continuity of shear stress at the surface
generated by wind stress, with a velocity-dependent
drag coefficient, and vanishing shear stress at the dee-
pest point of the lake (Hutter, 1993).
Heat is produced by the absorption of short-wave
radiation (400–3000 nm), E(z,t), and vertically distribu-
ted by diffusion. Accordingly, dynamic changes in
temperature, T, can be described as
@T
@t
¼
1
A
@
@z
AD
h
þ
D
z
s
h
@T
@z
1
rðTÞcðTÞ
@E
@z
; ð2Þ
where D
h
is the molecular diffusivity of heat in water,
and D
z
/s
h
is the turbulent diffusivity of heat expressed
as the turbulent diffusivity of momentum divided by the
Prandtl-number for heat, s
h
. The Prandtl-number is
parameterized as a function of turbulent kinetic energy,
turbulent dissipation rate, and temperature gradient. The
bathymetry of the lake is implicitly incorporated via the
area–depth relation, A(z), in such a way that heat is
conserved and radiation reaching the lake bottom is
absorbed by the sediment. We neglect heat exchange
between the water column and sediments, which is a
common assumption for deep lakes (unless the water is
heated from below by geothermal activity; e.g. Imboden
&Wu
¨
est, 1995). The second term on the right-hand side
describes the absorption of short-wave radiation, where
r(T) is the density of water and c(T) is the specific heat of
water, with both being functions of temperature (Chen &
Millero, 1986). Short-wave radiation in the water column
depends on the incident radiation at the water surface,
and decreases exponentially with depth according to
Lambert–Beer’s law. The two boundary conditions for
the heat equation are the continuity of heat flux at the
surface, driven by meteorological conditions, and a
vanishing heat flux at the deepest point of the lake
(Henderson-Sellers, 1984; Hutter & Jo
¨
hnk, 2004).
The dynamical distributions of turbulent kinetic en-
ergy, k, and turbulent dissipation rate, e, are described
by the two equations (Jones & Launder, 1972; Rodi,
1993; Mohammadi & Pironneau, 1994; Hutter & Jo
¨
hnk,
2004):
@k
@t
¼
@
@z
D
m
þ
D
z
s
k
@k
@z
þ P þ G e
@e
@t
¼
@
@z
D
m
þ
D
z
s
e
@e
@z
þ c
1e
P þ c
3e
G c
2e
eðÞ
e
k
;
ð3Þ
where P and G describe the production and loss of
turbulent kinetic energy induced by shearing of the flow
and buoyancy, respectively. Furthermore, s
k
and s
e
are
the Prandtl-numbers for turbulent kinetic energy and
turbulent dissipation rate, both assumed to be constant,
and c
1e
, c
2e
,andc
3e
are constants. Production and loss of
turbulent kinetic energy are functions of the velocity
gradient and the temperature gradient, respectively,
P ¼ D
z
@u
@z
2
; G ¼ yðTÞg
D
z
s
h
@T
@z
; ð4Þ
where y(T) is the temperature-dependent thermal expan-
sion coefficient of water, and g 5 9.81m s
2
is the Earth’s
gravitational acceleration. At the boundaries, the turbu-
lent fields are determined by heat and momentum fluxes
at the interface (Wilcox, 1993), extended for buoyant
fluids (Svensson, 1978). The constants and the Prandtl-
numbers of the ke turbulence model take their usual
values (Launder & Spalding, 1974; Rodi, 1993; Moham-
madi & Pironneau, 1994; Hutter & Jo
¨
hnk, 2004).
Finally, the vertical profile of turbulent diffusivity
when artificial mixing is switched off is calculated from
the relation
D
z
¼ c
m
k
2
e
; ð5Þ
where the constant of proportionality, c
m
, takes the
usual value of 0.09 (Rodi, 1993; Mohammadi & Piron-
neau, 1994; Hutter & Jo
¨
hnk, 2004). When artificial mix-
ing is switched on, the turbulent diffusivity is much
higher. Based on SCAMP measurements in 2003,
we estimated that artificial mixing contributed a
background turbulent diffusivity of 1.6 10
3
m
2
s
1
.
This value is added to Eqn (5) when artificial mixing
is on.
498 K. D. JO
¨
HNK et al.
r 2008 The Authors
Journal compilation r 2008 Blackwell Publishing Ltd, Global Change Biology, 14, 495–512
The full set of equations was integrated in time using
an implicit method on an equidistant space–time grid
with a spatial resolution of 0.1 m and a time resolution
of 4 min. Initial values of the hydrodynamic model were
a homogeneous temperature of 4 1C and zero values for
all other profiles on 1 January. For each simulation, the
hydrodynamic model was run over 3 years of meteor-
ological forcing to equilibrate the hydrodynamics with
the prevailing climate conditions. Only the last year of
each run was coupled to the phytoplankton competition
model.
Phytoplankton model
Simulations of the phytoplankton are based on a one-
dimensional phytoplankton competition model devel-
oped by Huisman and co-workers (Huisman et al.,
1999a, b, 2004, 2006). We consider a vertical water
column with a total number of n different phytoplank-
ton species. Because Lake Nieuwe Meer is a hyper-
trophic lake with very high nutrient concentrations
(Table 1), the model assumes that nutrient limita-
tion of the phytoplankton does not play a role. The
population dynamics of the phytoplankton species are
driven by light availability, temperature, and mixing
processes.
More precisely, let N
i
(z,t) denote the population den-
sity of phytoplankton species i at depth z and time t,
let I(z,t) denote the local light intensity in the PAR
range (400–700 nm), and let T(z,t) denote the local
temperature (in 1C). The population dynamics of
the phytoplankton species are described by a system
of reaction–advection–diffusion equations (Huisman
et al., 2004, 2006):
@N
i
@t
¼ m
i
ðI; TÞN
i
m
i
ðTÞN
i
þ
@
@z
v
i
ðTÞN
i
ðÞ
þ
@
@z
D
z
@N
i
@z
i ¼ 1; ...; n: ð6Þ
Here, m
i
(I,T) is the specific growth rate of species i as a
function of light intensity and temperature. The term
m
i
(T) is the specific loss rate of species i as a function of
temperature, and includes losses due to natural mor-
tality, zooplankton grazing, and virus attack. Further-
more, v
i
(T) is the vertical velocity of species i (with v
i
40
for buoyant species, and v
i
o0 for sinking species),
which depends on the viscosity of the water and, there-
by, on temperature. Finally, D
z
is the vertical turbulent
diffusivity, which is space- and time-dependent. Zero-
flux boundary conditions at the top and the bottom of
the water column assure that the phytoplankton species
neither enter nor leave the water column.
The underwater light conditions may change. For
instance, the water column will become more turbid
with increasing phytoplankton population densities.
More specifically, according to Lambert–Beer’s law,
the underwater light gradient can be described as
(Huisman & Sommeijer, 2002; Huisman et al., 2004):
Iðz; tÞ¼I
in
ðtÞð1 rÞ exp
Z
z
0
X
n
i¼1
k
i
N
i
ðs; tÞ
"#
ds K
bg
z
0
@
1
A
;
ð7Þ
where I
in
(t) is the incident light intensity (PAR) at the
water surface, r is a reflection coefficient to correct for
reflection losses at the water surface, k
i
is the specific light
attenuation coefficient of phytoplankton species i, K
bg
is
the background attenuation coefficient caused by all
nonphytoplankton components in the water column,
and s is an integration variable accounting for the non-
uniform phytoplankton population density distributions.
The model assumes that the specific growth rate is an
increasing saturating function of light intensity, as de-
scribed by a simple Monod equation (Huisman et al.,
1999a; Passarge et al., 2006):
m
i
ðI; TÞ¼
m
max;i
ðTÞI
m
max;i
ðTÞ=a
i
þ I
; ð8Þ
where m
max,i
(T) is the maximum specific growth rate of
species i at light-saturating conditions as a function of
temperature, and a
i
is the initial slope of the growth
function under light-limited conditions. The slope a
i
is
mainly determined by temperature-independent pro-
cesses of the light reaction of photosynthesis (like light
absorption), whereas m
max,i
is mainly determined by the
temperature-dependent processes of the dark reaction
of photosynthesis (Raven & Geider, 1988; Coles & Jones,
2000). We, therefore, assume that m
max,i
varies with
temperature, while a
i
is independent of temperature.
The temperature dependence of the maximum spe-
cific growth rate is described by an optimum curve,
which increases with temperature according to an Ar-
rhenius-type relationship, but decreases with tempera-
ture when the temperature optimum, T
opt,i
, is exceeded.
Our functional relationship is based on Robson &
Hamilton (2004), but reformulated to better describe
the point of the temperature optimum:
m
max;i
ðTÞ¼m
max;i
ðT
opt;i
Þ 1 þ b
i
R
TT
opt;i
1i
1
h
lnðR
1i
Þ
lnðR
2i
Þ
ðR
TT
opt;i
2i
1Þ
i
; ð9Þ
where m
max,i
(T
opt,i
) is the maximum specific growth rate
at the optimum temperature, and the parameters R
1i
,
R
2i
, and b
i
describe the form of the optimum curve for
HEATWAVES PROMOTE HARMFUL CYANOBACTERIA 499
r 2008 The Authors
Journal compilation r 2008 Blackwell Publishing Ltd, Global Change Biology, 14, 495–512
