The SFB is a large, tidal, semi-enclosed estuary located on the United States West Coast with river inflow on the order of 10³ m³ s⁻¹ 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⁻¹ (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.

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 m³ s⁻¹ 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).

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 m³ s⁻¹ (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⁻³, 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.

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.

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 pCO₂ relative to its initial state. If pCO_2ocean > pCO_2atmosphere after addition, the OAE will reduce the efflux of CO₂ to the atmosphere, conferring the same benefit as a change to or enhancement of an oceanic influx of CO₂. If the ocean becomes undersaturated with CO₂ relative to the atmosphere, there will be a CO₂ 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 CO₂ uptake. This tracer is added to seawater in the model, similar to alkalinity, but it also accounts for air-sea gas exchange, mimicking CO₂ equilibration with the atmosphere. Effectively, the tracer represents the DIC deficit in seawater resulting from the alkalinity addition and its subsequent replenishment by CO₂ ingassing. In the real ocean, the DIC deficit generated by alkalinity addition would be replenished by an air-sea flux of CO₂ from the atmosphere (or reduced outgassing from the ocean) until surface water CO₂ 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 CO₂ ingassing. By tracking the disappearance of the tracer c in the model, we can estimate the amount of CO₂ 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 CO₂, 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 CO₂ equilibration (i.e., that the buffer factor does not change as CO₂ 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, k_w, for the tracer, c, as 2.78 × 10⁻⁶ m s⁻¹ (0.01 m h⁻¹) at the surface, assuming a constant atmospheric CO₂ reservoir. This value is derived from a typical, global CO₂ gas transfer velocity of 5.56 × 10⁻⁵ m s⁻¹ (0.2 m h⁻¹), reduced by a factor of ∂DIC∂[H2CO3*]=20 (Sarmiento and Gruber, 2006, Ch. 8). This factor accounts for the fact that only a fraction of added alkalinity directly affects pCO₂ 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⁻¹ corresponds to the 10 m wind speed u_10m = 8 m s⁻¹, using standard parameterizations of air-sea CO₂ transfer (Wanninkhof, 1992): kw=au10m2(Sc/660)-1/2, with a = 0.31 cm hr⁻¹ s² m⁻², and Sc the Schmidt number dependent on sea surface temperature (Sc ~ 700). The kw=2.78×10-6 m s⁻¹ gas transfer velocity is implemented in the model as a surface tracer flux (sink) term, F_c [mass m⁻² s⁻¹],

calculated online in ROMS (i.e., at every time step). The fixed k_w 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 CO₂ does not change as a result of the air-sea flux of the tracer. Because the atmospheric CO₂ 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 CO₂. 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 CO₂ 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⁻³ (Supplementary material).

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⁻¹ 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 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⁻³] 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⁻³, 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⁻³. 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 m³ s⁻¹ for all rivers in the South Bay). The RMSE between the transect in Figures 2D and F is 3.05 kg m⁻³ and Figures 2E and F is 0.76 kg m⁻³. The RMSE between the average of the two transects and the predicted is 1.53 kg m⁻³.

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).

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.

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.

Tides drive the large-volume transport through the Golden Gate Bridge Strait, reaching over 1,000 m³ s⁻¹ (Figure 4E). Averaging the volume transport in time and summing in depth gives a residual of –1,432 m³ s⁻¹, 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⁻¹ outward, while the velocity near the bottom reaches up to 0.5 m s⁻¹ 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.

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 Fr₁ = 1 (Horner-Devine et al., 2015), is clearly seen in Figure 5E. The Fr₁ is defined as

where u_u is the upper layer velocity, h is the upper layer depth (depth of the buoyant plume), and g′=g(Δρ/ρ0) is the reduced gravity, where g is the gravitational acceleration, and Δρ is the difference in density relative to the ocean density ρ₀. Evaluating this at the liftoff point gives Fr₁ = 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 Fr₁ = 1 is difficult to identify because Fr₁ 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 (Fr₁ < 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⁻¹. The u velocity is strong behind the front and has magnitudes between 0.5 and 0.6 m s⁻¹, 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.

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 CO₂ 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).

The tracer equilibration at the surface determines the rate at which the perturbation to pCO₂ caused by OAE relaxes back to the equilibrium state, or, in other words, the rate of ingassing of CO₂. This tracer flux is equivalent to the net CO₂ uptake (i.e., the flux of CO₂ relative to the original condition). To quantify this rate, the area-integrated, normalized tracer flux F~c, expressed as a percentage, is calculated as follows:

where A is the horizontal area of integration and c_total is the total tracer released up to that point in the simulation. Figure 7A shows a map of tracer equilibration (percentage day⁻¹) 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.

Area-integrated tracer flux reveals a steady equilibration rate with inflections based on tides. Figure 7C shows F~c (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 Fc~ 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, CO₂ ingassing, on the shelf and in the open ocean.

Our study demonstrates the potential of ocean alkalinity enhancement in estuaries to draw down CO₂. 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 CO₂ 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 CO₂ 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.

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.

To obtain a crude estimate of the amount of atmospheric CO₂ 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⁻³ (100 mg NaOH L⁻¹). With a discharge of 2.19 m³ s⁻¹ (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 CO₂ 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 CO₂ 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 CO₂ year⁻¹. To put this in a global context, the IPCC target is to remove 5.5 Gt CO₂ year⁻¹ 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 CO₂ atmospheric extraction target and could be considered one implementation among several to address CO₂ 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⁻³ (equivalent to 2,500 mEq m⁻³) would, after initial dilution, yield a concentration of approximately 20 mEq m⁻³. 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⁻¹ 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).

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 CO₂ 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 k_w 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 CO₂ flux will be governed by the difference between CO_2air and CO_2sea, 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 k_w (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 CO₂ reservoir and excludes ocean-atmosphere feedback associated with CO₂ perturbations (Oschlies, 2009); however, such effects are unlikely to be significant at the scales of CO₂ removal considered here. Additionally, the assumption that the carbonate system behaves approximately linearly constrains the buffering effect of OAE to stay fixed even as CO₂ is taken up by seawater. Furthermore, the factor that reduces the base gas transfer velocity, ∂DIC∂[H2CO3*], 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 CO₂ 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 CO₂ 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 CO₂ 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 CO₂ uptake rate. The results presented here should therefore be considered conservative with respect to both gas exchange efficiency and the influence of physical transport.

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 CaCO₃ precipitation that reduce the capacity for CO₂ 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⁻³). 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 CO₂ 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.