Connections to Tidal Marsh and Restored Salt Ponds Drive Seasonal and Spatial Variability in Ecosystem Metabolic Rates in Lower South San Francisco Bay

We characterized five years of seasonal metabolic patterns in two tidal sloughs in Lower South San Francisco Bay, one surrounded by a broad tidal marsh and the other connected to restored salt ponds (ponds). Despite comparable seasonal dissolved oxygen (DO) patterns, the two sloughs exhibited markedly different seasonal metabolic rates. The in-depth analysis of these small intertidal systems reveals the complexity of processes at work in the fringing habitats of larger estuaries. In the marsh-connected slough, respiration rates peaked in summer (~10 g O2/m2/d\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{O}}_{2}/{\mathrm{m}}^{2}/\mathrm{d}$$\end{document}), coincident with highly dynamic net ecosystem production rates (NEP). The peaks in DO consumption were correlated with temperature and tidal elevation. The most negative NEP rates were associated with the highest nighttime tides, when strong over-marsh respiration rates cannot be balanced by production. Thus, the combination of warmer water and the coincidence of higher-high tides with dark hours during summer may largely explain seasonal DO patterns in the marsh-connected slough. In the pond-connected slough, strong net heterotrophy (NEP < -20 g O2/m2/d\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{O}}_{2}/{\mathrm{m}}^{2}/\mathrm{d}$$\end{document}) correlated with peak phytoplankton production and export from the adjacent ponds in spring, consistent with the hypothesis that organic matter exported from the pond exerts oxygen demand in the turbid slough. The impact of the elevated spring respiration rates appears to be offset by the super-oxygen-saturated water entering the slough from the highly productive pond. Oxygen minima in the slough (< 5 g/m3)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{g}/{\mathrm{m}}^{3})$$\end{document} occurred in summer when outflow from the pond was lower in DO, coinciding with weaker, but still very high, in-slough respiration (~ 10 g O2/m2/d\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm{O}}_{2}/{\mathrm{m}}^{2}/\mathrm{d}$$\end{document}).


Introduction
Nutrient concentrations play a direct role in the ecosystem metabolism of estuaries (Caffrey 2004;Hoellein et al. 2013). Dissolved oxygen (DO) concentrations reflect the balance of metabolic processes and, thus, can serve as indicators of the impacts of anthropogenic nutrient loading. In the absence of pre-urbanization water quality observations, attributing DO patterns to anthropogenic nutrient loading requires careful evaluation of the superimposed factors driving the production and consumption of dissolved oxygen. This can be particularly challenging in estuarine environments where tidal advection may expose water to a diversity of habitat types across short time scales (Nidzieko et al. 2014). Disentangling natural and anthropogenic drivers of oxygen condition in urbanized estuaries requires detailed characterization of metabolic rates across habitat types and time scales, but the majority of estuarine metabolism studies capture only timelimited snapshots of oxygen balances.
Direct rate measurements (e.g., bottle or sediment core incubations) can offer reasonable approximations of oxygen consumption or production rates, but linking these discretein-time-and-place measurements to ecosystem-scale oxygen balances requires broad extrapolation (e.g., Cai et al. 1999). In situ DO sensors enable characterization of oxygen dynamics at sub-tidal time scales. However, the use of these free-water measurements to estimate metabolic rates, a well-established technique in lakes (Staehr et al. 2010), has proven challenging in estuaries, primarily due to difficulties controlling for complex physical processes. Caffrey (2004) employed a free-water approach for characterizing ecosystem metabolism that was broadly applicable across diverse estuarine settings, with the important caveat that the technique is only applicable where tides minimally influence Eulerian observations of DO (see related analysis of advective effects in Howarth et al. 2014). A short-term manipulation experiment enabled Nidzieko et al. (2014) to control for tidal advection in estimating metabolic rates from free-water DO measurements in a coastal-California slough. But the time-limited nature of such an experiment limits exploration of metabolic balances under seasonally varying forcing conditions. Beck et al. (2015) leveraged the known periodicities of tidal components and the solar cycle to reduce advective effects on estuarine metabolism estimates. However, this method yields limited results when tides and the solar cycle are highly correlated (as is generally the case in estuaries with mixed semi-diurnal tides).
We present metabolic rate estimates from five years of freewater DO measurements in two San Francisco Bay (SFB; Fig. 1) tidal sloughs. We account for advective effects by leveraging monotonic salinity gradients to pair water-parcel measurements on flood and subsequent ebb tides, enabling quantification of net ecosystem production (NEP) and respiration (R) rates over five years (2015-2019) of nearcontinuous 15-min DO measurements. Our unique temporal characterization of metabolic balances, paired with co-located measurements of other water quality parameters, allows us to detail correlates with metabolic rates that help explain the occurrence of seasonal DO minima. Differentiating these patterns between the proximal but unique study sloughs (one surrounded by marshland and the other surrounded by restored salt ponds) highlights the complex interplay of habitat heterogeneity and seasonality on estuarine oxygen condition.
SFB presents a valuable case study for understanding the superposition of natural and anthropogenic water quality processes. Despite severe nutrient enrichment, DO concentrations in most of SFB rarely drop below regional indicator thresholds (San Francisco Bay Basin Plan 2019) for sustained periods (Cloern et al. 2020). The Lower South Bay (LSB; Fig. 1b) sub-embayment is a notable exception. The tidal sloughs and restored salt ponds (RSPs) surrounding LSB can exhibit sharp and sustained drops in DO (MacVean et al. 2018). Observations of depleted DO at the Dumbarton Narrows (where LSB connects to the rest of SFB; Fig. 1b-c) during low tides have guided a conceptual model wherein oxygen concentrations are drawn down in waters forced into the LSB periphery on flood tides, and these oxygen-depleted waters are advected bayward on ebb tides (Crauder et al. 2016). However, patterns in DO dynamics in the LSB can be highly variable and inconsistent, making a holistic understanding of the metabolic rates and balances of LSB difficult to pin down.
Previous studies provide limited insight on linkages between nutrient loading and DO dynamics in SFB. Sutula et al. (2017) evaluated the use of chlorophyll-a concentrations as thresholds for broader eutrophication symptoms in SFB. They cite uncertainty around chlorophyll-a as an indicator of oxygen impairment and specifically point to the need for quantifying metabolic contributions from the sloughs and RSPs surrounding LSB. Previous quantification of metabolic balances in this area has been limited. Grenz et al. (2000) linked elevated sediment oxygen demand in south SFB to settled organic matter from a spring algal bloom. Thébault et al. (2008) quantified the "delicately poised" balance between high rates of oxygen production and consumption in an RSP connected to LSB. However, like most studies of aquatic metabolic rates, these investigations explored only limited time periods (up to a few months) and did not include slough habitats.
This study was motivated by a bay-wide effort to understand and respond to the potential impacts of nutrient loading from wastewater treatment plants (WWTPs), the SFB Nutrient Management Strategy (NMS). The NMS has been tasked with quantifying the role of nutrients in SFB DO dynamics, including addressing: (1) what processes drive the observed DO patterns? (2) how and to what extent are these processes modified by anthropogenic nutrient loading? (3) what nutrient load reductions might be necessary to mitigate periodic hypoxia? Such questions are relevant to urbanized coastal ecosystems throughout the world (Boicourt et al. 2012). In addition to broadly informing nutrient management for urbanized estuaries, this analysis highlights the degree to which metabolic rates can vary within estuaries, in both space and time, due to the interplay of natural and anthropogenic processes.
Finally, while the direct aims of the study were motivated by local needs, the results, analyses, and discussion presented herein are relevant to estuaries throughout the world. Many active-margin coastlines (typified here by the open coastline of California) are dotted with small, highly dynamic estuaries that are similar in character to South San Francisco Bay; examples include Elkhorn Slough (Nidzieko and Monismith 2013), the numerous estuaries in the Southern California Bight (Harvey et al. 2020), and Pescadero Creek (Williams and Stacey 2016). In addition to revealing the complex linkages between physical and biogeochemical dynamics of these smaller systems, this study contributes understanding of the integrated processes of larger estuaries such as Chesapeake Bay or the Salish Sea, systems that are fringed by tidal marshes and flats, just as Newark and Alviso Sloughs fringe San Francisco Bay.

Methods
The NMS began continuous water quality monitoring in and around LSB in 2014. Prior to these efforts, discrete monthly channel measurements and occasional short-term studies represented the only observations related to DO dynamics and associated driving variables in LSB. The methodology described below is intended to leverage this near-continuous, multi-year dataset (2015-2019) to elucidate patterns in LSB oxygen balances across seasonal time scales.

Site Description
San Francisco Bay is a large estuary situated along California's central/north Pacific coast. The local climate exhibits a classic Mediterranean pattern, with a warm and very dry season (May to October) and a mild wet season (November-April) punctuated by "atmospheric river" precipitation events (most common in January-March). This study focuses on the southern region of SFB (South Bay; Fig. 1a), which is largely cut complex with monitoring station locations in Alviso Slough and at the Pond A8 outlet. In (b), (c) and (d), red stars and yellow arrows denote EXO2 water quality sensors and exchange gates/breaches into restored salt ponds respectively. Only the Alviso and Newark EXO2 stations were used for metabolism estimates. Note that bathymetry in (a) is referenced to the inset colorbar and bathymetry in (b-d) is referenced to the colorbar at right off from the major riverine freshwater inputs that flow through the North Bay and out to the Pacific Ocean. The South Bay (Fig. 1a) generally exhibits decreasing north-to-south salinity gradients that are driven by freshwater inputs from the surrounding urbanized watershed and WWTP effluent; the magnitude of these gradients increases markedly following winter precipitation events. Combined with mixed semi-diurnal tides (2.5 m between MLLW and MHHW), these longitudinal gradients cause rapid salinity fluctuations across tidal cycles.
San Francisco Bay's southernmost sub-embayment, the Lower South Bay (LSB; Fig. 1b), sees some of the highest nutrient concentrations and lowest DO concentrations in SFB. Constricted from the broader South Bay at the Dumbarton Narrows and ringed by interwoven tidal marsh, creeks, and sloughs, and by both commercially active and restored salt ponds, the LSB exhibits particularly dynamic biogeochemical signals. Water quality parameters, including DO, have strong spatial gradients and exhibit large swings over tidal to seasonal time scales. We focus this study on two LSB tidal sloughs whose distinct connections to the surrounding habitat broadly represent a range of settings characteristic of other LSB sloughs.
Newark Slough (Fig. 1c) extends from the northern edge of LSB into a broad tidal marsh; channel depth varies between roughly 1.5 and 4 m with the tides. Active commercial salt ponds (shown as light green in Fig. 1c) are separated from the slough and surrounding marsh by earthen levees. During higher-high tides, marsh flooding dominates the wetted subembayment habitat, both areally and volumetrically (Fig. 2a,  b). At low tides, waters retreat from the marsh into the slough channel. Newark Slough receives no external freshwater inputs beyond direct precipitation onto the surrounding marsh.
Alviso Slough extends into the leveed network of RSPs at the southeastern end of LSB ( Fig. 1d) and has channel depths comparable to those in Newark Slough (about 1.5 to 4 m). As part of the South Bay Salt Pond Restoration Project, the RSPs have been restored to tidal action via gated weirs and levee breaches; connections to Alviso Slough are shown as yellow arrows in Fig. 1d. Due to the surrounding levees, tidal-marsh flooding around Alviso is limited compared to that around Newark (Fig. 2c, d). During the wet months, inflows from the Guadalupe River freshen the salinity in Alviso Slough; inflows are minimal during the dry season. Coyote Creek, the deeper channel at Alviso Slough's outlet (Fig. 1b, d), receives loads of nitrogen and phosphorous from an upstream WWTP (dissolved inorganic nitrogen: ~ 5000 kg day −1 ; phosphate: ~ 200 kg day −1 ; Crauder et al. 2016). LSB's strong tides pump nutrient-enriched waters (100 Mμ DIN) up Alviso Slough and into adjacent RSPs during flood tides (SFEI 2017).
Tidal currents in Alviso Slough can exceed 50 cm/s during spring tides (Shellenberger et al. 2015). In situ velocity measurements were not available for Newark Slough, but comparable channel geometry and unpublished model results support estimates of similar-magnitude flow speeds.

Field Data Collection-Water Quality Monitoring Stations
Multi-parameter water quality sondes (YSI EXO2; www. ysi. com) were deployed at fixed elevations above bottom in Newark Slough (90 cm above bottom) and Alviso Slough (45 cm above bottom), at the Dumbarton Narrows (12.65 m above bottom; average depth 3 m), and also in the outlet channel (45 cm above bottom) connecting the RSP A8 outlet gate to Alviso Slough (Fig. 1c, d). The instruments collected burst averages (10 s at 1 Hz) every 15 min. Measured parameters included conductivity (with derived specific conductivity, SpC, and salinity, S), pressure/depth, temperature, chlorophyll-a fluorescence (fChla), and dissolved The vertical limits of the plots show approximate minimum and maximum observed tidal height during the five-year study period oxygen saturation and concentration. fChla was converted to chlorophyll-a concentration as 4 mg/L per relative fluorescence unit following an SFB-specific calibration ; reported concentrations should be considered approximate.
Instruments were field-serviced (cleaning and calibration checks) every 3-4 weeks and fully bench serviced (deep cleaning and re-calibration) every 2-3 months. Mechanical wipers and copper sensor guards minimized biofouling between servicing trips. All data were passed through automated filters and were manually reviewed ; invalid data due to biofouling and/or instrument drift/ malfunction were removed (see gaps in Fig. 3).

Field Data Collection-Tidal Elevation and Bathymetry Data
Six-minute tidal elevation data from a National Oceanic and Atmospheric Association station (NOAA #9414523; https:// tides andcu rrents. noaa. gov/ water levels. html? id= 94145 23) in the South Bay were used to estimate high-tide tidal magnitude in LSB. The depth signals from the EXO2 instruments offered reliable temporal signals for approximate local water column depth but, due to shifts in exact mounting height (order of 10 0 cm) , were less reliable for estimating datumreferenced high-tide magnitude over the 5-year study period.
The hypsopgraphic curves shown in Fig. 2 are based on a 2-m resolution bathymetric/topographic digital elevation model generated by the United States Geological Survey Earth Resources Observation and Science Center (www. usgs. gov/ core-scien ce-syste ms/ eros/ coned). Delineation of habitat type was generated by the San Francisco Estuary Institute as part of the San Francisco Bay Shore Inventory project (SFEI 2016).

Field Data Collection-Meteorological Data
An hourly wind speed time-series (normalized to 10-m elevation) was compiled for the LSB (122.0653 W, 37.4755 N) based on the bay-wide wind model developed by King (2019). The King (2019) model generates an hourly 10-m-elevation-normalized wind field by interpolating wind speed and direction measurements from 52 stations around SFB. Hourly incoming shortwave radiation data were compiled from the California Irrigation Management Information System (CIMIS;.cimis.water.ca.gov) station 171 (Union

Analytical Techniques
Free-water metabolism techniques were applied to the long-term dissolved oxygen data from Newark and Alviso sloughs to estimate NEP and R. These rates were compared to a range of water quality variables to evaluate drivers of metabolic variability. Estimates of in-slough gross primary production (GPP) were calculated using an empirical model (Cole and Cloern 1987; based on chlorophyll-a concentrations and light conditions) to support interpretation of the NEP and R signals.

Analytical Techniques-Lagrangian Water Mass Tracking and Free-Water Metabolism Estimates
In a Lagrangian frame, the change in dissolved oxygen concentration of a parcel of water can be described as a balance between surface-gas exchange and volumetric gross primary production ( GPP V ) and respiration ( R V ) (Odum 1956): C represents the mass concentration of DO; k is the surface-gas-exchange piston velocity. C sat is the saturation DO concentration for a given temperature and salinity. H is the mean depth of the vertically mixed water parcel.
Assuming a well-mixed water column, the combined benthic-water column respiration and production effects can be combined into aggregate areal terms R and GPP (with units of g O 2 ∕m 2 ∕d) , yielding Herein, GPP, R, and NEP should be taken to represent the areal ecosystem metabolic rates. Net ecosystem production is simply the balance between production and respiration,NEP = GPP − R . Thus, we can simplify Eq. (2) into an explicit Lagrangian estimate for NEP: If we isolate our analyses to dark (nighttime) hours, Eq. (3) expresses a relationship for ecosystem respiration: In the preceding equations, C and H are directly measured, and C sat easily derived. We follow the standard convention of estimating k as a function of flow speed, water-column depth, and wind speed (see Supplementary Material Section A for details).
The key challenge in applying these equations to estimate metabolic rates is the Lagrangian nature of the formulation; multiple measurements must be made on the same parcel of water. If we assume that advective transport dominates over dispersive mixing on the time scales between measurements (here, minutes to hours), then non-reactive tracers (e.g., salinity) can be used to identify the passage of parcels of water. This assumption is readily justified by estimating a Peclet number, Pe = uL K x , which compares advective versus dispersive rates over a distance L; u is velocity and K x is horizontal dispersion. Taking L = 2350 m and u = 0.2 m/s (see Supplemental Material B) and K x = 10 m 2 /s as an upper bound in tidal channels with significant trapping (Ralston and Stacey 2005), we find that Pe = 50, supporting the assertion that advection is much faster than dispersion with regard to the signals of interest.
In tidal sloughs, reversing flows allow fixed sensors to measure the same parcel of water on flood tides and again on subsequent ebb tides. When longitudinal salinity gradients are monotonic, parcels can be uniquely identified, and their water quality characteristics paired, from flood to subsequent ebb, enabling two time-lagged measurements per parcel. Note that this approach requires that slough channels be reasonably one-dimensional (cross-sectionally well mixed and with dominant longitudinal flow); see Nidzieko et al. (2014) for an example of this approach and Supplementary Section B for justification of cross-sectionally well-mixed conditions in Newark and Alviso sloughs for the present study. When these assumptions are justified, the change in DO concentration for a parcel of water is easily estimated from the flood (1) and subsequent ebb (2) measurements as and metabolic rates can be estimated from Eqs. (3) and (4) (standard oxygen balance methods). See Supplementary Materials Section C for technical details on the application of the Lagrangian water mass tracking technique for controlling for advective effects in our oxygen balance calculations. Metabolic estimates at the Pond A8 station were not possible due to multi-dimensionality of in-pond flow (and exchange among ponds) immediately upstream from the sensor.
Note that this technique limits rate estimates to time windows centered on high-slack tides. For our analysis, we collapse rate estimates from individual water-parcel pairings into averages associated with each high-slack-centered time window; we herein refer to averages within the high-slackcentered time windows as "high-tide" or "window" averages. These discrete rate measurements represent snapshots influenced by the solar insolation (and other processes changing on short time scales) during the high-slack-centered window rather than representing full-day-average rates. Interpretation of NEP estimates should consider time scales that are long enough to capture a random distribution of hour-of-day of high-slack tide timing. Interpretation of R estimates must consider covariation of potential metabolic drivers with the solar cycle (since R estimates are limited to dark hours). Further, the established approach of adding R (nighttime NEP) to NEP to estimate daily mean GPP (Staehr et al. 2010) is not applicable. To support interpretation of our calculated R and NEP rates, we approximate in-slough GPP using an established primary production model for SFB (Cole and Cloern 1987; see Supplementary Material D for details).

Seasonal Water Quality Patterns
The salinity signals at both Alviso and Newark sloughs follow patterns that reflect seasonal hydrology. Salinity decreases during the wet winter months (December-March) and increases through the dry late spring and summer (April-August) before leveling off in the fall (Fig. 3a). The magnitude of the observed fluctuations is consistent with inter-annual hydrologic patterns. Salinity levels reached study-period minima during the extremely wet winters in 2017 and 2019, and the winter minima were muted during the comparatively dry 2015, 2016, and 2018 rainy seasons. Differences in salinity between the sites reflect proximity to the open Bay and the effect of local freshwater inputs ( Fig. 1f shows representative LSB hydrologic conditions from Coyote Creek). Newark Slough is consistently more saline than Alviso Slough, with fall salinity levels approaching ocean concentrations. Alviso sees much larger salinity fluctuations across tidal cycles (shown by the thickness of the light-colored variance bands in Fig. 3a), reflecting the steep salinity gradients generated by freshwater inputs from the Guadalupe River.
Dissolved oxygen levels generally peak in winter and early spring at both sites, with peaks at Alviso regularly exceeding saturation (Fig. 3b). In drier years (2015,2016,2018), Alviso Slough reached DO minima in spring (< 50% of saturation), immediately following the winter peak. Following the wet 2017 and 2019 winters, the minima were delayed into summer. Summer DO in Newark Slough dropped to about 70-80% of saturation.
Chlorophyll-a peaks at Alviso Slough coincide with the dissolved oxygen maxima (Fig. 3b, d). The winter/spring peak was most pronounced during abnormally dry 2015. Extremely wet conditions in 2017 appear to have washed out or diluted the typical winter/spring bloom signal. However, this was not the case during the similarly wet 2019 winter/spring. Compared to Alviso Slough, Newark Slough chlorophylla levels were only rarely elevated. Modest chlorophyll-a increases, accompanied by increased DO concentrations, occurred in spring 2018 and 2019, coinciding with typical seasonal timing of increased production in the deep subtidal habitats of the open bay. Brief early-fall (September/October) spikes in chlorophyll-a do not appear to be reflected in the dissolved oxygen signal. Figure 4 illustrates three example tidal cycles among the hundreds that were used to quantify metabolic rates during the 5-year study period (1058 and 874 for Newark and Alviso, respectively). To estimate metabolic rates, we incorporated reaeration estimates into the individual ΔC∕ΔT values for each parcel (each dot in Fig. 4a), then averaged these values into the high-slack-centered window-average NEP rates. Estimates of R are simply the nighttime NEP values (tidal cycles with max(SW) = 0 in the high-slack-centered window). During nighttime tidal cycles, we see the expected negative relationship between time lag and change in DO concentration. During the daytime tidal cycle shown in Fig. 4a (green), this pattern is reversed, showing the positive relationship expected for a net productive time interval.

Quantification of Seasonal Metabolic Rates
Distinct within-slough seasonal patterns emerge when these rate estimates are aggregated into monthly distributions (Fig. 5). Note that seasonal covariation between dark hours and tidal elevation (Supplemental Figure S.C1) necessitates specifying that we are reporting nighttime R (herein R night ), which may not be representative of full-day-average R. That is, if tidal elevation has an effect on ecosystem respiration rates-which we suspect it might, given the dependence of habitat distribution on water level (Fig. 2)-distributions of R night would carry a seasonal bias (Supplemental Figure S. C1b). NEP is observed at all hours of the day and is therefore not biased by seasonal covariation between dark hours and tidal elevation (e.g. spring tide high-tide timing; Supplemental Figure S.C1a). The similar seasonal patterns between NEP (unbiased distributions) and R night (potentially seasonally biased distributions) point to a real metabolic impact associated with seasonally shifting phasing between the mixed tides and the solar cycle.
Newark Slough showed strongly seasonal metabolic rates. Median R night dropped to about 2 g O 2 ∕m 2 ∕d in winter, then increased beginning in March before hitting a June peak of 14.1 g O 2 ∕m 2 ∕d (Fig. 5a). Nighttime respiration rates remained elevated through August and September before dropping back to winter levels in October/November.
Lower-quartile NEP rates at Newark mirror the median R night pattern, with strong net heterotrophy in summer (NEP < -10 g O 2 ∕m 2 ∕d ; Fig. 5c). Median NEP rates were much closer to metabolic balance, ranging from -1.4 g O 2 ∕m 2 ∕d in January to -3.4 g O 2 ∕m 2 ∕d in August. Upper-quartile NEP values were near zero in winter but are consistently positive through the spring and summer seasons (varying between 0 and 3 g O 2 ∕m 2 ∕d ). Estimates of in-slough GPP were insufficient to account for the offsets differences between monthly distributions of NEP and −R night , pointing to GPP effects not accounted for in the simple slough-channel water-column GPP model. Seasonal metabolic patterns in Alviso Slough are markedly different from those in Newark Slough. Nighttime respiration rates increase into winter before peaking in spring, reaching a median of 27.2 g O 2 ∕m 2 ∕d in April (Fig. 5b). Median R night rapidly tapers in late spring, dropping to an annual minimum in July (8.1 g O 2 ∕m 2 ∕d ). Elevated winter/ spring chlorophyll-a concentrations (Fig. 3c) may represent some in-slough GPP (Fig. 5b). However, winter offsets between NEP and −R night exceed the in-slough GPP predictions, pointing to potential over-marsh production effects not accounted for in the slough-channel GPP model. In any case, primary production at Alviso, whether in the slough channel or over-marsh, was never sufficient to overcome the exceedingly high respiration rates. The seasonal NEP pattern at Alviso is dominated by seasonal respiration effects; almost no positive NEP values were recorded in the 5-year study period (Fig. 5d).
Seasonal patterns of NEP at the two sloughs were apparent across the five study years, but the timing and magnitude of the seasonal minima and maxima varied (Fig. 6). At Newark Slough, winter NEP rates were consistently near metabolic balance, if slightly net heterotrophic, but summer NEP showed pronounced variability between years. The most negative NEP rates were observed between June and August. With the exception of 2016, each year saw at least one summer month with median NEP rates below -5 g O 2 ∕m 2 ∕d . Newark Slough was consistently net productive during summer 2016, with median NEP reaching 2.6 g O 2 ∕m 2 ∕d in July. At Alviso Slough, median NEP rates consistently decreased from January until April/May/June and then increased from June to August before reaching peaks in late-summer to early-fall. In 2016 and 2019, Alviso NEP rates rapidly decreased from September into winter; this pattern was less pronounced in the other study years. Alviso exhibited stronger net heterotrophy than Newark across every month in each of the five study years.
Seasonal metabolic patterns were insensitive to uncertainty associated with reaeration estimates. Because oxygen saturation is directly measured, uncertainty in reaeration estimates is primarily driven by imperfect piston velocity estimates. Flow effects were more important than wind effects in setting the piston velocity (k) estimates (Supp. Figure S.A2). We tested the sensitivity of our metabolic rate estimates to a range of assumed flow conditions. For piston velocity estimates based on flow speeds scaled by factors of 0.5 and 2, monthly median NEP estimates for Alviso varied, on average, by 0.22 and 0.35 g O 2 ∕m 2 ∕d , representing 1.6% and 2.3% changes from the base case rate estimates. R night  estimates for the same sensitivity test yielded average differences of 0.19 and 0.27 g O 2 ∕m 2 ∕d, 1.4% and 1.8% changes from the base case rates respectively. For Newark Slough, monthly median NEP estimates varied on average by 0.17 and 0.23 g O 2 ∕m 2 ∕d, 8.5% and 9.2%, and R night estimates by 0.18 and 0.19 g O 2 ∕m 2 ∕d, 4.6% and 4.3%. Overall, the sensitivity of surface-gas exchange rates to a wide range of flow speeds, and therefore piston velocity estimates, is limited by the fact that the oxygen saturation conditions (which are directly measured) play a leading-order role in the exchange-rate magnitude and are the only term that can change the sign of the surface flux correction term in Eq. (3) (Supp. Figure S.A3). Our metabolic rate estimates are robust to the inherent uncertainty in estimating surface gas fluxes.

Drivers of Seasonal Metabolic Patterns
Metabolic rates were driven by different factors at Newark Slough compared to Alviso Slough. We used ordinary least squares linear regression to explore the univariate relationships of both R night and NEP with: (1) window-average temperature; (2) tidal-cycle high-tide elevation; (3) window-average chlorophyll-a concentration; (4) long-term chlorophyll-a concentration; (5) areal ratio of wetted habitat type at high tide (marsh:channel); and (6) ratio of water volume overlying habitat type at high tide (marsh:channel) (Fig. 7). Long-term chlorophyll-a concentration is calculated as the average chlorophyll-a concentration over the preceding 60 days in the continuous 15-min dataset. This value is intended to capture seasonal biomass in the sloughs that may be masked in our high-slack-centered windowing. Elevated chlorophyll-a concentrations are generally observed during low tides in the sloughs (up-slough or pond-sourced material), but these concentrations may impact high-tide metabolic rates (when water-column chlorophyll-a concentrations tend to be lower) via bed deposition of labile organic matter.
On the surface, there is an explicit lack of independence between our rate estimates (calculated partially as a function of water-column depth; Eq. 3) and the high-tide-elevation predictor variable. However, the use of the depth term in our rate estimates is to isolate a biological signal in a physically dynamic environment (water-column depth cannot be excluded from rate estimates). Despite the lack of independence between rate estimates and high-tide elevation, regressions between metabolic rates and high-tide elevation still offer insight into relationships between tidal dynamics (e.g., spring/neap) and ecosystem metabolism.
Newark R night had strong, positive correlations with temperature and tidal elevation, but very weak relationships with both short-and long-term chlorophyll-a concentrations (Fig. 7a-d). The tidal elevation relationship was particularly strong; at tidal elevations greater than 2.6 m (MLLW), R night was consistently elevated (Fig. 7b). This threshold is consistent with the approximate inflection depth in the Newark hypsographic curves (Fig. 2a, b) at which the surface area of flooded marsh increases nearly fourfold; this relationship is reflected in Fig. 7e,f, where R night shows similarly strong correlation with increased marsh:channel habitat ratio. Interestingly, Alviso Slough had negative relationships between R night and both temperature and tidal elevation. Alviso respiration rates were positively correlated with chlorophyll-a concentrations, more strongly so with averaged concentrations over the preceding 60 days (Fig. 7j).
The strong correlations between R night and high-tide elevation point to seasonal biasing in the monthly distributions of R night in Fig. 5a, b. Accordingly, care must be taken in interpreting relationships between R night and co-varying season, temperature, and tidal elevation. Relationships between NEP and potential metabolic drivers improve confidence in the relationships with R night (Fig. 7m-x). Observations of NEP do not carry a seasonal bias, yet they still exhibit relationships that match those observed in the R night data. Newark NEP is negatively correlated with temperature and tidal elevation (Fig. 7m, n), consistent with the positive correlations with R night . Alviso NEP is positively correlated with tidal elevation (Fig. 7t and negatively correlated with chlorophyll-a concentrations (Fig. 7u, v), consistent with decreasing and increasing respiration rates for these drivers respectively. Weaker coefficients of determination for NEP relationships may reflect scatter associated with variable light conditions for the discrete NEP measurements.
Seasonal covariation between nighttime tidal elevation and temperature obfuscates interpretation of their individual contributions. However, stronger correlations with tidal elevation, and the manifestation of the tidal-magnitude-biased seasonal R night pattern in the non-biased seasonal NEP pattern, allow us to confidently infer the importance of tidal elevation and its phasing with the solar cycle, on slough metabolism. We interpret these relationships in the context of the surrounding habitat and seasonal metabolic patterns in the discussion section that follows. Fig. 7 Ordinary least squares regressions of nighttime respiration rates and net ecosystem production (NEP) rates at Newark (a-f, m-r) and Alviso (g-l, s-x) versus potential drivers: (a/g/m/s) window-average temperature (℃); (b/h/n/t) high-tide elevation (meters above MLLW); (c/i/o/u) window-average chlorophyll-a concentration (μg/L); (d/j/p/v) 60-day backwards-looking average chlorophyll-a (μg/L); (e/k/q/w) ratio of surface area of wetted habitat type, marsh:channel, at tidal-cycle high tide; (f/l/r/x) ratio of volume overlying habitat type, marsh:chnnel: at tidl-cycle high tide. Goodness-offit (r 2 ) is displayed only where p < 10 −3 . At Newark, n = 1058 and 412 for NEP and R respectively. At Alviso, n = 874 and 268 for NEP and R respectively ◂

Evaluation of Rate Estimates
Estimates of metabolic rates using free-water dissolved oxygen methods are inherently noisy, particularly in physically dynamic environments like LSB (Staehr et al. 2010). While many tidal cycles exhibit the coherent and expected patterns shown in Fig. 4a, some do not; estimates of dC/dt can be variable among parcels within an individual tidal cycle. Variability in dC/dt within tidal cycles can be partially explained by water quality heterogeneity but may also be due to various forms of error. The 15-min data-collection time step is somewhat coarse for tracer-based parcel pairing across 4-h time windows. This is partially handled by interpolating to a 5-min time step, but linear interpolation in time includes its own uncertainty. While assumptions of crosssectionally well-mixed sloughs may be reasonable (Supp. B), it is unlikely that a given salinity (i.e. tracer) condition perfectly represents an individual, homogeneous fluid parcel. Peclet number estimates support the assumption of limited dispersion at parcel-pairing time scales, but over-marsh dispersion may introduce uncertainty into our method. Finally, reaeration estimates make assumptions about flow and wind speed conditions at the water surface, and application of these estimates across 1-to-4-hour parcel-pairing windows requires further averaging.
The scale of the dataset helps overcome the uncertainty associated with the methodological approach. Our results represent metabolic estimates from almost two-thousand tidal cycles, each of which is based on 5 to 18 paired parcel measurements, and we focus our reporting on aggregate seasonal patterns. We assume that error associated with interpolation and mixing assumptions is randomly distributed. We have shown that error associated with surface-gas exchange estimates generally represents only a small proportion of the rate estimates. The effect of reaeration on observed rate estimates is particularly small for the largest magnitude rates, which are perhaps some of the most noteworthy values from a management perspective (Supp. Figure S.A4).
The discrete nature of the rate-estimate technique can bias rate distributions if the measurement timing is non-random relative to factors affecting metabolic rates. Time windows for examining distributions of NEP estimates should include a random sampling across the hour of the day. The hour-ofday of NEP measurements is random within the monthly distributions shown in Fig. 5c/d. Randomness in hour-of-day is not possible for respiration measurements. Thus, we are careful to specify that we are reporting R night . Single R measurements are not biased-the reported values are meaningful measurements-but care must be taken in interpreting distributions that can be biased by seasonally shifting resonances between measurement timing and forcing variables. Finally, we must interpret results from this method remembering that they only represent high-tide-centered rates.
In spite of various uncertainties, coherent patterns emerge that are apparent across study years and are largely explainable by associated conditions. Our approach of pairing measurements on flood and subsequent ebb tides, modified from a comparable method used by Nidzieko et al. (2014), appears to be a valuable technique for accounting for tidal advection in free-water metabolism estimates in tidal sloughs, an ongoing challenge in estuarine metabolism studies (Loken et al. 2021). Application of this method comes with the important caveat that the study sloughs meet the following criteria: (1) monotonic longitudinal gradients exist in a conserved tracer (e.g. salinity); (2) slough-channel cross-sections are reasonably well mixed; and (3) sufficient data are available (both sampling frequency and study period) to minimize biasing associated with discrete measurements and for a signal to emerge from the noise. We expect that (3) would partially depend on the magnitude of the metabolic rates in a given system. Given the apparently strong rates of oxygen consumption, LSB serves as an excellent test case for this method.

Interpretation of Seasonal Metabolic Patterns
At Newark Slough, relationships between metabolic rates and environmental conditions are consistent with simple conceptual models and with observations in other estuaries. The dependence of both R night and NEP on temperature is fundamental (Yvon-Durocher et al. 2012). The dependence on tidal elevation is consistent with elevated respiration rates associated with increased marsh flooding at higher tides. This hypothesis is supported by Fig. 7e, f; respiration rates are strongly correlated with tidally driven marsh:channel habitat ratio. At 2.6 m MLLW, an approximate threshold for consistently elevated R night (Fig. 7b), marsh habitat constitutes 35% and 71% of the sub-embayment volume and surface area respectively; these proportions rapidly increase at tidal elevations above MHHW (Fig. 2a, b). Elevated respiration rates due to marsh flooding are consistent with existing SFB studies that estimate significantly higher respiration rates over marshes compared to in open water; Pearcy and Ustin (1984) estimated respiration rates in northern SFB marshes of about 30 g O 2 ∕m 2 ∕d ; Caffrey et al. (1998) estimated respiration rates of about 2 g O 2 ∕m 2 ∕d on the openwater shoal in the South Bay. Nidzieko et al. (2014) noted a similarly strong dependence of the metabolic balance of Elkhorn Slough on marsh flooding and further noted that the alignment of higher high tides with nighttime hours led to strong net heterotrophy during summer spring tides. Our results show a similar pattern, with the strongest net heterotrophy emerging during higher nighttime high tides (Fig. 8a). Higher high tides that align with daylight hours produced more positive NEP rates than lower high tides at similar light levels (Fig. 8a), pointing to the highly dynamic nature of the intertidal marsh ecosystem. The dependence of NEP on the interplay of light and tidal elevation conditions likely explains why nighttime respiration rates, which may not represent day-average R, still manifest in the seasonal NEP patterns shown in Fig. 5. The convergence of warmer temperatures with increased nighttime tidal magnitude during summer may explain why the most negative NEP rates were observed during summer months at Newark slough (Figs. 5c and 6a). This interplay explains 53% of the variability in observed R night rates (Fig. 8b).
Alviso Slough exhibits more complex metabolic patterns than those in Newark Slough. The negative univariate correlation between R night and temperature (Fig. 7g), which contradicts fundamental models, may largely be due to the coincidence of elevated biomass with cooler winter/spring water temperatures (Fig. 3d, e). The fact that this relationship drops out in the full set of NEP data (Fig. 7s) may be further indication that the effect shown in Fig. 7g is simply a relic of collinearity between temperature and tidal elevation in the R night data. The statistical association and seasonal covariation between oxygen demand and long-term chlorophyll-a concentrations (Fig. 7j, v) may be explained by the complex pond-slough exchange dynamics apparent during bloom events (Fig. 9).
The elevated mean chlorophyll-a concentrations in Alviso Slough during the winter/spring 2016 bloom (Fig. 3c) appear to be driven by outflow from Pond A8 that only arrives at the Alviso sensor well into the ebb (4-day example shown in Fig. 9b); these elevated concentrations are not apparent in our high-slack centered rate estimate windows. The lagged and muted DO and chlorophyll-a signals at Alviso relative to Pond A8, and the ebb-centered timing of their peaks, point to biomass export from Pond A8 rather than high rates of in-slough production (i.e., production in the main channel of Alviso Slough proper). This is consistent with our estimates of very limited in-slough GPP (Fig. 5b); with Thébault et al. (2008), who found that LSB RSPs can be highly productive; and with the fact that oxygen and chlorophyll-a peaks in Fig. 9 are more pronounced following daylight high tides (Fig. 9a). We consider the example shown in Fig. 9, during four days in February 2016, to be qualitatively representative of typical dynamics during winter/spring blooms in the Alviso complex.
By iterating the backwards-looking window size for calculating the long-term chlorophyll-a concentration, we found that a 60-day window yielded the strongest predictive power for Alviso metabolic rates. This relatively large time window, and the manifestation of elevated low-tide chlorophyll-a on high-tide-centered rate estimates, is consistent with the hypothesis that Alviso Slough respiration is driven, in part, by the buildup of labile organic material on the slough bed.
The apparently strong effect of slough-bed oxygen demand on Alviso metabolic dynamics is consistent with the negative relationship between R night and tidal elevation (Figs. 7h and 10b), and with the association of the most negative NEP rates with lower high tides (Fig. 10a). Metabolic relationships with depth can be interpreted, in part, as representing the effect of varying relative habitat type with water level (Fig. 7k, l, w, x). Decreasing oxygen demand (1) it is the opposite of what we observe in Newark Slough (just a few kilometers away); and (2) we expect that marsh oxygen demand should be quite high (based on both the literature and our Newark Slough results). The combination of long-term chlorophyll-a concentration-which we consider an indicator of RSP biomass deposition onto the slough bed-and tidal elevation-which we use to represent the shifting proportion of marsh:slough metabolic contributions-explains 27% of the variability in observed respiration rates at Alviso Slough (Fig. 10b).
The metabolic rates at Newark and Alviso Sloughs are roughly an order of magnitude greater than what might be anticipated, on average, for tidal systems of their size (Nidzieko 2018), though the high degree of variability evident in Figs. 8b and 10b shows that these systems do have rates that are close to the regression line in that paper. The presence of the restored salt ponds could be one explanation for this pattern; the influence of primary production in sub-habitats on respiration rates in other nearby sub-habitats has been observed elsewhere in the San Francisco estuary system (Loken et al. 2021).
The distinct seasonal metabolic patterns in Newark and Alviso sloughs, and the pronounced differences between the drivers of these patterns, point to the outsize impact of RSP connections on in-slough oxygen dynamics. NEP in RSPconnected Alviso Slough was consistently more negative than in Newark Slough (Fig. 6), even though annual peak respiration rates in Newark coincide with annual respiration minima in Alviso (Fig. 5). The extremely high respiration rates in Alviso Slough are comparable to observations from warm-water estuaries hundreds of kilometers south along North America's Pacific Coast, and from the sub-tropical southeastern USA (Caffrey 2004). We use our understanding of the complex seasonal metabolic patterns to guide simple models of seasonal oxygen patterns.

Metabolic Impacts on Aggregate Dissolved Oxygen Dynamics
The drivers of metabolic patterns at Newark and Alviso sloughs readily translate to predictions of daily dissolved oxygen minima. In Newark Slough, daily DO minima, like metabolic rates, can largely be explained by water temperature and tidal elevation (Fig. 11a). Warm-water spring tides lead to the lowest observed oxygen levels in Newark Slough. At Alviso Slough, daily dissolved oxygen minima are primarily driven by RSP effluent (Fig. 11b), quantified as the mean effluent conditions from Pond A8 (we assume pond-to-slough outflow when dH∕dt < 0 and instrument depth < 1.7 m). Elevated biomass in the RSP effluent appears to draw down oxygen concentrations in Alviso Slough, consistent with the relationships with R night and NEP (Figs. 7i,j,u,v and 10b). But this effect is easily offset, or exacerbated, by oxygen concentrations in the pond; chlorophyll-a and dissolved oxygen concentrations can be highly correlated in the RSPs (Thébault et al. 2008), and we see the manifestation of this pattern in Alviso Slough (Fig. 3b, c, d). A comprehensive understanding of oxygen dynamics in RSP-connected sloughs may require a more detailed understanding of seasonal RSP biogeochemical dynamics. Salt ponds in other regions of the world exhibit similarly dynamic metabolic behavior (Spivak et al. 2020) and can display a high degree of spatial heterogeneity (Spivak et al. 2017). Natural and man-made salt ponds are common features of urbanized estuaries; characterizing the connectivity to and biogeochemistry of these complex habitats is essential to understanding and managing eutrophication risk.
While the interplay between physical and biogeochemical processes in these systems is interesting in its own right, this study illustrates an important aspect of scale and connectivity that is relevant to all estuaries where fringing habitats contribute to the dynamics in a larger embayment. In this case, the metabolic patterns elucidated by our results help to explain broader, system-scale patterns in oxygen dynamics. During summer spring tides, lower-low tides carry water with distinctively low DO through the mouth of LSB at the Dumbarton Narrows (Fig. 12). Assuming that Newark and Alviso are representative of at least parts of the complex lattice of marshland, tidal sloughs, and restored salt ponds at the LSB periphery, we can use our results to frame potential explanations for the summer spring tide pattern. Most simply, spring-tide higher-high tides occur at night during the summer and precede the lower-low tides during which we observe the DO minima (Fig. 12). During these spring-tide higher highs (about 3 m above MLLW), tidal marsh makes up 37% of the LSB surface area and 11% of the LSB volume (compared to 7.3% and 2.8% respectively at mean sea level). Thus, summer spring tides constitute a convergence of warm temperatures and extensive dark-hours marsh flooding that combine, in part, to drive system-scale LSB DO minima. While oxygen demand in Alviso Slough may be at a minimum during summer months, this RSP-connected sub-embayment still contributes high rates of respiration and strong net heterotrophy during the summer months (Figs. 5 and 6). The strong net heterotrophy in the RSP-slough complexes may set a lower oxygen baseline that ultimately contributes to the magnitude of the marsh-driven fortnightly pattern in low DO at the Dumbarton Narrows.
Both the LSB-and slough-scale DO patterns appear to be driven by a combination of natural and anthropogenic factors. Oxygen drawdown during nighttime marsh flooding is intuitively a natural pattern, but we do not know the extent to which WWTP nutrient loads may be contributing to over-marsh metabolism and, thus, magnifying this pattern. Distinct summer oxygen minima may be partially driven by natural marsh flooding patterns, but low baseline oxygen conditions in RSP-connected sloughs almost certainly contribute to the magnitude of these minima. Nutrient-enriched waters, driven into RSPs during flood tides (SFEI 2017), likely catalyze the high rates of in-pond biomass production that ultimately contribute to strong net heterotrophy and, thus, the lower baseline DO concentrations in the surrounding sloughs.
Disentangling the full picture of LSB oxygen dynamics will require simultaneously zooming out and zooming in. Habitat-specific hypsographic information could be used to inform a full-scale LSB oxygen balance. Such an approach would enable exploration of the impacts of shifting LSB habitat distribution on dissolved oxygen condition. Doing justice to this approach would require a more detailed characterization of the diverse metabolic patterns in the tidal marshes, sloughs, mud flats, and particularly restored salt ponds of the LSB, all of which are likely to be highly variable at even short time scales (Ganju et al. 2020), and may be further influenced by such nuanced factors as pore-water exchange.

Conclusions and Recommendations for Future Work
This study quantified metabolic rates in an estuarine system defined by strong tides and complex biogeochemical heterogeneity. The marsh-dominated slough, Newark, exhibited peak respiration rates and the most dynamic net ecosystem production rates during summer. The combination of warmer water and the coincidence of high-high tides with dark hours during summer may largely explain seasonal DO minima at this site. Alviso Slough, which is connected to restored salt ponds, exhibited distinctly different seasonal metabolic patterns. Unlike Newark Slough, which was regularly in metabolic balance, or even net autotrophic, Alviso Slough was consistently net heterotrophic, with NEP rates rarely exceeding -10 g O 2 ∕m 2 ∕d . Peak respiration rates, and peak heterotrophy, were associated with the spring bloom period when biomass outflow from the highly productive RSPs may severely increase in-slough benthic oxygen demand. Our results point to highly elevated benthic oxygen demand in Alviso Slough that may exceed that in the surrounding tidal marsh. However, the impact of these elevated rates on winter/spring oxygen concentrations was generally offset by the super-saturated DO concentrations in RSP effluent during the bloom periods. Alviso exhibited DO minima in summer, when RSP effluent was lower in DO, and when in-slough respiration rates were at an annual minimum but were still high (~10 g O 2 ∕m 2 ∕d).
Tidal marsh and RSPs dominate the landscape surrounding LSB's peripheral sloughs and are common features of estuaries globally. The metabolic dynamics in each of these habitat types offer insight into the mechanisms underlying dissolved oxygen patterns. This study highlights the high degree of seasonal variability in estuarine metabolic rates, the striking differences in seasonal metabolic patterns between proximal estuarine habitat areas, and the importance of considering the superposition of these seasonal and spatial effects in disentangling drivers of broader dissolved oxygen condition. Fully characterizing system-scale Fig. 12 Example tidal and dissolved oxygen patterns at the Dumbarton Narrows, where LSB connects to south SFB (summer 2017). Note that the distinct oxygen minima shown during spring tides represent every other tidal cycle (the lower-low tides) oxygen patterns, in the LSB and in any urbanized estuary, requires detailing the complex and superimposed biogeochemical mechanisms within individual pelagic and intertidal habitats. For the LSB in particular, we recommend two pointed lines of inquiry: (1) detailed characterization of the seasonal metabolic patterns inside of the restored salt ponds; and (2) direct rate estimates (core incubations and/ or eddy flux measurements) of benthic oxygen demand from RSP-connected sloughs, non-RSP connected sloughs, and tidal marshes. If benthic respiration rates in RSP-connected sloughs are shown to exceed those in the surrounding tidal marsh (as we have hypothesized from our results), that would point to anthropogenic drivers of depleted oxygen in the LSB and could have important implications for the management of WWTP nutrient loads and of the massive South Bay Salt Pond Restoration Project.