|Home | About | Journals | Submit | Contact Us | Français|
Several theories have been proposed to explain the development of harmful algal blooms (HABs) produced by the toxic dinoflagellate Karenia brevis on the West Florida Shelf. However, because the early stages of HAB development are usually not detected, these theories have been so far very difficult to verify. In this paper we employ simulated Lagrangian coherent structures (LCSs) to trace potential early locations of the development of a HAB in late 2004 before it was transported to a region where it could be detected by satellite imagery. The LCSs, which are extracted from surface ocean currents produced by a data-assimilative HYCOM (HYbrid-Coordinate Ocean Model) simulation, constitute material fluid barriers that demarcate potential pathways for HAB evolution. Using a simplified population dynamics model we infer the factors that could possibly lead to the development of the HAB in question. The population dynamics model determines nitrogen in two components, nutrients and phytoplankton, which are assumed to be passively advected by surface ocean currents produced by the above HYCOM simulation. Two nutrient sources are inferred for the HAB whose evolution is found to be strongly tied to the simulated LCSs. These nutrient sources are found to be located nearshore and possibly due to land runoff.
The largest and most frequent harmful algal blooms (HABs; also known as “red tides”) caused by the toxic dinoflagellate Karenia brevis tend to occur along the southern portion of the West Florida Shelf (WFS) [Steidinger and Haddad, 1981; Kusek et al., 1999; Brand and Compton, 2007]. Associated with typical HAB events is widespread death of sea life resulting from the brevetoxins produced by K. brevis [Bossart et al., 1998; Landsberg and Steidinger, 1998; Landsberg, 2002; Shumway et al., 2003; Flewelling et al., 2005]. Brevetoxins are also known to cause injury to humans through ingestion of contaminated seafood, skin contact, or inhalation of aerosolized brevetoxins in coastal regions [Backer et al., 2003; Kirkpatrick et al., 2004].
Dinoflagellates have slow growth rates [Brand and Guillard, 1981]. As a result, a necessary condition for achieving a high dinoflagellate population density is that mixing does not act to significantly reduce the population density during the growth phase. Indeed, dinoflagellate blooms often occur in enclosed basins. The tendency of K. brevis to cause large HABs along the southern WFS may then be explained by the presence of a cross-shelf transport barrier that has been revealed in the analysis of drifter trajectories [Yang et al, 1999] and in an appropriate synthesis of surface ocean currents produced by a HYCOM (HYbrid-Coordinate Ocean Model) simulation of the WFS [Olascoaga et al., 2006]. In Olascoaga et al.  it was hypothesized that the presence of such a transport barrier, which was characterized as a Lagrangian coherent structure (LCS), provides favorable conditions for the development of HABs in the southern WFS by inhibiting cross-shelf transport and thereby allowing for a substantial nutrient build up and dinoflagellate concentration on the shoreside of the barrier.
Satellite ocean color sensors have demonstrated ability to provide a very useful HAB monitoring capability [Tester et al., 1998; Stumpf et al., 2003; Tomlinson et al., 2004; Hu et al., 2005; de-Araújo-Carvalho et al., 2007]. To be detectable by current satellite sensors, however, the K. brevis cell concentration must be of the order of 105 cell l−1 (equivalently 1 mg Chl m−3) [Tester et al., 1998], which may take a period of about one month from background concentration to be attained [Tester and Steidinger, 1997; Steidinger et al., 1998]. As a consequence, satellite ocean color imagery by itself is not capable of detecting HAB initiation, and hence cannot be used to determine what factors favor HAB development, which is a subject of debate.
Earlier work [Gunter et al., 1947; Rounsefell and Nelson, 1966; Dixon and Steidinger, 2004] has associated the development of HABs on the WFS with river runoff, pollution, and other nutrient enriching phenomena produced by human activities. Consistent with this earlier work and the increase in the frequency of HABs, a recent study [Brand and Compton, 2007] has indicated that K. brevis concentrations are indeed higher near the shoreline than offshore, and that nutrient-rich freshwater runoff from land is likely to be a major source of nutrients for the development of HABs. Other investigators [Hu et al., 2006] have argued that, in addition to river runoff, submarine groundwater discharge, particularly when enhanced by an active hurricane season, should play a role in stimulating HABs on the WFS. It has been also suggested [Tester and Steidinger, 1997] that some HABs initiate several tens of kilometers offshore, near frontal regions, and that upwelling events provide a mechanism for supplying the required nutrients for their stimulation. Another, more involved, theory has been proposed [Lenes et al., 2001; Walsh and Steidinger, 2001; Walsh et al., 2003; Walsh et al., 2006] which explains HAB development as a sequence of events preceded by iron-rich Saharan dust deposition. In short, the theory assumes that the Saharan dust stimulates Trichodesmium blooms, which fix nitrogen. The blooms decompose on sinking and release dissolved organic nitrogen, which stimulates a K. brevis seed population near the bottom. Coastal upwelling is then required to help this population reaching the sea surface, where light inhibition is alleviated and a HAB finally develops. One further theory has been recently proposed by Stumpf et al.  who argued that the Mississippi River may constitute an important source of subsurface nutrients that, once upwelled, may stimulate HAB development on the WFS. However, as pointed out above, the early stages of the development of HABs are usually not detected, which makes it very difficult to test the aforementioned theories.
In this paper we demonstrate: (i) the utility of simulated LCSs in tracing potential early locations of the development of HABs on the WFS; and (ii) the feasibility of acquiring a better understanding of the environmental factors that could possibly lead to their development by incorporating the inferred early location information in a population dynamics model. This is done by focusing on a HAB that was observed on the southern portion of the WFS during October–December 2004 (Figure 1) [Hu et al., 2005]. This HAB was detected offshore using MODIS (Moderate Resolution Imaging Spectroradiometer) fluorescence line height imagery, and in situ water sampling and microscopic counts confirmed that it contained mostly K. brevis
In this section we discuss details of the two main mathematical tools used in the paper. Namely, LCSs, which are used to trace potential early locations of the development of the abovementioned HAB, and a population dynamics model, which is used to infer the factors that could possibly lead to its development. Both LCS computation and the population dynamics model make use of a surface ocean flow description provided by a HYCOM simulation of the WFS, and thus is described first.
Numerical model output provides a surface ocean flow description u(x, t), where x denotes position and t time, that is suitable for use in surface ocean LCS computation and plankton advection. Also, it has the advantage of allowing for a spatiotemporal coverage that is impossible to be attained with existing observational systems. In this paper we have chosen to consider daily surface ocean velocity fields extracted in a domain containing the WFS from a data-assimilative 0.04°-resolution HYCOM simulation of the Gulf of Mexico, itself nested within a 0.08°-resolution Atlantic basin data-assimilative HYCOM simulation [Hogan et al., 2007]. The model assimilates satellite-derived sea surface height anomaly and temperature data as well as available vertical profile data, and is driven by realistic high-frequency forcing obtained from the U.S. Navy NOGAPS operational atmospheric model. The topography used is derived from the ETOPO5 dataset, with the coastline following the 2-m isobath. The model includes river runoff.
LCSs are distinguished material fluid curves in unsteady flows [Haller and Yuan, 2000]. Two types of LCSs, repelling and attracting, can be identified depending on the behavior of initially nearby fluid particles straddling the LCS. Specifically, neighboring fluid particles straddling a(n) repelling (attracting) LCS diverge (converge) at an exponential rate in forward time. As a consequence, repelling and attracting LCSs delineate the boundary between immiscible fluid domains where stretching rates are low. This property of LCSs, particularly of attracting LCSs, is exploited here as it allows to trace low mixing regions within which HAB development is favored. The LCSs cannot be guessed visually by inspecting velocity field snapshots. Rather, an appropriate tool is required to extract them from a given velocity field.
Computation of finite-time Lyapunov exponents (FTLEs) constitutes a practical and robust tool for identifying LCSs [Haller, 2001, 2002; Shadden et al., 2005; Lekien et al., 2005; Shadden et al., 2006; Mathur et al., 2007; Lekien et al., 2007; Green et al., 2007; Lekien and Coulliette, 2007]. The FTLE is a measure of the finite-time averaged separation rate of initially nearby fluid particle pairs, and is defined as
where denotes spectral norm and x(t0 +τ; t0, x0) is the position at time t = t0 +τ of a fluid particle that at time t = t0 was located at x0. The latter is obtained by integrating the particle trajectory equation regarded as a nonautonomous dynamical system obeying
where the overdot stands for time differentiation. Regions of maximum separation rates produce local maximizing curves or “ridges” in the FTLE field, which approximate LCSs. When the FTLE is computed by integrating particle trajectories forward (backward) in time, τ > 0 (τ < 0), a ridge in the FTLE field corresponds to a(n) repelling (attracting) LCS.
The FTLE computations reported in this paper are carried out numerically using the software package MANGEN, a dynamical systems toolkit designed by F. Lekien that is available at http://www.lekien.com/~francois/software/mangen. The general FTLE computation scheme is as follows: (i) the particle trajectory equation (2) is integrated for a grid of particles at time t0 to get their positions at time t0 + τ; (ii) the spatial derivatives in (1) are computed at each point in the initial grid by directly differencing with neighboring grid points; and (iii) the FTLE is computed at each point in the initial grid by evaluating (1). The previous three steps may be repeated for a range of t0-values to produce a time series of FTLE field. For the FTLE computations reported below we employ a fourth-order Runge–Kutta integrator with a fixed 1-h time step and a bilinear method for the required spatiotemporal interpolation of the numerically generated velocity; central differences are considered for the evaluation of spatial derivatives. FTLE fields are evaluated on a grid of roughly 4×104 points regularly distributed over a portion of the WFS extending from 27°N to 24°N, and 83.5°W to 81°W or the shoreline (Figures 2--33).
A number of biological and physical factors contribute to the initiation, growth, maintenance, and demise of HABs produced by K. brevis. We consider nutrient supply and surface ocean circulation constraints as the most important factors.
Consistent with the above assumptions, we formulate a population dynamics model in which the nitrogen concentration is determined in two compartments: nutrients, N(x, t), and phytoplankton, P (x, t), which is partly justified by the lack of K. brevis grazers [Steidinger, 1973]. We further assume that these densities are passively advected by surface ocean currents. The equations describing the evolution of N and P take the form [Riley, 1963; Wroblewski and O'Brien, 1976; Wroblewski et al., 1988]
In the above equations nutrient-limited phytoplankton growth is assumed to have a maximum rate μ = 0.4 d−1 and a half-saturation constant N0 = 0.5 mmol N m−3. The linear phytoplankton loss is parametrized by γ = 0.1 d−1. The assumed proportion of recylable nutrients in dead phytoplankton ε = 0.5. These parameter values are similar to those considered in earlier works [Steidinger et al., 1998; Walsh and Steidinger, 2001; Walsh et al., 2001, 2002, 2003], and result in a net population growth rate of the order of 0.25 d−1, in agreement with field observations [van Dolah and Leighfield, 1999]. No attempt is made to optimize them in order to improve the likeness between model and observations.
Equations (3) are solved using the method of the characteristics by tracing 104 passively advected particles carrying P and N. These passively advected but biologically active particles are released every other day on a uniformly distributed lattice extending from 27°N to 26°N and from the shoreline to roughly the 10-m isobath (Figure 4). Trajectory and population dynamics equations for each particle are integrated using a fixed 1-h time step fourth-order Runge–Kutta scheme; integration of trajectory equations for each particle additionally involves interpolation of the numerically generated velocity in both time and space, which is accomplished by using a bilinear method. For the simulation reported below we assume that each released particle initially carries a P-value of 0.01 mg Chl m−3 (the conversion equivalence for K. brevis is taken to be 1 mg Chl m−3 = 0:38 mmol N m−3 [Shanley and Vargo, 1993]). This P-value corresponds to a K. brevis cell concentration of 103 cell l−1, which lies within the range of typical background K. brevis cell concentrations found on the WFS [Geesey and Tester, 1993; Tester and Steidinger, 1997]. The initial N-value distribution on each set of released particles in the above lattice is nonuniform, and is chosen consistent with the results of the LCS analysis as explained below.
Figures 2 and and33 show FTLE field snapshots computed backward in time by evaluating (1) with τ = −60d. The regions of most intense red tones indicate attracting LCSs which conform a portion of the cross-shelf transport barrier described in Olascoaga et al. . The regions of less intense red through most blue tones also indicate attracting LCSs, albeit of smaller scale, which constitute the focus of this paper. Note the peculiar presence of a banana-shaped region demarcated by LCSs of the latter type, which is most clearly evident on mid December 2004 but can be traced in time from late October 2004. This fluid region is isolated from the surrounding fluid by material barriers, and is characterized by low stretching rates, which implies low mixing and thus favorable conditions for bioaccumulation within the region. Consequently, it constitutes a potential pathway for HAB development. Note that the region in question coincides remarkably well with the area spanned by the HAB described in Hu et al.  from mid November 2004 onward (Figure 1, panels e–i). The noted agreement between simulated LCSs and the boundary of the observed HAB distribution is exploited below to trace locations of potential early development of the observed HAB. This involves performing backward- and forward-time integrations of fluid particle positions, which are overlaid on the FTLE fields shown in Figures 2 and and33.
Consider first the particle positions lying within the aforementioned banana-shaped region of low mixing that resembles the area occupied by the HAB in question on 13 December 2004 (Figure 2, right panel). Immersed in the cloud of particles, which are indicated with black dots, is one particle indicated with a bold magenta dot that is used for reference in tracing the early location of the observed HAB. The position of this fiducial particle is chosen to lie approximately centered within the banana-shaped low mixing region on 13 December 2004. Integrating backward in time (2) during 30 d with this initial condition, it is found that the the early location of the HAB distribution in question could have restricted to a nearshore area between Estero Bay and Cape Romano (Figure 2, left panel). The integration time corresponds to the time it typically takes for K. brevis to develop a HAB from background concentration. Note that the fiducial particle remains for all dates shown in Figure 2 within the LCSs that delimit the banana-shaped low mixing region. The particle positions indicated with black dots are the result of forward-time integration of (2) during 30 d with initial positions surrounding that of the fiducial particle on 13 November 2004 in a nearshore region between Estero Bay and Cape Romano, which lies within the LCSs that bound the banana-shaped low mixing region on that date. Note that as time progresses from 13 November 2004 the particles fill the banana-shaped low mixing region while it stretches and bends. Note also that during the relevant time interval the particles approach, but do not traverse, the FTLE ridges that bound the banana-shaped low mixing region, which provides strong support for our assertion that the simulated LCSs are material fluid barriers.
Consider now the fiducial particle shown in the right panel of Figure 3. This fiducial particle is located in the center of the observed HAB distribution on 13 November 2004 (Figure 1, panel d), which is situated north of the banana-shaped low mixing region on that date. Integrating (2) backward in time during 30 d with this initial position, it is found that the early location of the HAB distribution in question could have restricted to an area near the mouth of the Caloosahatchee River (Figure 3, left panel). Integrating (2) forward in time during 30 d for a cloud of particles with initial positions surrounding that of the fiducial particle on 14 October 2004, it is found that these particles approximately cover the area spanned by the observed HAB on 13 November 2004. We emphasize that the motion of the particles is completely determined by the underlying LCSs. Note, in particular, the LCS on 14 October 2004 which appears to originate in the Caloosahatchee River mouth area, right before the aforementioned banana-shaped low mixing region starts to develop. Clearly, particles initially located near the shoreline north of the Caloosahatchee River mouth will never enter the banana-shaped low mixing region. Contrarily, these particles will disperse over a region characterized by a highly intricate underlaying LCS. It must be noted, however, that unlike the banana-shaped observed HAB distribution, which can be explained as a result of bioaccumulation within a confined low mixing fluid region, the HAB distribution in question can be explained as a result of folding of fluid elements in a region characterized by high FTLE values, which can effectively contribute to bioaccumulation but only during short intervals of time.
We now turn to the description of the results of a simulation based on the population dynamics model (3).
Recall that the population dynamics model equations (3) are solved by tracing particles. These particles are released every two days on a lattice spanning the domain shown in Figure 4. Recall also that it is assumed that the particles carry initially a unique P-value representing a background K. brevis concentration on the WFS. Nutrient sources are devised by attaching nonoligotrophic N-values only to particles of the above lattice in the vicinity of the two nearshore locations identified in the previous subsection using simulated LCSs. More precisely, to each particle released near the shoreline in the vicinity of the Caloosahatchee River (Figure 4, red subdomain) and between Estero Bay and Cape Romano (Figure 4, green subdomain) we initially attach an N-value selected according to
where t ranges from 1 October 2004 (start date of the simulation) through 13 December 2004 (end date of the simulation). Specific parameter values for particles released in the vicinity of the Caloosahatchee River are chosen to be Nmax = 13 mmol N m−3, tmax = 1 October 2004, and α = 0.02 d−2. Specific parameter values for particles released nearshore between Estero Bay and Cape Romano are assumed to be given by Nmax = 2 mmol N m−3, tmax = 13 November 2004, and α = 0.02 d−2. The basis for these parameter value choices is observational, and is explained below. To every particle released in locations different than the above two locations (Figure 4, gray domain) we initially attach an N-value of 0.1 mmol N m−3, which corresponds to normal oligotrophic conditions found on the WFS [Masserini and Fanning, 2000].
Figure 5 shows a sequence of snapshots of the evolution of phytoplankton distribution as simulated according to the population dynamics model (3). Note that the intense observed HAB distribution on 13 November 2004 (Figure 1d) is reasonably well-reproduced by the simulation (Figure 5, upper-right panel). Note also that the peculiar banana-shaped distribution acquired by the observed HAB from mid November 2004 onward (Figure 1, panels e–i) is well-captured by the population dynamics model (Figure 5, lower panels).
The noted good agreement between simulated HAB and observed HAB most salient features demonstrates that the observed HAB is quite strongly tied to the simulated LCSs. The results of the population dynamics model simulation are not found to be sensitive to small variations of model parameter values or the addition of a small amount of diffusion in (3). Rather, they are found to be sensitive to the location of the nutrient sources. This is consistent with results of Haller , who showed that LCSs are robust under small perturbations to numerically generated velocity data, and Lekien and Coulliette , who showed that tracer dynamics can be governed by advective LCSs rather than diffusion in quasiturbulent flows. Only microscopic tracer sources, such as those induced by fast chemical reactions [Giona et al., 2002; Károlyi and Tél, 2005], may overcome LCS action. To infer the factors that probably led to the development of the observed HAB, an interpretation of the characteristics of the nutrient sources considered in the population dynamics simulation is required.
For particles released in the vicinity of the Caloosahatchee River, the attached nonoligotrophic initial nutrient concentrations can be directly linked with nutrient-rich water from land runoff. The assumed nutrient concentration values are roughly adhere to measurements of nitrate-plus-nitrite concentration taken near the mouth of the Caloosahatchee River, which are fully utilizable (Figure 6, left panel). The measured nitrate-plus-nitrite concentration values are well-correlated with measured Caloosahatchee River flow discharge rate values (Figure 6, right panel). Note the high nitrate-plus-nitrite concentrations and the corresponding high flow discharge rates in late 2004. The latter are attributed to high precipitation rates associated with an active hurricane season.
For particles released between Estero Bay and Cape Romano, association of the assumed nonoligotrophic initial nutrient concentrations with nutrient loading into coastal waters as a result of land runoff is less direct, but cannot be ruled out. Indeed, images of colored dissolved organic matter (CDOM) appear to indicate that the nearshore region between Estero Bay and Cape Romano is reached by the Caloosahatchee River plume on the dates of the population dynamics simulation (see panel c of Figure 4 in Hu et al. ). An alternative hypothesis is that the Estero Bay–Cape Romano nutrient source could be associated with an upwelled intrusion of nutrient-rich Gulf of Mexico slope water on the WFS. This hypothesis is marginally compatible with the HYCOM simulation, consistent with some evidence of the occurrence of alongshore upwelling-favorable winds on the dates of the population dynamics simulation. On the other hand, although not conclusive, some evidence of upwelling can be found in satellite-derived sea surface temperature maps. However, there also exists in situ observational evidence indicating that salinity is low on the coastal area between Estero Bay and Cape Romano, which is consistent with the Caloosahatchee River plume hypothesis. It must be nonetheless noted that these hypotheses are not necessarily mutually excluding [Hu et al., 2004], and that the nature of their observational foundation is only mostly qualitative, which cannot be used to provide a strong support to the particular choice (4). In the absence of direct nutrient concentration measurements, the assumed maximum nutrient concentration value for the Estero Bay–Cape Romano source (2 mmol N m−3) corresponds to a diluted fraction of that measured near the Caloosahatchee River mouth (approximately 13 mmol N m−3) roughly estimated using salinity measurements.
Finally, it must be mentioned that although the main characteristics of the observed HAB distribution are well-reproduced by the population dynamics simulation, some discrepancies between simulated and observed HAB distributions are seen. For instance, note the nearshore simulated HAB distribution from 22 November 2004 onward (Figure 5, lower panels) for which there is no observational support (cf. Figure 1, panels f–i). Possible causes of these discrepancies include the following. First, the HYCOM simulation considered is not expected to accurately describe details of the coastal surface ocean circulation, and consequently also details of the associated LCS pattern, which strongly constrain plankton dispersal. Second, the population dynamics model (3) simulates the evolution of K. brevis only without contemplating the possibility that other phytoplankton species outcompete K. brevis for the same nutrient resources. Third, possible inhomogeneities in the initial K. brevis distribution are not considered in the simulation presented. Forth, three-dimensional effects are not simulated by population dynamics equations (3), which may play a role in determining the surface ocean distribution of the observed HAB.
In this paper we have traced the early stages of the development of a harmful algal bloom (HAB) produced by Karenia brevis which was observed on the southern portion of the West Florida Shelf (WFS) during OctoberDecember 2004. This HAB was detected offshore using ocean color imagery and field sampling. Potential early development locations of the observed HAB were traced with the aid of simulated Lagrangian coherent structures (LCSs). The LCSs determined fluid particle pathways which highly constrained the evolution of the HAB in question. The factors that could probably lead to the development of this HAB were inferred using a simplified nutrient–phytoplankton population dynamics model. Both LCS computation and population dynamics modeling were carried out using a surface ocean velocity description provided by a data-assimilative HYCOM (HYbrid-Coordinate Ocean Model) simulation. Two nearshore nutrient sources were inferred as responsible for the development of a simulated HAB with qualitatively similar characteristics as those of the event observed on the WFS in late 2004. One source could be directly associated with nutrient-rich water of land runoff origin as it was located at the mouth of the Caloosahatchee River. The other source was located near the shoreline between Estero Bay and Cape Romano, but it could be possibly remotely linked with the Caloosahatchee River, which calls for a detailed survey of its plume. Our results suggest that land runoff is an important environmental factor contributing to HAB development on the WFS. However, a deeper understanding of the development of HABs on the WFS demands that the importance of this factor be evaluated against that of other factors, such as Saharan dust deposition or groundwater discharge, which were not analyzed in this paper.
We are grateful to G. Halliwell and O. Smedstad, for proving the HYCOM model output, M. Brown and S. Smith, for critically reading an early version of the manuscript, I. Rypina and I. Udovydchenkov, for the benefit of discussions on dynamical systems theory, J. Triñanes and M. Framiñán for advice on the interpretation of satellite images, and three anonymous reviewers, for constructive criticism. MJO was supported by NSF grant CMG0417425, PARADIGM NSFONR NOPP grant N0000140210370, NSF grant OCE0432368, and NIEHS grant P50ES12736. FJBV and HK were supported by NSF grant CMG0417425. LEB was supported by NSF grant OCE0432368 and NIEHS grant P50ES12736.