Physical and Biological Controls on the Annual CO₂ Cycle in Agua Hedionda Lagoon, Carlsbad, CA

Abstract

Agua Hedionda Lagoon (AHL), a tidal estuary located on the southern California coast, supports a diverse ecosystem while serving numerous recreation activities, a marine fish hatchery, a shellfish hatchery, and the largest desalination plant in the western hemisphere. In this work, a 1-year time series of carbon dioxide data is used to establish baseline average dissolved inorganic carbon conditions in AHL. Based on a mass balance model of the outer basin of the lagoon, we propose that AHL is a source of inorganic carbon to the adjacent ocean, through advective export, at a rate of 5.9 × 10⁶ mol C year⁻¹, and a source of CO₂ to the atmosphere of 0.21 × 10⁶ mol C year⁻¹ (1 mol C m⁻² year⁻¹), implying a net heterotrophic system on the order of 6.0 × 10⁶ mol C year⁻¹ (30 mol C m⁻² year⁻¹). Although variable with a range throughout the year of 80% about the mean, the ecosystem remained persistently heterotrophic, reaching peak rates during the summer season. Using results from the mass balance, the annual cycle of selected properties of the aqueous CO₂ system (pH, pCO₂, and CaCO₃ saturation state) were mathematically decomposed in order to examine the relative contribution of drivers including advection, ecosystem metabolism, and temperature that act to balance their observed annual cycle. Important findings of this study include the identification of advection as a prime driver of biogeochemical variability and the establishment of a data-based estimate of mean flushing time for AHL.

Introduction

Semi-enclosed coastal systems (e.g., estuaries, lagoons, intertidal zones, wetlands) are important to a wide number of natural processes in addition to societal and commercial uses (Ramesh et al. 2015). In many cases, the proximity to civilization makes the ecosystems operating within these coastal systems highly influenced by anthropogenic effects such as eutrophication from nutrient loading and elevated organic and sedimentary inputs and habitat loss due to land use change (Bauer et al. 2013; Howarth et al. 2011; Windham-Myers et al. 2018). In addition to these long-recognized issues, climate change (warming, acidification, deoxygenation) in the coastal ocean (Feely et al. 2008; Hauri et al. 2009; Gruber et al. 2012; Kessouri et al. 2021) and shallow coastal systems (Feely et al. 2010; Waldbusser and Salisbury 2014; Cai et al. 2021) is now a well-established area of research.

Over the past decade, the US Pacific Northwest has been recognized as a bellwether in the study of coastal acidification (Hales et al. 2005, 2017; Feely et al. 2008; Evans et al. 2011; Barton et al. 2012; Fairchild and Hales 2021). Of particular note, Barton et al. (2012) provided an unambiguous demonstration that oyster larval development was negatively impacted by low aragonite (the calcium carbonate mineral of which larval shell is made) saturation state of the intake water at a commercial shellfish hatchery. This finding, along with other growing evidence of organismal sensitivity to increasing dissolved CO₂ and acidification (Doney et al. 2020), has led to the widespread recognition of the need to better observe the coastal ocean carbonate system globally. There is a need to understand the carbonate system within semi-enclosed coastal lagoons and estuaries and their net fluxes to the coastal ocean and to the atmosphere. Existing studies have a wide range of flux estimates and are not always well constrained on seasonal to annual timescales (e.g., Borges 2005; Borges et al. 2006; Crosswell, et al. 2017; Cai 2011; Wang et al. 2016; Paulsen et al. 2018; Windham-Myers et al. 2018; Cai et al. 2021).

In parallel with OA research, coastal ocean observing systems are steadily growing in breadth and autonomy (Tilbrook et al. 2019; Barth et al. 2019); and coastal management and stewardship programs such as the National Estuarine Research Reserves (NERR) increasingly rely on observing system data to develop their strategies (NOAA Office for Coastal Management 2017). Among these is an initiative by NOAA to partner with a selection of shellfish growers in an effort to better understand the baseline conditions of the growers’ local lagoons and estuarine systems (Barton et al. 2015; Hales et al. 2017). A central piece of equipment that has been deployed in this effort is a continuous-flow multiparameter instrument developed at Oregon State University (OSU) (Hales et al. 2004; Bandstra et al. 2006; Fairchild and Hales 2021). This system characterizes the full suite of CO₂ chemistry parameters (partial pressure of CO₂—pCO₂; total dissolved inorganic carbon—DIC; pH; total alkalinity—A_T; and aragonite saturation state—Ω_Ar) either directly or by derivation using well-established thermodynamic relationships.

Agua Hedionda Lagoon (AHL) is one of six tidal estuaries located along the northern San Diego County coast (Beller et al. 2014), AHL is characterized as a low-inflow estuary (LIE). LIEs are estuaries where the total freshwater inflow is small, episodic, and/or seasonal. LIEs are found throughout Southern California (e.g., Southern California Wetlands Recovery Project 2018; Doughty et al. 2018), but also worldwide in regions with steep watersheds and/or Mediterranean climates (Largier et al. 1997; Largier 2010).

In this work, we analyzed and interpreted a full year (2018) of carbonate system observations at AHL, collected using the OSU continuous-flow system. Combined with salinity observations and calculated flushing times, carbonate chemistry data were used in a mass balance analysis of the outer basin to estimate the rate of inorganic carbon export and net ecosystem metabolism (NEM). In addition to presenting the variability of key CO₂ parameters, the results of the mass balance were used to perform a mathematical decomposition of pH, pCO₂, and Ω_Ar. We found the decomposition to be an instructive step in visualizing how the effect of each driver (temperature, advection, NEM, gas exchange) evolves throughout the year to determine the composite average and seasonal cycle of CO₂ in AHL. Given the need to understand semi-enclosed coastal carbonate chemistry, this analysis contributes to the growing body of knowledge on estuarine carbonate chemistry, but also adds specifically to understanding these dynamics within LIEs.

Materials and Methods

Study Site

AHL, Carlsbad, CA, comprises three interconnected basins, commonly referred to as the outer (connected to the ocean), middle, and inner basins (Fig. 1). The original wetland was converted into the present lagoon structure in 1954 by the Encina Power Station (EPS) and is maintained in its present form by semi-annual dredging. AHL consists of > 75% open water with the remainder being marsh and mudflats (Beller et al. 2014). Water depths range from very shallow (< 1 m) up to 14 m in certain areas, with an average depth of 8 m (for bathymetry, see Figs. 2–1 to 2–4 in Elwany et al. 2005). The ocean, connected by an inlet on the western side of the outer basin, dominates physical forcing in all three basins, with tidal lags of up to 4 h at the head of the lagoon where Agua Hedionda Creek enters (Jenkins and Wasyl 2006). The inner basin receives episodic freshwater input from Agua Hedionda Creek (and its effective catchment area; see below), predominantly in winter and spring during rain events. During the rest of the year, the creek is dry, and there is no appreciable freshwater entering AHL.

Fig. 1

Map of Agua Hedionda Lagoon showing the instrument location (red X) and additional features of interest including the desalination plant (DP), Encina Power Station (EPS), the Carlsbad Aquafarm (CAF), Hubbs-Seaworld Fish Hatchery (HFH), strawberry fields (SF), and recreational sports and boating activities (RS). The 3 basins (inner, middle, and outer) used in the box model are labeled. The inset shows California with the red box indicating San Diego County, where the lagoon is located. Image from Google Earth

Fig. 2

Observations from the shore station system including pCO₂ (a), DIC (b), temperature and salinity (c), pH (d), A_T (e), and Ω_Ar (f). The data shown are low pass filtered at 24 h (7 days for temp, thin lines) and 30 days (thick lines) for each parameter. Horizontal lines represent ocean averages for temperature 18.47 ± 3.0 C and salinity 33.55 ± 0.2 (panel c)

Located in an urbanized area, AHL is a highly utilized and popular destination for the Carlsbad community and tourists and provides a thriving ecosystem for diverse species of plants and animals. The two primary industrial features include the EPS and the Carlsbad Desalination Plant, both of which rely on water intake from the outer basin for once-through cooling of the power plant boilers and as desalination source water. The EPS was decommissioned Jan. 2019 but may have operated intermittently during the course of this study, drawing water from the outer basin. The desalination plant diverts from the EPS intake and its byproduct (brine) outflow is released through the EPS discharge channel into a small basin connected directly to the ocean, and adjacent to the outer basin (Fig. 1). The total intake volume of the EPS and desalination facilities is variable and at peak periods may represent a significant fraction of the daily tidal prism (Elwany et al. 2005; City of Carlsbad 2005). Other features of AHL include agriculture (primarily strawberry fields bordering the inner basin), the Hubbs Marine Fish Hatchery, and the Carlsbad Aquafarm (CAF)—a sustainable mussel and oyster farm, which operates in the outer basin. Both the fish hatchery and aquafarm (which grow calcifying organisms) rely on adequate flushing by the ambient ocean in order to maintain oxygen, pH, and calcium carbonate saturation state (Ω) above thresholds critical to growth. The CAF was one of the sites along the Pacific West Coast provided with the automated instrument (described above) through NOAA’s Integrated Ocean Observing System (IOOS) program for real-time monitoring of the CO₂ system by local stakeholders and shellfish growers. In collaboration with the IOOS regional associations (NANOOS, AOOS, CeNCOOS, and SCCOOS), the common goal of the instrument was to provide real-time information on the CO₂ content of the waters along the coastal ocean and its impact on carbonate organisms at the aquafarm intake.

Sampling

Continuous pCO₂ and DIC measurements were carried out using the automated instrument supplied by OSU. The instrument’s water intake line was positioned near the CAF facility docks (Fig. 1) approximately 1 m below the surface. The lagoon water was pumped through about 25 m of PVC piping to the instrument, which was housed in a small lab at the CAF. The incoming water is first filtered through a nylon screen T-strainer to remove large debris from the water, which then passes through a line containing a Honeywell 4905 conductivity sensor and a Honeywell Durafet pH sensor (both connected to a UDA2182) to obtain salinity, temperature, and pH measurements. Water then enters an enclosed headspace showerhead equilibrator, which contains a bubbling tube to facilitate equilibration of the CO₂ between the headspace and water (Fairchild and Hales 2021; Hales et al. 2004). The CO₂ gas is circulated to a LI-820A non-dispersive infrared (NDIR) detector to measure the mole fraction of the CO₂ (xCO₂) gas. At hourly intervals, the instrument switches to DIC mode and the sample water flows through a separate DIC sample line where it passes through a stainless-steel tangential flow filter to remove micron-sized particles from the sample water. The filtered sample water is then acidified with a 10% v:v dilution of concentrated hydrochloric acid and is passed through a mixing coil. CO₂ is extracted in a hydrophobic gas permeable membrane contactor (Liqui-Cel, G543), where the evolved CO₂ gas stream is dried and directed to the NDIR detector to determine the xCO₂ of the exiting gas stream, which is directly linked to the inflowing seawater DIC via a mass balance (Bandstra et al. 2006). Additional CO₂ system calculations for seawater pH, A_T, and Ω for aragonite and calcite are performed within the LabView program which also collects and stores all instrument data. The sample frequency (1 Hz) and data reduction interval (15 s) are user-specified, as is the interval for DIC analysis (hourly, in this work). This study reflects data collected during a 365-day period of minimal instrument interruption from Dec 6, 2017, to Dec 5, 2018.

Instrument Calibration

All CO₂ system calculations discussed below were performed using CO2SYS in Matlab (van Heuven et al. 2011).

The instrument’s automated calibration sequence using CO₂ gas and DIC liquid standards occurs every 6 h. The set of gas standards were purchased from Scott-Marin at nominal concentrations of 200, 800, and 1500 ppm (reported accuracy of 1% relative to NIST) stored in gas cylinders. The set of DIC liquid standards contained a solution mixture of ultrapure water (> 18 MΩ resistivity), NaHCO₃, oven-dried Na₂CO₃, and 0.1 M HCl, targeted to maintain near-ambient pCO₂ for the varying DIC. The standards were gravimetrically prepared at the Scripps Institution of Oceanography (SIO) in concentrations of 1900, 2100, and 2300 µmol kg⁻¹ every 12–14 days and stored in custom-made gas impermeable Mylar bags (IMPAK P75C0919). A density correction is performed to account for differences between the low salt background of the liquid standards and the seawater sample. As further described in Fairchild and Hales (2021), during a gas calibration sequence, a linear regression is performed in real time to verify the accuracy of the NDIR and apply a correction to the respective pCO₂ data to account for any offsets in the raw xCO₂ measurements. A stopped-flow ambient atmospheric pressure reading allowed the final correction of xCO₂ to pCO₂. During the gas standard sequence, an atmospheric CO₂ measurement was made from a separate line that extends outside (Fig. 2a), serving as an additional validation for the calibration procedure. The liquid standard sequence immediately followed the gas standards and also produced a linear regression in real time. While the gas standard regressions typically achieved high linearity with R² = 0.999 or better, the liquid standards experienced intermittent dropouts indicated by unreasonably low R² (i.e., < 0.98). This issue was traced to running liquid standard bags to near-empty, resulting in the introduction of residual air bubbles, trapped in the bags during filling. To address this issue, R² ≤ 0.997 were rejected and the remaining liquid calibrations were averaged monthly. The monthly averaged liquid standards were applied to the raw data in postprocessing to produce corrected DIC data. While the monthly average liquid standard calibration is not ideal and may add several tenths of a percent to the overall uncertainty of the DIC value reported, this approach eliminated ostensible errors on the order of > 5%, which otherwise would have been flagged as bad data and resulted in a loss of ~ 25% of the time series.

Bottle samples were collected from the common line of the instrument for validation of pCO₂, pH, and salinity. A check on salinity (accuracy ~ 0.03 salinity) was performed monthly using a benchtop density meter (Mettler-Toledo DM45). The recorded instrument salinity and salinity derived from density were compared and the difference between the two was applied to the instrument reading in software if the difference exceeded ± 0.05 in Practical Salinity. Validation samples for CO₂ were collected and analyzed for pH and DIC three times throughout the study (glass stoppered borosilicate bottles were poisoned with mercuric chloride and analyzed against CO₂ Certified Reference Materials following best practices from Dickson et al. 2007). Corrections to the continuous DIC and pH data were performed somewhat differently. Using what was deemed the most reliable pH bottle sample (due to a problem later identified in the benchtop spectrophotometric pH measuring system), a pH offset of − 0.18 was applied based on an April 2018 bottle sample comparison. The majority of this offset (0.13) is expected, due to the factory calibration of the pH sensor on the NBS scale, which must be adjusted external to the UDA2182 pH readout to obtain pH on the total hydrogen ion scale used in the CO₂ system calculations. Samples analyzed at the end of 2018 determined an insignificant drift in pH (< 0.005) and no additional pH offsets were applied. DIC laboratory analysis of three bottle samples did reveal an offset in the instrument DIC of 29 ± 16 µmol kg⁻¹ (n = 3), which was added to the final DIC dataset to account for the difference. The initial instrument pCO₂ measurements were determined to be compromised due to an air leak; and so, the final pCO₂ dataset analyzed in this study was derived in CO2SYS from QC’d DIC and pH data. The additional calculated parameters including A_T and Ω_Ar were recalculated in CO2SYS using corrected pH and corrected DIC. We acknowledge that some of the data are not of the quality we had aimed for and, accordingly, address the resulting limitations of the dataset using a sensitivity analysis recognizing significant uncertainties associated with both measurement and model. Further details of the QC of this dataset are provided in Shipley (2022).

Data Processing

The flow-through instrument DIC, pCO₂, pH, salinity, and temperature data were combined with other publicly available data including wind speed (adjusted from 4.9 m height of the anemometer to 10 m assuming neutral stability conditions, U₁₀ = U_4.9 × 1.1 m s⁻¹) and rain (q_r converted from cm⁻¹ h⁻¹ to m⁻¹ h⁻¹), obtained from the NOAA climate database for the weather station at McClellan-Palomar airport in Carlsbad, CA, approximately four miles from AHL (Vose et al. 2014). The high-frequency data (15 s time step pCO₂, pH, salinity, temperature) were first transformed using a 1-h moving average. All data were interpolated onto a common 1-h (on the hour) time stamp (corresponding to the frequency of the DIC measurements). The interpolation of the 1-h averaged high-frequency data is effectively a 1-h mean rather than a downsampling of raw data by decimation.

Data were processed using 30-day and 24-h low-pass filters (LPF, Matlab filtfilt zero phase filter) to isolate the annual cycle and examine higher frequency extremes, respectively (Fig. 2). The hourly time series undoubtedly contains additional information (see Fairchild and Hales (2021), e.g.), yet we chose to focus this analysis on the LPF data because (1) mixing in tidal systems becomes increasingly difficult to parameterize as the model approaches tidal timescales; and data at these timescales (e.g., flow rates, bathymetry) were not available, and (2) some of the data were determined to be compromised in terms of intermittent sensor failure and/or noise and the LPF is used to filter out these signals. In regard to the latter, instrument salinity exhibited pronounced amplitudes and short-term offsets in the absence of rain events during the first half of the time series. The 30-day LPF averages over these errors but consequently removes real signals. In addition, diel temperature amplitude measured by the instrument was unduly high. The source of this problem was a delivery line from the lagoon to the instrument, extended across the roof of a building where the line underwent excessive heating and cooling due to its contact with the atmosphere and direct radiative effects. Based on in situ measurements from SeapHOx sensors in the lagoon (Shipley 2022), we determined that the 24-h LPF instrument temperature agrees with the average in situ lagoon temperature, making the 24-h LPF temperature a reasonable approximation of the daily mean lagoon temperature. Because the seawater in the delivery line was never in contact with the atmosphere, the carbon parameters should not have been affected by the delivery line heating. Correcting the hourly instrument data to reflect in situ temperature is possible, given a continuous measurement of in situ temperature. A validated in situ temperature was not available throughout the full span of the 365 d time series presented here, and though beyond the scope of this analysis, such data, if available, may provide additional information for further analysis of this dataset.

In the following sections, we describe a multi-box model that uses the 24-h LPF data to estimate water exchange rate with the ocean and a single-box model that uses the 30-day LPF data to evaluate biogeochemical rates. While the middle and inner lagoons are not constrained by measurements, it is feasible to include their volumes in the mixing model used to estimate the water exchange rate. Due to a lack of spatial carbonate system data, it is not possible to evaluate inter-basin carbon exchange; and for this reason, the biogeochemical model is set up as a single box reflecting the outer basin only. When reporting results on a m⁻² basis, this is effectively the same as treating the entire lagoon as a single box, constrained by data measured in only one basin. The underlying assumptions of the single-box biogeochemical model are that (1) biogeochemical rates are similar across all three basins and/or (2) flushing time is sufficient across the entire lagoon such that inter-basin gradients approach zero over the 30-day LPF time frame. When these conditions are met, the dominant advective term is between the ocean and the outer basin. Both models are executed on a 1-h time step, corresponding to the common time stamp mentioned above. The mixing model uses the 24-h LPF time series and the biogeochemical model uses the 30-day LPF and the 24-h LPF the gas exchange term.

Mixing Model

To model the annual CO₂ cycle based on the 24-h LPF time series, it is necessary to first approximate the average rate of water exchange between the lagoon and the ocean. In order to accomplish this, we developed a mixing model (Fig. 3) to interpret observations of salinity during a select period of high freshwater input where the instrument captured multiple cycles of freshening followed by seawater dilution. This analysis focuses on a 1-month period during Dec 2018 where two significant rain events occurred while the salinity measurement was functional and accurate. Rain events also took place during the previous winter of 2017/2018 (Fig. 2c), but, as noted above, this period also coincided with noisy salinity data that would not facilitate a mixing model. Dimensions used in the mixing model are given in Table 1.

Fig. 3

(Top) Mixing model used to determine flushing time based on salinity observations during Nov 2018. Exchange between the ocean and three basins of the lagoon is driven by tidal mixing and basin geometry. Freshwater input is parameterized based on rain and effective catchment area. Flow represented by gray arrows was not parameterized, but its effect is included in observed data that are used to constrain the model (see text and Table 1). (Bottom) Biogeochemical model used to estimate NEM, based on dissolved inorganic carbon data. NEM is the balance of observed changes, gas flux, and advection

Table 1 Input and output constants and variables used in the 3-box mixing model. V represents average box volume (mean tide height added to areas and volumes at MLLW taken from Elwany et al. 2005)

During Dec 2018, mixing was approximated using the tidal prism equation (Monsen et al. 2002)

$${T}_{f}=\frac{{V}_{OB}\times {T}_{period}}{\left(1-R\right){V}_{prism}}$$

(1)

where T_f is flushing time, V_OB is the average volume of the outer basin, T_period is the tidal period (set to 24 h here as this is the largest component even though tides are mixed semidiurnal), V_prism is the tidal prism of the outer basin, and R is the return flow factor that accounts for several processes, some poorly constrained, including stratification, incomplete mixing of the tidal prism and basin volumes, and intake volumes of the EPS and desalination plant. V_prism is calculated as the mean daily tidal range (0.87 m) multiplied by the mean area of the outer basin (2 × 10⁵ m²). The return flow factor is the average of flood and ebb return factors: R = (R_F + R_E)/2, which are determined using a 3-box model of the lagoon and an ocean boundary condition. The salt budget is

$$\frac{d{V}_{b}{S}_{b}}{dt}={q}_{in}{S}_{in}-{q}_{out}{S}_{b}$$

(2)

where V_b and S_b are the respective volume and salinity of the basin (b) in question with b = OB (outer basin), MB (middle basin), or IB (inner basin). For the outer basin during a flood tide S_in = S_OC (ocean). Flow rates (q_in, q_out) and their corresponding salinity (S_in, S_b) are the tide-phase-dependent volumetric inflow (in) and outflow (out) rates and salinity of the inflowing and outflowing water, respectively. For rising tide, the inflowing conditions are those of the seaward reservoir; for ebbing tide, they are for the landward reservoir. Since q = dV/dt, from Eq. 2, the salt change for a given time interval, dt, is

$$d{V}_{b}{S}_{b}={dV}_{in}\left({S}_{in}\right)-{dV}_{out}\left({S}_{b}\right)$$

(3)

which can be further discretized as

$${V}_{b,t}{S}_{b,t}-{V}_{b,t-1}{S}_{b,t-1}={dV}_{in}\left({S}_{in}\right)-{dV}_{out}\left({S}_{b}\right)$$

(4)

where subscript “t” represents the time step. Rearranging Eq. 4 gives salinity at time t as

$${S}_{b,t}=\frac{{{V}_{b,t-1}{S}_{b,t-1}+dV}_{in}\left({S}_{in}\right)-{dV}_{out}\left({S}_{b}\right)}{{V}_{b,t}}$$

(5)

Thus, for the outer basin, including the return flow factor terms and expressing V_b,t as a summation of V_b,t-1 and volume gained (including a term for runoff, V_r) and lost, the flood and ebb tide for a given time step are

$${S}_{OB,t}\left(flood\right)=\frac{{S}_{OB,t-1}{V}_{OB,t-1}+{dV}_{OB,t}\left(1-{R}_{F}\right){S}_{oc,t-1}+{dV}_{ob,t}{R}_{F}{S}_{ob,t-1}-{dV}_{mb,t}{S}_{mb,t-1}}{{V}_{ob,t-1}^{1}+{dV}_{ob,t}{- dV}_{mb,t}+{V}_{r}}$$

(6a)

$${S}_{OB,t}\left(ebb\right)=\frac{{S}_{OB,t-1}{V}_{ob, t-1}-{dV}_{ob,t}{S}_{ob,t-1}+{dV}_{mb,t}{(1-R}_{E}){S}_{mb,t-1}+{dV}_{mb,t}{R}_{E}{S}_{ob,t-1}}{{V}_{ob,t-1}-{dV}_{ob,t}{+ dV}_{mb,t}+{V}_{r}}$$

(6b)

where dV terms are calculated from changes in tide height and areas given in Table 1 and tidal lag is not considered. For clarity, the sign convention in E. 6 uses |dV|, resulting in a positive sign preceding a gain and a negative sign preceding a loss. Similar equations (not shown) are used for the middle and inner lagoons. Note that because we used the mean basin areas, we have simplified the model to boxes, not accounting for sloping side walls, which is an appropriate simplification to determine a mean flushing time. Return flow factors are only considered for the outer basin (i.e., for the middle and inner lagoon boxes R_F = R_E = 0). We found that tuning additional return flow factors for the middle and inner basin did not provide improvement, and furthermore, the observations being limited to the outer basin provides little constraint on adjacent lagoon return flow parameterizations.

Freshwater input (V_r, Eq. 6a, b) is based on local precipitation data. Though it is not a gauged stream and groundwater inputs were not quantified, we believe, based on empirical observations (e.g., hypersaline summer conditions in the inner basin, lack of visible creek flow during summertime visits) and adjacent watersheds, that freshwater input from Agua Hedionda Creek is mainly associated with rain events, and therefore, runoff is computed from hourly rainfall data with freshwater input otherwise set to zero. Rainfall (q_r) is multiplied by the effective catchment area (A_catch) and given a delay time of 1 day (established by matching modeled results to the phase of the freshwater spikes in observed salinity) to allow for runoff.

$${V}_{\mathrm{r}}={q}_{\mathrm{r},\mathrm{t}-1\mathrm{d}}\times {A}_{\mathrm{catch}}$$

(7)

Since the effective catchment area may depend on factors such as storm patchiness and soil permeability, V_r is treated as a tunable parameter (through CF below) along with R_F and R_E. The total watershed of Agua Hedionda Creek is 76 km², representing approximately 80 × the area of Agua Hedionda Lagoon and ~ 14% of the greater Carlsbad Hydrologic Unit which is divided among four separate lagoons (City of Carlsbad 2021). The effective catchment area A_catch was allowed to vary as a multiple (CF or catchment factor) of the area of AHL where A_catch must represent less than the 80 × watershed/lagoon area.

$${A}_{\mathrm{catch}}=\mathrm{CF}\left({A}_{\mathrm{OB}}+{A}_{\mathrm{MB}}+{A}_{\mathrm{IB}}\right)$$

(8)

Tuned output values (Table 1) and evaluation of uncertainty in T_f from ancillary data are discussed in “Results” and “Discussion.”

Biogeochemical Model

In the biogeochemical model of the outer basin (Fig. 3), a mass balance is estimated for the annual time series using

$$\frac{\mathrm{dDIC}}{\mathrm{dt}}h={F}_{\mathrm{adv}}+{F}_{\mathrm{gas}}+\mathrm{NEM}+\mathrm{NEC}$$

(9)

where F_adv is the advective flux, F_gas is the flux due to gas exchange with the atmosphere, NEM is the net ecosystem metabolism (i.e., gross primary production minus autotrophic plus heterotrophic respiration), NEC is the net ecosystem calcification (i.e., gross calcification minus dissolution), and the term on the left is the finite difference in observed DIC for each time step multiplied by water depth (h).

The 30-day LPF A_T exhibits no distinct annual signal and the range of 46 µmol kg⁻¹ corresponds to a range of 23 µmol kg⁻¹ in DIC (less than 20% of the observed variability), which is similar to the uncertainties discussed below. Normalizing the 30-day LPF A_T (Fig. 2e) with 30-day LPF Salinity (Fig. 2c) (nA_T = A_T × $\overline{\mathrm{S} }$/S) significantly increases (roughly doubles) the range of the nA_T compared to A_T while adding questionable signals to the nA_T, likely due to limitations in the salinity data discussed above. Due to the absence of a significant annual cycle in total alkalinity attributed to NEC, this term was excluded from the biogeochemical model. Due to the interpolation interval of 1 h (see above), all of the terms used in Eq. 9 are on an hourly time step even though the observations and advection terms have been passed through a low-pass digital filter. Because the mass balance is framed in carbon units from the perspective of bulk seawater composition, a gain (loss) of carbon is reflected as positive (negative) such that ingassing is positive, outgassing is negative, net heterotrophic NEM is positive, and net autotrophic NEM is negative. Air-sea CO₂ gas exchange (F_gas, mol m⁻² h⁻¹) is calculated using

$${\mathrm{F}}_{\mathrm{gas}} ={k}{\prime}\times {\mathrm{U}}_{10}^{2}\times \Delta {\mathrm{pCO}}_{2}/(365\times 24)$$

(10)

where U₁₀ is wind speed (adjusted to 10 m height) and the coefficient k′ = 7.7 × 10⁻⁴ (mol µatm year⁻¹) (m s⁻¹)⁻² incorporates multiple gas transfer terms (Schmidt number and CO₂-gas solubility) representative of constant ocean T and S (Wanninkhof 2014). The ΔpCO₂ is the difference between the measured atmospheric pCO₂ and pCO₂ derived from pH and DIC measurements. Though near-negligible in comparison to F_adv, the error in ΔpCO₂ and k′ (including variability in temperature and salinity) is assessed in the uncertainty analysis below.

The advective flux of inorganic carbon, F_adv, is calculated from the gradient between the outer basin (OB) and ocean using

$${\mathrm{F}}_{\mathrm{adv}}=\mathrm{Q}\times \left({\mathrm{DIC}}_{\mathrm{OB}}-{\mathrm{DIC}}_{\mathrm{ocean}}\right)$$

(11)

where the average flow rate (Q) is the flushing time (T_f) determined by the mixing model combined with the average volume of the outer basin:

$$Q=\frac{{V}_{\mathrm{outer}}}{{T}_{\mathrm{f}}}$$

(12)

DIC_OB is directly measured by the instrument in the outer basin and the ocean endmember, DIC_ocean, is derived from 2 years (2017–2018) of data from the La Jolla Shores Time Series (LJTS) measured at the surface in shallow near-shore waters (7–40 m isobath) 35 km south of AHL (Kekuewa et al. 2022). The 2 years of data collected by Kekuewa et al. (2022) were used to compile a monthly average climatology of DIC and A_T, and based on agreement with DIC estimated from a pH sensor deployed in the near-shore environment directly adjacent to AHL, an uncertainty in DIC_ocean of ± 21 µmol kg⁻¹ was assigned for propagation of error calculations (Fig. S1). Based on the absence of a distinguishable annual cycle in the 30-day LPF A_T, and no significant difference between the LJTS surface A_T (2239 ± 11 µmol kg⁻¹) and 30-day LPF A_T in AHL (2232 ± 12 µmol kg⁻¹, see “Results” below), advection of A_T is not parameterized in the biogeochemical model.

Excluding NEC and rearranging the mass balance,

$$\mathrm{NEM}=\frac{{\mathrm{dDIC}}_{\mathrm{OB}}}{\mathrm{dt}}h-{F}_{\mathrm{adv}}-{F}_{\mathrm{gas}}$$

(13)

A sensitivity analysis was employed to estimate uncertainty following the approximate method of numerical differentiation (Kragten 1994; Ellison and Williams 2012). In this analysis, errors for five input terms (Table 2) are estimated and propagated through Eqs. 10–13 using the approximation.

Table 2 Biogeochemical model inputs and outputs. U₁₀ and ΔpCO₂ are hourly measurements, DIC_OB is 30-day LPF. Negative output values represent a DIC loss from the lagoon. Range indicated by ± represents estimated uncertainty for inputs and combined standard uncertainty for outputs

$$u(y,x_i)\approx y\left(x_l,x_2...(x_i+u\left(x_i\right)\right)...x_n)-y\left(x_l,x_2,...x_l...x_n\right)$$

(14)

Propagated errors are then used to compute combined standard uncertainty for F_gas, F_adv, and NEM, respectively, using

$$u\left({F}_{gas}\right)=\sqrt{{\left(u\left({F}_{gas},\Delta {\mathrm{pCO}}_{2}\right)\right)}^{2}+{\left(u\left({F}_{gas},{k}{\prime}\right)\right)}^{2}}$$

(15)

$$u\left({F}_{adv}\right)=\sqrt{\begin{array}{c}{\left(u\left({F}_{adv},{T}_{f}\right)\right)}^{2}+{\left(u\left({F}_{adv},{DIC}_{OB}\right)\right)}^{2}\\ +{\left(u\left({F}_{adv},{DIC}_{ocean}\right)\right)}^{2}\end{array}}$$

(16)

$$u\left(NEM\right)=\sqrt{\begin{array}{c}{\left(u\left(NEM,{T}_{f}\right)\right)}^{2}+{\left(u\left(NEM,{DIC}_{OB}\right)\right)}^{2}+{\left(u\left(NEM,{DIC}_{ocean}\right)\right)}^{2}\\ +{\left(u\left(NEM,\Delta {\mathrm{pCO}}_{2}\right)\right)}^{2}+{\left(u\left(NEM,{k}{\prime}\right)\right)}^{2}\end{array}}$$

(17)

and reported as ± values in the output rows of Table 2. The gas exchange uncertainty was based on a 30 µatm error in ΔpCO₂ (corresponding to a 0.02 error in the pH sensor) and a 25% error in k′. Advection error was estimated from a combined error of T_f ± 2 days (see “Results” below) and errors of 16 and 21 µmol kg⁻¹ in DIC_OB and DIC_ocean, respectively. The error estimate for DIC_OB encompasses spatiotemporal sampling error between bottle and sensor in addition to errors due to the accuracy of the liquid standard which normally are lower than 10 µmol kg⁻¹ (Fairchild and Hales 2021). The error in DIC_ocean corresponds to a standard deviation reported by Takeshita et al. (2015).

Driver Decomposition

The results from the biogeochemical model were used to assess the change in observed parameters (Fig. 2) due to each “driver”: temperature, advection, gas exchange, and NEM. The effect of each driver on CO₂ parameters throughout the year was assessed by holding all other variables constant while allowing one to vary. Williams et al. (2018) applied this approach using the equations (with pH as the example):

$${\Delta \mathrm{pH}}_{\mathrm{t}}^{\mathrm{temp}}=\mathrm{pH}\left({{A}_{T}}_{\mathrm{t}0},{\mathrm{DIC}}_{\mathrm{t}0},{\mathrm{S}}_{\mathrm{t}0},{\mathrm{T}}_{\mathrm{t}}\right)-{\mathrm{pH}}_{\mathrm{t}0}$$

(18)

$${\Delta \mathrm{pH}}_{\mathrm{t}}^{\mathrm{adv}}=\mathrm{pH}\left({{A}_{T}}_{\mathrm{t}0},{\mathrm{DIC}}_{\mathrm{t}}^{\mathrm{adv}},{\mathrm{S}}_{\mathrm{t}0},{\mathrm{T}}_{\mathrm{t}0}\right)-{\mathrm{pH}}_{\mathrm{t}0}$$

(19)

$${\Delta \mathrm{pH}}_{\mathrm{t}}^{\mathrm{gas}}=\mathrm{pH}\left({{A}_{T}}_{\mathrm{t}0},{\mathrm{DIC}}_{\mathrm{t}}^{\mathrm{gas}},{\mathrm{S}}_{\mathrm{t}0},{\mathrm{T}}_{\mathrm{t}0}\right)-{\mathrm{pH}}_{\mathrm{t}0}$$

(20)

$${\Delta \mathrm{pH}}_{\mathrm{t}}^{\mathrm{NEM}}=\mathrm{pH}\left({{A}_{T}}_{\mathrm{t}0},{\mathrm{DIC}}_{\mathrm{t}}^{\mathrm{NEM}},{\mathrm{S}}_{\mathrm{t}0},{\mathrm{T}}_{\mathrm{t}0}\right)-{\mathrm{pH}}_{\mathrm{t}0}$$

(21)

where the driver in each equation is represented as an array of time, t, while all other terms are scalars equal to the first point in the time series, t0. The integrated effect of advection, gas exchange, and NEM on DIC (${DIC}_{t}^{adv}$, ${DIC}_{t}^{gas}$, and ${DIC}_{t}^{NEM}$) have been determined previously in the mass balance. Treating the time series as a closed system and using the driver terms as shown in Eqs. 18–21 works reasonably well when a given DIC driver: (1) exhibits a cycle resulting in little net change annually, and (2) varies on the same order as the observed DIC, in this case ~ 130 µmol kg⁻¹ (ca. 2100 to 1970, Fig. 2b). In a system such as AHL, where large uni-directional terms (i.e., advection and NEM) balance a much smaller observed annual signal, the decomposition calculation cannot be treated as a closed system as in Williams et al. (2018). This would effectively allow Eqs. 19 and 21 to freely drive DIC outside of the observed range, which is problematic due to the nonlinear nature of the CO₂ system. For example, converting integrated F_adv into an accumulated annual change in DIC results in a DIC decrease of about 3600 µmol kg⁻¹ which represents more DIC than present in the system and would reduce the initial DIC_t0 of 2048 µmol kg⁻¹ to a near-zero value after ~ 6 mo; whereas the actual observed DIC never falls below 1970 µmol kg⁻¹. Partial derivatives (aka Revelle Factors, Egleston et al. 2010; Frankignoulle 1994) of the CO₂ system are far different at near zero vs. the observed range in this work (DIC = 2100 – 1970) resulting in unrealistic results when using Eqs. 19 and 21. A more realistic approach to driver decomposition in a system such as AHL requires linearizing Eqs. 19 and 21 from partial derivatives representative of the local temperature, salinity, alkalinity, and observed range of DIC. Here, the linearization is performed over a 100 µmol kg⁻¹ range of DIC (2000–2100) at temperature and salinity at time = 0 to obtain

$$\frac{\partial \mathrm{pH}}{\partial \mathrm{DIC}}=0.0023; \frac{\partial {\mathrm{pCO}}_{2}}{\partial \mathrm{DIC}}=3.1; \frac{\partial {\Omega }_{\mathrm{Ar}}}{\partial \mathrm{DIC}}=0.0094;$$

and Eqs. 19 and 21 become

$${\Delta \mathrm{pH}}_{\mathrm{t}}^{\mathrm{adv}}={\frac{\partial \mathrm{pH}}{\partial \mathrm{DIC}}\mathrm{DIC}}_{\mathrm{t}}^{\mathrm{adv}}-{\mathrm{pH}}_{\mathrm{t}0}$$

(22)

$${\Delta \mathrm{pH}}_{\mathrm{t}}^{\mathrm{NEM}}={\frac{\partial \mathrm{pH}}{\partial \mathrm{DIC}}\mathrm{DIC}}_{\mathrm{t}}^{\mathrm{NEM}}-{\mathrm{pH}}_{\mathrm{t}0}$$

(23)

Results

Observational Data

The measured and derived physical (temperature and salinity) and CO₂ parameters from December 2017 to December 2018 are shown in Fig. 2. Each panel shows 24-h and 30-day LPF processed data in addition to the atmospheric pCO₂ reading by the instrument (panel a, a 1 min average of 1 Hz pCO₂). Temperature exhibits an annual amplitude of ~ 8 °C with a peak in Aug; while salinity remains relatively stable over the annual cycle, except during major rain events, in particular during January and December (Fig. 2c). The LPF pCO₂ ranges from ~ 400 to 800 µatm with a maximum in the spring and a minimum in the summer (Fig. 2a). The pH ranges from 7.7 to 8.1 annually, showing an inverse correlation with pCO₂ (Fig. 2d). The annual cycle of both properties is convoluted with the competing effects of transport, NEM, and temperature. The local ambient atmospheric CO₂ signal fluctuates from 400 to 500 μatm over the year and does not surpass AHL pCO₂, except during very short-lived spikes where a local boundary-layer contamination (e.g., from vehicle exhaust) may interact with the system’s intake port. DIC peaks in early spring to a maximum of 2100 μmol kg⁻¹ and reaches a minimum of ~ 1975 μmol kg⁻¹ by summer, resulting in an annual amplitude of ~ 125 μmol kg⁻¹. Similar to salinity, A_T (Fig. 2e), changes are restricted primarily to major rain events (Fig. 2c, e). During a pronounced rain, salinity may decrease to as low as 30 over a timescale of 1–2 days, resulting in only small decreases in the 30-day LPF salinity. Seawater saturation state with respect to aragonite (Ω_Ar) (Fig. 2f) ranges between roughly 2 and 3, with short-term excursions that may briefly approach 1.5 and > 3. An Ω_Ar < 1 is indicative of undersaturated waters and is considered unfavorable to calcifying (aragonite-forming) organisms (Harris et al. 2013), and while AHL surface waters do not go below this threshold, there are periods, especially during the winter months where Ω_Ar briefly reaches levels close to undersaturation.

While the LJTS A_T is indistinguishable from A_T observed in AHL (see above), the DIC_ocean (derived from the LJTS) is persistently lower than DIC_OB with a similar or slightly greater annual range, which is driven primarily by a bloom triggered by seasonal upwelling (Fig. 2b, see Kekuewa et al. 2022). Use of the LJTS climatology to represent DIC_ocean at a location 35 km north is justified as follows. On the one hand, this approach is supported by a concurrent 3.5 mo. time series of SeapHOx data (Fig. S1) where pH has been converted to DIC using the average alkalinity from the LJTS. The two DIC values clearly track (with an average deviation of 21 μmol kg⁻¹) during a period overlapping the 2018 study period discussed here. However, the sensor time series does not capture the most significant month of the LJTS (July) where in both 2017 and 2018, DIC was drawn down to its lowest levels (Fig. S1). We therefore emphasize both here and in the discussion below that the DIC_ocean term may contain additional errors not captured by the propagation of the ± 21 µmol kg⁻¹ uncertainty used in our analysis.

AHL-Ocean Exchange Rate

The three adjustable mixing model parameters CF, R_F, and R_E were tuned using the MATLAB function “fmincon” (interior-point method), resulting in values of CF = 21, RF = 0.0, and RE = 0.95 and RMSE = 0.14 in salinity (Fig. 4). Based on these results, an average R = 0.48 was used in Eq. 1 to compute a mean flushing time of T_f = 7.3 days. As explained in the “Discussion,” the biogeochemical model used T_f = 7 ± 2 days.

Fig. 4

Results of the mixing model compared to salinity observed in the outer basin during two major rain events in Dec 2018. Mixing parameters are tuned using a 24-h LPF of both model and data. Model results here show the best fit using T_f = 7.3 days, R = 0.48

Biogeochemical Rates

Areal flux and the associated depth + time integrated changes in each driver are shown in Fig. 5a and b, respectively, and reported in Table 2. Although the observed changes in DIC (ΔDIC) are close to balanced over the 12-month period of this study, the most significant drivers are far from balanced (Fig. 5b). The negative flux due to gas exchange represents a continuous source of CO₂ to the atmosphere throughout the year (on average 2.3 mmol C m⁻² day⁻¹, Table 2), which is evident in Fig. 2a where lagoon pCO₂ never falls below atmospheric levels. Advection represents a significant loss of DIC annually, exporting 5.9 × 10⁶ mol C year⁻¹ from the lagoon to the ocean. Similar to the gas flux, the characteristics of advective flux are also evident in Fig. 2b, where DIC_OB is persistently higher than DIC_ocean. The defined mass balance (Eq. 9) is closed by the NEM (Eq. 13), which counters the other directly calculated terms. As the advective term is the largest, NEM must necessarily counter that loss, resulting in a calculated NEM that largely mirrors F_adv (Fig. 5). The resulting NEM exhibits net heterotrophy and net advective export to the ocean during all 12 months, leading to an annualized NEM of 82 ± 35 mmol C m⁻² day⁻¹ (Table 2). The results of the sensitivity analysis (Table 2) highlight the importance of errors in DIC_OB, DIC_ocean, and T_f, which make up 31%, 40%, and 30% of the error in F_adv and 30%, 39%, and 29% of the error in NEM. Uncertainty in F_gas contributes less than 3% to the uncertainty in NEM.

Fig. 5

a Areal flux from mass balance (Eq. 13) and b the time and depth-integrated changes with error shown for each driver. NEM and advection (blue, left axes), and gas exchange and observed change (black, right axes). Note that the right axes cover a significantly smaller range of values than the left on both a and b

Driver Decomposition

The decomposition results (Fig. 6a–c) represent the accumulated effect on pH, pCO₂, and Ω_Ar throughout the year due to the “drivers” temperature, advection, gas exchange, and NEM. Because accumulated values shown in Fig. 6 are tallying changes of a near-stationary time series of DIC with large but balanced underlying drivers, the results may appear abnormal at first. Advection and NEM drove opposing accumulated changes of roughly ± 7 in pH, ± 10,000 in pCO₂, and ± 30 in Ω_Ar, representing a roughly tenfold greater effect than gas exchange and temperature (compare left vs. right axis ranges).

Fig. 6

Driver decomposition for pCO₂ a, pH b, and Ω_Ar c, following Eqs. 18–21 for each driver. Blue lines correspond to the left y-axis showing NEM and Advection. Black lines correspond to the right y-axis, showing the observed, gas exchange, and temperature data, which is an order of magnitude smaller than the left y-axis

Discussion

Physical Controls and AHL-Ocean Exchange Rate

Although the simple 3-box mixing model based on Eq. 6a, b (Fig. 3) does not adequately capture the tidal range in salinity, the underlying mixing trend (Fig. 4) representing the flushing time is the most important result for operating a biogeochemical model of the 1-year time series of LPF data. In tuning the mixing model, we used a 24-h LPF of salinity observed in the outer basin during a 1-month period where several rain events occurred (Nov 24–Dec 24); extending just beyond the 365 d DIC time series. During this period, the CO₂ instrument experienced a failure, but the salinity sensor continued to operate, capturing valuable information on the physical mixing of the lagoon. The multiplier for runoff (CF) represents an area 21 × greater than the lagoon itself, which seems somewhat reasonable given a watershed area of 80 × greater than the lagoon. The flood return flow factor (R_F) of 0 indicates that the outgoing plume from the lagoon is swept away from the mouth by the open ocean before the tide comes back in. This result also seems appropriate as the inlet of the outer basin is directly connected to the ocean where an alongshore flow often persists (Nam and Send 2011). The very high ebb tide return flow factor of 0.95 is not surprising, as the much smaller middle basin receives all of the outer basin flood tide (Fig. 1). The difference between a zero flood return flow factor and a non-zero ebb factor is also consistent with a well-mixed box model potentially being more representative during the flood tide, consistent with flood tides typically exhibiting more vertical mixing and less stratification in many estuarine systems (e.g., Simpson et al. 1990; Geyer and MacCready 2014).

The flushing time determined by the mixing model is nearly twice what one would estimate based purely on the tidal prism. For example, the average 24-h tidal range throughout the time series is 1.6 m. Combining this value with the area of the outer basin and the average volume of the outer basin (Table 1) gives a flushing time of only 3.8 days. Similar approaches based strictly on tidal prism, but accounting for temporal variability in the tidal stage, estimate flushing times across the different boxes ranging from 1 to 5 days but with roughly the same average of ~ 3 days (Elwany et al. 2005). As discussed by Monsen et al. (2002), it is not uncommon for flushing time approximated by basin geometry to differ from that derived from actual measurements of bulk properties such as salinity (or other closely related concepts such as residence time and age) by a factor of two or more. Moreover, recent studies utilizing high-resolution numerical models have shown that flushing time estimates can vary dramatically depending upon the method (Lemagie and Lerczak 2015) and that simple geometric flushing time estimates typically underpredict true residence time due to incomplete vertical mixing and return flows (MacCready et al. 2021).

Evidence suggests that flushing times determined from bulk properties during rain events (as done in Fig. 4) may underestimate average flushing time (Alber and Sheldon 1999). Exchange flow calculations from velocity measurements conducted during Spring 2016 (up to 0.1 m s⁻¹ during periods with freshwater input), combined with the dimensions of the full lagoon (Table 1) and an approximate cross-sectional area of ~ 45 m² at the mouth, suggest T_f in the range of ~ 6–14 days with shorter values during periods with freshwater input and longer values during extended periods without rainfall. Thus, the flushing time determined from salinity observations during a rain event in 2018 (Fig. 4) is in good agreement with that determined from velocity observations during a longer period. Combining the two observation-based estimates, we selected the range T_f = 5–9 days for use in the biogeochemical model, with 7 days (corresponding to a flow rate Q = 7.7 × 10³ m³ h⁻¹ in Eq. 12) as the central value used to report averages. However, we note that in extended periods without freshwater input, different mechanisms dominate exchange and flushing times could be longer. In summary, flushing time could be characterized in a number of different ways for this particular time series, but based on observations of salinity and velocity, we suspect that our estimate represents a more accurate value than that obtained based on geometry alone (Lemagie and Lerczak 2015; MacCready et al. 2021).

As noted above, a single-box treatment of the biogeochemical data requires a flushing time of the full lagoon sufficient to justify the assumption that inter-basin gradients approach zero over the averaging period of 30 days. The exchange flow calculations suggest that total lagoon flushing may meet this requirement. Although pH sensors were not located in all lagoons throughout this study, oxygen data measured during a 3-week period of this work support the assumption that inter-basin gradients between the middle and outer lagoon are small relative to the gradient between the ocean and outer basin (Fig. S2).

Biogeochemical Controls

Following its peak in March, the observed DIC begins a 6-month decrease until August (Figs. 2b and 5b; driving a change in sign of the dDIC/dt in Fig. 5a). This decrease might be separated into two periods, the first consisting of a gradual drop from March to late May followed by a second period of rapid decrease from June to Aug (Fig. 5b). The first period corresponds to the highest rates of gas exchange (Fig. 5a) under the largest air-sea gradient in pCO₂ (Fig. 2a) and moderately high lateral gradient in DIC (Fig. 2b). Other features, such as a drop in DIC in Jan, are likely associated with rain events rather than the annual cycle of the lagoon. Based on the mass balance, the lagoon remains net heterotrophic throughout the year, with maximum heterotrophic rates in the summer and minimum in the winter (Fig. 5a).

The advective export of DIC from the lagoon (Fig. 5b) results in a very large forcing to drive pH and Ω higher and pCO₂ lower (Fig. 6), which must be balanced by equally high respiration through the NEM. Because NEM is the mass balance constrained counter to the other processes, and advection is by far the greatest driver, NEM mirrors advection in all cases. It is important to note that the drivers are not mutually exclusive and may affect each other in complex interactions making it difficult to project the effect of a single driver in the absence of others. For example, advection has been linked mechanistically to rates of benthic processes including metabolism and CaCO₃ dissolution (Eyre et al. 2014). Similarly (but likely for different reasons), our biogeochemical model implies a strong link between NEM and advection, but as noted above, observations of A_T indicate little effect on NEC, possibly due to a small scale of calcification/dissolution (limited to CAF activities and not ecosystem-wide processes) relative to the exchange rate with the open ocean. Thus, extending these results to exclude one driver while assuming the others remain unchanged is speculative.

Similar to most of the open ocean, the system described in this study experiences an annual CO₂ cycle (Takahashi et al. 2002). The annual cycle of CO₂ in AHL appears to be almost entirely influenced by seasonality in advective flux (in our model due only to changes in DIC_OB – DIC_ocean, since T_f is held constant) and NEM (Fig. 5). Temperature contributes to a seasonal cycle, but the effect is masked by the overlapping and much larger advection signal, resulting in observed changes that bear little resemblance to the expected influence of temperature alone. From the perspective of organisms sensitive to Ω_Ar, pH, or pCO₂, it is clear from Fig. 6 that these parameters are primarily influenced by flushing and metabolic rate such that any perturbation to the advection (flushing of the lagoon) or respiration rate (organic and nutrient loads) may result in dramatic changes in these CO₂ system parameters. Thus, stakeholders (shellfish growers, for example) should be aware of changes that might impact flushing time, in particular.

Based on the estimated fluxes and their effects on DIC concentration (Fig. 5b), in a mass balance of AHL, the magnitude of advection will ultimately determine the estimate of NEM. Over the course of the 1-year time series, the integrated NEM was 82 ± 35 mmol C m⁻² day⁻¹ (positive = net heterotrophic). As a simple arithmetic check on this number, one can take a steady-state approximation where dDIC/dt = 0 and, ignoring F_gas, arrive at the simplified approximation:

$$NEM=\frac{\left({DIC}_{OB}-{DIC}_{ocean}\right)}{{T}_{f}}h$$

(24)

It is evident that, from the average lateral DIC anomaly between the ocean and the lagoon (71 mmol m⁻³), the average depth of 8 m, and the average flushing time of 7 days, the resulting average NEM (81 mmol C m⁻² day⁻¹) is essentially identical to that obtained by a more rigorous balancing of terms in Eq. 13. Nevertheless, it is instructive to use a temporally evolving mass balance rather than a steady-state simplification because it enables an examination of variability and driver decompositions as presented here. Furthermore, it is only possible to accurately assess the average lateral anomaly (Fig. 2b) after collecting a near-continuous annual time series, as sporadic point sampling throughout the year could lead to significant misrepresentation of a gradient that might be used in a steady-state approximation. Similar to DIC, including a temporally evolving flushing time in the biogeochemical model would certainly provide a more realistic time-dependent mass balance.

A NEM of 82 mmol C m⁻² day⁻¹ falls near the average of 55 mmol C m⁻² day⁻¹ (range of − 21 to 180) reported in the compilation of many US estuaries (42 NERRs sites) assembled by Caffrey (2004). While not a focus of this study, sources of organic material including the natural watershed along with a number of human-made installations (see “Study Site” above) likely contribute excess organic material fueling net heterotrophy in AHL. The name Agua Hedionda, or “foul-smelling water,” suggests either a high rate of respiration and/or a long flushing time. As we have seen, the heterotrophic rate in AHL is not abnormally high, nor is the flushing time unusually long. When given its name, AHL was only periodically connected to the ocean (e.g., during heavy runoff or seasonal high tides). The dredging and opening of AHL in 1954 radically altered the flushing (increasing it substantially) resulting in its anachronistic name. For comparison to other local data, in a recent study, under similar conditions in AHL but over sub-annual periods (n days = 52, 74 days), we estimated NEM in the range − 25 to 25 mmol C m⁻² day⁻¹ (Shipley et al. 2022). Interestingly, the Shipley et al. (2022) NEM estimates are based on the analysis of daily amplitudes from a different instrument using an approach that is independent of flushing time. Due to the differences in observational data and approach to NEM computation, Shipley et al. (2022) NEM values could be considered as mostly independent from the model approach used in this study. Of particular importance is the 1 month of corresponding data between the analysis of Shipley et al. (2022) and this work. During the month of Nov 2018, the in situ sensor used by Shipley et al. (2022) overlapped in time and was located within ~ 100 m of the instrument used in this study. For Nov 2018, the NEM derived from the daily amplitude of the in situ sensor was 20 mmol C m⁻² day⁻¹, while the LPF-based mass balance used in this study returned a value of 28 mmol C m⁻² day⁻¹. This level of agreement is quite encouraging, given the different approaches used between the two studies.

Since the outer basin is the only part of AHL connected directly to the ocean, and the biogeochemical model estimates advection based on measurements in the outer basin, the advective transport reported in Table 2 should reflect the export of the entire AHL. The most noteworthy comparison is the work of Paulsen et al. (2018) in the nearby San Dieguito Lagoon (SDL). The SDL provides an excellent comparison because it shares the neighboring watershed to AHL and, thus, a highly similar climate, ecosystem, and geomorphology. Using a classic approach of point sampling during selected periods of high/low tide, Paulsen et al. (2018) developed an estimate of the annual carbon export for the SDL. In their work, the combined inorganic and organic (DIC + TOC) export was estimated to be 10 to 25 × 10⁶ mol C year⁻¹. As pointed out by Paulsen et al. (2018), the inorganic to organic composition of carbon export is sensitive to runoff and frequency of storm events but, on average, total carbon export consisted of ~ 80% DIC, suggesting a DIC export from the SDL of ~ 8 to 20 × 10⁶ mol C year⁻¹ or roughly 2–3 times higher than that estimated for AHL in this study. It stands to reason that the much larger natural watershed area of the SDL compared to AHL (more than 10 ×) should lead to a greater export flux. On the other hand, the vast majority of river flow in all of San Diego County is captured by reclamation, leading to tidal lagoons with near-zero river input during much of the year and therefore making the natural watershed area less important than a more immediate or “effective” catchment area as defined above in the mixing model used here. In short, human-driven changes in land use may help to explain how a > 10 × difference in the natural watershed may translate into a 2–3 × difference in DIC export between the SDL and AHL. Also important is the frequency of storm events, which may lead to differences in the order of 2 × in export (Paulsen et al. 2018). Although storm events played a crucial role in our estimate of flushing time, the LPF data used in the biogeochemical model smoothed over these events. Incorporating the effect of storm events into a temporally resolved annual model remains a future goal and would be a logical next step beyond this work.

One of the most counterintuitive properties of coastal systems is the distinction between pCO₂ saturation and trophic status (net autotrophic vs. heterotrophic). Many coastal embayments, lagoons, and estuaries exist in a net heterotrophic state that sustains a pCO₂ higher than that of the atmosphere. However, since NEM (not pCO₂) defines trophic state, it is insufficient and incorrect to judge trophic status based on pCO₂ (or similarly O₂) alone. In AHL, pCO₂ is supersaturated with respect to the atmosphere, acting as a consistent source of CO₂ throughout the year 2018 (Figs. 2a and 5a). AHL appears to be net heterotrophic throughout the year, yet this (or any) system could transition to net autotrophy while continuing to act as a net source of CO₂ to the atmosphere, given a shift in controls such as the ratio of organic matter to nutrient inputs coupled with the ocean exchange rate.

Interesting in its own right is a comparison of the gas exchange flux between 2018 and 2020 (Fig. 7). As discussed by Shipley et al. (2022), a red tide in early 2020 resulted in an extreme event that drove AHL into a brief state of hypoxia. During a period of only a few weeks in May 2020, the outer basin reached peak heterotrophic rates of 140 mmol C m⁻² day⁻¹ and off-gassed the equivalent of > 50% of the full year of CO₂ flux observed in 2018. However, peak NEM rates observed in 2020 during the hypoxic event were no higher than the peak rates observed in 2018 (~ 140 mmol C m⁻² day⁻¹ during summer 2018, Fig. 5a). Based on nearby data at the Scripps Pier (Clements et al. 2020), the widespread nature of the 2020 red tide likely led to a high DIC in the adjacent ocean, which would result in diminished advective export, driving pCO₂ very high in the lagoon without the requirement for a massive jump in NEM. In summary, we hope that the simple comparison in Fig. 7 provides some insight into the potential for shallow coastal environments to exhibit abrupt changes influenced by the adjacent coastal ocean. This observation, in addition to the importance of storm events mentioned above, highlights the need for improved spatiotemporal observations in estuaries to better assess their contributions to the global carbon budget.

Fig. 7

Accumulated CO₂ air-sea flux during the 365-day period of 2017/2018 and a 4-month period of 2020. The rapid decrease in accumulated flux in 2020 was in response to a coastal red tide event that was also observed inside several estuaries (Shipley et al. 2022)

Conclusions

Based on a mass balance model of the outer basin of AHL, we found the system to be net heterotrophic, a source of inorganic carbon to the adjacent ocean, and a source of CO₂ to the atmosphere on an annual basis during 2018.

In this work, we have demonstrated the temporal variability of carbonate species and attempted to quantify the annual inorganic carbon budget of AHL. In this first-of-its-kind effort to quantify carbon fluxes in AHL, we identified the importance of collecting a temporally resolved dataset along with various challenges associated with interpreting data collected in a highly dynamic environment. While the average results of our simple biogeochemical model are in line with expectations based on similar lagoon systems, the rather large estimates of uncertainty associated with the reported annual budget reflect room for improvement in observing strategy.

The model used in this study is part mechanistic and part inverse. Physical processes (mixing, gas flux) have been parameterized mechanistically while the biological component (NEM) is determined as the balance of these physical processes along with observations. A full mechanistic model of the system could provide additional insight. For example, in forecasting AHL conditions under a new scenario such as enhanced circulation associated with the desalination plant and dredging, NEM might be assumed to follow the average pattern determined in this study, combined with a new flushing time. This approach will only provide satisfactory estimates of lagoon conditions if metabolic rates truly remain unchanged by flow. In the case where flow rates affect metabolic rates (e.g., by modulating nutrient and organic matter inputs), parameterization of NEM from variables such as nutrient and organic matter concentration, irradiance, and chlorophyll may prove useful. Similarly, a more complex treatment of physics through a hydrodynamic model coupled with a DIC mass balance could provide a more realistic description of the lagoon at the shorter time scales that were intentionally removed by the LPF in our analysis and allow for the incorporation of a possible seasonal variation in flushing times. Based on the demonstrated importance of advection, perhaps the most meaningful improvement to the model developed in this study would involve a better constraint on the temporal changes in DIC in the adjacent ocean, a more comprehensive set of total alkalinity data, and concurrent velocity observations. In addition, the important role of organic carbon should be addressed in the planning of any future study.