- 1Department of Atmospheric and Oceanic Sciences, University of California, Los Angeles, Los Angeles, CA, United States
- 2Southern California Coastal Water Research Project Authority, Costa Mesa, CA, United States
- 3[C]Worthy, LLC, Boulder, CO, United States
A high-resolution numerical ocean model is used to assess ocean alkalinity enhancement (OAE) in the San Francisco Bay (SFB) estuary. A novel tracer-based approach is introduced to simulate alkalinity release and the subsequent CO2 ingassing. The model is run for 6 days and accurately reproduces observational data of currents, density, and tides. Estuarine dynamics induce mixing, advect buoyant water out of the Bay, and transport the added alkalinity from deep in the estuary to the surface of the open ocean. Over 80% of the tracer remains in the upper 15 m throughout the simulation. The estimated air-sea equilibration rate of the added alkalinity is approximately 2% per day. Alkalinity exported to the open ocean plays a disproportionately large role in increasing the CO2 ingassing rate compared to that in the estuary. This rate is relatively fast compared to open-ocean OAE studies due to the San Francisco Bay buoyant plume, which confines the released alkalinity to the surface mixed layer. While estuaries offer many benefits for OAE releases, further studies are needed to quantify their biogeochemical and ecosystem impacts.
1 Introduction
Negative-carbon-emission technologies are needed to avoid 2 °C warming by the end of the century (Gasser et al., 2015). Ocean alkalinity enhancement (OAE) is one of several marine carbon dioxide removal (mCDR) strategies to remove CO2 from the atmosphere (National Academies of Sciences, Engineering, and Medicine, 2022). The addition of alkaline compounds to the ocean shifts the dissolved inorganic carbon (DIC) speciation to bicarbonate and carbonate. This lowers the partial pressure of CO2 (pCO2), thereby enhancing the ocean's capacity to absorb atmospheric CO2 (Renforth and Henderson, 2017). However, for this chemical shift to translate into actual atmospheric CO2 removal, the added alkalinity must remain near the surface long enough to drive air-sea gas exchange.
To estimate the spatiotemporal evolution of alkalized waters on scales relevant to the experimental and operational deployment of OAE, regional modeling is needed. Modeling can support field experiments, quantify OAE effects, and simulate when and where air-sea CO2 exchange occurs (Fennel et al., 2023). Critically, identifying the residence time, location, and rate of equilibration of an alkalinity anomaly is necessary to determine its suitability for field deployment. The effectiveness of the OAE release primarily depends on whether the alkalinity addition occurs in the mixed layer or the surface ocean for CO2 ingassing (Jones et al., 2014). Several modeling studies of OAE have been conducted in large rivers (Ou et al., 2025; Mu et al., 2023), the open ocean (Burt et al., 2021; Zhou et al., 2025; Köhler et al., 2013), open and closed seas (Wang et al., 2023; Yamamoto et al., 2024; Butenschön et al., 2021), near coast (He and Tyka, 2023), a coastal channel (Anderson et al., 2025), and in large estuaries (Li et al., 2025; Khangaonkar et al., 2024).
Estuaries have been identified as suitable locations for OAE releases. They are the confluence of riverine freshwater and seawater and are highly dynamic systems. They are characterized by strong tidal influence, energetic turbulence and mixing, river inflow, and fresh-saltwater exchange flow (MacCready and Geyer, 2010). Many estuaries are shallow, ranging from tens of meters deep to hundreds of meters deep in channels. Shallow bathymetry and vigorous mixing from strong tidal currents can weaken stratification depending on the vertical salinity gradient. Alkalinity added at the bottom of an estuary (e.g., through a wastewater outfall) can readily be transported through the water column. The buoyancy of estuarine water masses may facilitate the delivery of alkalinity to the ocean surface. The fine-scale dynamics relevant to estuarine OAE can only be investigated with a high-resolution numerical model.
In this study, we present two findings: (1) we develop a model simulation of the San Francisco Bay (SFB) using the Regional Oceanic Modeling System (ROMS) suitable for investigating regional OAE, and (2) we use it to examine the dynamics and equilibration timescales for a modeled OAE release from a wastewater outfall in this buoyant, shallow estuary with strong tidal mixing at high spatial and temporal resolution.
2 Materials and methods
We use ROMS to simulate the San Francisco Bay and its offshore shelf at high resolution, assess model physical fields against observational datasets, analyze how estuarine dynamics can affect OAE, and evaluate a one-day pulse of a tracer that mimics an alkalinity release from a wastewater treatment plant in San Francisco Bay.
2.1 Study site
The SFB is a large, tidal, semi-enclosed estuary located on the United States West Coast with river inflow on the order of 103 m3 s−1 primarily from the Sacramento-San Joaquin Delta situated in the northern part of the Bay. This inflow drives the horizontal and vertical gradients in salinity (Cloern et al., 2017; Geyer, 2010) and 90% of the freshwater to the SFB (Smith, 1987). The bathymetry inside the SFB and its offshore shelf is exceptionally shallow. The depth inside the Bay ranges from 2 m to about 100 m in the Golden Gate Strait and shoals to 15 m outside the Channel because of a sediment sill. The gently sloping shelf extends approximately 50 km offshore to a depth of 150 m. Because of the shallow depth, tides are barotropic (Monismith et al., 1996) and drivers of extreme water level changes and energetic currents (Nederhoff et al., 2021). The tidal range at the SFB mouth is approximately 1.8 m (https://tidesandcurrents.noaa.gov/stationhome.html?id=9414290) and current speeds can reach 200 cm s−1 (Downing-Kunz et al., 2021). A number of studies have used ROMS to model the region just outside the San Francisco Bay (Zhou et al., 2023; Chao et al., 2017; Liu et al., 2018; Sandoval-Belmar, in revision; Fiechter and Moore, 2024), but none have ventured inside the Bay using ROMS; hence, the need for a model evaluation.
2.2 Numerical model
ROMS is a free-surface, terrain-following coordinate model with 3-D curvilinear coordinates that solves the primitive equations using split-explicit time steps (Shchepetkin and McWilliams, 2005). The surface and bottom boundary layers are parameterized with the K-profile parameterization (Large et al., 1994). This configuration is part of a set of nested simulations beginning at a Pacific-wide 6 km resolution solution forced by the GLORYS reanalysis data set (Verezemskaya et al., 2021) that includes tides and surface geopotential forcing (McWilliams et al., 2024). The ROMS solution is downscaled to a 2 km resolution grid of the western United States (Figure 1A), a 600 m resolution grid covering central and northern California, a 200 m grid of the San Francisco and Monterey Coasts, and finally to a 60 m configuration of the SFB and its offshore region (Figure 1B). The model is initialized and forced at the boundaries by the 200 m solution, with a 6-h spin-up, and is run for 6 days between June 14 and June 20, 2017. This period falls during neap tide and corresponds to the 70th percentile of riverine discharge during 2016–2024. 2017 marks the wettest year in this period, with discharge reaching 10,000 m3 s−1 in February. The grid's inland limit is a closed boundary. A new open-boundary scheme is used to more accurately transmit high-frequency motions (e.g., internal waves and internal tides) from the parent to the child grids during the downscaling process (Molemaker et al., in prep).
Figure 1. Study domain and configuration. (A) Depth and domain of the U.S. West Coast-wide dx = 2 km ROMS solution and domains with dx = 600 m, dx = 200 m, and dx = 60 m. (B) The SFB domain and depth at dx = 60 m were used in this study. The 150 m isobath is contoured in gray. Markers indicate the locations of observational data used for model evaluation. SF = San Francisco. (C) Daily volume discharge from the Sacramento-San Joaquin Delta for the month of June 2017. (D) The approximate location of the wastewater outfall diffuser is shown in pink, and the pipeline leading to it is shown in black.
The numerical grid has 2,000 × 2,000 × 100 grid points. The horizontal resolution is dx = 60 m, and the vertical resolution ranges from dz < 0.5 m inside the SFB to dz = 1.5 m at the edge of the shelf at 150 m to dz ~ 20 m at the western edge of the domain. This high spatial resolution is required to accurately capture the complex, shallow bathymetry that shapes estuarine dynamics. A short time step of 5 s is used to satisfy the Courant-Friedrichs-Lewy (CFL) condition at this resolution. The atmospheric forcing is coupled one-way from the National Oceanic and Atmospheric Administration's (NOAA) High-Resolution Rapid Refresh (HRRR) forecast model, which has a horizontal resolution of dx = 3 km and is forced at an hourly frequency (Dowell et al., 2022). The use of this high-resolution atmospheric product is important for capturing the fast topographic winds in the Golden Gate Channel (Supplementary Figure S1). The bathymetry inside SFB is generated from a 10 m horizontal resolution digital elevation model from the United States Geological Survey (USGS) (Fregoso et al., 2017) combined with offshore bathymetry from Shuttle Radar Topography Mission (SRTM) 15 version 2.5 with a spatial resolution of 15 arc sec (approximately 365 m horizontal resolution for San Francisco) (Tozer et al., 2019). Flow gauge data from USGS are used to force 15 rivers in this domain, nine of which are within the SFB, including the Sacramento and San Joaquin Rivers (Figure 1C). Model output is saved at half-hourly intervals.
A tracer, acting as an alkalinity enhancement, is discharged through the diffuser of the East Bay Municipal Utility District (EBMUD) wastewater treatment plant, adjacent to Treasure Island (Figure 1D), for one day to estimate the effect of a single pulse of OAE. The outfall discharges approximately 1.7 km offshore to the center of the SFB near Treasure Island at a depth of 14 m, based on EBMUD outfall specifications. The effluent is injected at the bottom of the domain at three grid points (180 m long) with a volume discharge of 2.19 m3 s−1 (50 million gallons per day), a temperature of 17 °C, and a salinity of 1 PSU. This input has a small upward vertical velocity equal to the volume discharge divided by the grid cell area, which is much smaller than the buoyancy flux. The discharge volume is based on data reported by the treatment plant to the U.S. Environmental Protection Agency database, while the temperature and salinity are estimated from similar treatment plants in the region. A tracer (c), with a concentration of 1 tracer mass m−3, is added to the discharge to represent a hypothetical alkalinity addition. The tracer intentionally has no particular mass units to allow the results to be scaled to any alkalinity species.
2.3 Observational comparison
The objective of the model evaluation is to verify that the simulated fields accurately represent the key physical parameters relevant to this OAE study. We validate the model for the tidal frequency, currents, and density structure in and out of the SFB. Table 1 summarizes the observational dataset used to evaluate ROMS, and Figure 1B shows the location of each source. All comparisons are made at the same location in the model as the observational data and the same time period, if available.
Within the SFB, currents are compared to National Oceanic and Atmospheric Administration (NOAA) Center for Operational Oceanographic Products and Services (CO-OPS) and Physical Oceanographic Real-Time System (PORTS) stations s09010 Oakland Outer Harbor LB3 during June 2021 and s10010 Oakland Inner Harbor LB4 during June 2024, as these were the only periods close to 2017 with a full month of June data. The NOAA tide station 9414290 (San Francisco), just inside SFB past the Golden Gate Bridge, is used to validate sea surface height. Density is compared with stations encompassed in the model domain, stations 11–36, of the Water Quality of San Francisco Bay Research and Monitoring Project run by the USGS (Schraga et al., 2018) during June 2017.
Outside the SFB, density is compared to two transects from the Applied California Current Ecosystem Studies (ACCESS), Lines 4 and 6. (https://accessoceans.org, https://data.caloos.org). The averages for June, July, August, and September from 2016 to 2018 are used because the transects did not align with the model period. Surface currents are compared with high-frequency (HF) radar data from the NOAA Integrated Ocean Observing System (IOOS). Hourly HF radar data were averaged from June 14 to 20, 2017, at a 2 km resolution. The NOAA tide station 9415020 (Point Reyes) is used to compare sea surface height anomalies.
2.4 OAE anomaly tracer formulation
The goal of this experiment is to investigate the efficacy of OAE applied at an existing outfall. Comprehensively addressing OAE requires tracking both the added alkalinity and the corresponding DIC uptake by the ocean, which are governed by air-sea gas exchange and carbonate system chemistry. Alkalinity addition induces a depression in oceanic pCO2 relative to its initial state. If pCO2ocean > pCO2atmosphere after addition, the OAE will reduce the efflux of CO2 to the atmosphere, conferring the same benefit as a change to or enhancement of an oceanic influx of CO2. If the ocean becomes undersaturated with CO2 relative to the atmosphere, there will be a CO2 flux from the atmosphere into the ocean until the surface waters are re-equilibrated.
We use a single idealized tracer, c, to represent the effects of an alkalinity addition on CO2 uptake. This tracer is added to seawater in the model, similar to alkalinity, but it also accounts for air-sea gas exchange, mimicking CO2 equilibration with the atmosphere. Effectively, the tracer represents the DIC deficit in seawater resulting from the alkalinity addition and its subsequent replenishment by CO2 ingassing. In the real ocean, the DIC deficit generated by alkalinity addition would be replenished by an air-sea flux of CO2 from the atmosphere (or reduced outgassing from the ocean) until surface water CO2 re-equilibrates with the atmosphere. In the model, the tracer outgasses to the atmosphere, with an air-sea flux equal in magnitude but opposite in direction to the alkalinity-induced CO2 ingassing. By tracking the disappearance of the tracer c in the model, we can estimate the amount of CO2 taken up by surface waters. The tracer's outgassing rate follows a modified air-sea gas exchange model that reflects the slower equilibration rate of CO2, under the assumption of a linear C system with fixed chemistry and a constant gas transfer velocity. This simplified representation does not fully account for carbonate chemistry and assumes linearity in the relationship between alkalinity addition and CO2 equilibration (i.e., that the buffer factor does not change as CO2 is taken up by seawater), which may underestimate or overestimate equilibration rates under varying physical and biogeochemical conditions.
To mimic the alkalinity addition, we impose a gas exchange velocity, kw, for the tracer, c, as 2.78 × 10−6 m s−1 (0.01 m h−1) at the surface, assuming a constant atmospheric CO2 reservoir. This value is derived from a typical, global CO2 gas transfer velocity of 5.56 × 10−5 m s−1 (0.2 m h−1), reduced by a factor of (Sarmiento and Gruber, 2006, Ch. 8). This factor accounts for the fact that only a fraction of added alkalinity directly affects pCO2 on short timescales due to carbonate system buffering, and this value is typical of oceanic water. The base gas exchange velocity of 0.2 m h−1 corresponds to the 10 m wind speed u10m = 8 m s−1, using standard parameterizations of air-sea CO2 transfer (Wanninkhof, 1992): , with a = 0.31 cm hr−1 s2 m−2, and Sc the Schmidt number dependent on sea surface temperature (Sc ~ 700). The m s−1 gas transfer velocity is implemented in the model as a surface tracer flux (sink) term, Fc [mass m−2 s−1],
calculated online in ROMS (i.e., at every time step). The fixed kw in this gas exchange parameterization disentangles the direct effect of alkalinity addition from time-varying air-sea exchange parameters (e.g., wind, temperature) and provides a replicable benchmark of OAE efficiency.
Equation 1 implicitly assumes that the atmosphere does not contain any [c], allowing the added tracer c to escape completely to the atmosphere, thus simulating the anomaly response owing to OAE. It also assumes that atmospheric CO2 does not change as a result of the air-sea flux of the tracer. Because the atmospheric CO2 concentration is fixed, when the tracer disappears, there is no atmosphere-ocean feedback from the ingassing, and thus no change in the saturation concentration of CO2. Tests using a simple box model that includes both full C chemistry and the idealized tracer c show that the tracer approach provides an accurate representation of the CO2 uptake that would be obtained with full C system calculations and of its transient dynamics for small alkalinity perturbations similar to those typically involved in OAE and up to several mEq m−3 (Supplementary material).
3 Results
3.1 Model evaluation
Overall, the ROMS model captures observed currents, tides, and density well. Figure 2 shows the model comparison to observational data inside the SFB. The mean profiles of the u component of the horizontal currents match skillfully, with root-mean-square errors (RMSE) of 0.003 and 0.005 m s−1 for LB3 and LB4 (Figure 1B, white diamonds), respectively (Figures 2A, B). LB3 has a smaller standard deviation, whereas LB4 has a larger one than ROMS, likely due to differences in the time period and spring/neap cycles between the observational data and the model.
Figure 2. Model evaluation inside San Francisco Bay. (A, B) Vertical profiles of horizontal velocity component u from buoys (white diamonds in Figure 1B) and model. The red solid line is the mean over the modeled period, and the blue solid line is the mean for the same length of time in the observation. The blue-filled region and the dashed red lines are one standard deviation for the respective color. (C) Observed (blue line, southern purple triangle in Figure 1B) and modeled (dashed red line) sea surface height plotted for the duration of the simulation in station local time. (D, E) Observed density [σt = ρ(T, S)−1, 000 kg m−3, red circles in Figure 1B] and (F) modeled as a function of station and depth during June 2017.
Figure 2C shows a time series of the sea surface height during the modeled period compared to the observational data at the tide gauge next to the Golden Gate Bridge (Figure 1B, purple triangle near SF). The sea surface anomalies closely match the data, with an RMSE of 0.061 m.
Figures 2D–F shows the density [σt = ρ(T, S)−1, 000 kg m−3] as a function of station and depth, with stations 11–15 encompassing the North Bay (San Pablo Bay), 16–25 the Central Bay, and 26–36 in South Bay (Figure 1B, red dots, Schraga et al., 2018). The two USGS transects that occurred in the same month as the simulation are shown. Generally, ROMS captures the north-to-south density gradient and the weak stratification. The predicted data is generally less dense in the North Bay by approximately 4 kg m−3, which is expected based on the river inflow during June 16 compared to June 22 (Figure 1C). The Central Bay is denser than the rest of the SFB in both the model and observations due to the flood tide bringing in ocean water and has values comparable to both transects. The South Bay in ROMS is denser than the observational data by 5 kg m−3. Because the South Bay has a residence time of several months during the summer (Walters et al., 1985), this density difference may be a result of a lack of sufficient spin-up time that would prevent meridional water mass exchange coupled with weak river inflow during the modeled period (less than 2 m3 s−1 for all rivers in the South Bay). The RMSE between the transect in Figures 2D and F is 3.05 kg m−3 and Figures 2E and F is 0.76 kg m−3. The RMSE between the average of the two transects and the predicted is 1.53 kg m−3.
Outside the SFB, the model performs quite well, although some discrepancies between the model and observations exist (Figure 3). The simulated surface currents have the same qualitative pattern as the HF radar, with the southward alongshore current that becomes offshore at the southern and western ends of the domain (Figures 3A, B). The magnitudes of the currents at the western and central parts of the domain are similar. The largest discrepancy is at the Point Reyes headland, where the observed currents exhibit a strong offshore component, suggesting the presence of an eddy or Ekman transport that the model does not capture. Because this is a simulation without data assimilation, the model is not constrained by real-time observations to accurately capture eddy locations. While the model may not replicate the exact positions of eddies, it should capture their overall statistics and energy distribution. The model also does not include wave-current interactions on the shelf that can influence currents (McWilliams, 2022). Similar to Figure 2C, the sea surface height matches well with the tide gauge at Point Reyes (Figure 1B, purple triangle at the tip of Point Reyes), with an RMSE of 0.059 m (Figure 3C).
Figure 3. Model evaluation outside San Francisco Bay. Comparison of averaged surface currents during June 14–17, 2017, from (A) hourly high-frequency radar at 2 km resolution and (B) the model at the same resolution. (C) As Figure 2C for the NOAA tide station near Point Reyes, the northern purple triangle in Figure 1B. (D–G) Cross-shore sections of density (σt = ρ(T, S)−1, 000 kg m−3) from (D) ACCESS Line 4 [blue stars in (A)] and (F) Line 6 [pink squares in (A)] compared to (E, G) modeled. Note the colorbar range is 2 kg m−3.
Figures 3D–G shows the comparison of the modeled and observed density from ACCESS Lines 4 and 6. The nearshore-to-offshore and the surface-to-deep density gradients compare reasonably well. ROMS nearshore has a lower density compared to ACCESS Line 4 as a result of a rightward gravity current running out of the Bay (Fong and Geyer, 2002) that is not present in the observations, likely due to the mismatched time period (Figures 3D, E). ACCESS Line 6 is closer to the mouth of the Bay and reflects the lower nearshore density, which the model captures (Figures 3F, G), though the extent is larger in the model than in the observations. Both of these nearshore density anomalies can be attributed to differences in river discharge between the modeled and observed periods, as the observational data did not capture the modeled period.
3.2 Estuarine dynamics and tracer transport
The effluent's positive buoyancy, the strong tidal advection and mixing typical of the site, and the shallow depth promote rapid dispersion and surface expression of the tracer. Figures 4A–D shows surface tracer at flood and ebb tide during early and late simulation periods. Inside the SFB, the tracer spreads into the North Bay during flood tide, with a small amount even reaching the narrow riverine section of the domain (Figures 4B–D). Meanwhile, the South Bay is nearly entirely devoid of tracer, highlighting the weak exchange (Walters et al., 1985). During low and ebb tide, tracer leaves the SFB through the Golden Gate Bridge Strait, where a portion of it follows a coastal gravity current that naturally tends to the right when the Rossby number is low once the plume has slowed and expanded (Fong and Geyer, 2002; Horner-Devine et al., 2009), traveling northwest toward Point Reyes (Figures 4C, D). The beginnings of this gravity current can be seen during flood tide in Figure 4B. Some tracer can escape offshore during ebb tide, but a majority of it recirculates in and out of the Bay, with the highest concentration in the southern part of the Central Bay.
Figure 4. (A–D) Maps of surface tracer during (A, D) ebb and (B, C) flood tide (A, B) 2 days after release and (C, D) 4–5 days after release, log10 transformed. The blue line in (A) demarcates the inside and outside of the SFB in this study. The orange line is the eastern boundary of the grid. (E) Horizontally integrated volume transport through the Golden Gate Bridge Strait, shown by the white line in (B), as a function of depth and time. Positive values denote transport into the Bay and negative values out.
Tides drive the large-volume transport through the Golden Gate Bridge Strait, reaching over 1,000 m3 s−1 (Figure 4E). Averaging the volume transport in time and summing in depth gives a residual of –1,432 m3 s−1, which indicates a net flux out of the Bay. This Bay efflux is expected because of the strong river discharge from the Sacramento-San Joaquin Delta. The residual is approximately equal to three days' worth of the Delta discharge (June 18–20, Figure 1C). Because this simulation takes place during neap tide, the volume fluxes shown here are likely to be smaller than during spring tide.
Ebb tide forces a strong tidal plume (Horner-Devine et al., 2009) out of the Bay, revealing a sharp front that brings low salinity and tracer with it. Figure 5 shows instantaneous surface and cross sections of u (cross-shore current), salinity, and tracer at ebb tide. The front is strikingly apparent in all three fields, both on the surface and at depth. Inside the channel, the surface velocity is 0.9 m s−1 outward, while the velocity near the bottom reaches up to 0.5 m s−1 into the Bay (Figures 5A, D). This creates intense vertical shear from the straining of the horizontal density gradient (Geyer and MacCready, 2014, Figure 2) and generates mixing within the channel, uniformly distributing the tracer (Figure 5F). This vigorous mixing keeps the water column weakly stratified in the estuary (Figure 2), allowing the buoyant wastewater plume with OAE to shoal to the surface.
Figure 5. Surface fields of (A) u, the cross-shore velocity, where positive is into the Bay and negative is out, (B) salinity, and (C) tracer during an ebb tide on June 19, 2017, 10:00 UTC. (D–F) Cross-sections of the same fields along the diagonal black line in (A–C). The purple arrow in (E) shows the liftoff point. (G) Close-up views of the edge of the front created by the ebb tide of u, and (H) salinity. Note the different colorbars and scales for (G, H).
The dynamics of liftoff allow tracers from within and deep in the estuary to propagate into the ocean's surface mixed layer. In estuarine dynamics, the liftoff point is the location near the estuary mouth where estuarine discharge separates from the bottom and forms the buoyant layer, carrying constituents into the surface layer of the open ocean. The liftoff point, the location where the bottom-attached salt front detaches from the bottom and the upper layer Froude number Fr1 = 1 (Horner-Devine et al., 2015), is clearly seen in Figure 5E. The Fr1 is defined as
where uu is the upper layer velocity, h is the upper layer depth (depth of the buoyant plume), and is the reduced gravity, where g is the gravitational acceleration, and Δρ is the difference in density relative to the ocean density ρ0. Evaluating this at the liftoff point gives Fr1 = 1.06, a supercritical value where momentum exceeds the buoyancy of the plume layer, indicating that advection and shear mixing dominate the dynamics and near-field plume mixing is occurring (Hetland, 2010). The exact Fr1 = 1 is difficult to identify because Fr1 transitions from supercritical to subcritical over a substantial portion of the near-field area (Horner-Devine et al., 2015). The liftoff point is at the seaward edge of the sill, where the depth begins to increase (purple arrow, Figure 5E). The liftoff is caused by a widening of the flow as it exits the Bay, shoaling the interface (Horner-Devine et al., 2015). Because the tracer has mixed effectively in the estuary and inside the channel, the salinity front transports the tracer with it into the open ocean, shoaling along with the salinity because of the plume expansion offshore (Figure 5F). Away from the channel, past the sill, the upper-layer Froude number becomes subcritical (Fr1 < 1), and gravitational acceleration takes over as the buoyant term becomes important.
The front acts as a gravity wave bore and continues to move offshore because of gravitational acceleration (McCabe et al., 2009). To verify that it is a gravity bore, the phase speed C of the gravity wave in a two-layer fluid is calculated from Figures 5G, H:
and compared to the velocity behind the front (Akan et al., 2018). The phase speed of the gravity wave is C = 0.5 m s−1. The u velocity is strong behind the front and has magnitudes between 0.5 and 0.6 m s−1, then abruptly goes to near zero ahead of it. The u velocity being faster than the gravity wave fits within the general understanding that this is a nonlinear gravity wave bore.
3.3 Ocean alkalinity enhancement equilibration time
The time for an OAE addition to equilibrate with the atmosphere depends on whether the tracer is in the mixed layer. Figure 6A shows the vertically integrated tracer below 15 m depth averaged over day 6 of the simulation. 15 m is chosen as the surface mixed layer based on the upper layer depth during ebb tide (Figures 5D–F). We find that more than 80% of the tracer stays in the surface mixed layer over the course of the 6-day simulation (Figure 6B), and almost all the tracer that is not in the upper 15 m is in the deep channels contoured in blue in Figure 6A. This fraction increases to 90% when evaluating the upper 20 m (Supplementary Figure S2). The fraction of tracer inside the Bay follows the tidal signal, varying by approximately 15% with each flood and ebb tide. With each ebb tide, vigorous mixing in the channels is expected to resurface the majority of the deep tracer and eventually contribute to CO2 ingassing. The portion that may equilibrate over a potentially long timescale is the small fraction that ends up below 15 m offshore. Tracer that leaves the Bay is shoaled by the sill at the mouth of the Bay through the previously described liftoff process, and the higher buoyancy of this outflow keeps it in the surface layer. These dynamics and estuarine properties elevate the fraction of seaward tracer in the mixed layer (92%) compared to inside the Bay (80%) (Figure 6B).
Figure 6. (A) Vertically integrated tracer that is below 15 m, averaged over the last day of the simulation, 6 days after release. The 15 m bathymetric contour is in blue. (B) Time series of the fraction of tracer above 15 m inside the SFB (orange), outside the SFB (blue), and over the full model domain (dashed purple).
The tracer equilibration at the surface determines the rate at which the perturbation to pCO2 caused by OAE relaxes back to the equilibrium state, or, in other words, the rate of ingassing of CO2. This tracer flux is equivalent to the net CO2 uptake (i.e., the flux of CO2 relative to the original condition). To quantify this rate, the area-integrated, normalized tracer flux , expressed as a percentage, is calculated as follows:
where A is the horizontal area of integration and ctotal is the total tracer released up to that point in the simulation. Figure 7A shows a map of tracer equilibration (percentage day−1) averaged over day 6 of the simulation (i.e., Equation 4 without the area integration). A majority of ingassing occurs within the Bay, with the greatest in the south-Central Bay, where the highest surface tracer concentration resides. A significant amount of tracer equilibration occurs at the mouth of the estuary in the shape of a round plume, reflecting the effect of ebb tide and liftoff on surface tracer transport. The coastally trapped gravity current that runs northwest out of the SFB, where the bathymetry of the coastline is shallow (Figure 6A), plays a significant role in equilibration, with tracer flux magnitudes approximately equal to those of the estuarine plume near the mouth.
Figure 7. (A) Map of tracer equilibration rate averaged over the last day of the simulation, 6 days after release. (B) Time series of the fraction of tracer inside the SFB. (C) Time series of , the normalized tracer equilibration rate integrated over the area inside the SFB (orange), outside the SFB (blue), and over the full model domain (dashed purple). The division between inside and outside the SFB is defined in Figure 4A. (D) Time series of the cumulative integral of with the colored lines as in (C). This shows the total fraction of tracer equilibration over time.
Area-integrated tracer flux reveals a steady equilibration rate with inflections based on tides. Figure 7C shows (Equation 4) as a function of time and region. Over the full model domain (dashed purple line), the tracer flux varies from 1 to 3.5% per day after the one-day pulse. The tracer equilibration rate then levels off at approximately 2% per day and remains stable over the last three days of the simulation. At the end of the 6-day simulation, the total tracer equilibration is 12% (Figure 7D).
When the rate is decoupled inside and outside the Bay, demarcated by the blue line in Figure 4A, the rate inside the Bay is nearly entirely equal to the total rate for the first two days (Figure 7C). As more tracer escapes into the open ocean, the rate outside the Bay gradually increases until it is equal to the rate inside the Bay after about 6 days. The equivalence of the two rates highlights the importance of the tracer that has escaped the estuary. After 6 days, approximately 70% of the tracer remains inside the estuary (Figure 7B), and 75% of this is in the upper 15 m (Figure 6B), with the transport regulated by the tides. The exchange flow causes an estimated 5% of tracer to leave the Bay to the open ocean per day, and approximately 1.5% per day is equilibrated with the atmosphere (Figures 7B, C). Despite the majority of the tracer remaining inside the Bay, the tracer equilibration rates outside and inside the Bay are the same (Figure 7C), revealing the disproportionate role that the tracer that has left the estuary plays in equilibration. As the tracer propagates out of the Bay, the tracer concentration and surface area for air-sea flux increase (Figure 7A), elevating in the open ocean. The estuarine dynamics discussed in the previous section shoal and horizontally spread the tracer, substantially increasing the rate of tracer flux, and by proxy, CO2 ingassing, on the shelf and in the open ocean.
4 Discussion
4.1 Estuarine dynamics drive efficient CO2 drawdown
Our study demonstrates the potential of ocean alkalinity enhancement in estuaries to draw down CO2. In particular, the San Francisco Bay has several unique features that make it a good candidate for OAE, including considerable buoyant river influx, shallow bathymetry, and strong tides. The large freshwater inflow to the North Bay drives lower overall salinity and establishes a north-south haline gradient in SFB, allowing buoyant water export into the open ocean and propelling gravity currents along the shallow coast. The shallow bathymetry and tidally driven destratification in SFB allow the alkalinity addition to readily reach the surface, where it interacts with barotropic tides. These tides promote mixing, reduce vertical stratification, and advect in-Bay materials into the open ocean. This advection is expected to be greater during a spring tide. Specifically, ebb tide causes intense vertical shear, the liftoff of buoyant water masses into the ocean, and strong low-salinity plume fronts out of the estuary. These dynamics carry added alkalinity to the surface and induce CO2 ingassing.
The DIC equilibration rate of the OAE found here is approximately 2% per day, a relatively fast rate compared to other studies that add alkalinity to the open ocean. Extrapolating this rate, the OAE addition could reach full equilibrium in approximately 50 days. However, this equilibration rate is likely to decrease over time given other simulated uptake models show a slowdown of equilibration at some point (e.g., Köhler et al., 2013; Burt et al., 2021; He and Tyka, 2023; Khangaonkar et al., 2024; Yamamoto et al., 2024; Zhou et al., 2025; Anderson et al., 2025; Guo et al., 2025). The 50-day equilibration timescale likely represents a lower bound estimate, as equilibration rates are expected to decline with time as the surface layer approaches DIC saturation. This rate may increase during periods of higher riverine flow as the net volume export from the Bay and plume area will increase (Mazzini et al., 2025), leading to greater alkalinity transport to the open ocean, where we see that the alkalinity that is exported from the estuary plays a disproportionately large role in enhancing the CO2 ingassing rate. We find that over 80% of the tracer remains in the top 15 m and 90% in the top 20 m for the entire 6-day simulation. The amount of tracer is sensitive to the selected upper layer, and the mixed layer depth on the shelf outside the Bay is heavily modulated by the tides and estuarine outflow (Supplementary Figures S3–S5). Other studies have reported equilibration times of several months to years, depending on the region (e.g., Zhou et al., 2025; Yamamoto et al., 2024; Köhler et al., 2013; He and Tyka, 2023; Burt et al., 2021), primarily due to alkalinized water masses that do not remain in contact with the surface. This region is not subject to subducting water during the model period, avoiding this issue. The buoyant water outflow from the San Francisco Bay keeps alkalinity release in the mixed layer; however, this could change if southerly (downwelling-favorable) winds occur. Buoyant water from estuaries is advantageous for OAE because it ensures the alkalinity anomaly is in the mixed layer.
4.2 Seasonal and regional oceanic processes
Identifying the optimal time for alkalinity releases is important for maximizing OAE efficiency and will be determined by regional and seasonal phenomena. Due to the high computational cost of running this model, only six days of model output were generated. In addition, the domain's spatial extent imposes a practical limitation: once the tracer reaches the domain boundaries, the analysis can become biased, and an upscaling problem arises. The simulation duration represents a trade-off between capturing the key dynamics within the domain large enough to characterize the plume dispersal and the computational resources required to resolve them adequately. While modeling a full spring-neap tidal cycle at different times of year would provide greater guidance on optimal release periods, this is beyond the scope of this study. Insights from other regional and OAE studies may help inform this. The factors that affect OAE efficiency include surface residence time, gas transfer velocity, and offshore buoyant plume transport. Whether the alkalinity addition stays at the surface will be partially modulated by upwelling. For northern and central California, García-Reyes and Largier (2012) defines upwelling season as April-June, when wind stress is strongest and during which upwelling could aid in keeping offshore OAE in the mixed layer. Meanwhile, Mazzini et al. (2025) determined that the offshore SFB plume extent is greatest in winter when river discharge is high. In contrast, surface OAE releases from (Zhou et al. 2025), (Anderson et al. 2025), and (Guo et al. 2025) found that OAE efficiency is highest in summer due to stratification keeping the OAE addition at the surface, despite higher gas transfer velocities in winter from stronger winds. However, because of the estuarine plume's low salinity (Monismith et al., 2002), the SFB outflow is likely confined to the surface layer regardless of the season (Supplementary Figure S3) and can be found offshore year-round (Mazzini et al., 2025), suggesting that seasonal stratification plays a minor role. Based on these studies, a potentially good release period may be during the spring tide in April, when strong tides enhance SFB and oceanic exchange, upwelling is active, wind stress is elevated, and riverine inflow is higher than in summer. Further investigation of suitable release times for this region is needed.
4.3 Real-world potential and impacts
To obtain a crude estimate of the amount of atmospheric CO2 removed by releasing OAE in this location, a few scaling arguments can be applied. Assume that the OAE process increases the alkalinity of the wastewater effluent to 2,500 mmol NaOH m−3 (100 mg NaOH L−1). With a discharge of 2.19 m3 s−1 (50 million gallons per day), the total alkalinity (TA) added in one day is 473,040 mol TA, or 18.9 metric tons of NaOH. Using a value of 0.78 of moles of removed CO2 to moles of added alkalinity for the San Francisco shelf (Zhou et al., 2025), and supposing that all the TA is used up at a rate of approximately 2% a day, then the total CO2 removed would be 16.2 tons over O(2 months). If the OAE deployment were continuous for a year, the removal would be estimated at 4950±1,000 tons of CO2 year−1. To put this in a global context, the IPCC target is to remove 5.5 Gt CO2 year−1 from the atmosphere to keep warming at 1.5 °C by the end of the century. The removal approximated here represents a millionth of the CO2 atmospheric extraction target and could be considered one implementation among several to address CO2 emissions. It should be noted that the actual ingassing rate and carbon dioxide removal are expected to be diminished due to a slowdown in equilibration.
One concern of OAE is how the added alkalinity could affect seawater chemistry and have downstream effects on biology, particularly near field of the release where changes in pH and Ωarag are expected to be greatest. For a wastewater release, the short-term (minute-scale) near-field chemistry depends on the initial dilution of the effluent, the distance from the diffuser, and local currents (Kitidis et al., 2024). EBMUD's National Pollutant Discharge Elimination System (NPDES) permit reports an initial dilution of 124:1 during average discharge at slack tide and allows pH changes of up to 0.5 units, provided the pH does not exceed 8.5. The alkalinity addition from the previous example of 2,500 mmol NaOH m−3 (equivalent to 2,500 mEq m−3) would, after initial dilution, yield a concentration of approximately 20 mEq m−3. Using CO2SYS (Humphreys et al., 2022) with SFB conditions, this could increase the local pH by about 0.05 to 8.06. However, under continuous release and limited circulation, the near-field pH increase may be higher as OAE accumulates (He and Tyka, 2023; Anderson et al., 2025; Wang et al., 2023). In the SFB, tidal currents can reach up to 2 m s−1 through the Golden Gate channel, promoting rapid dilution of the alkalinity addition to more than 1,000:1 (Figure 4). On the other hand, the tides can also re-entrain perturbed waters, as seen by the recirculating tracer at the mouth of the Bay (Figure 4), and slack tide can prolong exposure to perturbations for several hours. Careful monitoring, reporting, and verification are therefore necessary to avoid unintended chemical disturbances (Schulz et al., 2023; Gill et al., 2024).
4.4 OAE tracer method applications and caveats
The novel surface tracer flux used to mimic an alkalinity addition and equilibration in this study allows broad application. This method simulates the essential features of an OAE experiment without needing to implement full biogeochemistry. The results are likely to differ in detail from a fully explicit treatment but are sufficiently robust to bound the rough time and space scales of interest. Nonetheless, there are several caveats to this approach. Foremost, this study uses a physical numerical ocean model with a simplified, idealized CO2 air-sea gas exchange at the surface, which does not include biogeochemistry or the carbonate system. We use a fixed gas transfer velocity that is typical of oceanic waters. Although the model circulation is forced with a realistic atmosphere, the gas exchange parameterization uses a constant kw to provide a reproducible baseline estimate of OAE efficiency and to isolate the effect of alkalinity addition from time-varying air-sea exchange drivers. In reality, the CO2 flux will be governed by the difference between CO2air and CO2sea, the mechanical factors of wind speed (Wanninkhof, 1992, 2014) and waves (Reichl and Deike, 2020), and the seawater temperature and salinity (Weiss, 1974). Capturing these effects would require coupling the instantaneous model fields to a time-varying kw (He and Tyka, 2023; Yamamoto et al., 2024; Khangaonkar et al., 2024; Guo et al., 2025; Anderson et al., 2025). Sensitivity analyses quantifying how variability in gas exchange rates may alter the conclusions presented here are an important next step that could be incorporated in future studies.
Because of a lack of a carbonate system, we neglect the potential feedback of alkalinity addition on natural alkalinity cycling (Bach, 2024). This approach also assumes a constant atmospheric CO2 reservoir and excludes ocean-atmosphere feedback associated with CO2 perturbations (Oschlies, 2009); however, such effects are unlikely to be significant at the scales of CO2 removal considered here. Additionally, the assumption that the carbonate system behaves approximately linearly constrains the buffering effect of OAE to stay fixed even as CO2 is taken up by seawater. Furthermore, the factor that reduces the base gas transfer velocity, , should be calculated using temperature, salinity, DIC, and alkalinity values specific to the surface waters of the chosen region and will change over the course of the uptake. This factor represents the need for CO2 to equilibrate with the entire DIC pool in the surface ocean. Under typical San Francisco Bay conditions, this buffer factor is approximately 8 (see Supplementary material), whereas a value of 20, which is representative of open-ocean conditions (Sarmiento and Gruber, 2006), was used here. Because the base piston velocity is divided by this factor, applying the regional value would roughly double the effective gas transfer velocity, leading to an increase in the rate of CO2 ingassing. Box model experiments (Supplementary Text) indicate that even a factor-of-two change in this parameter alters only the transient equilibration rate, not the total potential CO2 uptake. In a realistic oceanic setting, however, slower equilibration allows surface currents and vertical mixing to redistribute the alkalinity anomaly before complete equilibration occurs, effectively reducing the local CO2 uptake rate. The results presented here should therefore be considered conservative with respect to both gas exchange efficiency and the influence of physical transport.
4.5 Uncertainties of OAE for estuarine ecosystems
The effects of alkalization on the health of estuarine and ocean ecosystems are largely unknown. The buffering effect of OAE on pH tends to benefit estuaries as ocean acidification is mitigated (Doney et al., 2009; Renforth and Henderson, 2017; Cai et al., 2021), which could have positive outcomes for shell-forming, ecologically- and commercially important species (Bednaršek et al., 2014, 2020). Conversely, too large an alkalinity perturbation can lead to significant increases in pH and CaCO3 precipitation that reduce the capacity for CO2 uptake, decrease total alkalinity, and may influence nutrient availability (Moras et al., 2022; Schulz et al., 2023; Bach, 2024; Hooper et al., 2025). For short-term exposure of OAE to calcifiers, Bednaršek et al. (2025) reported that 60% of species exhibited non-neutral responses to relatively low concentrations of NaOH (50–500 mmol m−3). In contrast, Oberlander et al. (2025) observed minimal impact on most nearshore phytoplankton under short-term, elevated pH exposure. Additionally, Goldenberg et al. (2024) found that larvae and juveniles of two Atlantic fish species were largely unaffected by both short- and long-term OAE exposure. One potential benefit is that the addition of alkalinity to wastewater effluent can counteract acidification and respiration-induced CO2 in estuaries (Liu et al., 2018; Khangaonkar et al., 2024; Li et al., 2025) and the open ocean (Cai and Jiao, 2022; Kessouri et al., 2021; Ho et al., 2023). Overall, the biological effects of OAE in estuaries remain poorly understood and uncertain, necessitating greater caution and thorough quantification before alkalization is applied.
5 Conclusion
We use ROMS to simulate a 1-day pulse of OAE released from a wastewater treatment plant in San Francisco Bay and conclude that this estuarine region could be a suitable candidate for alkalinity addition. The model evaluation demonstrates that ROMS accurately reproduces observational data on currents, density, and tides. Estuarine dynamics induce mixing, advect buoyant water out of the Bay, and transport the tracer (which mimics added alkalinity) to the surface, with over 80% of the tracer remaining in the top 15 m and over 90% in the top 20 m throughout the simulation. The estimated air-sea equilibration rate of the DIC deficit caused by the idealized alkalinity addition is 2% per day. Alkalinity exported from the estuary to the open ocean plays an outsized role in the CO2 ingassing rate compared to the in-Bay alkalinity. This rate is relatively fast compared to open-ocean OAE studies due to the buoyant plume from San Francisco Bay, which confines the released alkalinity to the surface mixed layer. In summary, estuaries offer many appealing features for OAE releases, although further studies are needed to quantify their biogeochemical and ecosystem impacts.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://doi.org/10.5281/zenodo.15742864, https://doi.org/10.5281/zenodo.15742818, and https://github.com/CESR-lab/ucla-roms.
Author contributions
MH: Writing – original draft, Writing – review & editing, Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization. JM: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Supervision, Writing – review & editing. PD: Conceptualization, Data curation, Methodology, Resources, Software, Writing – review & editing. ML: Conceptualization, Data curation, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing – review & editing. DB: Conceptualization, Investigation, Methodology, Project administration, Resources, Software, Visualization, Writing – review & editing. JM: Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Writing – review & editing.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. This study was supported by the ClimateWorks Foundation grant 23-2582 and the National Science Foundation Award OCE-2343297.
Acknowledgments
Computational resources were provided by the Center for Earth System Research (CESR) at the University of California, Los Angeles (UCLA).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Gen AI was used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fclim.2025.1665329/full#supplementary-material
References
Akan, Ç., McWilliams, J. C., Moghimi, S., and Özkan-Haller, H. T. (2018). Frontal dynamics at the edge of the Columbia River plume. Ocean Model 122, 1–12. doi: 10.1016/j.ocemod.2017.12.001
Anderson, H. J., Mongin, M., and Matear, R. J. (2025). Ocean alkalinity enhancement in a coastal channel: simulating localised dispersion, carbon sequestration and ecosystem impact. Environ. Res. Commun. 7:041012. doi: 10.1088/2515-7620/adce5a
Bach, L. T. (2024). The additionality problem of ocean alkalinity enhancement. Biogeosciences 21, 261–277. doi: 10.5194/bg-21-261-2024
Bednaršek, N., Feely, R. A., Beck, M. W., Alin, S. R., Siedlecki, S. A., Calosi, P., et al. (2020). Exoskeleton dissolution with mechanoreceptor damage in larval Dungeness crab related to severity of present-day ocean acidification vertical gradients. Sci. Total Environ. 716:136610. doi: 10.1016/j.scitotenv.2020.136610
Bednaršek, N., Feely, R. A., Reum, J. C. P., Peterson, B., Menkel, J., Alin, S. R., et al. (2014). Limacina helicina shell dissolution as an indicator of declining habitat suitability owing to ocean acidification in the California Current Ecosystem. Proc. R. Soc. B Biol. Sci. 281:20140123. doi: 10.1098/rspb.2014.0123
Bednaršek, N., van de Mortel, H., Pelletier, G., García-Reyes, M., Feely, R. A., Dickson, A. G., et al. (2025). Assessment framework to predict sensitivity of marine calcifiers to ocean alkalinity enhancement - identification of biological thresholds and importance of precautionary principle. Biogeosciences 22, 473–498. doi: 10.5194/bg-22-473-2025
Burt, D. J., Fröb, F., and Ilyina, T. (2021). The sensitivity of the marine carbonate system to regional ocean alkalinity enhancement. Front. Clim. 3:624075. doi: 10.3389/fclim.2021.624075
Butenschön, M., Lovato, T., Masina, S., Caserini, S., and Grosso, M. (2021). Alkalinization scenarios in the mediterranean sea for efficient removal of atmospheric CO2 and the mitigation of ocean acidification. Front. Clim. 3:614537. doi: 10.3389/fclim.2021.614537
Cai, W.-J., Feely, R. A., Testa, J. M., Li, M., Evans, W., Alin, S. R., et al. (2021). Natural and anthropogenic drivers of acidification in large estuaries. Ann. Rev. Mar. Sci. 13, 23–55. doi: 10.1146/annurev-marine-010419-011004
Cai, W.-J., and Jiao, N. (2022). Wastewater alkalinity addition as a novel approach for ocean negative carbon emissions. Innovation 3:100272. doi: 10.1016/j.xinn.2022.100272
Chao, Y., Farrara, J. D., Zhang, H., Zhang, Y. J., Ateljevich, E., Chai, F., et al. (2017). Development, implementation, and validation of a modeling system for the San Francisco Bay and Estuary. Estuar. Coast. Shelf Sci. 194, 40–56. doi: 10.1016/j.ecss.2017.06.005
Cloern, J. E., Jassby, A. D., Schraga, T. S., Nejad, E., and Martin, C. (2017). Ecosystem variability along the estuarine salinity gradient: examples from long-term study of San Francisco Bay. Limnol. Oceanogr. 62, S272–S291. doi: 10.1002/lno.10537
Doney, S. C., Fabry, V. J., Feely, R. A., and Kleypas, J. A. (2009). Ocean acidification: the other CO2 problem. Ann. Rev. Mar. Sci. 1, 169–192. doi: 10.1146/annurev.marine.010908.163834
Dowell, D. C., Alexander, C. R., James, E. P., Weygandt, S. S., Benjamin, S. G., Manikin, G. S., et al. (2022). The high-resolution rapid refresh (HRRR): an hourly updating convection-allowing forecast model. Part I: motivation and system description. Weather Forecast 37, 1371–1395. doi: 10.1175/WAF-D-21-0151.1
Downing-Kunz, M. A., Work, P. A., and Schoellhamer, D. H. (2021). Tidal asymmetry in ocean-boundary flux and in-estuary trapping of suspended sediment following watershed storms: San Francisco Estuary, California, USA. Estuaries Coasts 44, 2194–2211. doi: 10.1007/s12237-021-00929-y
Fennel, K., Long, M. C., Algar, C., Carter, B., Keller, D., Laurent, A., et al. (2023). Modelling considerations for research on ocean alkalinity enhancement (OAE). State Planet 2-oae2023, 1–29. doi: 10.5194/sp-2-oae2023-9-2023
Fiechter, J., and Moore, A. M. (2024). Physical and biogeochemical properties of california current upwelled source waters. J. Geophys. Res. Oceans 129:e2023JC020164. doi: 10.1029/2023JC020164
Fong, D. A., and Geyer, W. R. (2002). The alongshore transport of freshwater in a surface-trapped River Plume. J. Phys. Oceanogr. 32, 957–972. doi: 10.1175/1520-0485(2002)032<0957:TATOFI>2.0.CO;2
Fregoso, T. A., Wang, R.-F., Ateljevich, E., and Jaffe, B. E. (2017). “A new seamless, high-resolution digital elevation model of the San Francisco Bay-Delta Estuary, California,” in Technical Report 2017-1067, U.S. Geological Survey (Santa Cruz, CA). doi: 10.3133/ofr20171067
García-Reyes, M., and Largier, J. L. (2012). Seasonality of coastal upwelling off central and northern California: new insights, including temporal and spatial variability. J. Geophys. Res. Oceans 117, C03028. doi: 10.1029/2011JC007629
Gasser, T., Guivarch, C., Tachiiri, K., Jones, C. D., and Ciais, P. (2015). Negative emissions physically needed to keep global warming below 2 °C. Nat. Commun. 6:7958. doi: 10.1038/ncomms8958
Geyer, W. R. (2010). “Estuarine salinity structure and circulation," in Contemporary Issues in Estuarine Physics, ed. A. Valle-Levinson (Cambridge: Cambridge University Press), 12–26. doi: 10.1017/CBO9780511676567.003
Geyer, W. R., and MacCready, P. (2014). The estuarine circulation. Annu. Rev. Fluid Mech. 46, 175–197. doi: 10.1146/annurev-fluid-010313-141302
Gill, S., He, J., Yin, J., Sutherland, K., Marsland, E., and Matlin-Wainer, M. (2024). Ocean Alkalinity Enhancement From Coastal Outfalls - Isometric. Available online at: https://registry.isometric.com/protocol/ocean-alkalinity-enhancement (Accessed May 31, 2024).
Goldenberg, S. U., Riebesell, U., Brüggemann, D., Börner, G., Sswat, M., Folkvord, A., et al. (2024). Early life stages of fish under ocean alkalinity enhancement in coastal plankton communities. Biogeosciences 21, 4521–4532. doi: 10.5194/bg-21-4521-2024
Guo, Y., Chen, K., Subhas, A. V., Rheuban, J. E., Wang, Z. A., McCorkle, D. C., et al. (2025). Site selection for ocean alkalinity enhancement informed by passive tracer simulations. Commun. Earth Environ. 6:535. doi: 10.1038/s43247-025-02480-1
He, J., and Tyka, M. D. (2023). Limits and CO2 equilibration of near-coast alkalinity enhancement. Biogeosciences 20, 27–43. doi: 10.5194/bg-20-27-2023
Hetland, R. D. (2010). The effects of mixing and spreading on density in near-field river plumes. Dyn. Atmos. Oceans 49, 37–53. doi: 10.1016/j.dynatmoce.2008.11.003
Ho, M., Kessouri, F., Frieder, C. A., Sutula, M., Bianchi, D., McWilliams, J. C., et al. (2023). Effect of ocean outfall discharge volume and dissolved inorganic nitrogen load on urban eutrophication outcomes in the Southern California Bight. Sci. Rep. 13:22148. doi: 10.1038/s41598-023-48588-2
Hooper, G., Findlay, H. S., Bell, T. G., Wilson, R. W., and Halloran, P. R. (2025). Removal of dissolved inorganic carbon from seawater for climate mitigation: potential marine ecosystem impacts. Front. Clim. 7:1528951. doi: 10.3389/fclim.2025.1528951
Horner-Devine, A. R., Hetland, R. D., and MacDonald, D. G. (2015). Mixing and transport in coastal River Plumes. Annu. Rev. Fluid Mech. 47, 569–594. doi: 10.1146/annurev-fluid-010313-141408
Horner-Devine, A. R., Jay, D. A., Orton, P. M., and Spahn, E. Y. (2009). A conceptual model of the strongly tidal Columbia River plume. J. Mar. Syst. 78, 460–475. doi: 10.1016/j.jmarsys.2008.11.025
Humphreys, M. P., Lewis, E. R., Sharp, J. D., and Pierrot, D. (2022). Pyco2sys v1.8: marine carbonate system calculations in python. Geosci. Model Dev. 15, 15–43. doi: 10.5194/gmd-15-15-2022
Jones, D. C., Ito, T., Takano, Y., and Hsu, W.-C. (2014). Spatial and seasonal variability of the air-sea equilibration timescale of carbon dioxide. Glob. Biogeochem. Cycles 28, 1163–1178. doi: 10.1002/2014GB004813
Kessouri, F., McWilliams, J. C., Bianchi, D., Sutula, M., Renault, L., Deutsch, C., et al. (2021). Coastal eutrophication drives acidification, oxygen loss, and ecosystem change in a major oceanic upwelling system. Proc. Nat. Acad. Sci. USA 118:e2018856118. doi: 10.1073/pnas.2018856118
Khangaonkar, T., Carter, B. R., Premathilake, L., Yun, S. K., Ni, W., Stoll, M. M., et al. (2024). Mixing and dilution controls on marine CO2 removal using alkalinity enhancement. Environ. Res. Lett. 19:104039. doi: 10.1088/1748-9326/ad7521
Kitidis, V., Rackley, S. A., Burt, W. J., Rau, G. H., Fawcett, S., Taylor, M., et al. (2024). Magnesium hydroxide addition reduces aqueous carbon dioxide in wastewater discharged to the ocean. Commun. Earth Environ. 5, 1–10. doi: 10.1038/s43247-024-01506-4
Köhler, P., Abrams, J. F., Völker, C., Hauck, J., and Wolf-Gladrow, D. A. (2013). Geoengineering impact of open ocean dissolution of olivine on atmospheric CO2, surface ocean pH and marine biology. Environ. Res. Lett. 8:014009. doi: 10.1088/1748-9326/8/1/014009
Large, W. G., McWilliams, J. C., and Doney, S. C. (1994). Oceanic vertical mixing: a review and a model with a nonlocal boundary layer parameterization. Rev. Geophys. 32, 363–403. doi: 10.1029/94RG01872
Li, M., Chen, Y., Doyle, R., Testa, J. M., Gagnon, A., Bott, C., et al. (2025). Wastewater alkalinity enhancement for carbon emission reduction and marine CO2 removal. Environ. Res. Lett. 20:044041. doi: 10.1088/1748-9326/adc1e3
Liu, Q., Chai, F., Dugdale, R., Chao, Y., Xue, H., Rao, S., et al. (2018). San Francisco Bay nutrients and plankton dynamics as simulated by a coupled hydrodynamic-ecosystem model. Cont. Shelf Res. 161, 29–48. doi: 10.1016/j.csr.2018.03.008
MacCready, P., and Geyer, W. R. (2010). Advances in estuarine physics. Ann. Rev. Mar. Sci. 2, 35–58. doi: 10.1146/annurev-marine-120308-081015
Mazzini, P. L. F., Pianca, C., Pareja-Roman, L. F., Cole, K. L., Walter, R. K., Castelao, R. M., et al. (2025). Spatio-temporal variability of San Francisco Bay Plume from space. Front. Mar. Sci. 12:1588441. doi: 10.3389/fmars.2025.1588441
McCabe, R. M., MacCready, P., and Hickey, B. M. (2009). Ebb-tide dynamics and spreading of a large river plume. J. Phys. Oceanogr. 39, 2839–2856. doi: 10.1175/2009JPO4061.1
McWilliams, J. C. (2022). Quasi-linear Theory for Surface Wave-Current Interactions. Mathematics of Planet Earth. Singapore: Springer Nature Singapore. doi: 10.1007/978-981-19-2876-5
McWilliams, J. C., Molemaker, M. J., and Damien, P. (2024). Baroclinic sea level. J. Adv. Model. Earth Syst. 16:e2023MS003977. doi: 10.1029/2023MS003977
Molemaker M. J. Damien P. McWilliams J. C. Dollery D. (in prep.). Accurate forcing of high frequency baroclinic waves in regional ocean models. JAMES.
Monismith, S. G., Burau, J. R., and Stacey, M. (1996). Stratification dynamics and gravitational circulation in northern San Francisco Bay. San Francisco Bay: Ecosystem 123:153.
Monismith, S. G., Kimmerer, W., Burau, J. R., and Stacey, M. T. (2002). Structure and flow-induced variability of the subtidal salinity field in Northern San Francisco Bay. J. Phys. Oceanogr. 32, 3003–3019. doi: 10.1175/1520-0485(2002)032<3003:SAFIVO>2.0.CO;2
Moras, C. A., Bach, L. T., Cyronak, T., Joannes-Boyau, R., and Schulz, K. G. (2022). Ocean alkalinity enhancement - avoiding runaway CaCO3 precipitation during quick and hydrated lime dissolution. Biogeosciences 19, 3537–3557. doi: 10.5194/bg-19-3537-2022
Mu, L., Palter, J. B., and Wang, H. (2023). Considerations for hypothetical carbon dioxide removal via alkalinity addition in the Amazon River watershed. Biogeosciences 20, 1963–1977. doi: 10.5194/bg-20-1963-2023
National Academies of Sciences Engineering, and Medicine. (2022). A Research Strategy for Ocean-based Carbon Dioxide Removal and Sequestration. Washington, DC: National Academies Press.
Nederhoff, K., Saleh, R., Tehranirad, B., Herdman, L., Erikson, L., Barnard, P. L., et al. (2021). Drivers of extreme water levels in a large, urban, high-energy coastal estuary - a case study of the San Francisco Bay. Coast Eng. 170:103984. doi: 10.1016/j.coastaleng.2021.103984
Oberlander, J. L., Burke, M. E., London, C. A., and MacIntyre, H. L. (2025). Assessing the impacts of simulated ocean alkalinity enhancement on viability and growth of nearshore species of phytoplankton. Biogeosciences 22, 499–512. doi: 10.5194/bg-22-499-2025
Oschlies, A. (2009). Impact of atmospheric and terrestrial CO2 feedbacks on fertilization-induced marine carbon uptake. Biogeosciences 6, 1603–1613. doi: 10.5194/bg-6-1603-2009
Ou, Y., Xue, Z. G., and Hu, X. (2025). A numerical assessment of ocean alkalinity enhancement efficiency on a river-dominated continental shelf—a case study in the northern Gulf of Mexico. Environ. Res. Lett. 20:024031. doi: 10.1088/1748-9326/adaa8b
Reichl, B. G., and Deike, L. (2020). Contribution of sea-state dependent bubbles to air-sea carbon dioxide fluxes. Geophys. Res. Lett. 47:e2020GL087267. doi: 10.1029/2020GL087267
Renforth, P., and Henderson, G. (2017). Assessing ocean alkalinity for carbon sequestration. Rev. Geophys. 55, 636–674. doi: 10.1002/2016RG000533
Sandoval-Belmar M. Damien P. Sutula M. Kessouri F. McWilliams J. C. Ho M. (in revision). Biogeochemical effects of Golden Gate Strait exchange other land-based inputs to the San Francisco Monterey Bay Coasts. J. Geophys. Res. Oceans.
Sarmiento, J. L., and Gruber, N. (2006). Ocean Biogeochemical Dynamics. Princeton, NJ: Princeton University Press. doi: 10.1515/9781400849079
Schraga, T. S., Nejad, E. S., Martin, C. A., and Cloern, J. E. (2018). USGS Measurements of Water Quality in San Francisco Bay (CA), 2016-2021 (ver. 4.0, March 2023). U.S. Geological Survey Data Release. doi: 10.5066/F7D21WGF
Schulz, K. G., Bach, L. T., and Dickson, A. G. (2023). Seawater carbonate chemistry considerations for ocean alkalinity enhancement research: theory, measurements, and calculations. State Planet 2-oae2023, 1–14. doi: 10.5194/sp-2-oae2023-2-2023
Shchepetkin, A. F., and McWilliams, J. C. (2005). The regional oceanic modeling system (ROMS): a split-explicit, free-surface, topography-following-coordinate oceanic model. Ocean Model 9, 347–404. doi: 10.1016/j.ocemod.2004.08.002
Smith, L. H. (1987). A Review of Circulation and Mixing Studies of San Francisco Bay, California. Reston, VA: Department of the Interior, U.S. Geological Survey. doi: 10.3133/cir1015
Tozer, B., Sandwell, D. T., Smith, W. H. F., Olson, C., Beale, J. R., Wessel, P., et al. (2019). Global Bathymetry and Topography at 15 Arc Sec: SRTM15+. Earth Space Sci. 6, 1847–1864. doi: 10.1029/2019EA000658
Verezemskaya, P., Barnier, B., Gulev, S. K., Gladyshev, S., Molines, J.-M., Gladyshev, V., et al. (2021). Assessing eddying (1/12°) ocean reanalysis GLORYS12 using the 14-yr instrumental record from 59.5°N section in the Atlantic. J. Geophys. Res. Oceans 126:e2020JC016317. doi: 10.1029/2020JC016317
Walters, R. A., Cheng, R. T., and Conomos, T. J. (1985). “Time scales of circulation and mixing processes of San Francisco Bay waters," in Temporal Dynamics of an Estuary: San Francisco Bay, eds. J. E. Cloern, and F. H. Nichols (Dordrecht: Springer Netherlands), 13–36. doi: 10.1007/978-94-009-5528-8_2
Wang, H., Pilcher, D. J., Kearney, K. A., Cross, J. N., Shugart, O. M., Eisaman, M. D., et al. (2023). Simulated impact of ocean alkalinity enhancement on atmospheric CO2 removal in the Bering Sea. Earth's Future 11:e2022EF002816. doi: 10.1029/2022EF002816
Wanninkhof, R. (1992). Relationship between wind speed and gas exchange over the ocean. J. Geophys. Res. Oceans 97, 7373–7382. doi: 10.1029/92JC00188
Wanninkhof, R. (2014). Relationship between wind speed and gas exchange over the ocean revisited. Limnol. Oceanogr. Methods 12, 351–362. doi: 10.4319/lom.2014.12.351
Weiss, R. F. (1974). Carbon dioxide in water and seawater: the solubility of a non-ideal gas. Mar. Chem. 2, 203–215. doi: 10.1016/0304-4203(74)90015-2
Yamamoto, K., DeVries, T., and Siegel, D. A. (2024). Metrics for quantifying the efficiency of atmospheric CO2 reduction by marine carbon dioxide removal (mCDR). Environ. Res. Lett. 19:104053. doi: 10.1088/1748-9326/ad7477
Zhou, J., Izett, J. G., Edwards, C. A., Damien, P., Kessouri, F., McWilliams, J. C., et al. (2023). Modeling the dispersal of the San Francisco Bay plume over the northern and central California shelf. Estuar. Coast. Shelf Sci. 287:108336. doi: 10.1016/j.ecss.2023.108336
Keywords: ocean alkalinity enhancement, marine carbon dioxide removal, estuarine dynamics, ocean modeling, San Francisco Bay, wastewater discharge
Citation: Ho M, Molemaker J, Damien P, Long MC, Bianchi D and McWilliams JC (2026) Ocean alkalinity enhancement in an estuary. Front. Clim. 7:1665329. doi: 10.3389/fclim.2025.1665329
Received: 14 July 2025; Revised: 19 November 2025;
Accepted: 24 November 2025; Published: 12 January 2026.
Edited by:
Maribel I. García-Ibáñez, Spanish Institute of Oceanography (IEO), SpainReviewed by:
Martin Johnson, Ecodiversity Ltd, IrelandHarris Anderson, CSIRO Oceans and Atmosphere, Australia
Copyright © 2026 Ho, Molemaker, Damien, Long, Bianchi and McWilliams. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Minna Ho, bWlubmFob0BhdG1vcy51Y2xhLmVkdQ==; bWlubmFoQHNjY3dycC5vcmc=
Pierre Damien1