Effects of Wind-Driven Lateral Upwelling on Estuarine Carbonate Chemistry

Estuaries are productive ecosystems that support extensive vertebrate and invertebrate communities, but some have suffered from an accelerated pace of acidification in their bottom waters. A major challenge in the study of estuarine acidification is strong temporal and spatial variability of carbonate chemistry resulting from a wide array of physical forces such as winds, tides and river flows. Most past studies of carbonate system dynamics were limited to the along channel direction, while lateral dynamics received less attention. Recent observations in Chesapeake Bay showed strong lateral asymmetry in the partial pressure of carbon dioxide (pCO2) and air-sea CO2 flux during a single wind event, but comparable responses to different wind events has yet to be investigated. In this work, a coupled hydrodynamic-carbonate chemistry model is used to understand wind-driven variability in the estuarine carbonate system. It is found that wind-driven lateral upwelling ventilates high DIC (Dissolved Inorganic Carbon) and CO2 deep water and raises surface pCO2, thereby modifying the air-sea CO2 flux. The upwelling also advects low pH water onto the adjacent shoals and reduces the aragonite saturation state Ωarag in these shallow water environments, producing large temporal pH fluctuations and low pH events. Regime diagrams are constructed to summarize the effects of wind events on temporal pH and Ωarag fluctuations and the lateral gradients in DIC, pH, and pCO2 in the estuary. This modeling study provides a mechanistic explanation for the observed wind-driven lateral variability in DIC and pCO2 and reproduces large pH and Ωarag fluctuations that could be driven by physical forcing. Given that current and historic mainstem Bay oyster beds are located in shallow shoals affected by this upwelling, a large fraction of the oyster beds (100–300 km2) could be exposed to carbonate mineral under-saturated (Ω


INTRODUCTION
Riverine water typically has lower dissolved inorganic carbon (DIC) and total alkalinity (TA) values and a higher DIC/TA ratio than seawater (Salisbury et al., 2008;Huang et al., 2015;Cai et al., 2021). Because of this difference between the river and ocean end members, distributions of TA, DIC, pCO 2 (partial pressure of carbon dioxide), and pH in estuaries feature strong gradients in the along-channel direction (Cai and Wang, 1998;Borges and Gypens, 2010;Cai et al., 2011). In stratified estuaries, strong vertical gradients in DIC and pH also develop where phytoplankton photosynthesis in the surface euphotic layer consumes DIC and respiration of organic material in the bottom layer produces DIC (Feely et al., 2010;Cai et al., 2011Cai et al., , 2017. These strong horizontal and vertical gradients make estuarine carbonate chemistry susceptible to disruptions from physical forcing. Turbulent mixing generated by tides and winds may mix acidic bottom water upward, leading to large changes in the airsea CO 2 flux (Cai et al., 2017;Paerl et al., 2018). Advection by tidal and wind driven currents across the horizontal gradient may cause large short-term pH fluctuations at fixed locations, which may add further stress to organisms for coping with extreme low pH events (Ribas-Ribas et al., 2013;Saderne et al., 2013;Pacellaa et al., 2018). There is an increasing recognition for the large temporal pH variability and extreme pH values in estuarine and coastal environments (Hofmann et al., 2011;Baumann et al., 2015;Baumann and Smith, 2018;Carstensen and Duarte, 2019), but most research has focused on metabolically driven pH fluctuations such as the diel photosynthesis/respiration cycle (Saderne et al., 2013;Baumann and Smith, 2018;Pacellaa et al., 2018). Relatively few studies have addressed how physical processes affect high-frequency carbonate chemistry in estuaries.
Wind-driven lateral upwelling is a well-known physical phenomenon in estuaries. Along-channel winds drive crosschannel transports via lateral Ekman transport in the surface layer of estuaries, leading to upwelling or downwelling flows at side boundaries and clockwise or counter-clockwise circulations in cross-channel sections. These lateral flows can transport momentum (Lerczak and Geyer, 2004;Scully et al., 2009), alter stratification (Lacy et al., 2003;Li and Li, 2011), and transport sediment (Geyer et al., 2001;Chen and Sanford, 2009). Using a numerical model of Chesapeake Bay, Li and Li (2012) found that the clockwise circulation generated under up-estuary winds is much stronger than the counterclockwise circulation generated under down-estuary winds. Xie et al. (2017) analyzed longterm mooring data at a cross-channel section in Chesapeake Bay and found that the lateral circulation strength increased linearly with wind speeds under up-estuary winds but was a parabolic and weaker function of wind speeds under downestuary winds.
Wind-driven lateral circulation and upwelling can transport biologically important materials such as nutrients and oxygen between the deep channel and shallow shoals (Malone et al., 1986;Sanford et al., 1990;Reynolds-Fleming and Luettich, 2004;Wilson et al., 2008). In Chesapeake Bay, Malone et al. (1986) observed higher phytoplankton production over the flanks of the main channel relative to production over the channel and attributed it to nutrient exchanges by the lateral circulation. Episodic upwelling of hypoxic water from the deep channel to shallow shoals has been observed in several estuaries, including Mobile Bay (Loesch, 1960), the Neuse River (Eggleston et al., 2005), Long Island Sound (Wilson et al., 2008), and Chesapeake Bay (Sanford et al., 1990;Breitburg, 1992). These upwelling events have the potential to generate short-lived, but severely negative conditions for organisms residing in these otherwise favorable shallow habitats, where in extreme cases, flurries of organisms congregate near shorelines as the escape to the upwelled waters (e.g., a crab "jubilee"). In a modeling study of Chesapeake Bay, Scully (2010) showed that wind-driven lateral exchange of oxygen between well-oxygenated shallow shoals and hypoxic deep channel may be more important than direct turbulent mixing in supplying oxygen to the hypoxic deep channel.
CO 2 dynamics often mirror O 2 dynamics, since all chemical species are advected or diffused by the same physical processes (currents and mixing) while the production and consumption of DIC and O 2 are affected by common biological processes such as phytoplankton production and organic matter respiration. However, while surface-water O 2 equilibrates with the atmospheric concentration relatively quickly, surface-water pCO 2 is slow to adjust and rarely reaches equilibrium with respect to the atmospheric pCO 2 due to the buffering effect of a much greater DIC pool on aqueous CO 2 (Cai et al., , 2021. Thus the effects of short-term wind effects may be different between O 2 and CO 2 . Huang et al. (2019) made repeated measurements of temperature, salinity, and pCO 2 in surface waters at a cross-channel transect in the middle reach of Chesapeake Bay. They observed large temporal and spatial variations of pCO 2 during a northerly wind event. The air-sea CO 2 flux over the eastern one-third of this transect was 20-40% higher than the flux over the middle section whereas the CO 2 flux over the western one-third of this transect was 30% lower. Huang et al. (2019) hypothesized that northerly (down-estuary) winds generated upwelling on the eastern shore, enhancing the release of respired CO 2 from acidic subsurface water. This raises a question whether southerly (up-estuary) winds will generate higher pCO 2 on the western shoal and under what wind conditions this lateral upwelling mechanism is effective in generating the lateral asymmetry in the CO 2 flux. A related question is how this wind-driven lateral upwelling affects the overall carbonate chemistry in an estuary.
Our study site, Chesapeake Bay, is an ideal system to address these questions, given that it is a large eutrophic estuary with relatively low buffer capacity (Cai et al., 2017) and wind forcing is comparable to tidal forcing  Figure 1A). It stretches for about 320 km from the mouth of the Susquehanna River at Havre de Grace, Maryland to its seaward end at Cape Charles and Cape Henry, Virginia. It is shallow, with a mean water depth of 6.5 m. However, a deep paleochannel running in the north-south direction dominates the bathymetry in the middle reaches of the main Bay. Chesapeake Bay is a partially mixed estuary: vertical salinity differences of 2-8 stratify the water column (Carter and Pritchard, 1988). It features a twolayer circulation with speeds on the order of 0.1 m s −1 (e.g., Goodrich and Blumberg, 1991). Compared with other estuaries, tidal forcing in the Bay is relatively modest with tidal range rarely exceeding 1 m (Browne and Fisher, 1988). Winds are episodic with dominant periods of 2-7 days.
Recent observations have mapped out the distributions of DIC, TA, pCO 2 , and pH in the main stem of Chesapeake Bay (Brodeur et al., 2019;Friedman et al., 2020). DIC and TA increased from surface to bottom and from north to south over the course of 2016. Upper, mid-, and lower bay DIC and TA ranged from 1000 to 1300, 1300 to 1800, and 1700 to 1900 µmol kg −1 , respectively. The pH range was large, with the maximum value of 8.5 at the surface and the minimum as low as 7.1 in bottom water in the upper and mid-bay (Brodeur et al., 2019). Significant spatial variability of pCO 2 was observed throughout the mainstem, with net uptake of atmospheric CO 2 during each season in the middle Bay and weak seasonal outgassing of CO 2 near the outflow to the Atlantic Ocean (Chen et al., 2020;Friedman et al., 2020). High-frequency pH data from a moored sensor showed high frequency fluctuations of DIC, pH, and pCO 2 influenced by circulation and biological processes (Shadwick et al., 2019). Coupled hydrodynamic-biogeochemical-carbonate chemistry models (ROMS-RCA-CC) have been developed for Chesapeake Bay (Shen et al., 2019a). Shen et al. (2019a) focused their modeling analysis on the large-scale carbonate chemistry dynamics and validation against field observations. Furthermore, Shen et al. (2019b) and Shen et al. (2020) conducted 30-year hindcast simulations to examine ecosystem metabolism and carbon balance and investigated the anthropogenic impacts on pH and aragonite saturation state. This modeling study takes the next step to investigate how physical processes affect carbonate chemistry in this estuary, seeking to test the lateral upwelling mechanism over a range of wind conditions, in particular during the up-estuary wind condition not observed in Huang et al. (2019). It also builds upon previous studies of O 2 dynamics to further our understanding of CO 2 dynamics.

MATERIALS AND METHODS
To investigate how winds affected carbonate chemistry in Chesapeake Bay, we used coupled hydrodynamicbiogeochemical-carbonate chemistry models to conduct hindcast simulations. The coupled models have three sub-models. The hydrodynamic model, based on the Regional Ocean Modeling System (ROMS, Shchepetkin and McWilliams, 2005;Haidvogel et al., 2008), has 80 × 120 grid points (∼1 km resolution) in the horizontal direction and 20 vertical sigma layers, as shown in Figure 1B (Li et al., 2005). ROMS simulates water level, currents, temperature and salinity. The biogeochemical model is based on the Row-Column Aesop (RCA) model, and includes a watercolumn component (Isleib et al., 2007) and a sediment diagenesis component (Di Toro, 2001). RCA simulates pools of organic and inorganic nutrients, two phytoplankton groups, and dissolved oxygen concentrations (Testa et al., 2014). The carbonate chemistry (CC) model simulates DIC, TA, and mineral calcium carbonate (aragonite CaCO 3 ) (Shen et al., 2019a(Shen et al., ,b, 2020. DIC is consumed by phytoplankton growth/photosynthesis and calcium carbonate precipitation. The sources of DIC include air-sea CO 2 flux, phytoplankton respiration, oxidation of organic matter, calcium carbonate dissolution, sulfate reduction, and sedimentwater fluxes. Calcium carbonate dissolution and precipitation is the primary source/sink for TA, but the contributions of several other biogeochemical processes (e.g., nitrification and sulfate reduction) to TA are also modeled (Shen et al., 2019a). Other carbonate chemistry parameters such as pH and pCO 2 are calculated from the CC model outputs using CO2SYS program (Lewis and Wallace, 1998).
The ROMS hydrodynamic model is forced by river flows at eight major tributaries, by wind stress and heat fluxes across the sea surface, and by sea level and climatology of temperature and salinity at the open boundary (Li et al., 2005(Li et al., , 2015Cheng et al., 2013;Xie and Li, 2018). At the upstream boundary of the eight major tributaries, freshwater inflows measured at USGS gauging stations were prescribed. Tidal forcing at the open ocean boundary was decomposed into tidal constituents from TPXO7 (TOPEX/Poseidon) dataassimilative global tidal model (Egbert and Erofeeva, 2002). For non-tidal forcing at the offshore boundary, de-tided daily sea-level observations acquired at NOAA Duck station were used. Air-sea momentum and heat fluxes were computed by applying standard bulk formulae (Fairall et al., 2003) to NARR (North America Regional Reanalysis from National Center for Environmental Prediction) products.
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 carbon, phosphorus, and nitrogen, and inorganic nutrients (ammonium, nitrate, phosphate) were obtained or derived from Chesapeake Bay Program biweekly monitoring data 1 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 was not considered.
The CC carbonate chemistry model is forced by the atmospheric CO 2 , the riverine loads and offshore concentration of TA and DIC. Time series of TA measurements in riverine inputs were obtained from the USGS station 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). The computed DIC is in very good agreement with the observed DIC when direct DIC measurements were made during the field cruises of 2016 (Supplementary Figure S1). 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 and temperature at the ocean boundary (Cai et al., 2010). DIC at the offshore boundary was calculated with the available TA, f CO 2 from SOCAT (Bakker et al., 2016), salinity and temperature using CO2SYS. The atmosphere pCO 2 was set to be 400 ppm. Regional ocean modeling system was initialized with outputs from a hindcast simulation from 1 January to 31 December 2012. The initial conditions of RCA were based on Chesapeake Bay Program monitoring data in December 2012. Initial conditions for DIC and TA were calculated from the two-end member mixing model. Model hindcast simulation for this study lasted from 1 January to 31 December 2013.
The ROMS-RCA-CC models have been validated against a wide variety of observational data (e.g., Li et al., 2005Li et al., , 2006Li et al., , 2016Zhong and Li, 2006;Testa et al., 2014;Xie and Li, 2018;Shen et al., 2019aShen et al., ,b, 2020Ni et al., 2020). The ROMS hydrodynamic model has been validated against water level measurements at tidal gauge stations 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). In particular, Xie and Li (2018) showed that ROMS reproduced the wind-driven lateral circulation patterns, including the current field and salinity distribution at a cross-channel section in the mid-Bay. The RCA biogeochemical model has been validated against biogeochemical data at a number of stations in Chesapeake Bay (including NO 3 , PO 4 , NH 4 , chlorophyll-a, dissolved oxygen, and organic carbon, nitrogen, and phosphorus), integrated metrics of hypoxic volume, rates of water-column primary production and respiration, and nutrient fluxes across the sediment-water surface Testa et al., 2013Testa et al., , 2014Testa et al., , 2017Li et al., 2016;Ni et al., 2020). 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  measurements of pH at a number of monitoring stations (Shen et al., 2019b(Shen et al., , 2020. The model-predicted along-channel distribution of airsea CO 2 flux is also in good agreement with that calculated from the observational data (Shen et al., 2019a). This paper uses the well-validated coupled hydrodynamic-biogeochemicalcarbonate chemistry models to investigate how winds affect carbonate chemistry in this estuary.

RESULTS
We examined the effects of up-estuary and down-estuary winds on lateral upwelling and carbonate system variability. Detailed modeling analysis revealed that interactions between wind direction and bathymetry determined the strength of upwelling and the associated bottom-surface exchange of low-pH and high-pCO 2 water. Regime diagrams were constructed to summarize the wind effects on carbonate system in the estuary and support a discussion of the associated impacts on Maryland's oyster beds.
An example sequence of a southerly wind event followed by a northerly wind event illustrated the contrasting responses of pH to winds. Southerly (up-estuary) winds blew over Chesapeake Bay between 22.00 LST (Local Standard Time) 18 September and 6.00 LST 22 September 2013, with a maximum wind speed of ∼7.5 m s −1 (Figure 2A). This was followed by northerly (down-estuary) winds between LST 6.00 on 22 September and LST 0.   Figure 1A). The three locations are further marked as red dashed lines in Figure 4D. The solid lines are original time series including tidal fluctuations and the dashed lines are low-passed with a 34-h Butterworth filter.
deep channel, and the eastern shoal) in a mid-Bay cross section. During the up-estuary wind event, the surface/bottom pH on the western shoal decreased from 8.1/7.9 to 7.9/7.7 while pH on the eastern shoal increased from 8 to 8.2 (Figures 2B,D). During the following down-estuary wind event, the opposite trend appeared: the surface/bottom pH on the western shoal increased from 7.9/7.8 to 8.2/8.1 while pH on the eastern shoal decreased from 8.1 to 7.9. In the deep channel, the bottom pH was much lower than the surface pH and the vertical pH difference decreased from ∼0.6 to ∼0.4 during the strong winds (21-23 September) ( Figure 2C). The differences in the pH time series at the three locations highlight the short-term variations in pH and carbonate chemistry in the estuary in response to wind forcing. Thus, we conducted a series of model simulations to quantify the carbonate system response to (a) weak winds, (b) up-estuary winds, and (c) down-estuary winds.

Weak Wind
As a baseline for comparison, we plot the distributions of physical and carbonate chemistry fields during a period when wind speeds were below 2 m s −1 , e.g., 22.00 LST 3 October -10.00 LST 5 October (Figures 3, 4). The tidally averaged current displayed a typical two-layer gravitational circulation, with seaward flow in the surface layer down to ∼10 m and landward flow in the underlying bottom layer ( Figure 3A). The along-channel salinity distribution showed sloping isohalines typical of a partially mixed estuary, with the top-bottom salinity difference reaching ∼6 in the deep mid-Bay ( Figure 3B). DIC and TA showed a similar along-channel distribution (Figures 3C,D). However, DIC had a larger vertical gradient than TA since phytoplankton production consumed DIC in the surface euphotic layer and respiration of organic material produced DIC in the bottom layer. The relative impacts of primary production and respiration on TA are small as is expected. pH also showed a strong longitudinal gradient, with lower pH values in the upper Bay due to the influence of poorly buffered, high-CO 2 riverine water and higher pH values in the lower Bay due to well-buffered oceanic water ( Figure 3E). In addition, pH had a strong vertical gradient, with the top-to-bottom pH difference reaching 0.8 in the upper part of the mid-Bay. The aragonite saturation state arag showed a similar along-channel distribution as pH, falling below 1 in the upper Bay and bottom waters of the mid-Bay but reaching 2-3 in the surface waters of the mid-Bay and everywhere in the lower Bay ( Figure 3F). The surface pCO 2 exceeded atmospheric equilibrium value in the upper Bay (north of 39.2 • N), was undersaturated in the mid-Bay (37.5-39.2 • N), and reached a near-equilibrium status in the lower Bay (south of 37.5 • N) ( Figure 3G). Consequently, the air-sea CO 2 flux showed outgassing of ∼5 mmol C m −2 d −1 in the upper Bay, and ingassing of (3-4) mmol C m −2 d −1 in the mid-Bay and near-zero exchange in the lower Bay ( Figure 3H). When integrated over the main stem of the Bay, there was a net influx of −10.5×10 6 g C d −1 into the estuary.
Since deep water in the mid-Bay is most susceptible to low pH and low O 2 in the estuary (Brodeur et al., 2019;Shen et al., 2019a), we examined a representative cross-channel section in this region. Salinity was vertically stratified and tilted slightly downward on the western shore as the Coriolis force confined the outflowing brackish water to the west ( Figure 4A). DIC and TA also showed strong vertical gradients and a slight west-east tilt ( Figure 4B), with lower DIC and TA within the western surface waters. pH isolines were also slightly tilted, but the lateral gradient in pH (∼0.2) was much smaller than the vertical gradient (∼0.7; Figure 4D).

Up-Estuary Wind
Winds led to large changes in the carbonate system. The upestuary winds forced water in the shallow lower Bay to move  landward at all depths ( Figure 5A). In the deep mid-Bay, the currents moved landward in the surface layer whereas the pressure gradient due to sea level pileup at the head of the estuary drove the bottom water seaward, thereby reversing the gravitational circulation. This reversed circulation strained the salinity field toward the vertical direction, promoting the ventilation of deep water. Furthermore, wind-driven turbulence generated a surface mixed layer, with its depth reaching ∼10 m between 38.5 and 39 • N ( Figure 5B). The vertical mixing and along-channel straining also strongly modified DIC and TA fields (Figures 5C,D). Their vertical gradients were greatly reduced, with nearly a uniform distribution down to a depth of 10-15 m. In the mid and lower Bay, pH was separated into two regions in the vertical direction (8-8.2 in the surface layer and 7.4-7.6 in the bottom layer) and was uniform (7.4) at all depths in the upper Bay ( Figure 5E). The aragonite saturation state also showed a sharp contrast between the surface and bottom layers, with highly under-saturated conditions in the bottom layer of the mid-Bay. The seaward flows in the bottom layer ( Figure 5A) -which under these wind conditions were opposite of the direction of mean flows (landward) -advected the under-saturated bottom water further seaward in the mid-Bay ( Figure 5F).
The southerly winds also impacted pCO 2 and air-sea fluxes. The surface water pCO 2 increased significantly in the mid-Bay as mixing and straining during the wind event raised the surface DIC ( Figure 5G). Around 39 • N, where bathymetry deepens rapidly, ventilation of the bottom water occurred, mixing high DIC-water upward and raising the surface water pCO 2 above atmospheric equilibrium. As a result, the peak in CO 2 efflux at 39 • N rivaled the peak in the uppermost part of the Bay (>39.5 • N) (compare Figures 3H, 5H). The wind event thus drove outgassing of CO 2 over an extended region in the upper and mid-Bay (down to 38.8 • N), with a magnitude of 40 mmol C m −2 d −1 . Further south in the wide mid-Bay (between 38.8 and 37.5 • N), however, the estuary absorbed CO 2 at −20 mmol C m −2 d −1 . When integrated over the entire estuary, there was a net influx of −12.1×10 6 g C d −1 .
The up-estuary winds also drove a strong clockwise lateral circulation in the cross-channel section, with the eastward velocity in the near-surface layer reaching 0.1 m s −1 (Figure 6A). This resulted in upwelling on the western shore, with an  Frontiers in Marine Science | www.frontiersin.org upwelling velocity up to 0.5 mm s −1 . Consequently, isohalines were tilted upward on the western part of the cross-section ( Figure 6B). Wind-driven turbulent mixing also homogenized salinity in the top 5 m. DIC and TA showed a similar upward tilt on the western shore, with the upwelling of high DIC and TA bottom water onto the broad western shore (Figures 6C,D). The upwelling changed acidity on the shallow shoals, where pH on the western shore was lowered (Figure 6E), reaching minima of 7.5-7.7 during the up-estuary wind (on the sea bed).
The upwelling resulted in much higher pCO 2 value on the western shore (∼350-400 ppm) than on the eastern shore (∼250 ppm). This lateral difference in the surface-water partial pressure led to a strong lateral gradient in the CO 2 flux ( Figure 6F). The mid-Bay region has been identified as a carbon sink on an annual scale. During this up-estuary wind event, the surface pCO 2 on the western shore increased, resulting in weak CO 2 flux there. On the other hand, the CO 2 flux on the eastern part of the cross section increased to −30 mmol C m −2 d −1 , producing a strong carbon sink. No significant change (<2 mg/L in chlorophyll a) in phytoplankton biomass was seen during the wind event (Supplementary Figure S2).

Down-Estuary Wind
Following the up-estuary wind event, northerly winds blew down the estuary (Figure 2A), with the maximum wind speed of 7 m s −1 . The down-estuary winds amplified the two-layer gravitational estuarine circulation, with the landward bottom current reaching 0.2 m s −1 (Figure 7A). This strong two-layer current strained the salinity field to flatten the isohalines, while the direct wind-induced mixing homogenized salinity in the surface layer down to 10 m ( Figure 7B). The landward flow in the bottom layer brought in higher salinity oceanic water into the estuary. This intrusion also caused bottom DIC and TA to increase (Figures 7C,D). pH showed an almost two-layer structure except in the shallow upper Bay where waters were vertically well mixed. In the mid-Bay, pH was 8.0-8.2 in the surface layer down to a depth of 15 m and 7.4-7.6 in the deep water below (Figure 7E). arag showed a similar distribution as pH ( Figure 7F). Due to strong wind mixing, the surface pCO 2 exceeded the atmospheric value nearly everywhere along the center axis of the Bay (Figure 7G) with net estuarine release of CO 2 except in a couple of regions of the mid-Bay (Figure 7H). When integrated over the entire estuary and over the wind event, there was a small net influx of −4.0×10 6 g C d −1 . Because the mid-Bay is much wider than the upper Bay, the total CO 2 flux was slightly negative during this down-estuary wind event even though the upper bay was a strong source, as the broad shallow shoals of the mid-Bay were a carbon sink.
In the cross-channel section, the down-estuary winds drove a counter-clockwise lateral circulation, with upwelling on the eastern shore and downwelling on the western shore ( Figure 8A). Isohalines were tilted upward, with water in the intermediate depths (∼15 m) uplifted toward the steep eastern shoal, with the isohalines outcropping at a distance from the coastline  ( Figure 8B). Due to the steep bathymetry on the eastern flank of the deep channel, the upwelling was the most intense slightly away from the shore rather than at the coastline. Both DIC and TA were uplifted upward over the eastern shoal, with the highest values at about 1 km from the shoreline, as shown in the doming of DIC and TA isohalines (Figures 8C,D). In the surface layer down to 15 m depth, pH developed a strong lateral gradient, where pH reached a low value of 7.6-7.7 on the flank of the eastern shoal (5-10 m depth) and a high value of 8.0-8.1 on the western shoal ( Figure 8E). The ventilation of high DIC water at 1 km off the eastern shoreline caused the surface water pCO 2 there to exceed the atmospheric value. The lateral distribution of air-sea CO 2 flux showed ingassing up to 15 mmol C m −2 d −1 on the western shoal and outgassing of up to 7 mmol C m −2 d −1 on the eastern shoal ( Figure 8F).

Summary of Wind Effects
We have so far focused on the detailed analyses of two wind events. To summarize the overall effects of wind on the carbonate system, we investigated a number of weather events that moved across Chesapeake Bay during 2013. This analysis captured how the carbonate system evolved in different seasons due to changing seasonal prevailing wind directions and seasonal evolution of CO 2 uptake and respiration. Two regime diagrams were constructed to identify general relationships between winds and variations in pH, DIC and other carbonate parameters.
First we investigated wind-driven temporal changes in pH. At the three stations in the mid-Bay cross-section (the cross-section marked in Figure 1A and the three stations in this section marked in Figure 4D), the time series of surface and bottom water pH were filtered to remove tidal fluctuations and the pH change from the beginning (pH begin ) to the end (pH end ) of a wind event was calculated: δpH = pH end − pH begin .
We also calculated similar quantity for the aragonite saturation state: where arag begin and arag end are the values of arag the beginning and end of each wind event. The wind impulse, defined as is an integrated measure of wind effects over a wind event and provides a good descriptor of the lateral upwelling effect. Here τ x is the along-channel wind stress at a mid-Bay location and the integration is over the entire duration of a wind event. Figures 9A,C,D,F show that δpH and δ arag at the two shallow shoals were approximately linearly proportional to I wind , and that there was additional variability beyond that explained by I wind .
With δpH reaching up to ±0.3, pH fluctuated over a wide range during a wind event. Changes in arag were also large, reaching ±1.0. Under the down-estuary winds, δpH and δ arag were positive (increased over time) on the western shoal but negative (decreased) on the eastern shoal. Under the up-estuary winds, δpH and δ arag were negative (decreased over time) on the western shoal but positive (increased) on the eastern shoal. They varied out of phase because the lateral upwelling affected the two shoals in opposite ways. In contrast, δpH and δ arag in the bottom water of the deep channel showed small changes, indicating that neither the direct wind mixing nor the lateral upwelling had much effect there (Figures 9B,E). δpH in the surface water of the deep channel showed much larger changes over each wind event and was mostly negative, indicating wind mixing injected lower pH subsurface water upward and reduced surface pH. In some cases, it became positive, suggesting that horizontal advection and straining may have brought high pH water to the mid-Bay section. δ arag in the surface water of the deep channel showed changes compatible to δpH there ( Figure 9E). With the exception of a few wind events, δpH and δ arag in the deep channel showed much smaller fluctuations than those on the shallow shoals. The lateral differences y DIC, y pH , and y pCO 2 were calculated as the DIC, pH and pCO 2 differences between the eastern and western flanks of the estuarine channel (each covering ∼1/3 of the cross-sectional area), averaged over each wind event. Instead of looking at one cross-channel section, we averaged y DIC and y pH over the mid-Bay (marked as the light blue area in Figure 1B) and over the entire main stem of the Bay (all shaded areas in Figure 1B). Figure 10A shows how y DIC averaged over the mid-Bay varied with I w . y DICwas generally negative under the up-estuary winds (positive I w ) and positive under the down-estuary winds (negative I w ). This simply said that DIC was higher on the western shoal than on the eastern shoal under the up-estuary winds. The reverse was true under the down-estuary winds: DIC was higher on the eastern shoal than on the western shoal. Despite the scatter in the plot, y DIC generally increased with I w . The stronger the wind speed or the longer the wind duration or both, the stronger the lateral gradient in DIC. There was a bias toward the upper right corner because of the downward tilt of isohalines and DIC-isolines due to the geostrophic balance (see Figure 4). y DIC can reach ±100 µmol/kg. Because the lateral asymmetry was weaker in the narrow upper Bay and shallow lower Bay, DIC over the entire Bay was smaller. Due to the geostrophic control in the lower Bay, y DIC was biased toward the positive value even under many up-estuary wind events ( Figure 10D). The overall wind effects on the lateral asymmetry in the entire Bay were not as strong as in the mid-Bay.
pH averaged over the top 5 m in the water column also displayed the lateral asymmetry that extended to all seasons and all wind events (Figures 10B,E). In the mid-Bay, y pH > 0 for I w > 0 and y pH < 0 for I w < 0 ( Figure 10B). Notwithstanding the scatters, pH appears to be a linear function of I w for I w > 0. Therefore, pH was lower on the western shoal than on the eastern shoal under the up-estuary Frontiers in Marine Science | www.frontiersin.org winds and this lateral difference y pH increased with the wind impulse. Under the down-estuary winds, the reverse asymmetry was true: pH was lower on the eastern shoal than on the western shoal. At I w < 0 the magnitude of y pH was generally less than 0.1 (except two events) whereas y pH could reach 0.2 at I w > 0. This asymmetry in y pH is related to the asymmetric response to wind direction: the lateral circulation strength and lateral upwelling are weaker under the down-estuary winds than under the-up-estuary winds (Li and Li, 2012;Xie et al., 2017). When averaged over the entire Bay, y pH versus I w showed surprisingly tight relationship as seen in the mid-Bay ( Figure 10E). The magnitude was smaller, down to ±0.1, as opposed to ±0.2 in the mid-Bay.
The lateral upwelling also generated pCO 2 variations in the lateral direction (Figures 10C,F). Under the up-estuary winds, pCO 2 on the western shoal was larger than that on the eastern shoal, i.e., y pCO 2 < 0 and its magnitude increased with the wind impulse. Under the down-estuary winds, the reverse was generally true, with y pCO 2 > 0. Huang et al. (2019) observed y pCO 2 = 50 ppm during a down-estuary wind event with the maximum wind speed of 7 m/s and wind impulse of 4.84 × 10 3 Nsm −2 . The observed data fit well with the scatter plot in Figure 10C, indicating a general agreement between the model-predicted and observed pCO 2 lateral differences. The numerical model allowed us to explore a wide range of wind forcing conditions, and Figures 10C,F showed that winds can generate lateral pCO 2 differences in the range of ± 200 ppm in the mid-Bay and in the range of ± 100 ppm when averaged over the entire estuary.

Effects of Upwelling on Nearshore Habitats
Most of the current and historic mainstem eastern oyster (Crassostrea virginica) beds in Maryland are located on the shallow shoals in the mid-Bay between 37.9 and 39.1 • N, as shown in Figure 11A (Smith, 1997). Wind-driven lateral upwelling could expose these oysters to episodic acidic conditions, particularly during the summer months when pH and arag reach seasonal minima in bottom waters, and exposure of the eastern oyster to low pH and under-saturated waters could have potentially negative consequences for growth and calcification (Waldbusser and Salisbury, 2014;Waldbusser et al., 2015). Under the weak winds, the aragonite saturation state arag on a narrow strip along the western shoal and a broader region on the eastern shoal exceeded 1, suggesting limited influence of deep-channel bottom water on most of the historic oyster beds (Figure 11B). In contrast to the shallow shoals, bottom water was highly undersaturated, with arag around 0.5. Under the up-estuary winds, arag fell below 1 almost everywhere on the western shoal, with the minimum value dropping down to 0.6 ( Figure 11C). On the other hand, arag increased on the eastern shoal, suggesting that dissolution-favorable conditions would not occur. Under the down-estuary winds, the most intense upwelling occurred channel-ward of the eastern coastline ( Figure 8C) such that arag remained above 1 over the shallowest regions on the eastern shoal ( Figure 11D). During the same down-estuary winds, arag increased to 2-3 on the fringe along the western coastline.
The duration of exposure to undersaturated conditions is an important factor for potentially negative growth and calcification effects on bivalves. We calculated the time which the bottom water arag < 1 persisted during each wind event, using the two shallow-shoal locations in the mid-Bay (Figure 4B) as examples. During the spring and summer, the bottom water on the western shoal was exposed to undersaturated conditions for a period ranging from 5 to 60 h ( Figure 12A). The bottom water was rarely exposed to arag < 1 during the fall. Since the location on the eastern shoal is further inshore than the mostintense upwelling region, the bottom water there was exposed to fewer wind events with an extended period of undersaturated conditions ( Figure 12B).
To examine the regional spatial impacts of wind-driven lateral upwelling, we calculated the size A oyster of the sea bed areas in the mid-Bay during each wind event. We computed the area between 37.9 and 39.1 • N where depth <8 m and arag < 1 for three different seasons (spring, summer, and fall). A oyster in spring and summer was a U-shaped function of the wind impulse, implying that more areas of the oyster beds were exposed to undersaturated conditions at higher wind speed or longer wind duration on both the eastern and western shoals ( Figure 12C). About 100-300 km 2 of this area could be exposed to arag < 1 during the spring and summer seasons, which is about 25-75% of the total oyster bed area. During the fall, pH and arag were higher such that only a small area (generally less than 50 km 2 ) was affected.

DISCUSSION
This study showed that winds can drive pH fluctuations of up to 0.3 in Chesapeake Bay over a period of a few days. These pH fluctuations are comparable to or larger than the long-term pH (0.1-0.2) decline that has been observed in Chesapeake Bay over the past three decades (Waldbusser et al., 2011;Shen et al., 2020). This adds to a growing body of evidence for large pH variability in coastal and estuarine systems (Carstensen and Duarte, 2019). Hofmann et al. (2011) analyzed high resolution time series of upper ocean pH obtained using autonomous sensors, over a variety of systems from polar to tropical regions. In the four coastal and estuarine sites on the United States West Coast, the standard deviations from the monthly mean pH ranged from 0.07 to 0.13, of a similar magnitude as found here. In a seagrass habitat of Hat Island in Puget Sound, Pacellaa et al. (2018) observed the mean diel range of pH of 0.39, with a maximum observed diel range of 0.74. The diel photosynthesis/respiration cycle was a major driver of this large pH variability, but tidal advection of metabolically altered parcels and wind-driven mixing events also played important roles. Similarly, measurements in a macrophyte meadow of the Baltic Sea showed that pH can vary from 8.22 ± 0.1 in August to 7.83 ± 0.4 in September (Saderne et al., 2013). Daily variations of pCO 2 due to photosynthesis and respiration were a major cause for this pH variability, but local upwelling of elevated pCO 2 water masses with offshore winds drove the variations at a time scale of days to weeks. Thus, winds and local productivity interact to drive high frequency variations in estuaries.
Winds may affect the carbonate system in an estuary or a coastal ocean in a number of ways. Depending on the air-sea difference in pCO 2 , winds can drive a large efflux or influx of CO 2 and change surface water DIC and pH (Huang et al., 2015). Wind-driven mixing can mix high DIC, high CO 2 and low pH bottom waters to the surface layer and reduce surface pH (Cai et al., 2017;Moore-Maley et al., 2018). Wind-driven upwelling can deliver acidic and carbonate undersaturated waters onto continental shelves (Feely et al., 2008;Hauri et al., 2009). Our analysis showed that wind-driven lateral upwelling could cause larger pH fluctuations on shallow shoals than those caused by direct wind mixing in Chesapeake Bay (Figures 2, 9). This upwelling mechanism is similar to that in large-scale coastal upwelling systems (e.g., Pacific coast of North America), even though winds are weaker and more variable. The upwelling velocity in Chesapeake Bay is <0.5 mm s −1 (Figures 6A, 8A), only a fraction of the typical upwelling velocity (2 mm s −1 ) found in large coastal upwelling systems (Johnson, 1977). However, it can still deliver deep acidic water to the shallow shoals because water depths in the estuary are only a few to tens of meters. For a synoptic weather event with a mean wind speed of 5 m s −1 and a duration of 2 days, the wind impulse is 8.6 × 10 3 Nsm −2 , well within the range reported in Figures 9, 12 and typically observed in estuarine and coastal regions. Our result for the dominance of upwelling over mixing in driving pH fluctuations is consistent with a related hypoxia modeling study in which Scully (2010) showed that wind-driven lateral upwelling and ventilation of deep hypoxic water onto shallow shoals are more important than wind mixing in supplying dissolved oxygen to the bottom water. This is probably more so for CO 2 than for O 2 as the air-sea gas exchange flux is generally much smaller for CO 2 than for O 2 (Jiang et al., 2019).
The wind-driven lateral upwelling generated large lateral differences in DIC. The relationship between y DIC and I w holds better in the mid-Bay than in the entire estuary (Figure 10), as the wind-driven lateral circulation is strongest in the wide mid-Bay that features a deep channel and broad shoals. Moreover, the vertical gradient in DIC and pH there is strongest due to DIC consumption in the surface layer and DIC production in the bottom layer, such that upwelling of acidic deep water has the largest effects on the shallow shoals. In comparison, the lateral circulation is much weaker in the shallow and narrow upper Bay. Strong lateral circulation can develop in the lower Bay (Xie and Li, 2018), but the vertical DIC gradient is much smaller there than that in the mid-Bay such that the upwelling does not induce large changes in surface water DIC. When integrated over the entire Bay, the relationship between y DIC and I w (and the lateral DIC asymmetry) is less robust than in the mid-Bay. y DIC is approximately a linear function of I w under the up-estuary winds but displays more scatters under the down-estuary winds. In a numerical modeling study, Li and Li (2012) showed that at the same wind speed the lateral circulation strength is 2-3 times stronger under the up-estuary winds than under the down-estuary winds. Using velocity measurements collected at a cross-channel array of AD, Xie et al. (2017) observed that the lateral circulation on the western part of the estuarine channel was weakened considerably due to the effect of baroclinic forcing under the down-estuary winds. Therefore, the larger scatter in y DIC versus I w under the down-estuary winds is likely caused by the weaker lateral circulation. We note that the scatter seen in Figure 10 are comparable to that in the plot for the lateral circulation strength reported in Xie et al. (2017) because the wind-driven lateral circulation is not only determined by the wind speed but also affected by river flows, offshore salinity and lateral baroclinic forcing.
Unlike DIC and pH, the air-sea CO 2 flux does not demonstrate a simple dependence on the wind speed/direction or wind impulse (Supplementary Figures S3-S5). As Shen et al. (2019a) and Chen et al. (2020) showed recently, Chesapeake Bay is releasing CO 2 into the atmosphere in the upper Bay, absorbing it in the mid-Bay and near the equilibrium condition in the lower Bay. While wind-driven mixing and straining/ventilation can result in locally enhanced CO 2 flux (Figures 5H, 7H) and lateral upwelling can generate cross-channel gradients in CO 2 flux (Figures 6F, 8F), the overall wind effects on the air-sea CO 2 exchange over the entire estuary cannot be easily quantified. As discussed earlier, wind effects are most pronounced in the deep hypoxic mid-Bay where respiration of organic matter creates strong vertical gradients in DIC. Typically, the mid-Bay is a net sink of CO 2 due to strong phytoplankton production in the surface euphotic layer and the surface water pCO 2 is in the range of 200-300 ppm which is less than the atmospheric pCO 2 of 400 ppm. Wind-induced mixing and ventilation in the deep channel raises surface water pCO 2 , but the net wind effect on the CO 2 flux depended on whether the surface water pCO 2 remains below or rises above the atmospheric value. If the surface water pCO 2 < 400 ppm, the winds would pump more CO 2 into the estuary and this influx/sink increases with the wind speed. For example, the up-estuary wind event U17 caused a peak CO 2 influx of −65.3 × 10 6 mol C d −1 over the mid-Bay and −61.4 × 10 6 mol C d −1 over the entire estuary (Supplementary Figure S5). If the surface water pCO 2 > 400 ppm, however, the winds will release CO 2 into the atmosphere, turning Chesapeake Bay from a weak CO 2 sink to a source. For example, during the down-estuary wind event D8, the surface water pCO 2 increased to 672.9 ppm, such that there was a peak CO 2 efflux of 133.7 × 10 6 mol C d −1 over the mid-Bay and 340.6 × 10 6 mol C d −1 over the entire estuary (Supplementary Figure S4). Similarly, the wind-driven upwelling changes the lateral distribution of pCO 2 , but its net effect on the CO 2 flux depends on the balance between the excess and deficits in the CO 2 flux at the two shallow shoals flanking the deep channel. Overall, the complicated non-linear dependence of CO 2 flux on the wind speed and strong spatial variability make it difficult to identify a simple relationship between the CO 2 flux and wind forcing in this estuary. Although numerical model estimates of air-sea@ CO 2 flux can help overcome challenges in making similar estimates from limited and error-prone observations (Herrmann et al., 2020), assigning causation to these changes under dynamic physical forcing is difficult even with model simulations.
Lateral upwelling of oxygen poor, acidic water onto otherwise oxic and high-pH shallow shoals in Chesapeake Bay has been documented in previous studies (e.g., Malone et al., 1986;Sanford et al., 1990;Scully, 2010;Huang et al., 2019), but its potential impact on oyster beds has not been specifically investigated in this estuary. In contrast, upwelling of acidic waters and it impacts on bivalves have been clearly documented along the Pacific coast of North America (e.g., Feely et al., 2008). We found that southerly or southeasterly (up-estuary winds) during the summer led to upwelling of corrosive water onto the western shoal where oyster beds have previously been identified, and that the duration of these events in 2013 led to 1-3 day periods with water with arag < 1. Although Waldbusser et al. (2011) documented reduced biocalcification in juvenile eastern oysters as W calcite declined, they also found that net calcification could occur under moderately corrosive conditions (W calcite ∼0.7-1). However, net dissolution in the eastern oyster occurred below arag = 0.6, and studies focused on other juvenile bivalves found reduced growth and mortality under comparable conditions (e.g., Green et al., 2009). The relatively brief duration of upwelling events in Chesapeake Bay may not be sufficiently chronic to induce reduced growth, calcification, or other physiological consequences (Talmage and Gobler, 2010;Keppel et al., 2016), but repeated exposures to these conditions may eventually be detrimental.
To account for accumulated acidification stress on larval bivalves, Gimenez et al. (2018) recently defined a new metric: ocean acidification stress index for shellfish (OASIS): where thrsh is the selected acidification stress threshold and D pft is the day(s) post-fertilization (spawning day = 0). OASIS has the unit of min day −1 and provides a measure of accumulated stress by integrating only the area of the arag time series that falls below the designated threshold over the entire larval period, from fertilization to metamorphosis. Although the experimental study by Gimenez et al. (2018) focused on the Pacific oyster, it is instructive to calculate OASIS using the model results presented in Figure 12. We adopted the same conservative threshold of thrsh = 1.5. Figure 13 shows how the OASIS index at the two shallow shoal locations in the mid-Bay cross section varies with the wind impulse, and the OASIS response largely mirrors that of low arag duration and area (Figure 12). OASIS was small during the fall, indicating limited acidification stress on the oyster larvae, but fall conditions may not be as important because oysters typically spawn between June and August (Kennedy, 1996). However, OASIS increased rapidly with I w during the spring and summer, a time when eastern oyster gametes are beginning to ripen and spawning occurs, respectively (Kennedy, 1996). At the western shoal location, OASIS reached 500-1500 min day −1 in the summer, which could cause a significant reduction in oyster larval survival rate (Gimenez et al., 2018) at a time when spawning is actively occurring. Furthermore, the interaction of exposures to both corrosive and hypoxic waters during these upwelling events may also have synergistic consequences for organisms (e.g., Miller et al., 2016). Clearly, there remains much to be learned concerning the potential impacts of upwelling events on biological communities, and the model simulations presented here provide a quantitative tool for prediction of the severity, duration, and frequency of such events in response to seasonally varying winds.
In summary, we used a numerical model to quantify the impacts of along-axis wind stress on the lateral dynamics of the carbonate system in Chesapeake Bay, a large estuary with moderate tides. Although the upwelling of corrosive and acidified water in response to wind is widely recognized in coastal oceans, this study provides the first mechanistic explanation that wind-driven lateral upwelling can deliver acidified deep water to shallow shoals in estuaries, producing large temporal pH and arag fluctuations. Fluctuations in pH could reach ±0.3, depending on the combined effect of wind speed and wind duration as represented by the wind impulse. Changes in arag were also large, reaching ±1.0. Regime diagrams are constructed to summarize the effects of wind events on temporal pH and arag fluctuations and the lateral gradients in DIC, pH, and pCO 2 in the estuary. Given the potential for expanded low-pH bottom waters under warming-enhanced respiration, elevated riverine inputs, and increasing atmospheric CO 2 , Chesapeake Bay may be more vulnerable to more frequent or intense future upwelling events. Future studies of estuarine acidification need to consider the intensity and control of these events, as well as the effects of these physically driven pH fluctuations on organisms' physiological response.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: https://zenodo.org/record/3964478# .XyBNRxIh2Uk.

AUTHOR CONTRIBUTIONS
ML, W-JC, and JT conceived and designed the numerical experiments. RL conducted the numerical model simulations. CS assisted in the model configurations. RL and ML analyzed the model results. ML wrote the manuscript with editorial contributions from W-JC and JT. All authors contributed to the article and approved the submitted version.

FUNDING
We are grateful to the NOAA Ocean Acidification Program (NOAA-OAP; Awards NA15NOS4780184 and NA18NOS4780179) for the financial support.