pmc.ncbi.nlm.nih.gov

Temperature and residence time controls on an estuarine harmful algal bloom: Modeling hydrodynamics and Alexandrium fundyense in Nauset estuary

. Author manuscript; available in PMC: 2016 Nov 1.

Published in final edited form as: Estuaries Coast. 2015 Mar 17;38(6):2240–2258. doi: 10.1007/s12237-015-9949-z

Abstract

A highly resolved, 3-d model of hydrodynamics and Alexandrium fundyense in an estuarine embayment has been developed to investigate the physical and biological controls on a recurrent harmful algal bloom. Nauset estuary on Cape Cod (MA, USA) consists of three salt ponds connected to the ocean through a shallow marsh and network of tidal channels. The model is evaluated using quantitative skill metrics against observations of physical and biological conditions during three spring blooms. The A. fundyense model is based on prior model applications for the nearby Gulf of Maine, but notable modifications were made to be consistent with the Nauset observations. The dominant factors controlling the A. fundyense bloom in Nauset were the water temperature, which regulates organism growth rates, and the efficient retention of cells due to bathymetric constraints, stratification, and cell behavior (diel vertical migration). Spring-neap variability in exchange altered residence times, but for cell retention to be substantially longer than the cell doubling time required both active vertical migration and stratification that inhibits mixing of cells into the surface layer by wind and tidal currents. Unlike in the Gulf of Maine, the model results were relatively insensitive to cyst distributions or germination rates. Instead, in Nauset, high apparent rates of vegetative cell division by retained populations dictated bloom development. Cyst germination occurred earlier in the year than in the Gulf of Maine, suggesting that Nauset cysts have different controls on germination timing. The model results were relatively insensitive to nutrient concentrations, due to eutrophic conditions in the highly impacted estuary or due to limitations in the spatial and temporal resolution of nutrient sampling. Cell loss rates were inferred to be extremely low during the growth phase of the bloom, but increased rapidly during the final phase due to processes that remain uncertain. The validated model allows a quantitative assessment of the factors that contribute to the development of a recurrent harmful algal bloom and provides a framework for assessing similarly impacted coastal systems.

Keywords: harmful algal bloom, hydrodynamic-biological model, Alexandrium fundyense, residence time, cyst germination, growing degree day

1. INTRODUCTION

Harmful algal blooms (HABs) are a global issue, causing significant public health, economic, and ecological problems (Hallegraeff, 1993; Anderson et al., 2012). Paralytic shellfish poisoning (PSP), caused by blooms of the dinoflagellate Alexandrium fundyense, is one example of a HAB human poisoning syndrome that frequently impacts the northeastern United States. A. fundyense blooms regularly occur over large regions of the Gulf of Maine and separately in embayments and estuaries along the coast (Anderson, 1997). Coastal states maintain extensive monitoring networks to minimize public health impacts, but contaminated and quarantined shellfish resources have significant negative impacts on local and regional economies (Shumway, 1990; Hoagland and Scatasta, 2006). Understanding the fundamental processes governing the development of A. fundyense blooms will allow more informed management of the public health risks and better predictive capabilities for how HABs may respond under changing climate conditions and anthropogenic inputs.

Previous studies of the physical and biological processes controlling the A. fundyense bloom in the Gulf of Maine led to the development of a physical-biological model (Franks and Anderson, 1992; Anderson, Keafer, et al., 2005; McGillicuddy et al., 2005; Stock et al., 2005; He et al., 2008). In turn, this model has been used to evaluate processes that influence bloom development, including basin-scale circulation, wind and river discharge patterns, and cyst and nutrient distributions (McGillicuddy et al., 2005, 2011; Li et al., 2009). The model has undergone continuous refinement and skill assessment as part of an effort to produce weekly and annual forecasts of HAB intensity in the region, and planning is underway to make the Gulf of Maine forecast model operational (McGillicuddy et al., 2011; R. He pers. comm.).

In addition to the large-scale Gulf of Maine bloom, independent blooms of A. fundyense occur in estuaries and bays along the northeast coast (Anderson, 1997; Anderson and Rengefors, 2006; Hattenrath et al., 2010; Borkman et al., 2014). One example is Nauset estuary on Cape Cod, MA, an area that experiences recurrent A. fundyense blooms and PSP toxicity (Anderson and Stolzenbach, 1985; Crespo et al., 2011; Ralston et al., 2014). Major blooms in Nauset led to shellfishing closures in 21 of 23 years from 1992 to 2014 (Massachusetts Division of Marine Fisheries). These blooms tend to occur earlier in the year than those in the Gulf of Maine and are localized around three ponds that form the uppermost termini of the estuary (Anderson, 1997; Crespo et al., 2011). Thus, Nauset frequently hosts multiple, concurrent but independent A. fundyense blooms each year, making it a natural laboratory for the investigation of A. fundyense bloom dynamics. Observations and model development from Nauset can inform predictive models of similar blooms occurring in shallow, stratified coastal environments around the world (Ralston et al., 2014; Raine, 2014).

Accurate modeling of HABs remains a significant challenge. Some HAB forecast models have demonstrated skill at categorical predictions of impacts on human health (Stumpf et al., 2009), but predictive skill of organism abundance is difficult. Bloom development depends on both physical and biological processes, so model uncertainties in water properties or circulation limit the ability to predict cell concentrations. For example, the A. fundyense model for the Gulf of Maine demonstrated skill in hindcasts of several years of blooms (He et al., 2008; Li et al., 2009), but uncertainties in water properties at the open boundaries led to overprediction of a bloom in 2010, a year when no significant concentrations of A. fundyense were observed (McGillicuddy et al., 2011). Similarly, concentrations of A. tamarense in the St. Lawrence estuary were not simulated accurately because they depended on plume advection processes that were poorly constrained by the model (Fauchot et al., 2008).

In this work, we develop, evaluate, and analyze a physical-biological model of A. fundyense in Nauset estuary. The aims are to use the model to better understand processes contributing to bloom development in Nauset and to use observations to evaluate and modify the A. fundyense model formulation so that it may be useful in similar coastal settings. Observational studies in Nauset in recent years (Crespo et al., 2011; Ralston et al., 2014; D. Anderson, unpub. data) provide extensive data to help constrain uncertainties in both physical and biological parameters. The A. fundyense component of the model builds off work in the Gulf of Maine (Stock et al., 2005; He et al., 2008), but is modified in accordance with observations from Nauset. The approach develops a single set of model parameters that are applied to conditions in multiple years, thereby minimizing over-fitting of the model to observations. Similarly, we strive to reduce model complexity that it is not supported by observational data. The model results are used to quantify factors controlling bloom development, including retention of cells due to physical processes and organism behavior, dependence on the distribution and emergence rate of cysts, and scaling with loss rates at bloom termination using a growing degree-day approach.

2. METHODS

2.1 Study location: Nauset estuary

Nauset estuary is a complex of marshes and submerged tidal kettle ponds on Cape Cod, Massachusetts connected by an inlet through a dynamic barrier beach to the Atlantic Ocean (Fig. 1). The central estuary has a network of tidal channels through a vegetated marsh platform. Three tidal kettle ponds (or salt ponds) - Salt Pond, Town Cove, and Mill Pond (north-to-south) - form the uppermost termini of the estuary and are connected to the central marsh by tidal channels. The maximum depths of the salt ponds range from 6 to 11 m, while the marsh tidal channels are typically 1 to 3 m deep. Tidal forcing is mainly semidiurnal with a range from 1 to 2 m.

Figure 1.

Figure 1

Nauset estuary location, bathymetry, and model grid detail. Sampling stations from the marsh-wide surveys are shown on the main map (gray circles), with the larger dots showing the locations of thetime series in Fig. 7. A detail of Salt Pond illustrates the model grid resolution, and also shows thelocations of transects A-A’ and B-B’ for field observations (casts at gray circles) and model results (redline) shown in Figs. 5 and 8. Note that the model bathymetry extends above mean higher high water toallow wetting and drying of the marsh, while plots of interpolated observations (Fig. 11) are shown withthe coastline as the mean sea level contour to demark the navigable channels.

Nauset estuary has no river input, so the primary sources of freshwater are groundwater and direct precipitation. Groundwater discharge from the Nauset and Monomoy lenses also provides a major source of nutrient loading to the estuary from a densely residential watershed with septic wastewater inputs (Giblin and Gaines, 1990; Portnoy et al., 1998; Colman and Masterson, 2008). Mean monthly precipitation from February to May is about 10 cm (National Climate Data Center (NCDC), Chatham Municipal Airport). The freshwater input and solar heating are sufficient to create stratification in the salt ponds due to both salinity and temperature gradients. Surface-to-bottom salinity gradients in the ponds were typically ~1-2 psu, and thermal stratification ranged from about 1 °C at the beginning of the bloom to greater than 10 °C toward the end (Crespo et al., 2011; Ralston et al., 2014).

2.2 Physical-biological model

2.2.1 Hydrodynamic model

The hydrodynamic model was developed using the Finite Volume Coastal Ocean Model (FVCOM) (Chen et al., 2003). FVCOM is discretized horizontally with an unstructured grid that allows variable resolution through the domain to resolve bathymetric complexity. The Nauset grid had node spacing ranging from a minimum of less than 10 m in the estuary to 4 km on the open boundary (Fig. 1). The model grid extends about 25 km north and south of Nauset inlet and 25 km offshore. FVCOM uses sigma layers vertically. The results presented here use 21 sigma layers; a 31 layer model was also tested but did not yield substantively different results. FVCOM incorporates surface heat fluxes using the TOGA-COARE bulk air-sea flux algorithm (Fairall et al., 1996, 2003). Turbulence closure is with the General Ocean Turbulence Model (Umlauf and Burchard, 2005), in these simulations using the k-ε model (Rodi, 1987; Burchard and Baumert, 1995).

In shallow estuarine flows, the hydrodynamic processes are extremely sensitive to bathymetry. High-resolution bathymetry was incorporated into the Nauset model from lidar-derived topographic maps of Cape Cod National Seashore from the U.S. Geological Survey (USGS) (Brock et al., 2007). Bathymetry in subtidal regions too deep for lidar penetration was based on acoustic surveys during this study and from previous observations by investigators from the USGS (Cross et al., 2006) and Woods Hole Oceanographic Institution (WHOI) (Aubrey et al., 1997).

Model boundary conditions were based primarily on observations. Water level at the open boundary was taken from the NOAA station at Chatham (#8447435), located 15 km south of Nauset Inlet. Salinity and temperature in the coastal ocean were based on buoy data in Massachusetts Bay from the Gulf of Maine Ocean Observing System (part of the Northeastern Regional Association of Coastal Ocean Observing Systems). Meteorological conditions, including wind speed and direction, air temperature, irradiance, and relative humidity, were from the Massachusetts Department of Environmental Quality monitoring station in Truro (17 km north of Nauset) and, when Truro data were unavailable, the Martha's Vineyard Coastal Observatory (WHOI). Groundwater fluxes were from a groundwater model developed by the USGS that covered the northern half of the estuary (Colman and Masterson, 2008); fluxes into the southern half of the estuary were based on the average groundwater flux in the northern half but made proportional to watershed area. Direct precipitation was added at the surface based on observations from the Chatham Municipal Airport (NCDC, #725069, 3 km east of the tidal station).

2.2.2 A. fundyense model

Although the physical and biological models are described separately, they were run as a single, integrated system. The A. fundyense model was based substantially on models previously developed for the Gulf of Maine (Stock et al., 2005; He et al., 2008). The biological model calculates germination and growth rates of A. fundyense as a function of environmental conditions including temperature, nutrients, and light based on laboratory experiments (Stock et al., 2005). To adapt the Gulf of Maine model framework to Nauset, the basic approach was to simplify and adjust the model to the extent consistent with observations. Model parameters from the calibration to Nauset observations are listed in Table 1, and others not listed are as described in Stock et al. (2005).

Table 1. A. fundysense model parameters.
Parameter Units Range of values Description
Gmax d−1 1.0-1.6 Maximum instantaneous growth rate
KN μM 0-0.5 Half saturation constant for nitrate
m d−1 0.2-0.5 Mortality
ws m d−1 10 Maximum vertical swimming speed
α g d−1 W−1 m2 0.055 Growth efficiency
kw m−1 0.2-0.5 Diffuse attenuation of light in water

Several important modifications were made to the A. fundyense model. For cyst germination, an endogenous clock places a rigid constraint on the timing of bloom initiation in the Gulf of Maine (Anderson and Keafer, 1987; Matrai et al., 2005). Peak germination rates in the Gulf of Maine occur April through June, but in Nauset the blooms are often well underway by mid-March (Crespo et al., 2011; Ralston et al., 2014). Observations of cysts from Perch Pond, a coastal embayment similar to Nauset, found that cyst germination was not controlled by an endogenous clock, and that instead cysts could germinate any time of year provided a suitable temperature regime (Anderson and Keafer, 1985; Anderson, 1997). Similarly, excystment experiments on A. minutum and A. tamarense cysts from Cork Harbor found no evidence of an endogenous clock, but that temperature did affect germination rate (Ní Rathaille and Raine, 2011). Given the absence of evidence of an endogenous clock for cysts in inshore systems like Nauset, that component was removed from the model. Temperature dependence of cyst germination was retained as in the Gulf of Maine model (Stock et al., 2005). Cyst distributions from benthic surveys of Nauset estuary were mapped in the falls of 2008, 2009, and 2011 (Crespo et al., 2011) and used to force the bottom boundary condition. Model sensitivities to the cyst germination rate and benthic cyst distribution are addressed in the results section.

As the water column is seeded by germinating cysts, additional A. fundyense cell concentration increases occur through vegetative cell division. Maximum growth rates depend on water temperature, irradiance, and nutrient availability, as in the Gulf of Maine model. The nutrient dependence was assumed to be only on the combined nitrate and nitrite concentration, and concentrations were specified based on observations from the weekly surveys. Nauset estuary is highly impacted by anthropogenic inputs, so nutrient concentrations are generally high (Giblin and Gaines, 1990; Colman and Masterson, 2008; S. Fox unpublished). Simulations were run using a range of half-saturation constants for nitrate (kN) from the literature. Model results using relatively low kN (< 0.5 M) were most consistent with observed cell concentrations, but skill did not differ significantly between low kN and removing the nutrient dependence entirely by setting kN=0.

Incorporating the daily variation in irradiance and associated diel vertical migration of A. fundyense was another important modification from the Gulf of Maine model. Observations have found that diel vertical migration substantially alters the vertical distribution of cells in Salt Pond (Anderson and Stolzenbach, 1985; Crespo et al., 2011). Vertical migration is incorporated through a time-varying vertical advection of cells, with upward velocities during the early morning and downward velocities in the afternoon, as in observations. During daylight hours, maximum cell concentrations typically were 2-3 m below the surface, while at night the maximum concentrations descended to 5-6 m. The environmental cues that are driving the observed vertical migration remain uncertain, but in the model we assume that it depends on solar irradiance. Specifically, cells in the model swim up during the day (I0 > 75 W m2) and down at night, with the top of the ambit at 1/kw and the bottom at 2/kw. The values for kw were based on fits to measured profiles of photosynthetically active radiation (PAR), described below. Typical values for kw were 0.3-0.5 m−1, leading to upper and lower migration limits consistent with the observed cell distributions (Anderson and Stolzenbach, 1985; Crespo et al., 2011). The migration in the model is in phase with the solar cycle rather than beginning ascent before sunrise and descent before sunset (Kamykowski, 1981), a simplification that could be refined. The daily variation in irradiance is also used to calculate cell growth rates, so maximum instantaneous growth rates (μmax) in the Nauset model are greater than growth rates in the Gulf of Maine model that use continuous, daily averaged irradiance. Averaged over the diel cycle, the maximum growth rates in the two models are similar.

The decline of A. fundyense blooms in the Gulf of Maine has been found to best match observations using a Q10 formulation that imposes a strong temperature dependence on cell loss (He et al., 2008). The Q10 formulation increases mortality exponentially with temperature and encompasses many potential loss terms including predation, parasitism, and encystment. A similar Q10 approach can be used to match observations in Nauset in any single year, but no one formulation simulated the decline of the bloom over the multiple years with different temperature histories. Instead of depending on water temperature, the decline tended to occur at a consistent phase of the bloom that could be described by a degree-day calculation (Ralston et al., 2014). The mechanisms behind the bloom decline in Nauset remain unresolved, as it may represent loss due to encystment, parasitism by Amoebophrya (Velo-Suárez et al., 2013), or some other process. In the results presented here, the mortality rate increased from 0 to the maximum value over degree days 500 to 550. The approach is highly empirical, and uses a local, Eulerian calculation of cumulative degree day (at each grid cell) to represent Lagrangian processes associated with advecting cells (Ralston et al., 2014). However, as shown in the results, this empirical approach more effectively represented the decline of blooms over multiple years than any single Q10 formulation. The degree day approach is not causatively linked to specific mechanisms driving the decline of the bloom.

2.2.3 Model skill assessment

To quantify model performance, we compared results with observations of physical variables from moored sensors (water level, temperature, salinity) and of A. fundyense concentrations from weekly surveys. For metrics, we focus on the correlation coefficient (r) and a skill score (SS). The correlation coefficient is the covariance between the model and the observations, while the skill score is the mean square error normalized by the standard deviation of the observations:

SS=1−σo−21N∑i=1N(Xm−Xo)2 (1)

where X is the variable of interest, δ is the standard deviation, N is the number of samples, and subscripts m and o represent the model and observations (Murphy, 1988). The skill score compares the model prediction to the mean of the observations, and is also described as a modeling efficiency (Stow et al., 2009). The skill score can be related to the correlation coefficient and two additional terms:

SS=r2−(r−σmσo)2−(Xm¯−Xo¯σo)2 (2)

where an overbar represents a time average. The second term is the difference in variance between the model and the observations, which vanishes when the slope of the regression line is equal to 1. The last term is the mismatch of the means, and is equal to the intercept of the linear regression. The maximum skill score is 1, and a skill of 0 represents a mean square error equal to the variance of the observations.

Most of the model-data comparisons shown here are from the spring of 2011, which has the most extensive set of physical and biological observations. Similar, less extensive, observations were made in 2009 and 2012. For consistency, model parameters were developed for the 2011 case and applied to all 3 years rather than tuning coefficients to match each year separately.

2.3 Observations

2.3.1 Large-scale weekly and high-resolution tidal-cycle surveys

Approximately weekly during the spring months of 2009, 2011, and 2012, large-scale surveys were made of Nauset estuary, with sampling at about 30 stations throughout the estuary (Fig. 1): 12 surveys from 24 March to 17 June 2009, 13 surveys from 23 March to 16 June 2011, and 11 surveys from 15 February to 8 May 2012 (Crespo et al., 2011; Ralston et al., 2014). Surveys occurred around daytime high tides to maximize navigability of the central marsh. Continuous vertical profiles of salinity and temperature were measured with a conductivity-temperature-depth (CTD, SeaBird Electronics SBE 19plus) sensor at each station, and in 2011 and 2012, instruments also measured chlorophyll fluorescence (470 nm excitation, 685 nm emission, SeaPoint Sensors) and photosynthetic active radiation (PAR, 400-700 nm, scalar detector, Biospherical Instruments). Water samples for A. fundyense cell counts and nutrient concentrations were collected at each station with 2.5 L Niskin bottles. Samples were taken at the surface and, where possible, at 3, 5, 7, and 10 m depths. Details on the sample processing are provided in Crespo et al. (2011) and Ralston et al. (2014).

The weekly surveys provided snapshots of Nauset estuary and captured large-scale spatial variability and seasonal trends. To complement the weekly surveys, a series of higher-resolution, higher-frequency surveys around the peak of the A. fundyense bloom were conducted in Mill Pond and Salt Pond in 2011 and 2012. These surveys had additional sampling stations to increase the horizontal resolution, and were repeated every 30 to 60 minutes through a tidal cycle. Sampling from a small-boat included the CTD and fluorescence profiles as above, with Niskin bottle samples collected for a subset of the stations at regular depths (surface, 2-3 m, 5 m, 7 m and 10 m, depending on water depth). Representative data from a high-resolution tidal-cycle survey in Salt Pond on 9 May 2011 are used below to illustrate the spatial structure of water properties and A. fundyense in the ponds.

The chlorophyll fluorometer with the CTD was calibrated in the laboratory using serial dilutions of A. fundyense in culture from 16 to 0.5 ×106 cell L−1. Regressions to convert from fluorometer voltage to cell concentration were developed based on the lab tests. Water column profiles from the fluorometer were also compared with cell counts from Niskin bottle samples (Fig. 2). Fluorometer values were extracted from the depth of the water sample, averaging over a depth range equal to the length of the bottle. The fluorometer compared well with the bottle samples when concentrations of A. fundyense were greater than about 104 cell L−1; at lower concentrations the chlorophyll signal likely was dominated by other phytoplankton species. The vertical structure in the fluorometer data when A. fundyense concentrations were high corresponded with the vertical structure from the bottle sampling profiles, albeit at much higher resolution. Concentrations derived from fluorometer readings tended to be greater than the cell counts by about a factor of 2, perhaps due to the presence of other phytoplankton, although no whole phytoplankton community counts were conducted to assess this. Cell counts from bottle samples were used for model evaluation, but the fluorometer data were used to characterize the spatial distribution of A. fundyense.

Figure 2.

Figure 2

Comparison of A. fundyense concentrations derived from bottle samples and from chlorophyll fluorescence. Data are from large-scale, weekly surveys in 2011 and 2012 and from high-resolution, tidalcycle surveys in Mill Pond and Salt Pond in 2011. In-situ fluorometer profiles were converted to cell concentrations using a laboratory calibration to cultured A. fundyesne, and extracted at depths corresponding with the sample bottles. The vertical dashed line is the concentration at which the fluorometer saturated on the CTD used for the weekly surveys, and lines with slopes of 1, 0.5, and 0.1 areshown for reference.

2.3.2 Moored instruments

Moored, internally recording instruments were deployed during the blooms to provide continuous records of water properties. The configuration of moored sensors varied among years, but typically water level and near-surface and near-bottom temperature and salinity were measured in each pond (Mill Pond, Salt Pond and Town Cove) and in the central marsh (Hemenway) (Fig. 1). In 2009, water level and temperature sensors (Onset HOBO) were deployed at these four locations. In 2010 and 2011, surface and bottom CTDs (Richard Brancker Research) were deployed in the deep parts of Salt Pond and Mill Pond, while Town Cove and Hemenway had shallower water level and temperature sensors. In 2012, surface and bottom CTDs were deployed in the centers of Town Cove, Salt Pond, and Mill Pond, and a water level and temperature sensor was at Hemenway. Sampling intervals for all the moored sensors were 10 minutes or less. The moored time-series were compared with the profiling CTD data from the weekly surveys to identify periods when fouling compromised moored data quality. Near-surface conductivity measurements were affected by fouling at several locations toward the end of the deployments, and these data were removed from the analysis.

2.3.3 Cyst mapping surveys

To map cyst distributions, sediment samples were collected in the falls of 2008, 2009, and 2011 at 73 stations around the estuary. Details on sample collection and cyst enumeration are given in Crespo et al. (2011). Cyst counts from the top (0-1 cm) sediment layer were mapped to the model grid to provide the bottom boundary condition for the A.fundyense model.

3. RESULTS

3.1 Hydrodynamics and physical conditions

Tidal forcing offshore of Nauset is predominantly semidiurnal with a tidal amplitude of 1 to 2 m. Tidal nonlinearities through the shallow central marsh make tidal water levels in the estuary highly asymmetric, with a brief flood period of elevated velocities (~3-4 h) and a longer, slower ebb (~8-9 h) (Aubrey and Speer, 1985). The hydrodynamic model reproduced the observed tides well, with skill scores typically greater than 0.95 (Fig. 3, Table 2). The nonlinearity of the tide was sensitive to the bottom roughness (z0) and to relatively modest changes in the bathymetry. Water level observations were used to calibrate the bottom roughness in the model, resulting in z0 = 0.024 m. In the absence of a spatially resolved map of bed composition, bottom roughness was assumed to be spatially uniform.

Figure 3.

Figure 3

Observations and model results from spring 2011. (a) Water level in Salt Pond. (b) Near surfaceand near bottom temperature in Salt Pond and (c) Mill Pond, (d) Town Cove, and (e) near-bottom temperature at Hemenway in the central marsh.

Table 2.

Model skill: correlations (r2) and skill scores (SS)

Parameter and location 2009 2011 2012
r2 SS r2 SS r2 SS
Water level
    Town Cove 0.98 0.97 0.96 0.94 0.97 0.96
    Mill Pond 0.97 0.96 0.94 0.93 0.97 0.90
    Salt Pond 0.95 0.94 0.97 0.97 0.97 0.95
    Hemmenway 0.96 0.95 0.97 0.97 0.97 0.95
Water temperature (surface/bottom, where available)
    Town Cove 0.86 0.79 0.78 0.72 0.95/0.91 0.92/0.90
    Mill Pond 0.85 0.62 0.00/0.25 0.00/−6.8 0.95/0.60 0.77/−2.8
    Salt Pond 0.83 0.35 0.01/0.80 −0.01/0.59 0.85/0.81 0.71/0.08
    Hemmenway 0.48 0.46 0.75 0.55 0.62 0.52
A. fundysense concentration (log transformed): degree-day mortality
    Town Cove 0.60 −0.76 0.87 0.62 0.46 −0.51
    Mill Pond 0.71 0.54 0.62 0.38 0.68 0.60
    Salt Pond 0.40 0.35 0.05 −0.13 0.55 0.40
    Hemmenway 0.29 −4.9 0.48 0.40 0.68 −1.2
A. fundysense concentration (log transformed): Q10 mortality
    Town Cove 0.35 0.28 0.52 0.39 0.20 −0.07
    Mill Pond 0.56 0.31 0.33 −0.20 0.08 −0.71
    Salt Pond 0.47 0.41 0.86 0.60 0.10 −0.42
    Hemmenway 0.15 −0.97 0.18 −0.83 0.79 0.57

Seasonal warming of water temperatures depended on several factors including the surface heat flux, tidal exchange, and bathymetry. Surface layers of the three ponds warmed faster than the central marsh due to the more limited influx of cooler water from the coastal ocean and the more limited exposure to wind forcing (Ralston et al., 2014). The deeper layers of the ponds remained cooler than the surface, as stratification limited the influence of surface heating at depth. Early in the spring, the salinity anomaly due to groundwater influx was the primary source of stratification, but as the surface layer warmed, temperature became an equal or greater contributor to stratification. The spatial and temporal distribution of seasonal warming in the model (Fig. 4) compared will with observations (Fig. 3 in Ralston et al., 2014). The ponds warmed faster than the central marsh, with temperature increases in Mill Pond and Town Cove leading Salt Pond.

Figure 4.

Figure 4

Maps of surface temperature from model results for spring 2011. Times shown are extracted to correspond with marsh-wide CTD surveys that occurred around high water, as in Fig. 3 of Ralston et al.(2014).

Variation in tidal energy affected stratification in the ponds at spring-neap time scales (Fig. 3). During neap tides, tidal velocities and exchange with the coastal ocean decreased, so the surface layers of the ponds warmed and stratification increased (e.g., Apr 24-29). During spring tides, stratification decreased. Wind events also affected stratification in the ponds, both through mixing of the surface layer and by driving coastal set-up that increased influx of cooler coastal ocean water. Diurnal solar heating of the surface layer and convective cooling at night altered water temperature at daily time scales. The model reproduced much of the daily to fortnightly variations in stratification. Stratification in Salt Pond (Fig. 3b) was well represented, including the warming of the surface layer during neap tides and the warming of the lower layer during spring tides. In Mill Pond the temperature variation in the surface layer was reasonably well represented, but the lower layer in the model became too warm, indicative of excessive mixing down from the surface. The extra vertical mixing was insensitive to the turbulence closure, and was apparently due to high rates of numerical mixing associated with steep topography, sigma coordinate systems, and lower order advection schemes (Burchard and Rennau, 2008). As shown later, numerical mixing also smoothed vertical profiles of A. fundyense. The skill metrics showed reasonable overall agreement with observed temperatures. For surface waters, most correlations (r2) were 0.8 or greater and skill scores were between 0.4 and 0.9 (Table 2).

Large differences in water properties were observed between the ponds and central marsh, as seen in transects through a tidal cycle in Salt Pond (Fig. 5). During the ebb (transects 1 and 2), the pond was strongly stratified and warm water ebbing out into the marsh formed a sharp temperature gradient with cooler water that entered from the coastal ocean during the previous flood (Fig. 5a,b). During the flood (transects 3 and 4), ocean water moved back into the marsh and into Salt Pond (Fig. 5c,d). The ocean water was cooler and saltier, and thus denser than the surface layer of the pond, leading to convective mixing downward near the pond inlet. The model captured the critical elements of this tidal pattern, including the maximum in stratification during the ebb, the front between pond and ocean water, and the downward mixing at the inlet during flood tides. More generally, both the survey and time series data indicate that temperatures in the model corresponded with observations at seasonal and tidal time scales, an important component for A. fundyense growth.

Figure 5.

Figure 5

Temperature sections from Salt Pond observations (left column) and model results (right) for 9May 2011. Sections are from the pond through the channel into the central marsh (A-A’, larger panels),and across the pond (B-B’, inset panels); section locations are shown on the Salt Pond inset of Fig. 1. Time of day is shown in the lower left of each panel, and tidal phases for the 4 sections are shown in the lower left panel with the water levels in Salt Pond (gray) and at Nauset Inlet (black). Salinity contours(0.2 psu interval) are overlaid on the model results.

3.2 A. fundyense population

The A. fundyense blooms in Nauset were highly localized within the ponds, with cell concentrations that were substantially greater than in the central marsh (Fig. 6), indicating little exchange of cells between the salt ponds (Crespo et al. 2011). The modeled spatial distribution and seasonal evolution of the bloom were consistent with observations in 2011 (Fig. 2 in Ralston et al., 2014) and other years (e.g., 2009 bloom in Fig. 7 of Crespo et al. 2011). As with temperature, the cell concentration maps were selected to correspond with observations that occurred near high-tide, and thus coastal ocean water in the central marsh had low cell concentrations.

Figure 6.

Figure 6

Maps of mean A. fundyense concentration from model results for spring 2011. Times shown are extracted to correspond with marsh-wide CTD surveys that occurred around high water, as in Fig. 2 of Ralston et al. (2014).

Figure 7.

Figure 7

Depth-averaged A. fundyense cell concentrations stations from weekly surveys and model results in spring 2011: (a) Town Cove, (b) Mill Pond, (c) Salt Pond, and (d) Hemmenway. Station locations are shown in Fig. 1. Vertical lines indicate the concentration range observed in samples collected at multiple depths

A comparison of the model results with mean cell concentrations measured in the weekly surveys finds that the model captured the seasonal trends at multiple stations, especially when variation with depth was considered (Fig. 7). However, the model did not reproduce the bloom phasing among the three ponds, a commonly observed pattern whereby cells accumulate first in Mill Pond, then Town Cove, and lastly Salt Pond. This may be due to discrepancies in water temperature, as the surface waters of Mill Pond and Upper Mill Pond were slightly cooler than observed, and thus growth rates were reduced. The daily variations in cell concentration in the central marsh were greater than at stations in the ponds, as tidal advection brought low concentrations from the coastal ocean during floods and higher concentrations from the ponds during ebbs. The weekly survey data occurred around high tide and therefore corresponded with the lower bound of the concentration envelope for the central marsh. Skill metrics for the cell time series varied by location and year (Table 2). Correlations (r2) were often greater than 0.5, and while the skill scores were lower, the A. fundyense model had positive skill scores in many cases, particularly in the ponds.

Tidal cycle transects in Salt Pond provide greater detail on the spatial and temporal variability in A. fundyense cell distributions around the peak of the bloom (Fig. 8). During the ebb, cells were concentrated in a relatively thin layer near the pycnocline, with much lower concentrations ebbing out near the surface (Fig. 8a). Late in the ebb as the water level dropped and surface layer thinned, export of cells from the pond increased, but concentrations in the marsh remained lower than in the pond (Fig. 8b). During the following flood, the incoming ocean water was cooler, saltier, and had lower concentrations (Fig. 8c,d). As it entered the pond and mixed down, the layer of cells near the pycnocline was mixed and displaced, reducing cell concentrations near the inlet. Late in the flood, maximum cell concentrations were deeper in the water column and nearer the periphery of the pond (Fig. 8d). The model reproduced the key elements of the observed A. fundyense distributions, including the local maxima at the pycnocline, the low concentrations exported late ebb, and the mixing and displacement of cells away from the inlet during the flood. The model had some discrepancies in phasing, for example, the flooding ocean water arrived slightly earlier in the observations than in the model, but overall the comparison with the high-resolution survey data suggests that the model is plausibly representing the hydrodynamics and A. fundyense in the pond. Subsequent sections use the model to analyze factors controlling bloom development, including the role of vertical migration for cell retention in the ponds, the dependence of the bloom on cyst distribution and germination rate, and the parameterization of bloom termination.

Figure 8.

Figure 8

Observed chlorophyll-a fluorescence (left column) and modeled A. fundyense concentration (right) from Salt Pond for May 9, 2011. Sections are the same as in Fig. 5. Time of day is shown in thelower left of each panel, and tidal phases for the 4 sections are shown in the lower left panel with the water levels in Salt Pond (gray) and at Nauset Inlet (black).

3.3 Vertical migration and residence time

Previous observational studies indicated that the vertical migration pattern of A. fundyense in Salt Pond was essential to the high rates of retention in the ponds and facilitated the intense blooms there (Anderson and Stolzenbach, 1985). Most of the observations in this study occurred during daylight hours, making it difficult to diagnose the diel variability in the cell population. Our tidal cycle survey in Salt Pond (Figs. 5 and 8) was also during the day, but additional profiles were collected in the center of the pond at midnight and 8 am the following morning (Fig. 9). The fluorometer profiles were consistent with the conceptual model that during the day cells are concentrated near the pycnocline, about 3 m below the surface and roughly corresponding with the 1/kw level of irradiance, and at night the center of mass is deeper, about 5 m below the surface. The following morning the cells appeared to be moving up toward the surface again. The movement of cells downward at 17:21 was likely due to a combination of vertical migration late in the day and convective mixing due to the jet of denser water flooding from the central marsh (Fig. 5d). The observed profile at 8:04 the following morning was shallower than in the model, which could be due to the phasing of the vertical migration leading the irradiance. Overall however, the observations were consistent with previous studies and with the simple parameterization of vertical migration in the model.

Figure 9.

Figure 9

Observed (upper panels) and modeled (lower panels) profiles of temperature and chlorophyll-a (observed) or A. fundyense concentration (model) in the center of Salt Pond on May 9-10, 2011. Time of day is in the upper right of each panel. Irradiance at the time of each profile is noted at the bottom of theupper panels, and tidal phase is shown with the red circle in the lower panels.

To quantify the role of vertical migration in retention of cells in the pond, we ran a series of model experiments. Cell growth and mortality were removed from the model to assess directly the dependence of advective losses on swimming behavior, tidal forcing, and stratification. Representative spring and neap tide cases from Salt Pond (April 17 and April 24, 2011, respectively) were selected as starting points for full model simulations, with growth and mortality turned off. The subsequent evolution of the bloom was due strictly to hydrodynamic processes and various swimming behaviors for A. fundyense, and cell concentrations decreased through cell export from the pond to the central marsh (Fig. 10a). Residence times were calculated by fitting an exponential decay curve to the total cell population in the pond, resulting in an e-folding time scale; the fitted curves are shown with the concentration time series for the spring tide cases (Fig. 10a).

Figure 10.

Figure 10

Residence time calculations from model results for different swimming and forcing cases. (a) Concentration of A. fundyense remaining in Salt Pond for spring tide cases: diel vertical migration up to 1/kw depth (swim), no vertical migration (don't swim), diel migration to the surface (swim to surface), anddiel vertical migration to 1/kw with barotropic physics and barotropic physics and no wind forcing. Exponential fits are shown for each case. (b) Residence time calculated from exponential fits for springtide (shown in (a)) and neap tide cases. For reference, the residence time for tidal exchange with a well-mixed pond (Vpond/Qtide) is shown for spring and neap tide cases.

Several swimming behaviors were examined. In addition to the base case where cells migrate between 1/kw and 2/kw (“swim”), cases were tested with no vertical migration (“don't swim”) and with vertical migration to the surface rather than to 1/kw (“swim to surface”) (Fig. 10b). All cases started with the same initial distribution of cells. In addition to the behavior of A. fundyense, stratification in the pond alters vertical mixing and thus cell retention, so a case was run with the default vertical migration strategy but barotropic hydrodynamics. The barotropic case assumes uniform water density, removing effects of stratification on turbulent mixing and increasing the vertical mixing of cells by wind and tidal forcing.

The numerical experiments showed that residence times for A. fundyense depended substantially on the tidal forcing, the swimming behavior, and the model hydrodynamics (Fig. 10b). For reference, the volumetric residence time, equal to the volume of the pond divided by the tidal volume exchange, is shown with dashed lines for spring and neap tides. The A. fundyense residence times were all greater than the volumetric exchange rate because cells were not completely, continuously mixed in the pond, as that analytical model assumes. The cases with the default swimming behavior still had longer residence times than the alternatives of not swimming or swimming to the surface. Residence times were longer during neap tides than during spring tides, particularly with vertical migration. During spring tides, not only did the volumetric exchange increase, but vertical mixing of cells into the surface layer also increased, accelerating export of A. fundyense from the pond.

Decreased residence times for the barotropic cases versus the “swim” and “don't swim” cases illustrated the importance of vertical mixing for cell retention. Stratification reduced vertical mixing at the pycnocline, which was typically near the upper limit of vertical migration. Turbulent mixing, particularly due to wind stress, was greater in the barotropic case, transporting cells into the surface layer for subsequent export. Wind mixing was the major source of the vertical flux, as a barotropic case without wind forcing had longer residence times than the full physics case (Fig. 10). The barotropic cases also did not include the energetic vertical mixing of cells downward by the dense jet of oceanic water entering the pond during spring flood tides (Fig. 8d), further reducing the tidal exchange of cells. Effective growth rates for these temperature and light conditions were 0.2-0.3 d−1 (Watras et al., 1982), equal to doubling times of 2.3 to 3.5 d. The doubling times were similar to or greater than the volumetric residence times, so cell accumulation is enhanced by retention due to active avoidance of the surface layer and the reduction in mixing by stratification.

3.4 Bloom initiation: Sensitivity to cyst distribution and germination rate

In the Gulf of Maine, the cyst bed distribution has leading order effects on the A. fundyense bloom (Anderson et al., 2014). In Nauset, cyst distributions varied from year-to-year, although cyst abundance was greater in the ponds than in the central marsh (Fig. 11, Table 3). Mean cyst concentrations in and near the ponds were a factor of 10 to 30 greater than mean concentrations in the central marsh, and mean values in the ponds varied by a factor of about 5 over the 3 years of sampling (Table 3).

Figure 11.

Figure 11

Maps of observed cyst distributions (0-1 cm sediment depth) in the falls of 2008, 2009, and 2011. Sample locations are shown with circles, and interpolated maps were used as bottom boundary conditions for the Alexandrium model. The coastline shown here is the mean sea level contour.

Table 3.

Observed A. fundysense cyst abundance in fall surveys (0-1 cm sediment depth)

Year Max (cysts cm−3) In and near ponds In central marsh
Samples Mean (cysts cm−3) Samples Mean (cysts cm−3)
2008 2,909 34 480 40 42
2009 4,965 34 762 38 47
2011 18,150 34 2,786 31 76

Ideally, the model would be forced with cyst data from surveys conducted the fall prior to each spring bloom, but surveys were only conducted in 2008, 2009, and 2011. The model was run with prior fall observations when available (2009 and 2012 blooms) and the 2011 bloom was initialized with the cyst distribution from the prior survey closest in time, from fall 2009. Additionally, simulated cyst fields were used to initialize the 2011 case to assess the sensitivity of the bloom to the cyst distribution, as described below. A formulation for temperature-dependent germination rates from the Gulf of Maine (Stock et al., 2005) was applied although the transferability of those rates to Nauset remains uncertain. The Nauset simulations did not apply the endogenous clock germination rhythm from the Gulf of Maine (Anderson and Keafer, 1987), and all cysts in the surface layer (0-1 cm) were assumed to be viable.

To test the sensitivity of the model to bloom inoculation by cysts, we altered both the germination rate and the spatial distribution. Observed cyst abundance near the sediment-water interface (0-1 cm depth in bed) had substantial heterogeneity around the estuary. To test sensitivity to this distribution, a case was run with a uniform cyst distribution equal to the mean of the mapped observations (Fig. 12). The uniform cyst case resulted in slightly higher cell concentrations early in the bloom for parts of the estuary (e.g. Town Cove and Salt Pond) and slightly lower in others (e.g., Mill Pond), but the differences were small compared with the concentrations at the peak of the bloom, and model skills were not significantly different from the base case. The result confirmed that the bloom is relatively insensitive to the cyst distribution and that vegetative cell division is the primary mechanism underlying bloom development, provided there is sufficient germination to initiate the bloom.

Figure 12.

Figure 12

Sensitivity to cyst distribution, comparing model results using the observed cyst distribution with a spatially uniform cyst distribution equal to the average cyst concentration and with cyst densities reduced by factors of 10 and 100. Cell concentrations from spring 2011 are shown for (a) Town Cove, (b) Mill Pond, (c), Salt Pond, and (d) Hemenway. Horizontal gray bar represents approximate range of concentrations in the ponds at the times that weekly toxicity sampling exceeded the regulatory threshold of 80 μg toxin per 100 g mussel tissue (MA Division of Marine Fisheries).

The sensitivity of the bloom to the number of cysts was examined with simulations using lower cyst abundance, reduced by factors of 10 and 100 from the observed distributions (Fig. 12). A factor of 10 reduction is similar to the observed range of spatial variability, while a factor of 100 reduction represents a significant decrease in cyst abundance, perhaps as might be associated with remedial action to reduce cyst concentrations. Blooms from the depleted cyst beds followed similar temporal development as the base case, but cell concentrations were reduced by factors similar to the 10 or 100 times reductions in cyst abundance. Cases were also run assessing the connectivity among the ponds by removing all cysts in the northern or southern part of the estuary. The ponds are isolated from each other by the hydrodynamics of the marsh, so cases with cysts only in the southern half had essentially no bloom in Salt Pond and cases with cysts only in the northern half did not have blooms in Mill Pond or Town Cove. Together, the cyst sensitivity results suggest that modifications to the cyst distribution in and near the ponds might alter the size or local geography of the blooms.

Along with the cyst distribution, bloom initiation depended on the rate at which cysts germinate into vegetative cells. In the base case, germination rates varied with temperature from about 1.5 % d−1 at less than 5°C to about 9 % d−1 at greater than 11°C (Anderson, Stock, et al., 2005). In 2012, January and February were mild and water temperatures were a few degrees warmer than in 2009 or 2011(Ralston et al., 2014). Also in 2012, the bloom occurred about 1 month earlier than in the other years. To test whether higher temperatures and increased germination rates may have accelerated bloom development, a case was run with the germination rate always equal to its maximum value, independent of temperature (Fig. 13). Faster germination early in the year produced greater cell concentrations in the first month of the bloom, but the discrepancy with the base case was only at relatively low concentrations (< 300 cells L−1) and was largely eliminated by late March. Early in the bloom, the total number of germinated cysts was large compared with the number of vegetative cells, but as the water warmed and growth rates increased, cell concentrations were several orders of magnitude greater than the number of germinated cysts.

Figure 13.

Figure 13

Sensitivity of model results to cyst emergence rate, comparing the base case of temperature-dependent cyst emergence with a case with cyst emergence at the maximum rate from lab data, independent of temperature. Cell concentrations in the water column and the total number of germinatedcysts are shown for each case; both are averaged over the entire Nauset marsh to remove effects of spatial variability in cyst concentration.

3.5 Bloom termination: Interannual variability in loss rates

Modeling the termination of A. fundyense blooms in Gulf of Maine has been challenging. A constant mortality rate does not account for variability associated with increased temperature and development of a thermocline toward the end of blooms, and it was difficult to separate increased mortality from decreased nutrient availability as causes of bloom decline (Stock et al., 2005). Subsequently, a temperature-dependent loss term (Q10) was introduced that better simulated bloom termination (He et al., 2008). This function represented a variety of mechanisms including predation, parasitism, mortality and encystment.

A similar Q10 formulation was initially adopted for the Nauset model. The approach gave reasonable results for both the 2009 and 2011 blooms using a single set of parameters (Fig. 14). However, as discussed above, the winter of 2012 was anomalously warm and the Nauset bloom occurred about 1 month earlier than in the previous years (Ralston et al., 2014). January and February were several degrees warmer in 2012 than usual, spurring rapid growth earlier in the year, but temperatures in March and April were similar to 2009 and 2011. The bloom duration and maximum concentrations were similar in the three years, but the 2012 bloom was shifted about 1 month earlier. As a result, water temperatures during the termination phase of the bloom were several degrees lower in 2012 than in the other years, and the Q10 formulation significantly underpredicted loss rates (Fig. 14c).

Figure 14.

Figure 14

A. fundyense concentration from observations and model results from Salt Pond for (a) 2009, (b) 2011, and (c) 2012. Observations are from weekly surveys each year, and vertical bars represent therange of concentrations from samples at multiple depths. Model results are shown using two different mortality formulations: based on temperature (Q10) and on growing degree days.

Alternatively, a simple degree-day model for cumulative, temperature-dependent growth can collapse temporal variability among the 3 years (Ralston et al., 2014). When loss rates were made a function of degree day, the termination phase of the bloom was reasonably represented by the same formulation in all 3 years (Fig. 14). Due to the model architecture, it was necessary to calculate degree days locally within each grid cell, an approach that was not optimal given the diel ambit of individual A. fundyense cells through temperature stratification. The degree day formulation is therefore entirely empirical and specific to the Nauset observations and model, but it does provide greater skill than a single Q10 approach over multiple years (Table 2).

4. DISCUSSION

4.1 Factors contributing to bloom development and structure

The coupled hydrodynamic-biological model compared well with the observed physical conditions and A. fundyense concentrations over multiple years, and the model can be used to assess the factors contributing to the development of the recurrent HAB in Nauset estuary. As noted from observational data, the effect of temperature on growth rate was a key factor controlling bloom development (Watras et al., 1982; Crespo et al., 2011; Ralston et al., 2014). Temperature in the ponds depended both on the seasonal increase in surface heat flux and on spring-neap variability in tidal exchange, with warming of the surface layer during neap tides and cooling from coastal ocean water during spring tides. The spring-neap variability in exchange affected the residence time of cells in the ponds, but for cell retention to be substantially longer than the cell doubling time required both active vertical migration and stratification that inhibits vertical mixing of cells into the surface layer. The model results were consistent with a previous observational study in Salt Pond using a dye tracer (Anderson and Stolzenbach, 1985), but emphasized that stratification is important in addition to migration . Similarly, a HAB in Cork Harbor (Ireland) has been found to develop when spring tides are relatively weak around the summer solstice, reducing the loss of cells to tidal mixing and advection (Raine, 2014). In Thau Lagoon (France), wind is the dominant energy source and control on retention, and Alexandrium blooms are restricted to periods of weak wind mixing (Laanaia et al., 2013).

Unlike A. fundyense blooms in the Gulf of Maine, interannual variability in cyst distribution does not appear to be the dominant control on the intensity or geography of the Nauset blooms. Instead the Nauset system appears to be more similar to Cork Harbor, a coastal embayment the intensities of Alexandrium blooms have been shown to be independent of pre-bloom cyst abundance (Cosgrove et al., 2014). In Nauset, comparisons between cyst abundance and peak toxicities or maximum cell concentrations were inconclusive due to the limited number of realizations (3 years). Maximum toxicity recorded at shellfish monitoring stations in Salt Pond and near Mill Pond were correlated with cyst abundances in the ponds the previous fall (n = 6, r2 = 0.80, p=0.02), but the correlation was dominated the elevated toxicity and cyst abundance recorded in Mill Pond for the 2012 bloom; without that data point the relationship fell apart (n=5, r2 = 0.01, p=0.86). The relationship between cyst abundance and bloom intensity was much stronger in the model results. Although cell concentrations varied little between cases with a realistic or a spatially uniform cyst distribution (Fig. 12), peak cell concentrations decreased nearly proportionately to 10- or 100-fold decreases in cyst abundance. An important caveat to the model is that bloom termination depends heavily on the degree-day heuristic. This relationship between heat accumulation and increasing cell loss is derived from observations, without an understanding of the underlying mechanism. Other mortality formulations such as a concentration-based approach could allow the bloom to continue to develop to similar maximum concentrations. .

Modifications to the biological model suggest differences in the A. fundyense population from the Gulf of Maine. An endogenous germination rhythm for cysts in Nauset, if it is present, is shifted earlier in the year. Population losses at the end of the bloom are highly nonlinear, but do not appear to be strictly a function of temperature. Instead, given that inferred losses are extremely low during early phases of the bloom and increase rapidly at similar cell concentrations across multiple years, termination may be concentration-dependent. The apparent skill of the degree day approach then relies on a linear relationship between temperature and growth rate until exceedance of a cell concentration threshold that triggers rapid losses to parasitism and/or sexual encystment. More observations are needed to isolate mechanisms of termination to provide the basis for a more mechanistic model.

4.2 Model attributes and limitations

The highly resolved 3-d hydrodynamic and A. fundyense model represents strong spatial gradients in water properties and cell concentrations, both spatially among the ponds and central marsh and vertically due to stratification. A coarser or more idealized model, such as depth-averaged or box model approach, would be unlikely to capture the dominant role of bathymetric features in determining the structure of the bloom. The model demonstrated robust skill for hydrodynamics and water properties, and provided order-of-magnitude predictive capability for cell concentrations across multiple years.

While the comparisons with observations were favorable, key uncertainties remain in the model formulation. The high loss rates at bloom termination may be due to multiple processes including encystment, parasitism, and grazing. The controlling factors are not obvious, but are likely to depend on the life history of the organism (sexual induction and encystment), infection by the parasite Amoebophrya (Velo-Suárez et al., 2013), grazing pressure, and physical properties of the environment. A Lagrangian modeling approach might better represent the declining phase of the bloom by tracking individual particles. Similarly, the approach in the current model to vertical migration would benefit from a mechanistic approach based on external cues (light, nutrients, density, velocity shear) and organism state (life stage, nutrient status; Ralston et al., 2007) for a more robust examination of how vertical migration influences residence time and bloom development.

The dependence of the Nauset bloom on nutrients remains unresolved. Observationally, growth rates based on cell concentrations from weekly surveys were not correlated with measured nutrient concentrations (Ralston et al., 2014). To force the model, samples from Niskin bottles at discrete stations and depths did not resolve well the strong spatial gradients in nutrient concentrations. High resolution surveys of nitrate concentration along the shoreline of Salt Pond detected spatial gradients in nitrate corresponding with groundwater seeps (J. Colman, pers. comm.). Late in the blooms, ammonium concentrations below the pycnocline in Mill and Salt Ponds were high enough to be toxic to organisms through inhibition of nitrate uptake (Dugdale et al., 2007). However, the swimming ability of A. fundyense may allow it to actively avoid inhibitory concentrations in the lower water column. Similarly, we could not assess whether limitations of other nutrients or altered nutrient stoichiometry associated with anthropogenic loading (Glibert et al., 2013) may affect bloom development.

The model presented here is consistent with the available observations from Nauset, but that does not exclude alternative model formulations. For example, similar skills can be obtained removing nutrients as a growth factor by setting kN = 0. Similarly, the mortality function assumes loss processes occur in the water column (grazing, encystment, parasitism), but benthic grazing may also be important. As in many ecosystem models, the justification for adding complexity is limited by the availability of field or lab observations to constrain the model. The results here are perhaps the most resolved and validated simulation of an HAB to date, but observations should continue to inform both its conceptual and numerical development.

Acknowledgements

This work was supported by the National Science Foundation (OCE-0430724, OCE-0911031, and OCE-1314642) and National Institutes of Health (NIEHS-1P50-ES021923-01) through the Woods Hole Centerfor Oceans and Human Health, and by National Park Service (NPS) Cooperative Agreement H238015504. Special thanks to B. Keafer, K. Norton, and other members of the Anderson lab for datacollection. We are grateful to J. Colman, V. Cross, A. Sorenson, and G. Voulgaris for providing data, andto G. Cowles and R. He for assistance in model development.

References

  1. Anderson DM. Bloom dynamics of toxic Alexandrium species in the northeastern US. Limnology and Oceanography. 1997;42:1009–1022. [Google Scholar]
  2. Anderson DM, Cembella AD, Hallegraeff GM. Progress in understanding harmful algal blooms: Paradigm shifts and new technologies for research, monitoring, and management. Annual Review of Marine Science. 2012;4:143–176. doi: 10.1146/annurev-marine-120308-081121. doi:10.1146/annurev-marine-120308-081121. [DOI] [PMC free article] [PubMed] [Google Scholar]
  3. Anderson DM, Keafer BA. Dinoflagellate cyst dynamics in coastal and estuarine waters. Toxic Dinoflagellates: Proc. 3rd Intl. Conf. Elsevier. 1985:219–224. [Google Scholar]
  4. Anderson DM, Keafer BA. An endogenous annual clock in the toxic marine dinoflagellate Gonyaulax tamarensis. Nature. 1987;325:616–617. doi: 10.1038/325616a0. doi:10.1038/325616a0. [DOI] [PubMed] [Google Scholar]
  5. Anderson DM, Keafer BA, Geyer WR, Signell RP, Loder TC. Toxic Alexandrium blooms in the Western Gulf of Maine: The plume advection hypothesis revisited. Limnology and Oceanography. 2005;50:328–345. doi:10.2307/3597905. [Google Scholar]
  6. Anderson DM, Keafer BA, Kleindinst JL, McGillicuddy DJ, Jr., Martin JL, Norton K, Pilskaln CH, Smith JL, Sherwood CR, Butman B. Alexandrium fundyense cysts in the Gulf of Maine: Long-term time series of abundance and distribution, and linkages to past and future blooms. Deep Sea Research Part II: Topical Studies in Oceanography, Harmful Algae in the Gulf of Maine: Oceanography, Population Dynamics, and Toxin Transfer in the Food Web. 2014;103:6–26. doi: 10.1016/j.dsr2.2013.10.002. doi:10.1016/j.dsr2.2013.10.002. [DOI] [PMC free article] [PubMed] [Google Scholar]
  7. Anderson DM, Rengefors K. Community assembly and seasonal succession of marine dinoflagellates in a temperate estuary: The importance of life cycle events. Limnology and Oceanography. 2006;51:860–873. [Google Scholar]
  8. Anderson DM, Stock CA, Keafer BA, Bronzino Nelson A, Thompson B, McGillicuddy DJ, Jr., Keller M, Matrai PA, Martin J. Alexandrium fundyense cyst dynamics in the Gulf of Maine. Deep Sea Research Part II: Topical Studies in Oceanography, The Ecology and Oceanography of Toxic Alexandrium fundyense Blooms in the Gulf of Maine The Ecology and Oceanography of Toxic Alexandrium fundyense Blooms in the Gulf of Maine. 2005;52:2522–2542. doi:10.1016/j.dsr2.2005.06.014. [Google Scholar]
  9. Anderson DM, Stolzenbach KD. Selective retention of two dinoflagellates in a well-mixed estuarine embayment: The importance of diel vertical migration and surface avoidance. Marine Ecology Progress Series. 1985;25:39–50. [Google Scholar]
  10. Aubrey DG, Speer PE. A study of non-linear tidal propagation in shallow inlet/estuarine systems Part I: Observations. Estuarine, Coastal and Shelf Science. 1985;21:185–205. doi:10.1016/0272-7714(85)90096-4. [Google Scholar]
  11. Aubrey DG, Voulgaris G, Spencer WD, O'Malley SP. Tidal Circulation and Flushing Characteristics of the Nauset Marsh System (Technical No. WHOI-97-11) Woods Hole Oce; Woods Hole, MA: 1997. p. 112. [Google Scholar]
  12. Borkman DG, Smayda TJ, Schwarz EN, Flewelling LJ, Tomas CR. Recurrent vernal presence of the toxic Alexandrium tamarense/Alexandrium fundyense (Dinoflagellata) species complex in Narragansett Bay, USA. Harmful Algae. 2014;32:73–80. doi:10.1016/j.hal.2013.12.005. [Google Scholar]
  13. Brock JC, Wright CW, Patterson M, Nayegandhi A, Travers LJ. EAARL Topography— Cape Cod National Seashore (Open-File Report No. 2007-1375) U.S. Geological Survey; St. Petersburg, FL.: 2007. [Google Scholar]
  14. Burchard H, Baumert H. On the performance of a mixed-layer model based on the κ-ε turbulence closure. Journal of Geophysical Research: Oceans. 1995;100:8523–8540. doi:10.1029/94JC03229. [Google Scholar]
  15. Burchard H, Rennau H. Comparative quantification of physically and numerically induced mixing in ocean models. Ocean Modelling. 2008;20:293–311. doi:10.1016/j.ocemod.2007.10.003. [Google Scholar]
  16. Chen C, Liu H, Beardsley RC. An unstructured grid, finite-volume, three-dimensional, primitive equations ocean model: Application to coastal ocean and estuaries. Journal of Atmospheric and Oceanic Technology. 2003;20:159–186. [Google Scholar]
  17. Colman JA, Masterson JP. Transient simulations of nitrogen load for a coastal aquifer and embayment, Cape Cod, MA. Environmental Science & Technology. 2008;42:207–213. doi: 10.1021/es070638b. [DOI] [PubMed] [Google Scholar]
  18. Cosgrove S, Rathaille AN, Raine R. The influence of bloom intensity on the encystment rate and persistence of Alexandrium minutum in Cork Harbor, Ireland. Harmful Algae. 2014;31:114–124. doi: 10.1016/j.hal.2013.10.015. doi:10.1016/j.hal.2013.10.015. [DOI] [PubMed] [Google Scholar]
  19. Crespo BG, Keafer BA, Ralston DK, Lind H, Farber D, Anderson DM. Dynamics of Alexandrium fundyense blooms and shellfish toxicity in the Nauset Marsh System of Cape Cod (Massachusetts, USA). Harmful Algae. 2011;12:26–38. doi: 10.1016/j.hal.2011.08.009. doi:10.1016/j.hal.2011.08.009. [DOI] [PMC free article] [PubMed] [Google Scholar]
  20. Cross VA, Bratton JF, Crusius J, Colman JA, McCobb T. Submarine Hydrogeological Data from Cape Cod National Seashore (Open-File Report No. 2006-1169) U.S. Geological Survey; Woods Hole, MA.: 2006. [Google Scholar]
  21. Dugdale RC, Wilkerson FP, Hogue VE, Marchi A. The role of ammonium and nitrate in spring bloom development in San Francisco Bay. Estuarine, Coastal and Shelf Science. 2007;73:17–29. doi:10.1016/j.ecss.2006.12.008. [Google Scholar]
  22. Fairall CW, Bradley EF, Hare JE, Grachev AA, Edson JB. Bulk Parameterization of Air– Sea Fluxes: Updates and Verification for the COARE Algorithm. Journal of Climate. 2003;16:571–591. [Google Scholar]
  23. Fairall CW, Bradley EF, Rogers DP, Edson JB, Young GS. Bulk parameterization of airsea fluxes for Tropical Ocean-Global Atmosphere Coupled-Ocean Atmosphere Response Experiment. Journal of Geophysical Research. 1996;101:3747–3764. [Google Scholar]
  24. Fauchot J, Saucier FJ, Levasseur M, Roy S, Zakardjian B. Wind-driven river plume dynamics and toxic Alexandrium tamarense blooms in the St. Lawrence estuary (Canada): A modeling study. Harmful Algae. 2008;7:214–227. doi:10.1016/j.hal.2007.08.002. [Google Scholar]
  25. Franks PJS, Anderson DM. Alongshore transport of a toxic phytoplankton bloom in a buoyancy current: Alexandrium tamarense in the Gulf of Maine. Marine Biology. 1992;112:153–164. doi:10.1007/BF00349739. [Google Scholar]
  26. Giblin A, Gaines A. Nitrogen inputs to a marine embayment: the importance of groundwater. Biogeochemistry. 1990;10:309–328. doi:10.1007/BF00003150. [Google Scholar]
  27. Glibert PM, Kana TM, Brown K. From limitation to excess: the consequences of substrate excess and stoichiometry for phytoplankton physiology, trophodynamics and biogeochemistry, and the implications for modeling. Journal of Marine Systems, Advances in Marine Ecosystem Modelling Research III. 2013;125:14–28. doi:10.1016/j.jmarsys.2012.10.004. [Google Scholar]
  28. Hallegraeff GM. A review of harmful algal blooms and their apparent global increase. Phycologia. 1993;32:79–99. [Google Scholar]
  29. Hattenrath TK, Anderson DM, Gobler CJ. The influence of anthropogenic nitrogen loading and meteorological conditions on the dynamics and toxicity of Alexandrium fundyense blooms in a New York (USA) estuary. Harmful Algae. 2010;9:402–412. doi:10.1016/j.hal.2010.02.003. [Google Scholar]
  30. He R, McGillicuddy DJ, Anderson DM, Keafer BA. Historic 2005 toxic bloom of Alexandrium fundyense in the western Gulf of Maine: 2. Coupled biophysical numerical modeling. Journal of Geophysical Research. 2008;113:C07040. doi: 10.1029/2007JC004602. doi:10.1029/2007JC004602. [DOI] [PMC free article] [PubMed] [Google Scholar]
  31. Hoagland P, Scatasta S. Granéli PDE, Turner PDJT, editors. The Economic Effects of Harmful Algal Blooms. Ecology of Harmful Algae, Ecological Studies. Springer Berlin Heidelberg. 2006:391–402. [Google Scholar]
  32. Kamykowski D. Dinoflagellate growth rate in water columns of varying turbidity as a function of migration phase with daylight. Journal of Plankton Research. 1981;3:357–367. doi:10.1093/plankt/3.3.357. [Google Scholar]
  33. Laanaia N, Vaquer A, Fiandrino A, Genovesi B, Pastoureaud A, Cecchi P, Collos Y. Wind and temperature controls on Alexandrium blooms (2000–2007) in Thau lagoon (Western Mediterranean). Harmful Algae. 2013;28:31–36. doi:10.1016/j.hal.2013.05.016. [Google Scholar]
  34. Li Y, He R, McGillicuddy DJ, Jr., Anderson DM, Keafer BA. Investigation of the 2006 Alexandrium fundyense bloom in the Gulf of Maine: In-situ observations and numerical modeling. Continental Shelf Research. 2009;29:2069–2082. doi: 10.1016/j.csr.2009.07.012. doi:10.1016/j.csr.2009.07.012. [DOI] [PMC free article] [PubMed] [Google Scholar]
  35. Matrai P, Thompson B, Keller M. Circannual excystment of resting cysts of Alexandrium spp. from eastern Gulf of Maine populations. Deep Sea Research Part II: Topical Studies in Oceanography. 2005;52:2560–2568. doi:10.1016/j.dsr2.2005.06.013. [Google Scholar]
  36. McGillicuddy DJ, Anderson DM, Lynch DR, Townsend DW. Mechanisms regulating large-scale seasonal fluctuations in Alexandrium fundyense populations in the Gulf of Maine: Results from a physical–biological model. Deep Sea Research Part II: Topical Studies in Oceanography. 2005;52:2698–2714. doi:10.1016/j.dsr2.2005.06.021. [Google Scholar]
  37. McGillicuddy DJ, Townsend DW, He R, Keafer BA, Kleindinst JL, Li Y, Manning JP, Mountain DG, Thomas MA, Anderson DM. Suppression of the 2010 Alexandrium fundyense bloom by changes in physical, biological, and chemical properties of the Gulf of Maine. Limnology and Oceanography. 2011;56:2411–2426. doi: 10.4319/lo.2011.56.6.2411. [DOI] [PMC free article] [PubMed] [Google Scholar]
  38. Murphy AH. Skill scores based on the mean square error and their relationships to the correlation coefficient. Monthly Weather Review. 1988;116:2417–2424. [Google Scholar]
  39. Ní Rathaille A, Raine R. Seasonality in the excystment of Alexandrium minutum and Alexandrium tamarense in Irish coastal waters. Harmful Algae. 2011;10:629–635. doi:10.1016/j.hal.2011.04.015. [Google Scholar]
  40. Portnoy JW, Nowicki BL, Roman CT, Urish DW. The discharge of nitrate-contaminated groundwater from developed shoreline to marsh-fringed estuary. Water Resources Research. 1998;34:3095–3104. [Google Scholar]
  41. Raine R. A review of the biophysical interactions relevant to the promotion of HABs in stratified systems: The case study of Ireland. Deep Sea Research Part II: Topical Studies in Oceanography. 2014 doi:10.1016/j.dsr2.2013.06.021. [Google Scholar]
  42. Ralston DK, Keafer BA, Brosnahan ML, Anderson DM. Temperature dependence of an estuarine harmful algal bloom: Resolving interannual variability in bloom dynamics using a degree-day approach. Limnology and Oceanography. 2014;59:1112–1126. doi: 10.4319/lo.2014.59.4.1112. doi:10.4319/lo.2014.59.4.1112. [DOI] [PMC free article] [PubMed] [Google Scholar]
  43. Ralston DK, McGillicuddy DJ, Townsend DW. Asynchronous vertical migration and bimodal distribution of motile phytoplankton. Journal of Plankton Research. 2007;29:803–821. doi:10.1093/plankt/fbm061. [Google Scholar]
  44. Rodi W. Examples of calculation methods for flow and mixing in stratified fluids. Journal of Geophysical Research. 1987;92:5305–5328. [Google Scholar]
  45. Shumway SE. A Review of the Effects of Algal Blooms on Shellfish and Aquaculture. Journal of the World Aquaculture Society. 1990;21:65–104. doi:10.1111/j.1749-7345.1990.tb00529.x. [Google Scholar]
  46. Stock CA, McGillicuddy DJ, Jr., Solow AR, Anderson DM. Evaluating hypotheses for the initiation and development of Alexandrium fundyense blooms in the western Gulf of Maine using a coupled physical–biological model. Deep Sea Research Part II: Topical Studies in Oceanography. 2005;52:2715–2744. doi:10.1016/j.dsr2.2005.06.022. [Google Scholar]
  47. Stow CA, Jolliff J, McGillicuddy DJ, Jr., Doney SC, Allen JI, Friedrichs MAM, Rose KA, Wallhead P. Skill assessment for coupled biological/physical models of marine systems. Journal of Marine Systems. 2009;76:4–15. doi: 10.1016/j.jmarsys.2008.03.011. doi:10.1016/j.jmarsys.2008.03.011. [DOI] [PMC free article] [PubMed] [Google Scholar]
  48. Stumpf RP, Tomlinson MC, Calkins JA, Kirkpatrick B, Fisher K, Nierenberg K, Currier R, Wynne TT. Skill assessment for an operational algal bloom forecast system. Journal of Marine Systems. 2009;76:151–161. doi: 10.1016/j.jmarsys.2008.05.016. doi:10.1016/j.jmarsys.2008.05.016. [DOI] [PMC free article] [PubMed] [Google Scholar]
  49. Umlauf L, Burchard H. Second-order turbulence closure models for geophysical boundary layers. A review of recent work. Continental Shelf Research. 2005;25:795–827. doi:10.1016/j.csr.2004.08.004. [Google Scholar]
  50. Velo-Suárez L, Brosnahan ML, Anderson DM, McGillicuddy DJ., Jr. A Quantitative Assessment of the Role of the Parasite Amoebophrya in the Termination of Alexandrium fundyense Blooms within a Small Coastal Embayment. PLoS ONE. 2013;8:e81150. doi: 10.1371/journal.pone.0081150. doi:10.1371/journal.pone.0081150. [DOI] [PMC free article] [PubMed] [Google Scholar]
  51. Watras CJ, Chisholm SW, Anderson DM. Regulation of growth in an estuarine clone of Gonyaulax tamarensis Lebour: Salinity-dependent temperature responses. Journal of Experimental Marine Biology and Ecology. 1982;62:25–37. doi:10.1016/0022-0981(82)90214-3. [Google Scholar]