Metabolic rates in the Agua Hedionda Lagoon during the 2020 Southern California red tide event

Estuarine environments are uniquely diverse ecosystems located at the land-river-ocean interface. While these systems exhibit great biogeochemical heterogeneity (Bauer et al., 2013; Paulsen et al., 2018), coastal ocean systems are affected by many of the same sets of local drivers, including the adjacent ocean and surrounding landscape through runoff and river discharge (Paerl et al., 2006; Howarth et al., 2011; Windham-Myers et al., 2018). In comparison to offshore ecosystems, near-shore ecosystems are subjected more often to anthropogenic stressors, which makes understanding the metabolic state of these systems of great importance (Hewitt et al., 2005). Foremost, the balance between production and respiration is necessary for estimating net ecosystem metabolism (NEM), which defines whether a system is heterotrophic or autotrophic (Caffrey, 2003). Many estuarine systems experience both trophic states in a single day, some shift between states on a seasonal basis, and the vast majority appear to exhibit net heterotrophy, on average, presumably due to organic loading from land runoff (Smith and Hollibaugh, 1993; Hopkinson and Vallino, 1995; Gattuso et al., 1998; Caffrey, 2004). Dissolved oxygen (DO) is the most common measurement used to assess trophic conditions, as the value provides immediate information about the various thresholds (e.g., hypoxia, anoxia). DO is also the variable most widely used to establish NEM rates for coastal ocean systems (Breitburg et al., 2008). Both oxic status and NEM rates are fundamental to understanding ecosystem health and, specifically, deleterious effects of eutrophication and red tides (Conley et al., 2009).

Though not fully predictable, red tides can develop whenever conditions favor a strong phytoplankton bloom and tend to correspond to both natural processes such as upwelling and human sources such as nutrient and wastewater runoff from land (Horner et al., 1997). A particularly intense red tide event occurred from the end of March 2020 through June 2020 off the southern California coast. Based on satellite data, the bloom extended from Baja California to Los Angeles, with peak chlorophyll-a (Chl-a) concentrations (sampled manually at the Scripps Pier at a depth of approximately 5 m) reaching 867 mg m^–3 on April 27, 2020 (Kahru et al., 2021). This event was driven by the bioluminescent phytoplankton Lingulodinium polyedra, a predominantly non-toxic species in the Southern California region, but known to cause discoloration and post-bloom oxygen depletion (Horner et al., 1997; Armstrong and Kudela, 2006). During the peak of the bloom at the Scripps Pier, L. polyedra counts were 9,170 cells mL^–1 (Kahru et al., 2021). Non-toxic red tide events can be dangerous to aquatic species, often resulting in adverse impacts to ecosystem health by creating anoxic conditions that lead to fish kills (Anderson et al., 2012).

The primary objective of this study was to analyze biogeochemical observations from two different periods at the Agua Hedionda Lagoon (AHL), Carlsbad, California: a “normal” period in 2018 serving as the baseline condition against which we compare the extreme conditions observed during 2020, when a heavy storm preceded a red tide event. The datasets used in this work derive from a Southern California Coastal Ocean Observing System (SCCOOS) monitoring site that was originally installed for the purpose of providing near real time water quality information to an aquaculture facility. Due to a continuous presence in the lagoon, the sensor-based time series captured an event that otherwise would have been missed due to its transient nature and the considerable time required (typically months) to plan, mobilize, and deploy in situ assets. As an opportunistic study of an acidification/hypoxic event (at the height of a global pandemic), the data are naturally limited to assets in the water at the time of the bloom. Because ancillary measurements (e.g., concentrations of various sulfur and nitrogen species) that would normally be employed in a planned study of hypoxia/anoxia are not available, the work we present here is limited to an assessment of biogeochemical rates that can be derived from pH and oxygen data. To achieve this assessment, a weighted regression model (WRM) was used to filter out the tidal influence on DO and dissolved inorganic carbon (DIC) in order to obtain the biological signal, expressed here as the metabolic rates of production, respiration, and the net (production – respiration). Following this effort, we offer a speculative discussion of the chemical pathways that could be invoked to explain the observed dynamics in oxygen and pH; however, we must stress that the limited data prevent a definitive statement about the dominant hypoxic/anoxic reactions present during the study.

The AHL, located in Carlsbad, CA, is comprised of three interconnected basins, including an outer (26.7 × 10⁴ m²), middle (10.9 × 10⁴ m²), and inner basin (1.2 × 10⁶ m²; Figure 1; 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 hours at the creek (Jenkins and Wasyl, 2006). The original wetland was converted into the present lagoon structure in 1954 by the Encina Power Station and is maintained in its present form by semi-annual dredging. The AHL consists of >75% open water with the remainder being marsh and mud flats (Beller et al., 2014). Water depths at the lagoon range from very shallow (i.e., < 1 m) up to approximately 14 m in certain areas, with an average depth of 8 m (Elwany et al., 2005). The inner basin receives freshwater input from Agua Hedionda Creek during rain events which occur primarily in winter and spring. During the rest of the year, the creek runs dry, and the lagoon is purely tidal.

The AHL is highly utilized and a popular destination for the Carlsbad community and tourists and provides a thriving ecosystem for many diverse species of plants and animals. The locations of the various businesses within and surrounding the lagoon are also shown in Figure 1. The two primary industrial features include the Encina Power Station (fully decommissioned by January 2019) and the Carlsbad Desalination Plant, both of which rely on water intake from the outer basin for once-through cooling of the power plant and as desalination source water. During peak operation (prior to decommission) the power plant intake, located in the outer basin, could divert greater than 90% of the tidal prism over a tidal cycle; on average, around 50% of the tidal prism was diverted (Elwany et al., 2005). The desalination plant, which began operating in 2016, continues to divert water from the outer basin to generate freshwater at a rate of 2.3 × 10⁸ L d^–1, which accounts for around 20% of the tidal prism being diverted on a daily basis. Other features of the 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 that operates in the outer basin. Both the fish hatchery and aquafarm (which grows calcifying organisms sensitive to pH) rely on adequate flushing of the lagoon by the ambient ocean to maintain oxygen and calcium carbonate saturation levels above thresholds critical to growth. The primary tool used to accomplish the desired state of the ecosystem in AHL is dredging, which began in 1954 for the power station and transformed the lagoon from a backwater slough into an ocean-connected tidal lagoon.

Continuous biogeochemical observations were made from a mooring deployed in the outer basin of the AHL (Figure 2). The 2018 mooring was deployed for 69 days, from November 2, 2018, to January 10, 2019. In 2020, the mooring was deployed for 145 days, from January 22 to June 15. Sensors deployed on the mooring included a SeapHOx (Bresnahan et al., 2014) at a fixed height of 1 m above bottom to measure pH (Honeywell Durafet), DO (Aanderaa optode), temperature, salinity (Sea-Bird SBE-37), and pressure (Honeywell M5200 series). Each sensor was set to record measurements at 30-minute intervals. A surface buoy housed a cellular modem and controller connected directly to the SeapHOx, allowing hourly, real-time data access (Bresnahan et al., 2020). The data were quality-controlled and made publicly available online through SCCOOS and ERDDAP (Simons, 2020). The mooring was deployed near the CAF mussel and oyster cages (see Figure 1) for accessibility and relatively gentle flow (to minimize mooring drag due to tidal current).

Raw, 30-minute data stored onboard were re-processed and quality-controlled following recovery of the mooring. Raw oxygen data were salinity-corrected, following the manufacturer’s manual, using data from the SBE-37. Validation samples for pH were taken during the recovery of the mooring at the beginning of 2019 but were not possible in 2020 due to restrictions associated with the pandemic. For a bottle validation sample, a Niskin bottle was filled near the sensor and then subsampled into a 250-ml borosilicate bottle and poisoned with mercuric chloride (HgCl₂). The sample was later analyzed for pH using an Agilent 8453 spectrophotometer with temperature-controlled 10-cm cell following conventional best practices (Dickson et al., 2007). A pH correction was applied to the 2018 dataset following Bresnahan et al. (2014). Because the pandemic did not allow for validation samples during the 2020 deployment, the pH sensor offset determined in 2019 was carried forward for 2020 (the same Durafet pH sensor was used in the SeapHOx). Based on our previous experience, we estimate the accuracy of the time series to be 0.05 pH units, and in line with the manufacturer’s stated accuracy for oxygen (±2 mmol m^–3) and salinity (±0.003). As discussed previously (Takeshita et al., 2014; Kapsenberg et al., 2017; Shangguan et al., 2022), the ability of a Durafet pH sensor to accurately resolve short-term (e.g., hourly, diel) changes is greater than 10× the stated accuracy because, once deployed and conditioned, the sensors have been shown to remain highly stable and exhibit ideal Nernstian behavior. Based on the sensor resolution of approximately 0.0005 pH, a conditioned sensor can reliably observe relative changes in pH on the order of 0.001, making this sensor highly capable for capturing the amplitude of diel pH cycles, which ranged from 0.6 to 0.006 in this study.

CO₂ system calculations were performed using CO2SYS for MATLAB (van Heuven et al., 2011), using equilibrium constants as recommended by(Dickson and Millero 1987). Phosphate and silicate as inputs to CO2SYS were set to zero because, although they were not measured in this study, a brief campaign in 2016 observed concentrations in these nutrients at levels (maximum observed: phosphate = 0.35 µmol L^–1; silicate = 6.7 µmol L^–1) that are near-negligible in their effect on pH (approximately 0.001 pH unit). In our analysis, the partial pressure of carbon dioxide (pCO₂) and dissolved inorganic carbon (DIC) were calculated from SeapHOx temperature, salinity, and pH (total proton scale) combined with total alkalinity (A_T), estimated from salinity. A local A_T-S relationship was established from data recorded during a runoff event in 2018 where five paired pH and DIC measurements were used to establish the relationship: A_T (µmol kg^–1) = 29.2 × S + 1238 (RMSE = 6.3, R² = 0.98, n = 5; Shipley, 2022). Based on this relationship, an error of 10 µmol kg^–1 in A_T was used in the sensitivity analysis below. In converting pH to DIC, we have treated A_T as a property conservative with salinity because the canonical factor of ΔA_T:ΔDIC = 0.16 associated with Redfield stoichiometry is small relative to the pH signals observed and errors introduced by estimating A_T from S. To account for gas exchange, the hourly CO₂ flux (outgassing = negative, ingassing = positive) was calculated using Equation 1 (Wanninkhof, 2014):

where hourly wind speed (U) was available through the NOAA climate database, measured at the McClellan-Palomar airport in Carlsbad, CA, approximately 4 miles from the AHL (Vose et al., 2014). The ΔpCO₂ is the difference between the measured atmospheric pCO₂ (2018) and an average atmospheric pCO₂ of 411 µatm for 2020 and pCO₂ calculated from pH and A_T(S) in CO2SYS. The sensitivity to atmospheric pCO₂ is very small compared to the uncertainty due to pH. A 0.05 pH error translates to a 60 µatm change in pCO₂, making the error in the ΔpCO₂ term completely dominated by errors in the pH sensor. Any error within a few µatm in our knowledge of atmospheric pCO₂ is not significant for the output of the WRM (see below). For the model (next section), the cumulative sum of the flux (F, µmol kg^–1 hr^–1) is calculated over the time period of each dataset and then subtracted from the derived DIC in µmol kg^–1 to obtain a DIC value that is corrected for gas exchange. A comprehensive analysis of the errors in gas exchange estimates was not performed as it was determined that this term is a relatively small factor in the WRM estimate of rates.

The weighted regression model developed for riverine systems by Robert Hirsch (USGS) and modified for tidal systems by Marcus Beck (Tampa Bay Estuary Program) was used to filter out variability in DO and DIC due to tidal effects (Hirsch et al., 2010; Beck et al., 2015). The model code (written in R) was accessed on May 20, 2020, through the Github repository (https://github.com/fawda123/WtRegDO). Briefly, the WRM includes a series of functions that remove tidal influence (detide), estimate NEM, and evaluate metrics associated with performance criteria. The detiding function (Equation 2) is a multiple linear regression of observed properties (DO or DIC) in time (t) and tidal height (H).

See below for determination of β. The detided output is used to estimate daily NEM independently for DO (Caffrey, 2003) and DIC using Equations 3 and 4, respectively.

Total respiration (R_T) is calculated as the average nighttime rate of change for DO or DIC, multiplied by 24. Gross production (P_g) represents the observed change in DO or DIC during the day minus R_T. The F term represents air-sea gas exchange. O₂ gas exchange is handled within the WRM as originally written by Beck et al. (2015). The O₂ gas flux is computed from the difference between the DO at saturation (calculated as a function of temperature and salinity) and the observed value, multiplied by a volumetric reaeration coefficient, k_a (Beck et al., 2015).

Because the original version of the model was developed for DO only, several minor modifications were necessary to process the DIC time series. The primary change was to simply omit the gas exchange portion of the code during a DIC run and feed the model with the CO₂ gas exchange–corrected DIC, applied outside of the model (see above). Some additional unit changes were made to the model for the purpose of running DIC, including conversion of initial input of DIC to mmol m^–3 and several DO-specific unit conversions were omitted. To compare the two NEM values, NEM_DO was converted to carbon units using the canonical Redfield ratio of 106 C:138 O₂. Using this convention, P_g in C units is negative and a negative NEM represents net autotrophy; and R_t in C units is positive, corresponding to net heterotrophy.

The WRM applies adjustable half-window width settings (day, hour of day, tide height) to tune the model by adjusting each of the β variables as a function of time. The window widths are applied to time-series data to achieve gradual weighting from the center of the window over a specified range, where values that fall outside of this range are given a weight of zero. Beck et al. (2015) extensively tested the WRM by using it to extract NEM from a simulated control case with various random and systematic sources of error. For a thorough discussion of these results, including uncertainty analysis and comparison of NEM derived from detided versus unfiltered time series, please refer to Beck et al. (2015).

In tuning the WRM, we followed the approach prescribed by Beck et al. (2015), which used a span of five different widths for each of the three half-window width settings (day, hour, and tidal height) resulting in a total of 125 unique combinations of window width settings (we used the exact same 125 combinations outlined in Beck et al., 2015, p. 739). The tuned model is then considered to be the case (of the 125) which results in the optimal set of performance metrics provided in the summary statistics function of the model and outlined as follows. The first metric includes correlations between tidal changes and DO_obs, P_g, and R_t, which should decrease after filtering. Standard deviation for P_g and R_t is the second metric, which should also show a maximum decrease in the filtered output. The final metric is the percent of anomalous results, which is defined as the percentage of points in the time series with a negative P_g or positive R_t (for DO units). This percentage should be maximally reduced and ideally would be zero after detiding. The goal in tuning the model was to find the combination of weights where all three criteria were met. Finding the maximal reductions for each of the three metrics for any given combination of weight settings was not possible; therefore, a combination in which all three metrics were reduced in comparison to the unfiltered data was considered the optimal result for the 2018 and 2020 datasets (Table 1). There is one additional metric described in Beck (2015) using the annual mean for P_g and R_t, which was not applied to the AHL datasets due to the shorter time series used. A sensitivity test was performed to determine potential error in the WRM NEM output by introducing the associated sensor errors of 0.05 unit for pH and 2 mmol m^–3 for DO. This test resulted in no change for NEM_DO and a 0.2% error for NEM_DIC. An additional test for A_T was performed where an error of 10 µmol kg^–1 was introduced, which also yielded a 0.2% error in NEM_DIC. The insensitivity of the WRM to systematic errors (representing sensor accuracy error) is not surprising, as the model depends on short term changes in pH and DO, both of which are retained quite accurately even when pH and DO sensors may be slightly out of calibration.

In Figure 3 we compare baseline conditions of 2018 (panels A–E) to the anomalous red tide conditions observed during 2020 (panels F–J). The 2-month period from November–December 2018 was selected because it covers both dry and wet (rain) conditions. Temperature during both 2018 and 2020 ranged between 14°C and 20°C, although, due to the difference in time of year observed, the lagoon was cooling and heating during the 2018 and 2020 datasets, respectively. In the absence of rain, salinity remains very close to 33.8 during both years. Depending on duration and intensity, rain events can drive salinity down by several units, such as during a strong storm in 2020 when salinity briefly fell below 30. Following a rain event, salinity fully recovers in roughly one week, reflecting the average residence time of the entire AHL. During long periods without rain, the inner reaches of the lagoon trend toward slightly hypersaline relative to the surrounding ocean salinity of 33.6 (i.e., behaving as an inverse estuary). Average pressure (Figure 3C and H) reflects a difference of several meters between deployment locations in 2018 and 2020, with the latter deployment at a shallower depth of approximately 5 m. In 2020 the pressure sensor appeared to have become obstructed in mid-March around the time of a rain event but resumed regular operation by April 1. The rainfall in 2020 was more intense than that observed in 2018, with measured rainfall rates occasionally reaching close to 1 cm hr^–1 (Figure 3B and G).

Biogeochemical observations show nominal pH and DO in late 2018 (Figure 3D and E) and early 2020 (Figure 3I and J), where pH ranged between 7.9 and 8.1 and DO ranged from 212 to 262.5 mmol m^–3 with percent O₂ saturation ranging from 87 to 110%. These ranges typify the lagoon as slightly lower in both pH and DO than the adjacent surface ocean (Shipley, 2022; unpublished data), which are commonly observed features in a nearly balanced but, on average, net heterotrophic system (Cai, 2011). Following intense rain events in March and April 2020, the observed range of both pH and DO increased dramatically with pH reaching a maximum of 8.5 and a minimum of 6.9 in May. In connection with these extreme pH conditions, DO climbed as high as 517 mmol m^–3, which corresponds to about 225% saturation, and fell to zero (anoxic) on several days (Figure 3).

Figure 4A–D shows the 2018 and 2020 estimated rates for P_g, R_t, and NEM, including DO output for each period (Figure 4A and C) and the DIC results (Figure 4B and D). The average values for each rate are listed in Table 2 (with 2020 divided into Phases I–IV, as discussed below), where 2018 DO and DIC averages for P_g (–150 and –144) and R_t (156 and 158) are quite similar, although the 2020 plots (Figure 4C and D) show a much larger range due to the red tide event. Comparing Phase-I 2020 conditions to 2018, P_g and R_t are somewhat smaller in magnitude and NEM appears to be in a net autotrophic state.

The 2018 model output for DO and DIC is shown in Figure 5. Detiding removes both high and low spikes associated with tidal mixing, resulting in an overall reduction in diel amplitude, sometimes of greater than 50% (Figure 5A and B). The detiding step does not remove variability on timescales longer than 24 hours, however. Implications of this preservation of low-frequency variability are addressed in the Discussion. NEM for DO and DIC generally run at positive values with day-to-day changes (Figure 5C), which reflect conditions typical of a coastal lagoon where factors such as precipitation, sunlight, and spring and neap tide introduce variability. All daily NEM estimates fall within reasonable bounds, making the assignment of outliers resulting from model noise difficult (e.g., insufficient removal of tides). However, questioning the validity of spikes or abrupt shifts of NEM between adjacent days is warranted. Accordingly, to address noise in model NEM, a 6-day low-pass filter (LPF) was used to generate a smoothed time series (Figures 5C and 6C). The NEM LPF is also helpful in defining the four phases discussed below for the 2020 data. On average, NEM_DO in November–December 2018 was 6.75 ± 9 mmol C m^–2 d^–1 (n = 52) and NEM_DIC was 14 ± 11 mmol C m^–2 d^–1 (n = 52). Although the NEM_DIC average is twice that of NEM_DO, the 2-month mean (or integrated daily) NEM is not statistically different between NEM_DO and NEM_DIC, and the daily and LPF time series are in good agreement (Figure 5C).

Observed (black lines) and detiding results (blue lines) for (A) dissolved oxygen (DO) and (B) dissolved inorganic carbon (DIC), and estimates of net ecosystem metabolism (NEM) for (C) DO in carbon units (black lines) and DIC (blue lines) after applying a 6-day, lowpass Butterworth filter function (LPF, solid lines) and unfiltered NEM (dotted lines). Positive values indicate heterotrophic NEM; negative values indicate autotrophic NEM.The early portion of the 2020 time series (Phase I; Figure 6) is quite similar to the “baseline” conditions observed during 2018 with a net autotrophic Phase-I NEM_DO average of –6.6 ± 9.6 and an NEM_DIC average of –8 ± 16.1 mmol C m^–2 d^–1. Somewhat problematic is that variability in daily NEM is large relative to the mean, even under baseline conditions in both years, leading to mean NEM values that are statistically indistinguishable from zero. Whether late 2018 could be considered slightly net heterotrophic and Phase-I 2020 slightly net autotrophic is not clear, as stated above. Regardless of the trophic labeling of these time periods, the result of the model is that both periods appear to be operating in a near steady state of balanced NEM.

The massive red tide that developed during the spring of 2020 along hundreds of miles of coastline began influencing the AHL around mid-April. We defined the bloom onset (corresponding to the transition between Phases I and II) as the point where both NEM LPF lines inflected toward net autotrophy. Phase II culminates in a period of bloom termination and “hyperventilation” (Beck and Bruland, 2000), characterized by pronounced daily amplitude in both DO and DIC. Phase III, defined as the point when NEM LPF values from both DO and DIC exceed the largest rates observed under baseline conditions (approximately 50 mmol C m^–2 d^–1), represents bloom demise and an abnormally high respiration rate, ending in a weeklong period of sustained hypoxia (DO ranging between 0 and 100 and on average 49 mmol m^–3) and sporadic, limited anoxia. Phase IV, defined by the re-emergence of greater oxygen represents a weeks-long recovery and return to baseline conditions. Similar to 2018, detiding reduces daily extreme values, while retaining the lower frequency shift to a very low average DO (<100 mmol m^–3) and higher average DIC (>2300 mmol m^–3) during the peak conditions of Phase III. Within 1 month of anoxic conditions, DO and DIC had fully recovered to baseline conditions (due to the temperature change between Phases I and IV, DO recovery is slightly lower, but O₂% saturation (Figure 3J) recovers to 100%).

In general, detided results are devoid of the occasional spikes and steps observed during Phase I in both DO and DIC (Figure 6). As noted by Beck et al. (2015), suboptimal conditions for estimating NEM with the WRM include high correlation between sun angle and tide, and (surprisingly) decreased magnitude of the advection signal. In our data, the sun angle to tidal correlation was verified to be within the acceptable range for both time series, albeit closer to the acceptable limit within ±0.2 in 2018. There may also be true metabolic variability as a result of tides (Nidzieko et al., 2014), which the WRM might incidentally remove. Moreover, tidal stage is an imperfect proxy for advection. The effect of advection or, rather, gradients in general, is more difficult to assess. In their error analysis using simulated data, Beck et al. (2015) treat advection as a tidal term, linked to tide height, which is almost certainly the primary driver of signals associated with gradients. However, signals associated with phenomena such as diffusion, internal tides, and intake by the power plant in 2018 and the desalination plant in 2020 may escape the detiding process of the WRM. Regarding this effect, perhaps most noteworthy are two periods during 2020 (April 26–28 and May 22–25) when NEM_DIC significantly overshoots NEM_DO in the negative (net autotrophic) direction (Figure 6C). Assuming no massive shift in C: O stoichiometry during these two oxygenated periods, the most likely source of the discrepancy is a physical process that escapes detiding. Around both anomalies, the unfiltered DIC and DO (black line in Figure 6A and B) both exhibit pronounced, low amplitude (> daily) variability; however, the detided DIC (blue line in Figure 6B) retains a greater slope (particularly in late April). While this subtle retention of low amplitude structure may be the source of discrepancy between NEM_DO and NEM_DIC at times, we do not believe that it explains the difference observed during peak respiration of Phase III, when NEM_DIC overshoots NEM_DO in the positive (net heterotrophic) direction to extreme values.

Performance of the WRM under differing conditions is captured in 7-day snapshots of three different periods (Figures 7–9). Both 2018 (Figure 7) and 2020 Phase I (Figure 8) demonstrate the remarkable capability of the WRM to isolate the biological signal. In both cases, detided DO and DIC are in phase with the day/night cycle (shown in vertical bars) as expected. Figures 9 and Figure 4C and D show the one anomalous period (inverted P_g or R_t, as defined above) during the hypoxic event in mid-May 2020. During this period, the tidal signal is observed in the detided time series and, as discussed next, leads to questionable results.

One of the most intriguing and perplexing features of the 2020 time series is a marked difference between the NEM_DO and NEM_DIC during Phase III (Figure 6C), indicating a DIC respiration signal much stronger than that of DO. Peak NEM_DIC reaches 140 mmol C m^–2 d^–1 compared to a peak of around 50 mmol C m^–2 d^–1 for DO. A check, independent of daily amplitude, was performed to verify the approximate magnitude of the daily NEM rates generated by the model (Figure 10A and B). This approach involved applying a simple regression to arbitrarily selected multi-day periods of the detided DIC time series. The difference between the multi-day regression and the model output daily NEM therefore fundamentally lies in the difference between high frequency (daily period, peak-trough diel cycles) and low frequency (trend over multiple days) in the detided data. In theory, if detiding is perfectly efficient and NEM is a constant value, then a linear trend over multiple days will exactly match that obtained from the model’s peak-to-trough calculation of daily NEM. Of course, neither of these conditions holds true in practice, but as seen in Figure 10, the four trendlines corresponding to Phases I–IV are in reasonably good agreement with the daily NEM output of the model. The trendline check was only performed for NEM_DIC because the filtered (detided) DIC is already corrected for gas exchange (in the modified version we fed the model gas exchange-corrected DIC), whereas the original model, as designed for DO, performs a gas exchange correction internally and does not supply a gas exchange-corrected detided DO value in the output file. In summary, the simplest interpretation of Phase III data shown in Figure 6C is an evaluation of the average (LPF) output keeping in mind that daily NEM clearly contains additional noise due to the tidal signal making it through the filtering process during the week of extreme conditions. While detiding may not wholly remove 100% of tidal effect, it offers significant improvement to the more commonly used methods for assessing estuarine trophic status (Odum, 1956; Caffrey, 2003), which do not account for tidal influences at all. If the model noise over multiple days is random, then some of the underlying trend in the detided data will still influence the daily calculated P_g, R_t, and NEM. In this case, the LPF value still reflects the correct number as it averages over the 6 days.

If the differences in Phase III between LPF NEM_DIC and NEM_DO reflect a real signal, what process is responsible? The simple observation that average pH continues to decrease during the hypoxic portion of Phase III (Figure 3I and J) suggests that oxic respiration remained the dominant process. To visualize the multiple processes that may be present, we illustrate the effect of selected oxic and anoxic reactions in pH-DO space (Figure 11). The vector origins for each condition (oxic, anoxic) are somewhat arbitrary. For example, the symbol representing oxic conditions at pH = 8, DO = 200 along with the four vectors representing photosynthesis, respiration, calcification, and calcite dissolution could be located anywhere in the range of observed data where DO > 0. Similarly, the anoxic condition could be selected over a wide range of pH where DO = 0 (note that the anoxic vectors in Figure 11 are offset from zero DO for clarity). Important features of Figure 11 include: 1) a clockwise hysteretic cycle from Phases I through IV; 2) a primary signal due to P-R along with, perhaps, a secondary signal due to either mixing or CaCO₃ precipitation and dissolution seen most clearly in Phases I and II as a decrease in slope of the subset; 3) a wide range of pH under near-zero oxygen concentrations.

Regarding the final feature of Figure 11, if, for example, denitrification became significant, the A_T increase relative to DIC increase associated with the reaction

would be 0.96 rather than the Redfield value of 0.16 (Tromp et al., 1995). Thus, under denitrification in a closed system, ΔDIC:ΔTA is close to 1, resulting in a near-zero pH change, counter to what we observe in Figure 3I and J during Phase III. Under similar conditions, in the case of sulfate reduction, following the reaction (Tromp et al., 1995),

the pH shows an overall increase, which is the opposite of what we observe during Phase III where pH continues to decrease. Additionally, sulfide (pyrite) precipitation and oxidation reactions (Łukawska-Matuszewska and Graca, 2018)

show that, during these processes, A_T can decrease resulting in a drop in pH. In a closed system, sulfide formation followed by precipitation or oxidation will balance, resulting in no changes in A_T or pH. However, if decoupling (e.g., export) occurs when sulfide is formed versus consumed, a pH change will result. In Figure 11, the sulfide arrow is shown in the direction of decreasing pH, which could occur from a pulse of sulfide from the upper lagoon or an imbalance in sediment flux. While determining the exact processes responsible for decoupling during the hypoxic/anoxic period in 2020 is difficult, part of the continued decrease in pH may have been driven by imbalances in the sulfur cycle within the AHL.

We hypothesize that, although several brief (1–6 hour) anoxic periods occurred (Figure 9), oxygen remained in sufficient supply (in part through gas exchange) to fuel continued respiration of the terminating bloom. To represent this process, we include an additional vector in Figure 11 for cryptic oxygen consumption (Garcia-Robledo et al., 2017) in the direction of negative pH. This vector represents a system poised at near-zero oxygen where oxygen supplied (in this case through gas exchange rather than photosynthesis as proposed for anoxic marine zone systems by Garcia-Robledo et al., 2017) is immediately consumed with little or no accumulation, driving pH lower through normal oxic respiration channels.

Average rates for P_g, R_t, and NEM are reported in Table 2, along with ranges from a compilation of US estuaries (Caffrey, 2003; 2004) for comparison. Caffrey’s NEM compilation (table 2 in Caffrey, 2004) for US estuaries ranges from –21 to 180 mmol C m^–2 d^–1, with an average of 55 mmol C m^–2 d^–1 (g O₂ reported by Caffrey was converted to mmol C using –106 C:138 O₂ and 32 g mol^–1 O₂). Recent estimates of the regional averages for the US East Coast fall between 5.5 and 20 mmol C m^–2 d^–1 (Najjar et al., 2018). These estimates are in line with our findings for the AHL, where we observed a range of short-term NEM of –50 to 150 mmol C m^–2 d^–1 and baseline averages of –8 to 14 mmol C m^–2 d^–1. The P_g and R_t time series (Figure 4A–D) generally fall within the ranges for US estuaries. The 2018 P_g and R_t rates for DO at the AHL (Figure 4A and B) fall within the lower half of the range reported by(Caffrey 2004) with P_g ranging from –62 to –236 mmol C m^–2 d^–1 and R_t rates from about 60 to 243 mmol C m^–2 d^–1. We see the P_g and R_t rates reaching the extremes reported by(Caffrey 2004) during Phase II of the 2020 red tide. Perhaps the most noteworthy feature across Table 2 (and Figures 4 and 6) is that the most extreme P_g and R_t occur during Phase II while the most extreme NEM occurs during Phase III. This apparent decoupling between gross and net rates is not a widely reported observation, likely due to the requirement of continuous observations to derive the rate terms.

Applying an open source WRM to observations before, during, and after a red tide event has provided insight on the temporal progression of metabolic rates in a tidal estuary. While physical forcing through tidal advection is a major driver of variability in the lagoon, the WRM successfully isolated the local biological signal and quantified metabolic rates that drove the lagoon into a state of hypoxia. Baseline NEM established for 2018 and Phase I of 2020 are in line with estimates from other estuaries (Caffrey, 2003; 2004; Najjar et al., 2018), providing some confidence in the rates determined by the same model under extreme conditions in 2020. The wide range of pH observed at near-zero DO suggests that anoxic reactions related to sulfur may have occurred. The low pH signal expected for denitrification makes stating whether this process occurred based solely on pH and DO data impossible. However, we can state that processes other than denitrification were at work during the anoxic period as denitrification alone cannot explain the pH variability. In addition to the sulfur cycle decoupling, we hypothesize that there may have been a brief occurrence of a cryptic oxygen cycle like those discovered recently in anoxic marine zones (Garcia-Robledo et al., 2017). As an opportunistic and unplanned study of hypoxia, the suite of measurements employed was limited for studying anoxic reactions. As such we would recommend that any future study planned expressly for this purpose include, at a minimum, an autonomous nitrate sensor (e.g., Sea-Bird SUNA) alongside the pH and DO sensors.

The WRM was shown to be a useful tool that required minimal inputs to obtain NEM estimates for both baseline and extreme conditions in a coastal ocean setting. The production and respiration rates estimated in the model show the occurrence of anomalous periods, which are likely a result of the WRM not being able to filter out completely the tidal signals under extreme conditions (particularly for DO near anoxia), stratification, or significant calcification (in the case of DIC as a model input). Though beyond the scope of this work, a multi-box model of the AHL, including features of the lagoon such as sediments and aquaculture, would be a useful tool to further investigate the rise in DIC (decrease in pH) during the hypoxic period. Such a model could also demonstrate the influence of the desalination plant on residence time and differences in residence time between DO and DIC that might lead to the decoupling observed in Phase III.

Given some of the concerns mentioned above, an alternative approach to attempting to draw meaningful results from short peak-to-trough diel signals would be to estimate NEM from trendlines (or secants, finite differences, etc.) in the detided time series over multiple days to weeks of LPF detided data. Another approach to interpret the data more efficiently and possibly yield better NEM estimates would be to focus on specific time periods in the data that exhibit the lowest sun angle correlation or when the correlation is zero, as suggested in Beck (2015). In summary, we found the WRM and the novelty of the DIC addition to be a great resource and hope that our contributions outlined here may encourage others to explore similar uses of this tool.

The AHL observational data for 2018 and 2020 are available at the UC San Diego Library Digital Collections. https://doi.org/10.6075/J0MK6D3G.

We thank Thomas Grimm and the Carlsbad Aquafarm for providing small boat access to the Lagoon, Matt Steinke who helped with boat operations, Vicky Rowley from SCCOOS for assistance with the real-time data portal. We also thank Maria Herrmann and an anonymous reviewer for providing comments that greatly improved the quality of the manuscript.

Observations were supported by NOAA IOOS Award NA16NOS0120022 and NSF Award 1736905.

The authors have declared that no competing interests exist.

Contributed to conception and design: TM, KS, PB, TW.

Contributed to acquisition of data: TM, KS, PB, TW.

Contributed to analysis and interpretation of data: TM, KS, PB, TW.

Drafted and/or revised the article: TM, KS, PB, TW.

Approved the submitted version for publication: TM, KS, PB, TW.

How to cite this article: Shipley, K, Martz, T, Bresnahan, P, Wirth, T. 2022. Metabolic rates in the Agua Hedionda Lagoon during the 2020 Southern California red tide event. Elementa: Science of the Anthropocene 10(1). DOI: https://doi.org/10.1525/elementa.2022.00018

Domain Editor-in-Chief: Jody W. Deming, University of Washington, Seattle, WA, USA

Associate Editor: Kevin R. Arrigo, Department of Earth System Science, Stanford University, Stanford, CA, USA

Knowledge Domain: Ocean Science

Part of an Elementa Special Feature: Red Tide: Multidisciplinary Studies of an Exceptional Algal Bloom in Southern California