Detailed description of SCHISM and ICM are available at https://www.schism.wiki. Only a short description is given here. SCHISM employs a flexible unstructured grid with a highly efficient semi-implicit finite-element Eulerian-Lagrangian algorithm to solve the physical (Equations 1 and 2) and transport equations (Equation 3; Zhang et al., 2016):

where ∇ is the differential operator del ( ∂ ∂ x , ∂ ∂ y ) ; η(x,y,t) is the surface elevation; h(x,y) is the bathymetry; u (z,y,z,t) is the horizontal velocity; w is the vertical velocity; F represents other forcing terms such as baroclinic gradient ( − g ρ 0 ∫ z η ∇ ρ d ζ ) , horizontal viscosity, Coriolis force, tide, and atmospheric pressure; C is the tracer concentration; ν is the vertical viscosity; κ is vertical eddy diffusivity; and F_h represents the horizontal diffusion and source/sink terms. Flexibility in space and time is the key feature in SCHISM for this study. The unstructured grid ensures that the model fits the complex shoreline in the Corsica River, and the terrain-following vertical grid guarantees fine vertical resolution without bathymetry smoothing. The semi-implicit feature enables the flexibility in the time step and model advancement during the simulation.

ICM has 36 state variables, including three phytoplankton groups; two zooplankton groups; four types of nutrients; labile and refractory dissolved organic carbon (DOC), nitrogen (DON), and phosphorus (DOP); labile (G1), refractory (G2), and inert (G3) particulate organic carbon (POC), nitrogen (PON), and phosphorus (POP); DO; chemical oxygen demand (COD); and five classes of sediments (sand, silt, clay, organic detritus, and total inorganic solids) ( Figure 2A ). The model structure of ICM is flexible, and the state variables can be turned on or off depending on the interest of each application. In this application, the zooplankton compartments, refractory dissolved organic matters, inert particulate organic matter, silicate, and sediment variables were turned off ( Figure 2B ). The model has four benthic modules: submerged aquatic vegetation (SAV), tidal wetland, benthic microalgae, and shellfish. In this application, only the benthic microalgae module was activated. Detailed kinetics and parameterization can be found in previous publications (Cerco and Noel, 2004; Tian et al., 2021; Cai et al., 2022). Briefly, the governing equation of phytoplankton (B) is:

where μ is the phytoplankton growth rate, α_m is the respiration loss term, α_p represents the predation loss and the last two terms are advection and turbulent diffusion, respectively (Equation 4). The phytoplankton growth rate (μ) is controlled by water temperature, photosynthetically active radiation (Jassby and Platt, 1976) and nutrient resource, which is parameterized using the Michaelis-Menten function:

where µ_max is the maximum growth rate; T_O is the reference temperature where phytoplankton growth rate reaches its maximum, K_T(1,2) is the coefficient determining the temperature control on phytoplankton growth with K_T(1) for temperature < T_O and K_T(2) for temperature > T_O , I is the photosynthetically active radiation (PAR), K_I is the the growth-radiation coefficient, N and P are nitrogen and phosphorus concentrations, and K_N and K_P are the half-saturation constants, respectively (Equation 5).

DO envolves most of the variables and its decription provides an overview of the ICM water quality model. DO is determined by photosynthesis production, respiration consumption, remineralization of DOC, nitrification, chemical oxygen demand (COD), aereation at the sea surface, and sediment oxygen demand (SOD) at the bottom:

where a_OC is the ratio between oxygen and carbon in phytoplankton and organic matters, p_NH is the preference of ammonium uptake by phytoplankton, a_m is the phytoplankton respiration coefficient, a_ON is the ratio between oxygen and nitrogen in nitrification, NT is nitrification, a_DOC is the DOC remineralization rate, k_OC is the half-saturation constant of DO for DOC remineralization, a_COD is the oxidation rate of COD, k_COD is the half-saturation constant for COD oxidation, DO_s is DO saturation at a given temperature and salinity, DZ_s is the thickness of the surface layer, DZ_B is bottom layer thickness (Equation 6). DO aeration is applied to the surface layer and SOD to the bottom layer. When carbon fixation is based on nitrate uptake, 30% more oxygen is released as compared to ammonium uptake (the 1.3 constant in the first term of Equation 6; Morel, 1983).

The Di Toro (2001) sediment diagenesis model was incorporated in ICM, which simulates fluxes and exchanges at the sediment-water interface. Brady et al. (2013) and Testa et al. (2013). have provided comprehensive description on the diagenesis model application in Chesapeake Bay. Basically, the diagenesis of organic matter deposited to the sediment yields sulfide, methane, and ammonium whose oxidation constitutes the sediment oxygen demand (SOD). Under anaerobic conditions, sulfide and methane can be released directly to the bottom water to support additional oxygen consumption (COD).

The benthic microalga model based on Cerco and Seitzinger (1997) is coupled with SCHISM with an unstructured grid for coastal shallow water systems. The mass balance of benthic microalgae is determined by growth, respiration, and predation loss:

where B is the benthic microalgal biomass in g C m⁻² and G, R, and P represent the growth, respiration, and predation terms, respectively (Equation 7). Benthic microalgal growth is controlled by light, temperature, and nutrient availability. Benthic microalgal self-shading is parameterized as a function of benthic microalgal biomass and its light attenuation coefficient:

where I_BA is light available to benthic microalgal photosynthesis, I_B is light at the sediment surface, k_S is the sediment solids light attenuation coefficient, and k_BI is benthic microalgal light attenuation for self-shading (Equation 8). The light–growth curve is formulated as the Jassby–Platt function (Jassby and Platt, 1976):

where the light half saturation h_BI ² is related to the ratio of the growth rate [G_Bf(T)] to the initial slope of the light growth curve α_B . f(T) is the temperature influence on benthic microalgal growth and is formulated as an exponential function:

where kb is the exponential coefficient and Tb_O is the reference temperature set at 20°C (Equation 10). Nutrient limitation on benthic microalgae is parameterized with the Michaelis–Menten function:

where h_N is the half-saturation constant and N is the available nutrient (nitrogen or phosphorus) to benthic microalgae (Equation 11). The available nutrient concentration is computed as the sum of the bottom water concentration and sediment flux divided by the bottom cell thickness. Both respiration and predation loss of benthic microalgae are also formulated as an exponential function:

where α m b ( or   α p b ) is the exponential coefficient for metabolism (or predation) (Equation 12). All the effects of benthic microalgae activities on solute constituents are added to sediment fluxes, including DO photosynthesis production and respiration consumption, nutrient uptake and respiration release, and DOC flux from benthic microalgal metabolism. On the other hand, solids production of benthic microalgae is added to the corresponding sediment components, essentially organic carbon, nitrogen, and phosphorus. Model parameter definition and values are listed in Table 1 .

The simulation domain covers the entire Corsica River ( Figure 1 ). Grid resolution is approximately 100 m at the river mouth to 20 m near the coastline, with 5,614 cells, 3,159 nodes, and five vertical sigma layers. The simulation time step was set at 120 s. The model was first calibrated with the observation of the entire year 2006 without benthic microalga simulation. This run was used as the benchmark for comparison and called the “control run.” Upon the control run, the benthic microalga simulation was activated and called the “scenario run,” the comparison of which with the control run allowed us to assess the impact of benthic microalgae on DO and nutrient fluxes at the sediment–water interface and primary production, chlorophyll concentration, and DO in the water column, and these for both near-field stations where benthic microalgae grew and far-field stations where benthic microalgae were prohibited by environmental conditions.

Short and long-wave radiation data were obtained from the Reanalysis V5 (ERA5) of the European Centre for Medium-Range Weather Forecasts (https://www.ecmwf.int/en/forecasts/dataset/ecmwf-reanalysis-v5). Air temperature, wind, precipitation, pressure, and specific humidity data were downloaded from the North American Regional Reanalysis domain (NARR; https://www.ncei.noaa.gov/products/weather-climate-models/north-american-regional). River discharge and nutrient loads were simulated by the Hydrological Simulation Program–FORTRAN (HSPF), calibrated with the USGS River Input Monitoring (RIM) stations of the Chesapeake Bay Program (Shenk et al., 2012; Shenk and Linker, 2013). Daily river discharge and nutrient loads were available at seven loading points in the Corsica River domain for the simulation year ( Figure 1 ). Open boundary conditions were based on the CH3D-ICM simulation in the entire Chesapeake Bay, calibrated against long-term monitoring data for regulatory purposes over the past 30 years (Cerco and Noel, 2019).

Twenty-one discrete sampling events were carried out from April through December by the Department of Natural Resources of Maryland (DNR-MD, USA). Four stations were occupied for measurement of water temperature, salinity, and chlorophyll during each sampling event: the tidal water head station COR005, the upper estuary station XHH3851, the mid-estuary station XHH4931 and lower estuary station XHH4916 ( Figure 1 ). Measurements of sediment oxygen demand (SOD) and sediment-water ammonium fluxes collected during previous projects were also used to validate the diagenesis model simulation (Boynton et al., 2009; Boynton et al., 2018).

Generalized Additive Model (GAM) was used to identify the major predictors and characterize the nonlinear relationships between the predictors and the dependent variable (Wood, 2004; Wood, 2006; Harding et al., 2016). GAM from the “mgcv” package in R was applied to the simulated timeseries data of benthic microalgae production and DO flux at the sediment water interface with cubic spline (Hastie and Tibshirani, 1986; Murphy et al., 2022). Spectral analysis was performed on the simulated DO flux to identify the major frequencies within the complex variations in the timeseries data (Olson, 1986; Sanford et al., 1990; Fleming et al., 2012). The contribution of each individual frequency is characterized in the periodogram and the variance of each signal is quantified by the spectrum power density, defined as the amplitude power of the signal. Taylor Diagram was used to compare the scenario run with benthic microalgae and the control run without benthic microalgae. Taylor Diagram compares between simulations and observations in terms of correlation coefficients, standard deviations of both simulation and observation and centered root mean squared errors (CRMSE) on the same diagram (Taylor, 2001; Tian et al., 2014). The correlation coefficient between simulation and observation is expressed as the angle from the y axis, the normalized standard deviation (std) of the simulated results (simulation std divided by observation std) is given by the distance from the origin, and the centered root mean squared error of the simulation is measured by the distance between the simulation point and the observation point (see illustration in the “Results” section).

There was a flushing event in the summer of 2006, with freshwater discharge reaching 27 m³ s^-1 and dissolved inorganic nitrogen (DIN) loads up to 4,200 kg N day^-1 ( Figure 3A ). River discharge and nutrient loads were also relatively elevated earlier in the year and in the fall. There were two dry periods between the flushing events (Day 60 to 170 and 210 to 240) when river discharge was below 2 m³ s^-1. Air temperature displayed a typical seasonal cycle, up to 36 °C in summer and as low as 0 °C in winter ( Figure 3B ). In addition to the seasonal cycle, higher frequency variations on the order of weeks to a month was observed. Wind also showed high frequency variability, with frequencies in the order of hours to days ( Figure 3B ). Stronger wind was observed in winter and fall than in summer. Both seasonal and diel cycles were observed in the timeseries of photosynthetically active radiation (PAR) ( Figure 3C ). PAR was about 200 W m^-2 during the day in winter and reached up to 450 W m^-2 in summer. On top of the seasonal and diel cycles, cloud coverage interrupted the continuous regular variation, leading to low radiation from time to time regardless of the season.

The control run without benthic microalgae was calibrated with the data and served as a benchmark for comparison with the scenario simulation with benthic microalgae ( Figure 4 ). The model mostly reproduced the observed features in a variety of variables, including temperature ( Figures 4A–D ), salinity ( Figures 4E–H ), DO ( Figures 4J–L ), chlorophyll ( Figures 4M–P ), SOD ( Figures 4Q–T ), and ammonium flux from the sediment ( Figures 4U–X ). Temperature was dominated by the seasonal cycle with low temperature down to 0°C in winter and high temperature up to 35°C in summer. Modeled results compared well with observations, including high-frequency variations on the order of weeks. The major feature in salinity was the freshening event in summer due to high precipitation and freshwater discharge shown in Figure 3 . Low salinity was also observed in late winter–early spring due to elevated discharge. The model reproduced these major features in both magnitude and duration. The DO simulation was dominated by a seasonal cycle with high values in spring and fall, low values in summer, and transitional during other periods of the year ( Figures 4J–L ). On top of the seasonal cycle, the model predicted diel high-frequency variability in DO, resulting from net photosynthetic production during daytime and net respiration consumption during the night. Field observation data were mostly within the range of model prediction. A similar pattern can be observed in the chlorophyll simulation ( Figures 4M–P ). The model predicted higher chlorophyll concentrations during spring and summer and lower concentrations in the fall and winter. The model also predicted diel frequency variability on top of the seasonal cycle of chlorophyll concentration, which is a phenomenon revealed in shallow water systems by electronic sensor-based continuous monitoring (Shen et al., 2008; Graziano and Jones, 2017; Duvall et al., 2022). The high chlorophyll production in summer was related to the precipitation event that brought high nutrient loads in summer. The chlorophyll data were scattered with large variations. The data tended to support high chlorophyll concentrations in summer at the two upper estuary stations ( Figures 4M, N ) but showed relatively lower chlorophyll concentrations at the lower estuary stations as compared to the model prediction ( Figures 4O,P ). These two low estuary stations were closer to the open boundary and more subject to the open boundary conditions. SOD is expressed as negative, implying DO flux from the water column to the sediment ( Figures 4Q–T ). Larger values were simulated during the summer season than the rest of the year. Variability with frequencies on the order of weeks was simulated in addition to the seasonal cycle. The SOD simulation and observation data were comparable, with data mostly within the range of model variations. Similarly, the ammonium fluxes were also comparable between the data and simulation ( Figures 4U–X ). Higher values were simulated and observed in summer than during other seasons. Overall, the control run provided a reasonable solution as compared to the observation.

The biomass of benthic microalgae stayed at a relatively elevated level of approximately 5 g C m⁻² early in the year in January and February at the tidal headwater station COR0056 but decreased considerably in March to 3 g C m⁻² ( Figure 5 ). Benthic microalgae bloomed up to 5 g C m⁻² around April, followed by a long period of limited abundance through the summer until October and November when the biomass increased up to 7 g C m⁻² and then decreased to a moderate level of approximately 5 g C m⁻² as at the beginning of the year. The annual average biomass of benthic microalgae was 3.2 g C m⁻², which is within the range of field observation. Gould and Gallagher (1990) reported abundance of benthic diatoms ranging from 2 g C m⁻² to 15 g C m⁻² in the coastal region of Massachusetts Bay. A DNR survey of the Maryland Atlantic coastal regions reported active benthic chlorophyll abundance ranging from 24 mg chlorophyll m⁻² to 52 mg chlorophyll m⁻² (DNR, 2016). Assuming a C:chlorophyll ratio of 50 (19–60; Gould and Gallagher, 1990), the benthic microalgal biomass would range from 1.2 g C m⁻² to 2.6 g C m⁻². Cahoon and Safi (2002) reported active benthic chlorophyll of up to 250 mg m⁻². Benthic microalgal primary production was generally lower during the summer months as compared to that in spring and fall based on remote sensing data and spectral information (Méléder et al., 2020; Jacobs et al., 2021).

Benthic microalgae extended to 2 m deep along the coast, but higher abundance was mostly within the 1-m isobath of bathymetry ( Figure 6 ). Significant abundance of benthic microalgae was simulated in all the tidal headwaters and in most of the coastal regions. Among the four observation stations, only the tidal headwater station COR0056 is located within the benthic microalga zone and all other stations are outside the areas with significant benthic microalga abundance. Consequently, station COR0056 is considered as a near-field station and other stations as far-field stations. Station XHH4931, the mid-estuary station, is in an area where benthic microalgae were particularly scarce so that it is used as an example of far-field stations in the following sections.

High-frequency variability of DO flux with large amplitudes was simulated in the tidal headwater station COR0056 under the influence of benthic microalgae ( Figure 7A , green line). The amplitudes reached up to 10 g O m⁻² day⁻¹, and the frequencies were within diel cycles. The control run without benthic microalgae did not generate high-frequency variability ( Figure 7A , blue line). The DO flux variability in the simulation without benthic microalgae was mostly within 2 g O m⁻² day⁻¹, and the frequencies were on the order of weeks on top of the seasonal cycle. Both simulations predicted larger DO flux in summer than the rest of the year (negative DO flux indicating DO flux from the water column to the sediment). At the far-field station XHH4931, the model did not predict high-frequency variations in DO flux ( Figure 7B ). However, the simulation with benthic microalgae ( Figure 7B , green line) predicted DO flux significantly larger than the control run without benthic microalgae during the summer season ( Figure 7B , blue line). Also, the simulation with benthic microalgae predicted DO flux variations on the order of weeks, whereas these variations were absent in the control run without benthic microalgae.

Patterns similar to that of the DO flux can be observed in the nutrient flux predictions ( Figures 7C–H ). High-frequency variability was predicted in the phosphate simulation in the scenario run with benthic microalgae, whereas the control run without benthic microalgae did not generate similar high-frequency variability ( Figure 7C ). In contrast to the DO flux for which the high-frequency variations were skewed to the positive side (indicating fluxes from the sediment to the water column), the high-frequency variations in phosphate simulation were more skewed to the negative side, indicating phosphate fluxes from the water column to the sediment. In the control run without benthic microalgae ( Figure 7C , blue line), the phosphate fluxes were mostly positive, i.e., phosphate release from the sediment to the water column. The seasonal cycle remained in both runs with and without benthic microalgae, with higher phosphate release from sediment in summer than in the rest of the year. However, no significant difference was predicted at the far-field station XHH3941, where the two simulations were practically identical ( Figure 7D ). Overall, phosphate release from the sediment tended to be lower at the far-field station in the mid-estuary than at the upper estuary.

The high-frequency variability of nitrate fluxes at the sediment–water interface was mostly negative, indicating absorption of nitrate from the water column to the sediment ( Figure 7E ). Only during summer was the nitrate flux positive, i.e., release from the sediment. It is interesting to note that during a short period of time around day 180, overall negative nitrate fluxes were simulated in both scenarios with and without benthic microalgae. Similar events occurred during other periods of time of the year but with shorter durations and lower amplitudes. The high-frequency variations in the ammonium fluxes were also skewed on the negative side, but significant positive fluxes were simulated during the summer season in the scenario with benthic microalgae ( Figure 7G , green line). On the other hand, ammonium fluxes were mostly positive in the control run without benthic microalgae ( Figure 7G , blue line). At the far-field station XHH4931, the two simulations were practically identical for nutrient fluxes at the sediment–water interface, meaning that benthic microalgae did not have notable influence on nutrient fluxes at the far-field stations ( Figures 7D, F, H ).

Benthic microalgae had a significant impact on chlorophyll concentration in the water column at the near-field station COR0056 ( Figures 8 , 9 ). The Taylor diagram shows an overall improvement in model–data comparison of DO and chlorophyll in the water column ( Figure 8 ). DO simulation in the scenario run with benthic microalgae (red dots and squares) are closer to the data point in the Taylor diagram as compared to the control run without benthic microalgae (blue dots and squares), indicating higher correlation coefficients and smaller root mean square errors with the observation. Most of the correlation coefficients between DO simulation and observation are higher than 0.85, and the centered root mean square errors are smaller than 0.5. The chlorophyll comparison was not as good as the DO simulation, with most of the correlation coefficients lower than 0.85. It is a challenge to compare discrete sampling data to a time series of data with high-frequency variability in which the timing of the variations can considerably degrade the comparison in terms of correlation coefficient and root mean square errors. Moreover, suspension of benthic microalgae can alter the chlorophyll concentration in the water column and affect the model–data comparison that the model does not have the parameterization at the current stage of model development.

The monthly average showed that chlorophyll concentration was lower in the scenario run with benthic microalgae than in the control run without benthic microalgae over all the months, but differences were particularly higher in spring and summer than during the rest of the year ( Figure 9A ). Chlorophyll concentration in the benthic microalga simulation was only 63% of that in the control run in April (38 versus 60 μg L⁻¹). But the two scenarios were more similar in January and December when the difference was within 10%. Chlorophyll concentration was also lower in the scenario run with benthic microalgae than in the control run at the far-field station XHH4931, with differences smaller than that at the near-field station ( Figure 9B ). Chlorophyll concentration was 81% in the scenario run as compared to the control run in April and at the same level between the two runs in January and December. The annual average chlorophyll concentration in the scenario run with benthic microalgae was 76% of the control run at the near-field station and 90% at the far-field station.

The monthly average DO concentration was slightly higher in the scenario run with benthic microalgae than in the control run during most of the time ( Figure 9C ). However, in spring from March through June, the monthly average DO concentration was slightly lower in the scenario simulation than in the control run. As the two extrema, the monthly average DO concentration was 19% higher in the scenario simulation than in the control run in July and 4% lower in May. On an annual basis, the monthly average DO concentration was 6% higher in the scenario simulation than in the control run at the near-field station. A similar pattern can be observed at the far-field station, but with reduced differences ( Figure 9D ). The monthly average DO concentration was slightly higher in the scenario simulation than in the control run during most of the months but slightly lower in March and April. The annual average DO concentration stood at 9.2 mg L⁻¹ in the scenario simulation and 8.8 mg L⁻¹ in the control run, i.e., a difference of 4%.

Benthic microalgae had a significant impact on the primary production in the water column ( Figure 10 ). At the near-field station COR0056, phytoplankton production was systematically lower in the scenario run with benthic microalgae than in the control run throughout the year ( Figure 10A ). Phytoplankton production in the scenario run was only 80% of that in the control run in February and 94% in average. Benthic microalgal production was lower in summer than in the rest of the year at the tidal headwater station COR0056, whereas phytoplankton production was the highest in the summer. The correlation coefficient between benthic microalgae and phytoplankton production was −0.86. Benthic microalgal production was higher than the phytoplankton production in November, December, and January but lower than the latter for the rest of the year. On an annual basis, benthic microalgal production was 22% of the phytoplankton production at the near-field station COR0056. Even though phytoplankton production was lower in the scenario run than in the control run, the total primary production with phytoplankton and benthic microalgae combined was higher in the scenario run than in the control run during most of the months. Only in May through July was the total primary production lower in the scenario run with benthic microalgae than in the control run. The annual total primary production was higher by 14% in the scenario run than in the control run at the near-field station.

Benthic microalgae did not grow at the far-field station, yet the phytoplankton production in the water column was lower in the scenario run than in the control run ( Figure 10B ). The difference was the highest in April when the phytoplankton production in the scenario run was only 74% of that in the control run. Annual average phytoplankton production in the scenario simulation with benthic microalgae was lower by 11% than that in the control run without benthic microalgae. The total benthic microalgal production in the entire estuary was 40% of phytoplankton production and 28% of the total primary production with benthic microalgae and phytoplankton production combined. Phytoplankton production integrated in the entire estuary in the scenario run with benthic microalgae was 95% of that in the control run without benthic microalgae, but the total primary production with phytoplankton and benthic microalgae combined was higher by 35% than the total phytoplankton production in the control run.

The fact that benthic microalgae are mostly limited within the 2-m isobath of water depth results from the controlling factors that impact benthic microalgal growth in deeper waters. Benthic microalgal growth is determined by light, nutrients, and water temperature. During the early months of the year from January through March, temperature was the limiting factor that restricted benthic microalgal growth at both the near-field and far-field stations ( Figure 11 ). Phosphorus started to be limiting in spring, followed by nitrogen limitation in summer at the near-field station COR0056 ( Figure 11A ). Light availability was also relatively lower in summer but sufficient for benthic microalgae to grow. Phosphorus became more limiting in the fall, followed by temperature restriction in December. At the far-field station, light availability was a limiting factor from spring through fall when the light limiting factor was mostly below 0.2. Consequently, light was the primary controlling factor in determining the spatial distribution of benthic microalgae that was mostly limited to the 2-m isobath of water depth. A non-linear GAM function was fitted between benthic microalgal production and each of the four controlling factors ( Figure 12 ). Light turned out to be the dominant predictor, explaining 65% of the benthic microalgal production variance, far larger than the rest of other factors. Temperature was the second predictor with 5% variance explained, followed by phosphorus, 3%, and nitrogen, 2%. Barranguet et al. (1998) and Blackford (2002) found that nutrient limitation played a minor role in controlling benthic microalgal production, and Bowman et al. (2007) reported that benthic microalgae were often controlled by factors other than nutrient availability. It is noteworthy that none of the relationship between the benthic microalgal production and predictors is linear. Benthic microalgal production reached the peak at a light limitation factor of approximately 0.85 and gradually decreased beyond ( Figure 12A ). There was a positive relationship between benthic microalgal growth and temperature limitation factor up to 0.6, followed by a decreasing trend ( Figure 12B ). Limitations from other factors and respiration acceleration are the possible causes for the decreasing trend with temperature. The relationship between benthic microalgal production and nitrogen and phosphorus limitation factors showed similarities with light and temperature in that the benthic microalgal growth rate decreased at the high end of the nutrient limitation factors ( Figures 12C, D ). The non-linearity between benthic microalgal production and its controlling factors reflects the complexity of coastal shallow water dynamics and the interactions among them.

Where present, benthic microalgae resulted in high-frequency variability in DO fluxes. Spectral analysis of the DO flux time-series data of the control run without benthic microalgae did not resolve any high-frequency variability ( Figure 13A ). An annual cycle with seasonal variation was the only spectral signal resolved in the control run. In the spectral analysis of the DO flux time-series data of the scenario run with benthic microalgae, several high-frequency signals were resolved ( Figure 13B ). First, the diel signal with a 24-h period constituted the most prominent spectrum in the time series. This was due to benthic microalgal production during the daytime, which produced oxygen that outgassed from the sediment to the water column. Secondly, there was a signal of 12.4 h, which was coherent with the M2 tide frequency. This means that there was a physical signal in the DO flux time series. Tian et al. (2022) reported significant tidal impact on water quality in shallow water systems, and Kwon et al. (2014) observed tidal signals in benthic microalgal production, which ultimately affected DO flux. Tide influences benthic microalgal production mostly through altering environmental conditions, such as temperature, nutrient abundance, water depth, and light availability. Certain species can migrate in response to the tidal phase (Mitbavkar and Anil, 2004). The third signal had a period of 12 h, which happened twice a day. DO flux rapidly increased from sunrise until 10:00 to 11:00 am and then decreased more gradually to sunset, constituting the major component of the diurnal cycle ( Figure 13C ). Occasionally the peak of DO flux was split, which can be caused by interaction with physics and tide. Unlike the daytime peak, DO flux was mostly unchanged during the night, forming two inflection points at sunrise and sunset. This can be part of the 12-h semidiurnal spectral signal.

GAM fitting showed that benthic microalgal production and bottom water temperature explained 99% of the DO flux variance ( Figure 14A ). Benthic microalgal production alone explained 65% of the DO flux variance, which accounted for most of the high-frequency variability. The relationship between DO flux and benthic microalgal production is practically linear, with high benthic microalgal production leading to high DO flux ( Figure 14B ). However, GAM fitting showed a non-linear relationship between DO flux and water temperature ( Figure 14C ). DO flux increased with water temperature until approximately 12°C and then decreased with increasing temperature. High temperature explained most of the extremely low DO flux in summer. Acceleration in mineralization and respiration are behind the negative relationship between DO flux and temperature at the high end.

Interactions between benthic microalgae and phytoplankton constitute a major dynamic in coastal and estuarine shallow water systems (Hope et al., 2020). The model revealed significant resource competition between pelagic phytoplankton and benthic microalgae. Phytoplankton constitute a limiting factor of light reaching the benthic microalgae on the bottom. PAR was high during summer, yet light availability to benthic microalgae was low during the same period. Phytoplankton light attenuation and absorption can considerably reduce light penetration through the water column that sustains benthic microalgal photosynthesis. Yamaguchi et al (2007) found a negative relationship between chlorophyll concentration in the water column and benthic microalgal production, and Darrow (2007) reported the shading effect of phytoplankton on benthic microalgae. On the other hand, benthic microalgae constitute a competitor of phytoplankton nutrient uptake. During the benthic microalga productive season (e.g., April and October), phytoplankton production was significantly reduced, due to the competition of nutrients from benthic microalgae. Overall, nutrient availability to phytoplankton was reduced by one-third over an annual cycle at the near-field station. Laboratory experiments showed a similar phenomenon that benthic microalgae significantly reduced nutrient fluxes from the sediment (Sundback and Graneli, 1988). Field measurement also showed that benthic microalgae sequestered nutrients from being released to the water column during the productive season (Webster et al., 2002).

Shallow water systems account for only approximately 7% of the world ocean surface area but contribute up to 30% of the total primary production (Andersson and Mackenzie, 2004). Shallow waters (≤2 m) occupy 23.7% of the Chesapeake Bay surface area. Being located at the land–ocean interface, shallow water systems are the primary receiver of nutrient loads from watersheds. Understanding benthic microalgal dynamics constitutes a significant element in management decision-making for coastal ecosystem restoration and conservation. The benthic microalga model will be applied to the entire Chesapeake Bay and other tributaries in the coming years. Wave-drive resuspension and climate change can potentially exacerbate benthic algal impact, which are not parameterized in the current microalga model and need further investigation in future applications.