The ROMS hydrodynamic model was configured for Chesapeake bay and its adjacent shelf (Figure 2), consisting of 80 × 120 grid points in the horizontal direction and 20 evenly distributed vertical sigma levels in the vertical direction (Li et al., 2005). An orthogonal curvilinear coordinate system was used to follow the general orientation of the deep channel and the coastlines of the main stem of the bay. Coastal boundaries were specified as a finite-discretized grid via land/sea masking.

ROMS is forced by freshwater discharge at river heads, water levels at the open boundary, and heat and momentum flux across the sea surface. The freshwater input was prescribed for the eight major tributaries of Chesapeake bay, based on measurements at US Geological Survey gaging stations (USGS). The offshore boundary water level consists of tidal and non-tidal components. The tidal component was provided by global tidal model TPXO7 (TOPEX/POSEIDON) (Egbert and Erofeeva, 2002), and the non-tidal component was extracted from daily sea level measured at Duck, North Carolina, by the National Oceanic and Atmospheric Administration (NOAA). The air-sea heat flux and momentum flux were computed by using the North America Regional Reanalysis (NARR) data. The vertical eddy viscosity and diffusivity were parameterized using the k-kl turbulence closure scheme with the background value of 1 × 10^-6 m² s^-1, and the horizontal eddy viscosity and diffusivity were set to be constant (1 m² s^-1). The ROMS model was initialized using climatological temperature and salinity conditions and run for a spin-up period of 2 years to get the initial condition for year 2006. A detailed description of the model configuration can be found in Li et al. (2005). This hydrodynamic model was previously validated against water level measurements at tidal gauge stations (Zhong and Li, 2006; Zhong et al., 2008), salinity and temperature time series at monitoring stations (Li et al., 2005; Ni et al., 2020), salinity distributions collected during hydrographic surveys, and current measurements (Li et al., 2005; Xie and Li, 2018; Xie and Li, 2019), including the year of 2006 (Ni et al., 2020).

The RCA biogeochemical model includes a water-column component (Isleib et al., 2007) and a sediment diagenesis component (Di Toro, 2001), coupled to the ROMS hydrodynamic model in an offline mode. RCA simulates pools of organic and inorganic nutrients, two phytoplankton groups (one representing winter-spring diatoms and one representing summer dinoflagellates), and dissolved oxygen concentrations (Testa et al., 2014). The RCA biogeochemical model is forced by loads of dissolved and particulate materials from the eight major rivers. Riverine constituent concentrations for phytoplankton, silica, particulate and dissolved organic C, phosphorus (P), and nitrogen (N), and inorganic nutrients [ NH 4 + , NO 2 − + NO 3 − hereafter  NO 3 − , PO 4 3 − ] were obtained or derived from Chesapeake bay Program biweekly monitoring data as described in Testa et al. (2014). The ocean boundary concentrations were acquired from the World Ocean Atlas 2013 and Filippino et al. (2011). Atmospheric deposition of nutrients was much smaller than the riverine nutrient loading and thus not considered, following the previous studies (Ni et al., 2020; Li et al., 2020b; Zhang et al., 2021). The initial conditions of RCA were based on Chesapeake bay Program monitoring data in December 2005. The RCA model has been validated against biogeochemical data at a number of stations in Chesapeake bay (including NO 3 − , PO 4 3 − , NH 4 + , chlorophyll-a, dissolved oxygen, and organic C, N and P), integrated metrics of hypoxic volume, rates of water-column primary production and respiration, and nutrient fluxes across the sediment-water surface (Brady et al., 2013; Testa et al., 2013; Testa et al., 2014; Li et al., 2016; Testa et al., 2017; Ni et al., 2020).

The CC model simulates Dissolved Inorganic Carbon (DIC), Total Alkalinity (TA), and mineral calcium carbonate (aragonite CaCO₃) and has previously been coupled to ROMS-RCA for Chesapeake bay (Shen et al., 2019a; Shen et al., 2019b; Shen et al., 2020). DIC is consumed by phytoplankton growth/photosynthesis and calcium carbonate precipitation. The sources of DIC include air-sea CO₂ flux, phytoplankton respiration, oxidation of organic matter, calcium carbonate dissolution, sulfate reduction, and sediment water fluxes. Calcium carbonate dissolution and precipitation are the primary source/sinks for TA, but the contributions of several other biogeochemical processes (e.g., nitrification and sulfate reduction) to TA are also modeled. Other carbonate chemistry parameters such as pH and pCO₂ are calculated from the CC model outputs using the CO2SYS program (Lewis and Wallace, 1998). A detailed description of the CC model and its coupling to RCA is described in Shen et al. (2019a). The CC model is forced by the atmospheric CO₂, the riverine loads and offshore concentration of TA and DIC. Time series of TA measurements in riverine inputs were obtained from the USGS stations in the Susquehanna and Potomac Rivers (Raymond et al., 2000). The riverine DIC concentrations were calculated through CO2SYS with the available TA and pH (Shen et al., 2020). Carbonate chemistry data for the other smaller tributaries were estimated using empirical relationships as functions of freshwater discharge (Shen et al., 2019a). TA at the ocean boundary was directly estimated with the empirical equation based upon salinity at the ocean boundary (Cai et al., 2010). DIC at the offshore boundary was calculated with the available TA, fCO₂ from SOCAT (Bakker et al., 2016), salinity and temperature using CO2SYS. The atmosphere pCO₂ was set to be 400 ppm in 2006, a year chosen for historical validation, according to the observation from NOAA-ESRL (https://www.esrl.noaa.gov/gmd/ccgg/trends; Li et al., 2020a). Initial conditions for DIC and TA were calculated from the two-end member mixing model. The CC model has been validated against extensive surveys of DIC, TA and pH collected during ten cruises in 2016 (Shen et al., 2019a) and long term (1985-2015) measurements of pH at a number of monitoring stations (Shen et al., 2019b; Shen et al., 2020).

The P. minimum model is a mechanistic HAB model that has previously been embedded within RCA for Chesapeake bay (Zhang et al., 2021). A rhomboid strategy was used: that is, P. minimum is modeled individually while the two other plankton populations (winter diatoms and summer dinoflagellates) are represented by the aggregate functional classes. The model parameters for P. minimum are given in Zhang et al. (2021) and the parameters for the winter-spring diatoms and summer dinoflagellates can be found in Testa et al. (2014).

For the P. minimum model, the growth rate depends on temperature, light and nutrient concentrations, while the mortality terms include both grazing and respiration, and the model parameters have been calibrated according to published physiological experiments on P. minimum and numerical sensitivity-analysis experiments. The equation for P. minimum biomass is given by

where proro is the biomass of P. minimum measured by carbon (in unit of mgC L^-1), G is the growth rate, R_res is the respiration rate, and R_gz is the grazing rate. To include the effects of pH on P. minimum growth, G is calculated as

where G_T is the specific growth rate depending on temperature, G_par represents the effects of light availability, G_N represents the effects of nutrient limitation on growth, and G_pH represents the effects of pH on growth. The formulae for G_T, G_par and G_N can be found in Zhang et al. (2021). In Eq. (2) G_T and G_pH were assumed to be independent of each other, but some previous laboratory experiments showed that elevated CO₂ alone led to a higher growth rate but elevated CO₂ in concert with temperature increase had no significant effect on P. minimum growth (Fu et al., 2008). Equation (2) can be readily modified to consider nonlinear interactions between higher CO₂ and warming when more experimental data are available.

G_pH was estimated by fitting an empirical relation to previously published experimental data on P. minimum culture grown under different pH conditions (Hansen, 2002; Berge et al., 2010; Figure 3). In those experiments, P. minimum growth rates were obtained in a laboratory setting under a range of pH conditions while other parameters such as temperature and salinity were held at fixed values. To capture the sole effect of pH on the growth rate, G_pH was normalized by the respective mean value in each laboratory experiment such that P. minimum growth is enhanced by OA when G_pH > 1 but suppressed by it when G_pH < 1. pH values in Hansen’s (2002) experiments ranged from 7 – 10, far beyond the pH range observed in Chesapeake bay, but only experimental data in the pH range of 7.2 – 8.5 are shown in Figure 3. G_pH increased by about 10% as pH decreased from 8.5 to 7.2, representing a modest enhancement of growth rate.

The boundary conditions for P. minimum at the river heads and continental shelf were set to 0. The initial condition of P. minimum for the entire estuary were interpolated using the distribution reported in the estuary-wide surveys reported in Tyler and Seliger (1978). Model sensitivity-analysis experiments in Zhang et al. (2021) showed that the prediction of P. minimum blooms is insensitive to the initial condition as long as a small seed population exists at the beginning of the year.

The new coupled models were first used to conduct a hindcast simulation for the year 2006. The results were used to compare with the simulations by the original ROMS-RCA-Prorocentrum model which did not consider pH effects (Zhang et al., 2021).

To project the future effects of OA on P. minimum blooms, two climate projection runs were conducted by increasing the atmospheric pCO₂ to 550 ppm for the mid-21^st century and 800 ppm for the late-21^st century (IPCC, 2021) and using the corresponding DIC at the oceanic boundary to drive the CC model. The oceanic boundary DIC was calculated by assuming that the pCO₂ difference between the atmosphere and the ocean surface is the same as year 2006. The initial DIC conditions for two projection runs were calculated from the two-end member mixing model with the new oceanic end DIC. Other boundary and initial conditions as well as hydrodynamic forcings were assumed to be unchanged.

To show the effects of OA on P. minimum blooms, the time series of the daily surface pH, G_pH and P. minimum biomass concentration at four stations along the Chesapeake bay mainstem (their locations marked in Figure 2A) are presented for the year of 2006 and the model runs with or without considering pH effects are compared.

At the upper bay station CB3.1, surface pH increased from ~7.5 to ~7.9 from winter to early spring (Figure 4A). It then decreased from the middle of April and reached a minimum (~7.0) in August, before recovering to higher values during the fall. Since pH at the upper bay station was relatively low, G_pH > 1 all year, yielding ~5% amplification in the modeled growth rate (Figure 4B). Accordingly, the P. minimum concentrations in the model run incorporating pH effects were moderately higher than those without considering pH effects (Figure 4C). The peak concentration in May reached 1.18 × 10⁶ cells L^-1 in the model run with pH effects, as compared to 1.13 × 10⁶ cells L^-1 in the model run without pH effects. This represented a 5-10% increase in the biomass during the primary May bloom (Figure 4D), although larger increases (up to 20%) were seen in the fall season during which a smaller bloom developed. At station CB3.3C further downstream, the Seasonal variation of surface pH was similar to station CB3.1 except pH increased dramatically to 8.5 in June (Figure 4E). Though the surface pH was about 0.15 larger, G_pH was still above 1 during the bloom period (Figure 4F). The peak concentration increased by 0.05 × 10⁶ cells L^-1 (Figure 4G), representing a 3-5% increase in the bloom size (Figure 4H). The model-predicted P. minimum cell density is in good agreement with the observed cell density at CB 3.3C as well as at the stations in the mid-bay (CB 4.3C) and lower-bay (CB 5.2) where the monitoring data were available (Figures 4G, K, O).

At the mid bay station CB4.3C, surface pH generally increased from 7.9 to 8.2 between January and July (Figure 4I). It then dropped and reached the minimum (~ 7.5) in August and subsequently increased from September to December. The corresponding G_pH was slightly >1 during most of the time from January to May and then decreased to <1 in June and July (Figure 4J). G_pH increased to >1 in August and then decreased to slightly <1 again from September to December. Thus, surface P. minimum concentrations considering pH effects were also slightly higher than those without pH effects (Figure 4K), with the peak concentration during the spring bloom increasing by <5% (Figure 4L). At the lower bay station CB5.2, surface pH averaged around 8 and showed a weaker seasonal variation (Figure 4M). G_pH was slightly less than 1 all year except in August when pH reached the minimum (Figure 4N). There was virtually no difference in the cell density between the two model runs except in the late fall (Figures 4O, P). It was surprising to see higher P. minimum concentration at CB4.3C (4.2% increase) and CB5.2 (3.0% increase) during late fall when pH was high and G_pH < 1 (Figures 4J, N). On the other hand, P. minimum concentrations at the upper bay station CB3.1 remained elevated due to consistently low pH (Figure 4D). Seaward estuarine outflow could advect the higher biomass downstream, raising the cell density in the mid and lower bay even though the local growth rate was lower.

Focusing on the primary spring bloom period, monthly averaged surface pH in April increased from 7.7 in the upper bay to 8.0 in the mid bay while values were slightly lower, 7.9, in the lower bay (Figure 5A). G_pH was > 1 in the upper bay and < 1 in the mid-bay (Figure 5B). In the lower bay, G_pH was almost equal 1. Consequently P. minimum concentrations were moderately higher in the model run considering the pH effects (Figures 5C, D). In May, surface pH values in the mid bay and lower bay were similar to April, while the minimum surface pH in the upper bay declined to ~ 7.5 (Figure 5E). G_pH increased by 5% in the upper bay while dipping slight less than 1 in the mid-bay where pH remained above 8 (Figure 5F). The P. minimum bloom in May was mostly confined to the mid and upper bay, covering the mainstem between 38 and 39.2°N, as well as the Potomac River (Figure 5G). P. minimum concentrations increased everywhere in the bay, with the largest increase in the area between 38.7 and 39.3 °N (Figure 5H). The maximum increase in the cell density was 7 × 10⁴ cells L^-1 in the upper bay.

Along-channel distributions of pH, G_pH, and cell concentrations provide further information on the pH effects on the cell distributions (Figure 6). pH showed an along-channel gradient but also a vertical gradient, with the top-to-bottom difference reaching ~0.8 in the mid-bay (Figures 6A, E). G_pH of the bottom water was >1 due to lower pH (Figures 6B, F). In April, P. minimum concentrations were highest in the bottom waters of the lower bay as cells were advected landward by the estuarine return flow (Figure 6C). Since P. minimum growth was severely limited by light availability in the bottom waters, P. minimum concentrations did not increase despite the lower pH values (Figure 6D). In May, the bloom developed in the surface waters of the mid bay and upper bay (Figure 6G). Low pH in the upper bay enhanced G_pH (Figure 6F) and amplified the bloom size (Figure 6H). Although G_pH < 1 in the surface waters of the mid-bay, seaward advection of higher biomass from the upper bay compensated for the lower growth rate such that the P. minimum concentration in the mid-bay did not decrease (Figure 6H).

To explore the impact of OA on P. minimum blooms in the future climate, the time series of the daily surface pH, pH decline, G_pH and P. minimum concentrations were compared for the years 2006, mid-21^st century 2050 and late-21^st century 2100 at the three stations (Figure 7). At the upper bay station CB3.1, pH is projected to decrease by ~0.1 in 2050 and ~0.2 in 2100 (Figures 7A, B). Correspondingly, G_pH time series shift upwards by ~0.007 in 2050 and ~0.015 in 2100 (Figure 7C). The averaged P. minimum concentrations in April and May increased 2.7% by 2050 and 5.5% by 2100 (Figure 7D). At the mid-bay station CB4.3C, pH shows a larger reduction, decreasing by ~0.15 in 2050 and ~0.3 in 2100 (Figures 7E, F). G_pH shifts upwards and remains > 1 all year in 2100 (Figure 7G). The average bloom size from April to June increases 2.0% in 2050 and 5.0% in 2100 (Figure 7H). pH reduction at the lower bay station CB5.2 is as large as that at CB4.3C (Figures 7I, J). In the mid- and late-21st century, G_pH is projected to stay above 1 during most of the year (Figure 7K). P. minimum concentrations during the spring blooms would increase 2.7% in 2050 and 5.3% in 2100 (Figure 7L). The time series of the concentration ratio between future scenarios and year 2006 show the same trend at the three stations. During the spring blooms, the percentage increase in P. minimum concentrations reaches a maximum in April, as pH decline is relatively larger (Figures 7B, F, G). It should also be noted that though the percentage increase in P. minimum concentrations during late fall is very large, the fall blooms are much smaller than the spring blooms.

Surface distributions of monthly averaged pH, G_pH, and P. minimum concentrations in May are altered in the future scenario relative to the year 2006 (Figure 8). In 2050, monthly averaged surface pH in May ranges from 7.4 to 7.9, with an averaged decrease of ~ 0.1 as compared with year 2006 (Figure 8A). The area where G_pH is >1 would cover the whole bay (Figure 8B). Thus, P. minimum concentrations increase almost everywhere in the bay (Figure 8D). In 2100, monthly averaged surface pH in May decreases to 7.3 – 7.7 (Figure 8E). With such low pH, G_pH is even larger, reaching 1.05 in most area of the bay (Figure 8F). This results in even larger increases in P. minimum concentrations (Figure 8H).

Five stations were chosen to further illustrate the effects of OA on the peak bloom size of P. minimum (Figure 9). At the upper bay station CB3.1, the peak concentration increases 1.47% in the mid-21^st century and 2.86% in the late-21^st century as compared with the year 2006. Stations CB3.3C, CB4.1C, and CB4.3C are located in a region (38.5 – 39 °N) where the P. minimum blooms typically occur (Figure 5G). The peak P. minimum concentration at these three stations is projected to increase by 1.51%, 0.87%, 2.42% in 2050 and by 2.59%, 2.68%, 4.37% in 2100, respectively. At the station CB5.2 at the lower bay, the peak P. minimum concentration increases by 2.04% in 2050 and 3.73% in 2100.