MurKy waters: Modeling the succession from r to K strategists (diatoms to dinoflagellates) following a nutrient release from a mining facility in Florida

The impacts of pulsed nutrient injections or extreme runoff events on marine ecosystems are far less studied than those associated with long‐term eutrophication, particularly in regard to mechanisms regulating the response of plankton community structure. Over 800 million liters of nutrient‐rich water from a fertilizer mine were discharged over a 2‐week period into Tampa Bay, Florida, in 2021, providing a unique opportunity to document the plankton response. A 3D‐coupled hydrodynamic biogeochemical model was developed to investigate this response and to understand the observed succession of a large, short diatom bloom followed by a secondary Karenia brevis bloom that lasted through the summer. The model reproduced the observed changes in nutrient concentration, total chlorophyll a, and diatom and K. brevis biomass in Tampa Bay. With a faster growth rate and spring temperature close to the optimal window of growth, diatoms had an initial competitive advantage, with 2/3 of the nutrient uptake due to ammonium and 1/3 due to nitrate. However, exhaustion of external nutrients led to the rapid decline of the diatom bloom, and the associated particular organic nitrogen sank onto the bay sediment. Enhanced sediment release of ammonium during the weeks following, and summer remineralization of dissolved organic nitrogen provided sufficient regenerated nitrogen to support slow‐growing K. brevis that could capitalize on low nutrient conditions. Modeling analysis largely confirmed Margalef's conceptual model of r to K‐selected species succession and provided additional insights into nutrient cycling supporting the initial diatom bloom and the subsequent bloom of a slow‐growing harmful algal species.

upwelling conditions), small diatoms cells should dominate the first phase of bloom response.Diatoms have a high surface/volume ratio and high potential growth rates.In the later stage of blooms, when nutrients have been depleted and stratification has been established, diatoms tend to sink out, and larger dinoflagellate species with lower growth rates predominate.In the field of plant ecology, Grime (1977) provided another species succession theory describing a three-way trade-off between organisms, known as the competitors-stress tolerators-ruderals (CSR) triangle, which described a succession from high production strategists (competitor) to stress-tolerant strategists to Ruderal species along the decrease in available nutrient resources.Since Margalef's papers, it has been shown diatoms are true photoautotrophs, while dinoflagellates have a more flexible nutrient acquisition strategy, strategies that affect biogeochemical cycling with feedback implications for growth that challenge the classic Margalef's successional models (reviewed by Glibert and Mitra 2022).Furthermore, it remains unclear if plankton response to a short pulse of intense nutrient loading can be interpreted within this theoretical framework, especially the nutrient recycling processes that could support a succession of plankton blooms over an extended period of time.
In 2021, a major nutrient injection into Tampa Bay, Florida, occurred over 2 weeks following the breach of a holding pond of wastewater from an abandoned phosphorus (P) mining plant at Piney Point.This spring nutrient pulse provided the opportunity to observe the succession of phytoplankton taxa and to use a modeling approach to understand the underlying processes that led to a several months-long successional sequence from diatoms to the dinoflagellate Karenia brevis that was able to be sustained and thrive through the summer, a time when K. brevis typically does not bloom.
Chemical spills/releases from mining facilities or fertilizer plants pose a major threat to coastal oceans.Over the past decade, about 55 major spills/emergence releases were reported around the world, discharging hundreds of million cubic meters of wastewater and nutrients into the ocean and sensitive nearshore habitats (WISE 2022).The released nutrients stimulate primary production and may lead to cascading effects on entire marine and coastal ecosystems (Ghosh et al. 1999;Yeager et al. 2005).Furthermore, disproportionate releases of certain nutrient types can change nutrient stoichiometry in the water, favoring harmful algal species, many of which thrive in unbalanced nutrients (Glibert 2016).For example, over 41 million liters of swine wastewater was discharged into the New River estuary after the accidental rupture of a waste-holding lagoon in 1995, triggering a noxious algal bloom in the following 2 weeks (Burkholder et al. 1997).An anomalous release of phosphogypsum into the Grand Bay National Estuarine Research Reserve in 2005 resulted in a 10-fold increase in PO 3À 4 concentration and an 18-56% increase in chlorophyll a (Chl a; Beck et al. 2018).
Similarly, stormwater runoff over polluted land or agricultural fields during flood events can bring huge amounts of pollutants and nutrients into rivers, bays, and coastal oceans over a short period of time (Paerl et al. 2006a,b;Chen et al. 2019).For example, Hurricanes Dennis, Floyd, and Irene passed through North Carolina in 1999, delivering over half of the typical annual nitrogen (N) loading in a very short time span and resulting in 3-5 times higher phytoplankton biomass over the next 6 months (Paerl et al. 2001).In particular, cryptophytes and chlorophytes responded rapidly to storm events as they have higher inherent growth rates than other taxa (Hall et al. 2013;Paerl et al. 2014).In comparison, dinoflagellates with low growth rates were susceptible to flushing during high flow conditions, but their biomass increased during moderately wet storms when higher nutrient loadings and limited flushing provided a favorable combination for their growth (Paerl et al. 2001;Hall et al. 2008).More recently, Hurricane Ian produced extensive watershed flooding and a large estuarine plume laden with nutrient and dissolved organic matter, precipitating an intense K. brevis bloom on the West Florida Shelf in 2022 (https:// www.flickr.com/photos/myfwc/sets/72157635398013168).
Despite the previous studies, the impacts of pulsed nutrient injections due to chemical spills/emergence releases or extreme runoff events on marine ecosystems are not as well understood as those related to long-term eutrophication (Rabalais et al. 2009;Harding et al. 2015;Song et al. 2022).Although Margalef's successional conceptual models would suggest that a strong pulse of nutrients should initially favor fast-growing phytoplankton species (Margalef 1967b;Reynolds 1987), much is unknown about how intensive nutrient pulses affect the plankton community composition over the longer term, including a delayed response of slowgrowing harmful algal species and potential long-lasting changes in marine ecosystems (Paerl et al. 2001;Valdes-Weaver et al. 2006).The fertilizer wastewater discharge into Tampa Bay in 2021 provided a unique opportunity to test ecological concepts of phytoplankton succession and biogeochemical cycling in response to a strong pulse of nutrient loading; such intense nutrient pulses are becoming more common in a warming climate featuring more extreme precipitation events (Knutson et al. 2019(Knutson et al. , 2020)).
With one of the largest storages of phosphatic ore in the United States, central Florida hosts numerous P mining facilities and related phosphate (PO 3À 4 ) processing plants.The Piney Point plant, located 3 km from the shore of Tampa Bay, left a legacy of over 1.9 billion liters of such nutrient-rich disposal water (Garrett et al. 2011), creating a perilous situation prone to accidental spills/emergence releases due to both human errors or structural breaches that can occur during hurricanes and tropical storms and increasing volumes from rainy season inputs.Leakages were detected in the plastic liner in the reservoir's walls in late March 2021 (Beck et al. 2022).To avoid a large-scale collapse of the reservoir structure, an emergency release of the stacked water was authorized by the Florida Department of Environmental Protection.Over a period of 10 days, approximately 814 million liters of this water was discharged into Tampa Bay through the outfall near Port Manatee, yielding a nutrient load in this period equivalent to that normally received by the lower Tampa Bay annually (Beck et al. 2022).
To understand how the phytoplankton community responded to the nutrient release from the Piney Point event and to interpret the observed succession of diatom and blooms following the nutrient release, a coupled hydrodynamic-biogeochemical model was developed for Tampa Bay and the West Florida Shelf, using the modeling framework that was previously used to study water quality and harmful algal bloom (HABs) in Chesapeake Bay (Testa et al. 2014;Li et al. 2016;Ni et al. 2019;Ni et al. 2020;Li et al. 2021;Zhang et al. 2021;Li et al. 2022;Ni and Li 2023).Specifically, we investigated whether the rapid nutrient pulse would result in a succession of r to K strategists, how and whether the influx of different forms and amounts of nutrients, through various recycling pathways, would support a succession of phytoplankton taxa, and the extent to which the unusual bloom of K. brevis that was subsequently sustained throughout the summer months may have been supported and enhanced by this unusual nutrient load.

Methods
A 3D coupled hydrodynamic-biogeochemical model was developed to investigate the plankton response to the nutrient release from the Piney Point plant.The hydrodynamic model is based on the Regional Ocean Modeling System (ROMS) (Shchepetkin and McWilliams 2005;Haidvogel et al. 2008) and the biogeochemical model is based on the Row Column Aesop (RCA) model (Di Toro 2001;HydroQual 2004;Isleib et al. 2007).

Hydrodynamic model (ROMS)
The ROMS model was configured to cover Tampa Bay and the West Florida Shelf, spanning 87.63 W to 80.20 W and 24.14 N to 31.95N.An orthogonal curvilinear coordinate is used to follow the general orientation of the West Florida Shelf.There are 456 Â 186 grid points in the horizontal directions, with the grid resolution ranging from $ 800 m inside Tampa Bay $ 3 km at the offshore boundary.A terrain-following coordinate is used in the vertical direction and consists of 21 layers.The vertical eddy viscosity and diffusivity are parameterized using the k À ε turbulence closure scheme with the background value set at 10 À6 m 2 s À1 (Lanerolle and Patchen 2011).The Smagorinsky parameterization is used to calculate the horizontal eddy viscosity and diffusivity (Smagorinsky 1963).
The ROMS model is driven by the atmospheric, oceanic, and riverine forcing at its boundaries.At the sea surface, the air-sea fluxes of heat and momentum are calculated with the standard bulk formula (Fairall et al. 1996;Fairall et al. 2003).Wind speed, air pressure, and heat fluxes were interpolated from the NCEP Climate Forecast System Reanalysis (CFSR) dataset (https://rda.ucar.edu/datasets/ds094.0/).River discharge, salinity, and temperature were prescribed at the boundaries of 11 major tributaries, including the Lake Tarpon Canal, Hillsborough River, Alafia River, Little Manatee River, and Manatee River in Tampa Bay.The flows were obtained from daily measurements at USGS gauging stations in the rivers (https://waterdata.usgs.gov/nwis)and the reported wastewater release (Beck et al. 2022), while daily riverine salinity and temperature data were obtained from the Coastal and Heartland National Estuary Partnership Water Atlas (CHNEP, https://chnep.wateratlas.usf.edu/).
At the offshore boundaries, the Chapman boundary condition was used for free surface elevation and the Flather boundary condition was used for the barotropic velocities (Flather 1976;Chapman 1985).The Orlanski radiation boundary condition was used for the baroclinic velocities, temperature, and salinity (Orlanski 1976).Tidal components, including the tidal elevation and barotropic velocity, were extracted from the Oregon State University global inverse tidal model TPXO8 (Egbert and Erofeeva 2002), while the subtidal surface elevation was interpolated from the fine-resolution HYCOM model for the entire Gulf of Mexico (https://www.hycom.org/data/gomu0pt04/expt-90pt1m000).Since the Loop Current may interact with the West Florida Shelf (Hetland et al. 2001;Liu et al. 2016), a buffer zone is introduced at the 10 outermost grid points near the offshore boundary, in which temperature, salinity, and baroclinic velocities in the ROMS are nudged to those produced from HYCOM.To ensure a smooth transition between the two models, the nudging coefficient increases exponentially from 0 at the innermost points in the buffer zone to 0.25 at the offshore boundary.
The ROMS model was initialized using climatological temperature and salinity values from the World Ocean Atlas (https://www.ncei.noaa.gov/products/world-ocean-atlas)and ran for a spin-up period of 1 year.The outputs from this model run were used to set the initial condition for the hindcast model run of the chemical release on 01 March 2021 and the model integration ran until 31 December 2021.

Biogeochemical HAM models (RCA-HAB)
The RCA biogeochemical module (HydroQual 2004) was hardcoupled to ROMS by importing all RCA codes and subroutines into the ROMS software system, thus enabling parallel computations (Fig. 1).RCA includes a water-column component (Isleib et al. 2007) and a two-layer sediment diagenesis model (Di Toro 2001;Brady et al. 2013).The water-column model has state variables describing dissolved inorganic nitrogen (DIN: NO À 3 and NH þ 4 ), dissolved inorganic phosphorus (PO 3À 4 ), dissolved silica (SiO 4À 4 ), dissolved and particulate organic forms of nitrogen and phosphorus (DON, PON, DOP, and POP) of varying levels of reactivity, three functional groups of phytoplankton, and dissolved O 2 (DO).The sediment diagenesis model has one aerobic layer and one anaerobic layer, simulating the cycling of carbon (C), N, P, Si, and DO equivalences (sulfur and methane).The water column biogeochemical model and sediment diagenesis model are coupled through the water-sediment fluxes, including particulate matter deposition from the water column and diffusive exchanges of dissolved inorganic nutrients.
In the RCA model configured for Tampa Bay and West Florida Shelf, the three major phytoplankton groups are diatoms, K. brevis (HAB) and cyanobacteria Synechococcus sp.The equations describing the temporal changes in the biomass of these plankton groups are shown below: where the subscripts (i = 1, 2, and 3) denote the biomass associated with diatoms, K. brevis, and Synechococcus, respectively.In Eq. ( 1), the temporal change ∂P i =∂t at a fixed location is on the left-hand side.Terms on the right-hand side can be divided into three groups.
, where the velocity fields (u, v, w) are obtained from the ROMS hydrodynamic model.The third group is the turbulent diffusion term where the horizontal diffusivity K h and vertical eddy diffusivity K v are also calculated in ROMS.μ max i is the maximum growth rate of phytoplankton while G T i describes the growth dependence on temperature, G par i is the P-I curve, and G N i represents the nutrient limitation.For K. brevis, G Sal is added to account for the salinity regulation of its growth rate according to published experimental data (Magaña and Villareal 2006).
Following the approach taken in the CoSiNE model (Chai et al. 2002;Liu et al. 2018), the nutrient limitation on phytoplankton growth is separated into two parts, one representing the NO ): where the half-saturation constants for NO À 3 , NH þ 4 , PO 3À 4 , and SiO 4À 4 are denoted by , k P i , and k DSii , respectively.
The equations for NH þ 4 and NO À 3 are as follows: where nutrient and chlorophyll a concentration were monitored, and the locations of four stations (B, C, D, and E) selected for the following modeling analysis.The yellow dashed line marks the mid-bay area near the outfall that is used for the regional integration shown in Figs.7 and 8.
where the subscripts (L and R) denote the labile and refractory form of ON. r and n refer to the remineralization and nitrification coefficients.In both equations, the temporal change at a fixed location is on the left-hand side.The source/sink terms of Eq. ( 8) includes the remineralization from DON (r L LDON,r R RDON) and regenerated NH þ 4 from phytoplankton mortality (R NH þ 4 ), while the loss terms are phytoplankton uptake (U NH þ 4 ) and nitrification (nNH þ 4 ).In Eq. ( 9), NO À 3 gains through nitrification (nNH þ 4 Þ but is lost due to phytoplankton uptake (U NO À 3 ) and denitrification (DNIT).In addition, NO À 3 and NH þ 4 are affected by advection and diffusion processes through the other terms on the right-hand side of Eqs. ( 8) and ( 9), such as the riverine inputs, benthic flux, and the upwelling of deep ocean nutrients.The remineralization rates of the labile and refraction dissolved organic matter (DOM) are given by and the nitrification rate is given by where r Lt , r Rt , and n t are the rates at 20 C, and θ rL , θ rR , and θ n are the temperature coefficients.TPHYT refers to the total phytoplankton biomass, k phyt is the corresponding half-saturation coefficient, and k mNIT is the half-saturation coefficient for nitrification.
Values of the parameters used in the RCA model were chosen according to the published literature (Testa et al. 2014;Li et al. 2016;Ni et al. 2019;Ni et al. 2020;Ni and Li 2023) as well as the previous studies on the West Florida Shelf (Walsh et al. 2003;Walsh et al. 2012) 4 , and dissolved Si concentrations on the shelf and their time series at the offshore boundaries.Meanwhile, the initial concentrations of DON on the shelf were estimated by cruise samples from 2007 to 2010 (Dixon et al. 2014).Inside Tampa Bay, the initial nutrient concentrations were interpolated from the CHNEP data, which almost cover the whole bay.Since observations of diatoms and K. brevis within Tampa Bay were limited during the spring of 2021, their initial concentrations were estimated by combining the cell count data collected by the Mote Marine Laboratory with those reported to the Fish and Wildlife Research Institute website (https://myfwc.com/research/redtide/).

Short-term response-A large diatom bloom
The discharge of wastewater from the Piney Point plant began on 30 March and ceased on 09 April 2021, releasing about 180 metric tons of total N, mostly in the form of NH þ 4 (15,000 μmol L À1 ) and little in the form of NO À 3 (0.029 μmol L À1 ), and over 110 metric tons of P (5100 μmol L À1 ) into lower Tampa Bay.The concentration of PO 3À 4 reached a maximum of $ 20 μmol L À1 near the release site (Fig. 2a).As the currents moved back and forth and dispersed the wastewater, PO 3À 4 was distributed throughout the estuary.NH þ 4 also spread across the estuary, but a highconcentration patch (20-40 μmol L À1 ) was mostly confined to the mid-estuary region surrounding the outfall site (Fig. 2e-h).

Concentrations of NH þ
4 showed a marked decline from 10 to 15 April that was not seen in the PO 3À 4 distributions (compare Fig. 2g,h with Fig. 2c,d).In comparison with NH þ 4 , NO À 3 was mostly confined around the release site (Fig. 2i-l).Its maximum concentration was much lower than that of NH þ 4 , $ 4 μmol L À1 near the outfall, but was 2-3 orders of magnitude higher than that in the wastewater.The co-locations of high NH þ 4 and NO À 3 patches near the outfall suggest that NO À 3 was generated by nitrification.
Since phytoplankton inside Tampa Bay are normally limited by N availability rather than P availability (Wang et al. 1999;Greening and Janicki 2006), the injection of high concentrations of NH þ 4 and NO À 3 resulted in the dominance of fast-growing diatoms (Fig. 2m-p).A large diatom bloom (with a concentration up to 1 mg C L À1 ) developed on 10 April and its biomass continued to accumulate through 12 April at its peak, while NH þ 4 and NO À 3 declined due to nutrient uptake.There was a phase lag in the temporal evolution of NH þ 4 , NO À 3 , and diatom biomass following the nutrient release (Fig. 3), but the magnitude differed by station location.Sta.A is located close to the outfall.NH þ 4 reached a peak of 24 μmol L À1 at the end of the controlled release on 09 April and decreased to a low background concentration ($ 0.2 μmol L À1 ) by 19 April (Fig. 3a).Although NO À 3 concentration in the wastewater was comparatively low, NO À 3 increased presumably due to nitrification, reaching its peak concentration around 10 April.Diatom biomass increased as NO À 3 and NH þ 4 became available and reached a peak bloom size of $ 1 mg C L À1 rapidly, and decreased gradually to a background concentration of 0.3 mg C L À1 on 25 April.A phase lag can be noticed between the peak nutrient concentration and peak diatom biomass, but there was virtually no lag between the peaks of NH þ 4 and NO À 3 concentration.Although NH þ 4 and NO À 3 at the downstream Sta.B were considerably lower, a diatom bloom still developed there, with the maximum biomass of 0.5 mg C L À1 on 14 April (Fig. 3b).The diatom biomass at Sta. C near the mouth of Tampa Bay was far lower than that at Stas.A and B and a temporal response to the nutrient pulse could not be detected there (Fig. 3c).
These model predictions are in reasonable agreement with the observed nutrient concentrations and Chl a, as shown in Fig. 4 for station A (denoted as Piney 6 station by the Florida Department of Environmental Protection).The observed NH þ 4 concentration rapidly increased from 30 March to 09 April, with a maximum concentration of 20 μmol L À1 .The agreement between the model-predicted NH þ 4 and observed values includes the months of May and June, well after the initial nutrient injection, when both the observed and predicted NH þ 4 fell onto 0.2 μmol L À1 (Fig. 4a).The predicted NO À 3 time a.Yoder (1979) b.Kaeriyama et al. (2011) c.Magaña andVillareal (2006) d. Vargo (2009) e. Tilney et al. (2019) f.Kana and Glibert (1987) g.Mackey et al. (2013) h.Walsh et al. (2012) i. Walsh et al. (2003) j.This study k.Jeong et al. (2005) l.Glibert et al. (2009) m.Procise ( 2012) n.Smayda (1971) o.Eppley et al. (1969) p. Romeo and Fisher (1982) q.Killberg-  series is also in good agreement with the observed NO À 3 , with a peak concentration of around 1.7 μmol L À1 (Fig. 4b).There was no phase lag between the peak NH þ 4 and peak NO À 3 , indicating that nitrification occurred at a fast rate.The surface NH þ 4 and NO À 3 concentrations dropped to very low levels after the diatom bloom in April, indicating that most of the nutrients discharged from the chemical release were exhausted by the beginning of May.The predicted and observed PO 3À 4 time series are also in very good agreement, rising from an initial concentration of 2 μmol L À1 to a maximum of 10 μmol L À1 on 09 April and declining to $ 2 μmol L À1 after early May (Fig. 4c).The observed Chl a showed a large peak of 18 μg L À1 around 12 April, followed by a second peak of $ 10 μg L À1 near the end of June (Fig. 4d).The predicted Chl a time series is in very good agreement with the observed one, not only in terms of the bloom magnitude but also in the timing of the two blooms: the first larger peak around 12 April and the second peak around 18 June.Diatom Chl a accounted for $ 99% of the total Chl a between 30 March and 30 April (Fig. 4d).In contrast, K. brevis Chl a accounted for $ 70% of the total during the summer (01 June to 31 August).A short transition period can then be noticed between 01 May and 15 May when diatom Chl a decreased and K. brevis Chl a started to increase (Fig. 4d).Long-term response-A summer K. brevis bloom A K. brevis bloom occurred $ 1 month after the diatom bloom had disappeared and the surface NO À 3 and NH þ 4 concentrations dropped to low levels ($ 0.01 μmol L À1 and $ 0.2 μmol L À1 ), as shown in Fig. 4. Most of NO À 3 was depleted during the diatom bloom in April and was at near-zero concentrations except on the shelf during May-August (Fig. 5a-d).The surface NH þ 4 concentration was maintained at a level of 0.1-0.7 μmol L À1 , in the range of the reported halfsaturation coefficient (0.4 μmol L À1 ) for the NH þ 4 uptake by K. brevis (Fig. 5e-h; Table 1).More tellingly, the bottom NH þ 4 concentration was higher and in the range of 0.5-1 μmol L À1 (Fig. 5i-l).Given that the water column was well mixed in the shallow waters of Tampa Bay, NH þ 4 from bottom waters was constantly injected into the surface euphotic layer, providing a source of N to support K. brevis growth.The K. brevis cell density began to increase in late April, reached a peak bloom density of 10 6 cells L À1 in early July, and lasted until the end of August (Fig. 5m-p).Its cell density was higher at the coast and in the lower Tampa Bay, in contrast to the diatom bloom which was centered in the middle of Tampa Bay (compare Fig. 2o,p with Fig. 5o,p).
The time series in Fig. 6 provides additional insights into how NH þ 4 was regenerated and how this regenerated NH þ  0.08 μmol L À1 d À1 around 15 May, and was maintained at this level until the end of September.The higher remineralization of DON on 1-15 May was likely associated with the decaying diatom bloom (compare Figs. 4d,6c).In contrast, NH þ 4 efflux from the sediment was lower in May and September but considerably higher in June-August (Fig. 6a).The sediment diagenesis rate increased with temperature and was most active during the hot summer months.It is interesting to note that NH þ 4 release from the sediment was 2-3 times higher than the integrated remineralization of DON in the water column at Sta.A (Fig. 6c).A similar story with regard to NH þ 4 regeneration can be seen at Sta. C located near the mouth of Tampa Bay (Fig. 6b).At this location the sediment-water NH þ 4 flux and the remineralization of DON in the water column were of similar magnitude, as a large portion of the PON derived from the fast-sinking diatoms had already been deposited onto the bay sediment further upstream (Fig. 6d).A comparison between the NH þ 4 and K. brevis cell density time series suggests the NH þ 4 uptake drove the K. brevis bloom (Fig. 6a,b).Indeed, at Sta. A, NH þ 4 uptake increased from 0.0 to 0.3 μmol L À1 d À1 and K. brevis cell density increased from 0.01 to 0.16 mg C L À1 , while NH þ 4 declined from 0.25 μmol L À1 on 01 May to 0.18 μmol L À1 on 30 June (Fig. 6a,b).At Sta. C near the mouth of Tampa Bay, declining NH þ 4 concentration peaked around mid-June, and then went down in August.Also, snapshots of deposition and benthic flux in the whole bay (Supporting Information Fig. S2) showed that areas with nutrient pulse and diatom bloom had larger benthic flux in the subsequent summer month.
In the current ROMS-RCA-HAB models a simple method was used to parameterize the mixotrophic feeding of K. brevis by adding a grazing term Rg kb P 2 = P 2 þ k kb ½ ð Þ P i 2 in the equation for the cyanobacteria Synechococcus (Eq.1).As K. brevis biomass increased, Synechococcus decreased (Supporting Information Fig. S3).A similar inverse relationship between the two taxa was previously shown on West Florida Shelf through signature pigment distributions (Heil et al. 2007).A more sophisticated description of mixotrophic feeding would consider the interactions of C, N, and P within a mixotrophic cell, accounting for photosynthesis, inorganic nutrient uptake and consumption and digestion of prey, as well as the potential contribution of photosynthesis from the incorporation of a kleptochloroplast (i.e., ingested photosystem; Flynn and Mitra 2009).Such a model has recently been developed for Karlodinum veneficum and its cryptophyte prey in the Chesapeake Bay (Li et al. 2022), and efforts are ongoing to develop such a model for K. brevis.
In this study dead fish killed by brevotoxin were not considered, but previous studies suggested they could generate an additional nutrient source to support K. brevis blooms (Heil et al. 2014).However, the close agreement between the predicted and observed NH þ 4 and Chl a concentration suggests that nutrient recycling due to dead fish was unlikely to be a major source of nutrients during the summer of 2021.Macroalgae is another ecosystem component that could play a role in the plankton response to the nutrient release from the Piney Point plant.Macroalgae, the cyanobacterium Dapis pleousa, bloomed between the diatom and K. brevis blooms and the decomposition of the macroalgae could provide a nutrient pathway connecting the two blooms since similar macroalgal blooms were observed following a nutrient release from Piney Point plant in 2003-2004(Switzer et al. 2011).In addition, Tropical Storm Elsa (2021) passed through Tampa Bay during the peak of the summer K. brevis bloom (https:// www.weather.gov/tae/2021_tropicalstorm_elsa).It could also have played a role in affecting or even promoting the summer K. brevis bloom and fish kills, either through enhanced watercolumn mixing and nutrient recycling or additional nutrient flows from the watershed.However, as an early-season Category 1 hurricane, the post-storm flood into Tampa Bay was relatively small.
The modeling study of the diatom and K. brevis response to the wastewater release in Tampa Bay offers a comparison with the extensive observational analysis of phytoplankton community response to tropical storms in the Neuse River estuary and Pamlico Sound (Paerl et al. 2018).The short-term response in both systems seems to fit well with Margalef's succession model that favors fast-growing species: diatoms in the Tampa Bay and cryptophytes and chlorophytes in the Neuse River.Although diatoms also have fast growth rates in the Neuse River-Pamlico Sound, they constitute less than 20% of phytoplankton biomass which may explain their relatively muted response to storms.Furthermore, there are differences between the two systems and two types of pulsed nutrient loadings that should be noted.Although the flushing time was a key parameter in regulating phytoplankton biomass accumulation in the Neuse River, it was not an important factor in Tampa Bay since the large nutrient loading was delivered by extremely high nutrient concentration and relatively low flow discharge rate.

Conclusions
In summary, the plankton response to the nutrient event from the Piney Point plant can be divided into three phases (Fig. 9).In the first phase, intensive effluent rapidly raised the Fig. 9.A schematic diagram of the phytoplankton group succession after a massive nutrient input event.
nutrient concentration near the outfall and triggered a bloom of fast-growing diatom species.As a response to this large NH þ 4 influx, local nitrification flux also increased and resulted in increasing NO À 3 , which would provide the required NO À 3 for diatom metabolism and growth.In the second phase, nutrient concentrations dropped as the diatom exhausted the external input of nutrients.Species with higher nutrient affinity and slower growth rate benefit from the resulting relatively oligotrophic environment.A more moderate but prolonged bloom was able to develop with a moderate amount of generated nutrients.POM stored in the sediments provided an additional source for the second bloom through remineralization and benthic flux.Compared with a previous study on the N cycle in the bay (Bronk et al. 2014), our model may actually underestimate the remineralization flux.
The succession of species herein followed the general conceptual model of r to K-selected species sequence predicted by Margalef (Margalef 1958(Margalef , 1967a,b),b).This study supports his descriptive sequence through a detailed analysis of the physiological, biogeochemical, and physical processes.The early descriptions of phytoplankton succession in relation to nutrient supply (Margalef 1958(Margalef , 1967a,b),b), which was initially developed in the context of upwelling systems may be as common in aquatic systems receiving intensive nutrient injections from wastewater releases or storm events.These results also underscore that from the perspective of management response to such events, the initial response should not be taken as the only ecosystem response, as lagged responses may be as important from a human health or ecosystem health perspective.Moreover, the enduring consequences of this incident can provide valuable lessons for similar facilities worldwide.Most retired phosphorus fertilizer plants pose significant risks to the surrounding environment due to toxic wastewater and radioactive byproducts, but the level of investment in research and management efforts remains insufficient (Nelson et al. 2021).As such, proactive measures and stronger oversight are necessary for businesses that even continue to operate as usual, to prevent taxpayers and the environment from bearing the external costs of industry.

Fig. 1 .
Fig. 1.(a) Schematic of the coupled hydrodynamic (ROMS)-biogeochemical (RCA)-HAB modeling framework that is one-way nested within the Hybrid Coordinate Ocean Model-Navy Coupled Ocean Data Assimilation global reanalysis (HYCOM-NCODA); (b) a map of the Tampa Bay (depths shown by the gray color scale), the location of the wastewater release outfall from the chemical plant (the orange filled circle), the location of Florida Department of Environmental Protection station Piney 6 (marked the blue dot with the label Sta.A)where nutrient and chlorophyll a concentration were monitored, and the locations of four stations (B, C, D, and E) selected for the following modeling analysis.The yellow dashed line marks the mid-bay area near the outfall that is used for the regional integration shown in Figs.7 and 8.

Fig. 2 .
Fig. 2. Snapshots of nutrients and diatom biomass after the release.(a-d) Surface PO 3À

4
stimulated the K. brevis bloom.As most of the NH þ 4 from the nutrient discharge was exhausted by 01 May, two potential sources of DON (the first two terms on the righthand side of Eq. 8), and the NH þ 4 release across the sedimentwater interface due to the sediment diagenesis of particulate organic matter (POM) that had settled down onto the bay bottom.At Sta.A close to the outfall site, the remineralization of DON was 0.14 μmol L À1 d À1 on 01 May, decreased to

Fig. 5 .
Fig. 5. Snapshots of nutrients and Karenia brevis biomass after the release.(a-f) Surface NO À (Walsh et al. 2012;Iversen and Ploug 2013)k terms, including the gross growth rate G i , the respiration rate Rr i , and mortality rate Rg i .The loss rate due to sinking is prescribed for the diatoms with a constant sinking velocity w(Walsh et al. 2012;Iversen and Ploug 2013), while an additional mortality rate on Synechococcus is included based on the K. brevis biomass P 2 .The second group includes the advective terms ∂uP . Table1lists the key parameters, their values, and units.

Table 1 .
Part of model parameters for Diatom, Karenia brevis, Synechococcus (i = 1, 2, 3), and nitrogen cycling.The references for the values are labelled with lower-case superscript letters.