Biophysical interactions control the progression of harmful algal blooms in Chesapeake Bay: A novel Lagrangian particle tracking model with mixotrophic growth and vertical migration

Climate change and nutrient pollution contribute to the expanding global footprint of harmful algal blooms. To better predict their spatial distributions and disentangle biophysical controls, a novel Lagrangian particle tracking and biological (LPT‐Bio) model was developed with a high‐resolution numerical model and remote sensing. The LPT‐Bio model integrates the advantages of Lagrangian and Eulerian approaches by explicitly simulating algal bloom dynamics, algal biomass change, and diel vertical migrations along predicted trajectories. The model successfully captured the intensity and extent of the 2020 Margalefidinium polykrikoides bloom in the lower Chesapeake Bay and resolved fine‐scale structures of bloom patchiness, demonstrating a reliable prediction skill for 7–10 d. The fully coupled LPT‐Bio model initialized/calibrated by remote sensing and controlled by ambient environmental conditions appeared to be a powerful approach to predicting transport pathways, identifying bloom hotspots, resolving concentration variations at subgrid scales, and investigating responses of HABs to changing environmental conditions and human interference.

tracking and biological (LPT-Bio) model was developed with a high-resolution numerical model and remote sensing. The LPT-Bio model integrates the advantages of Lagrangian and Eulerian approaches by explicitly simulating algal bloom dynamics, algal biomass change, and diel vertical migrations along predicted trajectories. The model successfully captured the intensity and extent of the 2020 Margalefidinium polykrikoides bloom in the lower Chesapeake Bay and resolved fine-scale structures of bloom patchiness, demonstrating a reliable prediction skill for 7-10 d. The fully coupled LPT-Bio model initialized/calibrated by remote sensing and controlled by ambient environmental conditions appeared to be a powerful approach to predicting transport pathways, identifying bloom hotspots, resolving concentration variations at subgrid scales, and investigating responses of HABs to changing environmental conditions and human interference.
Harmful algal blooms (HABs) are increasing globally, posing an increasing threat to public and ecosystem health and coastal economies Glibert 2020). Climate change and escalating nutrient pollution contribute to the expanding footprint of HABs (Glibert 2020). A better understanding of the fundamental processes causing blooms to form and more accurate predictions of HABs and their responses to changing ocean conditions depend on developing new monitoring and laboratory techniques and numerical modeling . Fueled by the increased frequency and geographic coverage of observational data, many approaches to formulating HAB models have been established (Franks 2018), including conceptual, statistical, and processbased/mechanistic models (McGillicuddy Jr 2010; Anderson et al. 2015;Ralston and Moore 2020).
Among techniques used to simulate HABs, Lagrangian particle tracking (LPT) is proving to be a powerful tool for identifying sites of HAB initiations, transport pathways, and dispersion patterns in coastal regions (Onitsuka et al. 2010;Velo-Su arez et al. 2010;Giddings et al. 2014;Gillibrand et al. 2016;Pinto et al. 2016;Weisberg et al. 2019;Lim et al. 2021). Combined with monitoring and hydrodynamic numerical modeling, LPT models are potential tools for HAB early warning (Wynne et al. 2013;Cusack et al. 2016;Davidson et al. 2021). LPT models are initiated by introducing numerical particles into a system and instructing them to be transported passively through physical processes. Most previous studies use these models to represent the transport of harmful algae without considering algal bloom dynamics, algal cell behaviors, or other environmental factors affecting algal growth. The prediction accuracy of these applications needs to be interpreted with caution since, while transported as particles, algal blooms are also highly regulated by ambient environmental conditions, that is, different algal species have growth optima and respond to environmental conditions favorable for their proliferation.
Several LPT models have imparted particles with behaviors/ responses consistent with the organism of interest's biology and behavioral characteristics, for example, temperature-, salinity-, and light-regulated growth (Aleynik et al. 2016;Gillibrand et al. 2016;Zhou et al. 2021), and vertical migrations (Henrichs et al. 2015). Yet, no LPT models to date have incorporated nutrient limitation and mixotrophic growth, factors that may contribute to the competitive success of dinoflagellate groups (Mulholland et al. 2018;Hofmann et al. 2021;Qin et al. 2021). To better predict distributions of harmful algae and their response to the changing environment as they are transported to connected waterways and to disentangle biophysical controls, we developed a novel LPT and biological (LPT-Bio) model that integrates the complex growth strategies of mixotrophic dinoflagellates with physicochemical environmental factors and coupled that with an advanced high-resolution numerical hydrodynamic model and remote sensing. The LPT-Bio model explicitly simulates algal growth, respiration, mortality, and diel vertical migration (DVM) patterns along the predicted transport trajectories.
Blooms of the harmful algal species Margalefidinium polykrikoides have become a global concern and the geographical distributions extend to coastal areas worldwide (Iwataki et al. 2008;Kudela and Gobler 2012;Al-Azri et al. 2014;L opez-Cortés et al. 2019;Roselli et al. 2020;Yñiguez et al. 2021). Nearly annual late summer blooms of M. polykrikoides in the Chesapeake Bay (CB) threaten many fisheries and aquaculture industries (Hudson 2018). Blooms of M. polykrikoides have expanded from those first-described in the York River (Mackiernan 1968) to those reported later in the mainstem and lower CB tributaries (Marshall and Egerton 2009;Mulholland et al. 2009;Morse et al. 2011Morse et al. , 2013. In the late summer of 2020, high chlorophyll a (Chl a) concentrations were evident in the lower CB and the adjacent Atlantic Ocean based on Ocean Land Color Imager Level 3 products. Several coastal and CB monitoring programs confirmed the bloom to be M. polykrikoides. Satellite images recorded the complete bloom cycle, including the uncommon phenomenon of M. polykrikoides export from the lower CB to coastal Atlantic waters following Tropical Storm (TS) Isaias (Xiong et al. 2022). A similar expansion was only reported in 2007 (Mulholland et al. 2009) although the CB outflow plume persistently delivers materials toward the Atlantic Ocean (Jiang and Xia 2016). Currently, no model exists to accurately predict the trajectory and extent of M. polykrikoides in the CB and the adjacent continental shelf once a bloom has initiated. The extensive long-term observations and recent research parameterizing the key physiological processes and behavioral strategies of the regional M. polykrikoides (Hofmann et al. 2021;Qin et al. 2021) make the lower CB an ideal test bed for validating the prediction ability of the LPT-Bio model. The advent of high-resolution daily ocean color imagery covering the CB region and the adjacent continental shelf (Wynne et al. 2021) further provides a unique opportunity to study HABs by merging the LPT-Bio model with satellite data. The newly developed LPT-Bio model and the successful simulation of the 2020 bloom demonstrate the predictive potential for HABs in the CB.

Bloom monitoring
Surveillance data for the 2020 M. polykrikoides bloom include dataflow data, fixed site observations, water sample collections, and satellite imagery. The dataflow mappings and continuous measurements of Chl a at one fixed site inside bloom-impacted regions were collected by the Hampton Roads Sanitation District. The cell counts and species identifications via water samplings were made using standard microscopic methods at Old Dominion University and Virginia Institute of Marine Science and were reported to the Virginia Department of Health (VDH) for inclusion on the algal bloom surveillance map (https://www.vdh.virginia.gov/waterbornehazards-control/algal-bloom-surveillance-map/).
The satellite-based daily monitoring, at 300 m resolution, was derived from Copernicus Sentinel-3 satellite data from the European Organisation for the Exploitation of Meteorological Satellites (EUMETSAT) and processed by NOAA's National Centers for Coastal Ocean Science (https://coastwatch.noaa. gov/cw_html/NCCOS.html). The satellite-derived Chl a were determined by a near-infrared to red ratio (Gilerson et al. 2010;Wynne et al. 2021) and the satellite-derived Red Band Difference (RBD) was used as a proxy of relative chlorophyll fluorescence to highlight areas with high algal biomass (Wolny et al. 2020;Jordan et al. 2021).

LPT-Bio model
Numerical particles were used to represent M. polykrikoides cells and were equipped with a virtual sensor to record varying cell densities driven by physical and biological processes. The horizontal and vertical locations of particles are (Chiu et al. 2018) where U, V, and W are velocities in Cartesian coordinates x, y, and z. R is a real random number with a zero mean and a uniform distribution between À1 and 1. K x , K y , and K z are turbulent diffusion coefficients. W swim is the daytime (06 : 00-18 : 00) ascent and nighttime (18 : 00-06 : 00) descent velocity (Park et al. 2001). We integrated algal bloom dynamics (i.e., mixotrophic growth, respiration, temperaturemodulated mortality, and density-related aggregation settling) into the LPT-Bio model to simulate cell density variations, which are dynamically driven by environmental parameters (i.e., temperature, salinity, irradiance, and nutrients). The cell density (C t i ) recorded by particle i is controlled by where G is the mixotrophic growth rate, including phototrophic (G P o ) and heterotrophic growth (G H o ). The phototrophic growth is where G P opt is the phototrophic growth rate at optimal conditions. f T ð Þ, f S ð Þ, f I ð Þ, and f DIN ð Þ are temperature-, salinity-, irradiance-, and DIN-limited functions, respectively. Details for limitation functions are given in Qin et al. (2021).
The heterotrophic growth is where G H opt is the maximum heterotrophic growth rate at optimal conditions. The growth-limiting function for the bioavailable organic matter (DOM and POM smaller than 12 μm, Jeong et al. 2004), f OM 12 ð Þ, is adjusted based on bloom intensities (Supporting Information Table S1).
Without considering interactions between phototrophy and heterotrophy, the mixotrophic growth rate equals the heterotrophic growth rate at night and consists of phototrophy and heterotrophy during the daytime and cannot exceed the maximum growth rate at a certain temperature and salinity (Qin et al. 2021). The mixotrophic growth rate is thus The cell density loss terms in Eq. (4) include respiration ( Res, Qin et al. 2021), temperature-modulated mortality (M T , Hofmann et al. 2021), and density-related aggregation settling is the aggregation parameter; Lima and Doney 2004). The high-resolution unstructured-grid model SCHISM (Zhang et al. 2016; Supporting Information Fig. S1) was utilized to provide 3D velocities, temperature, salinity, and nearsurface irradiance to drive particle transport and cell density variations. Blooms of M. polykrikoides are usually patchy with local accumulations due to algal behavior and physical transport. The superior boundary fitting and local refinement enable unstructured-grid models to perform better in regions with complex bathymetry and shoreline geometry (Nunez et al. 2020). The present hydrodynamic model did not activate the wave module, yet wave effects, for example, wave-induced sea surface roughness and Stokes drift advection, are suggested to be important factors impacting simulated particle trajectories, especially during wind events (Mao and Xia 2020). Moreover, nutrient and turbidity data associated with algal growth were from Chesapeake Bay Program (CBP) monitoring data (https://data.chesapeakebay.net/WaterQuality; Fig. 1).

LPT-Bio model simulations
The initial particle seeding locations and total particle numbers were determined by satellite-derived Chl a (mg Chl a m À3 ), which were converted into cell density (cells m À3 ) using a conversion factor of 1.69 Â 10 À8 mg Chl a cell À1 for M. polykrikoides (Hofmann et al. 2021). Satellite-derived RBD values identify locations with high algal biomass, and pixels with RBD greater than 2.00 Â 10 À4 in the lower CB were chosen as particle-releasing regions (Supporting Information Fig. S2). Particle numbers within selected pixels were determined as the ratio of satellite-derived cell density to the initial cell density of each particle, C 0 (cells m À3 particle À1 ). C 0 = 1.00 Â 10 7 cells m À3 particle À1 was chosen to balance the computational cost and model robustness. The sensitivity of simulated Chl a to C 0 is presented in Supporting Information Figs. S3, S4. Particles were then seeded randomly inside selected pixels.
The M. polykrikoides bloom was first observed on 22 July 2020, in the lower James River and then observed outside the James River mouth 5 d later (Supporting Information Figs. S5, S6). Based on available high-quality satellite images, three dates (22 July, 28 July, and 08 August) were selected to release particles (i.e., as reinitializations), with the tracking periods ranging from 7 to 10 d. The total released particles were 76,809, 120,448, and 325,841, respectively (Supporting Information Table S1). The satellite overpass time was around 14 : 30-15 : 30 and particles were released at 15 : 00 at a depth of 0.5 m on the three dates. Each particle is implemented with DVMs with a daytime ascent rate of 30 m d À1 and a nighttime descent rate of 70 m d À1 (Sohn et al. 2011). Sensitivity simulations were performed to assess the importance of DVMs on bloom development.
To compare model results with satellite data, the simulated cell densities recorded by all particles in the same pixel (coordinates given by satellite images) were integrated together, then converted to Chl a, and saved every half an hour. A best-match strategy was applied to assess model performance given the inherent uncertainties in satellite and model configurations.
The closest model-satellite pair with respect to Chl a for each pixel was searched between 10 : 00 and 18:00. The normalized Taylor diagram (Taylor 2001) was used to evaluate the model prediction skill, in terms of correlation coefficient, normalized centered root-mean-square difference, and normalized standard deviations of the model against satellite data (Hofmann et al. 2008).

HAB development and progression
Satellite images revealed the progression of M. polykrikoides blooms from the lower James River to the lower CB and then Fig. 2. Simulations with particles released on 28 July 2020. Column 1: simulated Chl a concentrations (Chl a, depth < 1 m); Column 2: satellite-derived Chl a; Column 3: satellite-derived RBD value; Column 4: best-match time determined daily between 10 : 00 and 18 : 00; Column 5: cumulative distribution function of simulated and satellite-derived Chl a. The black lines represent satellite data. The red lines are simulations of Chl a from the best-match searching, while faded red lines are simulations saved every half-hour between 10 : 00 and 18 : 00. Column 6: histogram of simulated (best-match) and satellite-derived Chl a.
coastal Atlantic waters (Fig. 2;Supporting Information Figs. S7,S8). The bloom first appeared in the lower James River on 22 July and subsequently propagated coastward, with a density of 7610 cells mL À1 (Supporting Information Fig. S6) measured on 27 July outside the James River mouth. A marked coastward progression and a massive bloom at Virginia Beach occurred after TS Isaias (Supporting Information Fig. S1) passed the CB watershed on 04 August. Cell densities reached 8640 cells mL À1 at the southern bank of the CB outside the James River mouth, and 6990 cells mL À1 were recorded at the Virginia Beach oceanfront on 06 August following the storm (Supporting Information Fig. S6), resulting in discolored and foul-smelling waters. The dataflow mappings and satellite images indicate the spatial patchiness of the bloom, while fixed site Chl a measurements show marked temporal variations in bloom intensity (Supporting Information Fig. S5).
Although M. polykrikoides blooms occur almost annually in the lower CB and its tributaries, export to the Atlantic Ocean has only been reported in 2007 (Mulholland et al. 2009(Mulholland et al. ), 2020.

LPT-Bio model results
Evolutions of cell concentrations and coastward transport of HABs from the James River to the lower CB and coastal Atlantic waters were well-captured by the LPT-Bio model (Fig. 2), despite some biases in bloom intensity and extent; a good model performance was achieved (on 05 and 06 August) with correct reproduction of the massive bloom at Virginia Beach (Fig. 2). The dense bloom streaks at the James River mouth were also well reproduced ( Fig. 2; Supporting Information Fig. S8). The LPT-Bio model generally predicts more Fig. 3. Simulated Chl a (depth < 1 m) with different DVM speeds. Particles were released on 28 July 2020. The vertical migration speeds were annotated above each column.

Xiong et al.
Harmful algal blooms in Chesapeake Bay accurately when transport processes dominate the bloom proliferation, that is, the best performance was achieved after the passage of TS Isaias when physical processes dominated bloom transport (Supporting Information Fig. S9). The quantitative comparisons between the LPT-Bio model and satellite images are shown in Supporting Information Fig. S9, where most data points cluster within the 0.5-1.0 normalized centered rootmean-square circle and the 0.5-1.0 normalized standard deviation circle. The average model-data correlation coefficient for all comparison dates is 0.58 (Supporting Information  Table S2). In terms of simulating Chl a in the CB, the LPT-Bio model skill is comparable to other 3D biogeochemical models (Yu and Shen 2021). Since ocean color images cannot differentiate the chlorophyll signal from M. polykrikoides and other co-occurring phytoplankton groups, the simulated Chl a specific to M. polykrikoides is generally smaller than the satellite-derived Chl a (Fig. 2, Columns 5-6; Supporting Information Figs. S7, S8). Overall, the model demonstrated a reliable predictive ability over the 7-to 10-d period when the 2020 M. polykrikoides bloom initiated and was transported to the coast. Once a bloom initiation has been identified through coastal monitoring programs, the LPT-Bio model, incorporated into a hydrodynamic and water quality forecast system, can predict bloom transport trajectories, intensities, and variations.
Effects of DVMs M. polykrikoides vertically migrates with daytime ascents and nighttime descents, as do other dinoflagellates (Park et al. 2001;Jeong et al. 2015;Noh et al. 2018). DVMs are thought to optimize cell exposure to high concentrations of near-surface sunlight during the day and high concentrations of near-bottom nutrients at night where they might also avoid predators (Jeong et al. 2015). Here, we suggest that DVMs may also allow fast-swimming algae to circumvent physical washout due to coastal ocean circulation and thereby achieve high cell densities as observed in the lower CB (Fig. 3). Particles released on 28 July were used to evaluate DVMs because the best model performance was achieved during this tracking period. Without DVMs, the passive particles were transported by currents, leading to a large-scale dispersal of bloom organisms outside the bay mouth (Column 1 in Fig. 3). Under this scenario, the predicted Chl a remained low and almost no blooms developed. The model better reproduced cell abundances in surface waters after applying a daytime ascent term with a swimming speed of 30 m d À1 (Column 2 in Fig. 3). Employing both terms for daytime ascents and nighttime descents, the bloom organisms accumulated in the nearshore region and the simulated bloom aggregations better matched observations (Columns 3-5 in Fig. 3). The vertical migration speed applied here is sufficient to create the bloom intensity in shallow regions, although algal cells might adjust their swimming behaviors based on phototaxis, geotaxis, temperature, turbulence intensity, stratification, internal biochemical state, chain length, life-cycle stage, and so on (Heaney and Eppley 1981;Park et al. 2001;Erga et al. 2015;Henrichs et al. 2015;Lovecchio et al. 2019;Brosnahan et al. 2020;Shikata et al. 2020;Lim et al. 2022). More field observations and laboratory experiments are expected to better constrain the DVMs of M. polykrikoides. It is clear, however, as seen in our findings and earlier studies (Gillibrand et al. 2016;Qin et al. 2021) that vertical migrations are critical for bloom development and accumulations, especially in nearshore regions with strong flushing (Shulman et al. 2012). Fig. 4. Simulated Chl a concentrations (depth < 1 m) from sensitivity tests of mixotrophy and algal bloom dynamics with particles released on 28 July 2020, by disabling (a1-a5) heterotrophy and (b1-b5) all growth and biomass loss processes.

Necessity of incorporating mixotrophy and algal bloom dynamics into LPT models
The LPT models are useful for predicting the transport and spatial distributions of HABs (Fernandes-Salvador et al. 2021) and detecting bloom initiation sites (Zhang et al. 2020;Clark et al. 2021). In principle, the tracked particles could spread over the whole model domain if they are allowed to disperse randomly. Yet, only locations with favorable conditions can trigger bloom initiation since biological and environmental factors come into play. Bloom intensities and locations are highly related to environmental conditions, including temperature, salinity, available light, nutrient concentrations, and local transport conditions (Qin and Shen 2019).
To demonstrate the necessity of incorporating mixotrophy and algal bloom dynamics into LPT models, we first disabled the heterotrophic growth of M. polykrikoides, which can assimilate organic nutrients to likely outcompete other algal species (Mulholland et al. 2018) and found that the phototrophic growth alone cannot sustain the bloom intensities (Fig. 4). Although particles reached Virginia Beach after the storm (Fig. 4), the simulated Chl a concentrations were about 50% of the base scenario (Supporting Information Table S3), consistent with previous findings that the degree of heterotrophy significantly impacted the timing, duration, distribution, and magnitude of the M. polykrikoides bloom (Hofmann et al. 2021;Qin et al. 2021). A further decrease in bloom intensity was predicted when all growth (heterotrophy and phototrophy) and loss terms related to the bloom dynamics were disabled, that is, particles only carried the initial cell density without being modified during the tracking (Fig. 4). Therefore, fully coupling particle tracking models with algal bloom dynamics regulated by surrounding environmental factors is a must to accurately predict HABs, especially under changing climate conditions.

Conclusions
Individual-based modeling utilizing particle tracking techniques with implemented algal bloom dynamics offers a more accurate approach to simulating HABs. This type of model specifies algal behavior patterns and examines interactions between physical transport and biological processes. Here, we developed a novel LPT-Bio model with algal bloom dynamics driven by surrounding environmental parameters. The model was coupled with a high-resolution numerical hydrodynamic model and remote sensing. The satellite-derived Chl a was used to initialize and calibrate the LPT-Bio model. The physical transport and DVM parameterizations used in the model well captured the bloom intensity/patchiness, coastward expansion, and penetration into the Atlantic Ocean.
The relatively good model predictions of the 2020 M. polykrikoides bloom over the 7-10 d of its development and the full incorporation of environmental factors suggest that the present model framework may serve as the foundation for a HAB forecast system for M. polykrikoides blooms in the CB, as well as investigate responses of HABs to the changing ocean environments. Further improvements can be made in model configurations and monitoring sources, for example, the inclusion of cyst dynamics (Hofmann et al. 2021), and the availability of satellite data with a higher temporal resolution. The present modeling technique can be transferred to other motile organisms such as raphidophytes (Handy et al. 2005) or nonmotile algal species and water bodies if remote sensing (or synoptic bloom mappings) and parameterizations of cell physiological processes are available.