Main publications: Renault et al., 2021, Progress in Oceanography; Deutsch et al., 2021, Progress in Oceanography.

With funding from NOAA and the California Ocean Protection Council, we developed a state-of-the-art Earth systems model for the California Current System (CCS), incorporating both physical and biogeochemical components. The physical model consists of the Regional Ocean Modeling System (ROMS), with a series of advances to improve representation of small-scale coastal circulation, including high-resolution atmospheric forcing, wind drop-off and current feedbacks on surface stress. The biogeochemical component employs the Biogeochemical Elemental Cycling (BEC) model, comprising multiple phytoplankton functional groups, zooplankton, and the cycles of carbon, alkalinity, nutrients (N, P, Si, Fe), and oxygen. As part of prior work, we built a series of nested domains, from coastwide mesoscale-resolving (4-km horizontal resolution), to nearshore submesoscale-permitting (1-km), to very nearshore submesoscale resolving (300-m) configurations that allows inclusion of terrestrial inputs from rivers and POTW. Prior work with ROMS-BEC built on a series of nested domains, from coastwide mesoscale-resolving (4-km horizontal resolution), to nearshore submesoscale-permitting (1-km), to very nearshore submesoscale resolving (300-m) configurations. Analysis of two decades of simulations (1997-2017) demonstrated excellent consistency with observations from large to local scales for both the physical and biogeochemical states (Renault et al., 2020; Deutsch et al., 2021). Submesoscale-permitting simulations also showed dramatic intensification and variability on the shelf (Kessouri et al., 2020), and emergence of local hotspots of enhanced OA and oxygen loss (Damien et al. 2023, see also Section 3.1).

Main publications: Sutula et al., 2021, Marine Pollution Bulletin; Kessouri et al., 2021a, Journal of Advances in Modeling Earth Systems.

We developed a 300-m ROMS-BEC configuration for the Southern California Bight (SCB) to study coastal phytoplankton and biogeochemical cycles and investigate the effects of terrestrial inputs from rivers and POTW, and their consequences for HABs and other stressors (acidification and O2 loss).

Terrestrial inputs (Sutula et al., 2021). This 300-m SCB configuration was forced with a 1997-2016 daily time series of spatially explicit terrestrial nutrient inputs, including freshwater flow, N, P, silica, and organic C representing natural and anthropogenic sources. These data were compiled for: 1) WWTP ocean outfalls, 2) riverine discharges, 3) atmospheric dry and wet deposition and 4) atmospheric CO2 air-sea exchange. WWTP effluent data were compiled from permit monitoring databases. Riverine runoff from California’s ~400 coastal confluences are from model simulations (Sengupta et al. 2013) and monitoring data. The baseline covers 1971—2017 for large WWTP discharging >50 million gallons per day (MGD) and 1997–2017 for small WWTP and rivers. PS are the dominant nitrogen source, with contributions of 70% of the total annual freshwater discharge and 95% of nitrogen loads. WWT upgrades have reduced organic nitrogen loads by 73% since 1971. Inorganic nitrogen loads are generally held constant (35–40 Gg y−1) for the large WWT plants. This baseline represents a period prior to extensive wastewater and stormwater recycling that is increasing in the region (Figure-terrestrial input compilation). 

Figure – Terrestrial input compilation. Daily time series of 1997–2017 of all rivers summed for (a) riverine discharge and (b) total nitrogen load in the Southern California Bight. After Sutula et al., 2021.

SCB ROMS-BEC validation and evaluation (Kessouri et al., 2021). The 300-m SCB ROMS-BEC simulation was evaluated against a broad suite of observational data throughout the SCB, showing realistic depiction of the mean state and its variability with satellite and in situ measurements of state variables and biogeochemical rates (Figure-SCB model). Specifically, this “focused validation” was carried out to assess the accuracy of simulations in replicating anthropogenic gradients within the SCB. This validation utilized a dataset specifically tailored to match the spatial and temporal scales relevant to these gradients. The evaluation was conducted across five coastal sub-regions from Santa Barbara to San Diego associated with ocean outfalls. Various statistical tests and metrics were employed to quantify the disparities between model simulations and observational data. These metrics included measures of correlation, magnitude, variability, and persistent offset, and the results were categorized into performance levels such as excellent, very good, reasonable, or poor based on established criteria. 

The focused validation showed that the 300-m SCB ROMS-BEC simulation reproduces the main structure of the seasonal upwelling front, the mean current patterns, the dispersion of wastewater plumes, as well as their seasonal variability. Furthermore, it reproduces the mean distributions of key biogeochemical and ecosystem properties and their variability. Biogeochemical rates reproduced by the model, such as primary production and nitrification, are also consistent with measured rates. Overall, the model strikes a balance of capturing the forcing by U.S. Pacific Coast-wide phenomena, while representing the small-scale features that affect the transport of nutrients from natural and human sources. Moreover, it allows simulations at time scales that approach the interannual frequencies of natural ocean variability (e.g., ENSO). Therefore, we conclude that the 300-m SCB ROMS-BEC model is an appropriate tool to identify, investigate, and communicate uncertainty to stakeholders to support management decisions on local anthropogenic nutrient discharges to coastal zones, and their effects on the coastal biogeochemistry and ecosystem.

Figure – SCB model. (a) Regional Oceanic Modeling System-Biogeochemical Elemental Cycling (ROMS-BEC) model configurations. dx = 4 km is the black box, dx = 1 km is the blue box, and dx = 0.3 km is the red box. Background color shading show the topography from dx = 4 km. (b) Schematic of the Biogeochemical Elemental Cycling (BEC) model. The schematic shows state variables (boxes) and biogeochemical rates and feedback (arrows). The inner red box shows the Southern California Bight domain. After Kessouri et al., 2021.

Main publications: Sandoval Belmar et al., in preparation-a.

Physical setup and evaluation. A similar model development as for the SCB was conducted in the San Francisco Monterey Bay (SDMB) region. With funding from San Francisco Estuary Institute through to Chris Edwards of UCSC we constructed a 300 m ROMS-BEC nest and developed a targeted physical simulation of the water exchange between the San Francisco Bay and the Gulf of Farallones through the Golden Gate Strait. A description of such configuration is presented in Zhou et al., 2023. Simulations of the dispersion patterns of the San Francisco Bay (SFB) plume across the northern-central California continental shelf from 2011 to 2012 show that the model effectively replicated surface current dynamics and state variables observed in the region (Zhou et al., 2023). Results show that, after entering the Pacific Ocean through the Golden Gate, the SFB plume followed three primary routes: southward towards Monterey Bay, northward along the northern coast towards Point Arena, and an offshore path constrained within the shelf break. The northward-dispersed plume had an impact zone along the coast about 1.5 times longer than that of the southern branch. The plume’s dispersal exhibited significant temporal variability, attributed to wind and surface-current forcings. Results suggest that the SFB plume could spend up to 50 days away from the San Francisco Bay, before being flushed away from the Gulf of the Farallones.

Biogeochemical model setup. To incorporate exchange of biogeochemical tracers through the Golden Gate Strait, we combine the SFB east boundary forcing described in Zhou et al. (2023) with output from the CoSiNE model of Wang et al. (2020), based on the SCHISM framework. Firstly, data from the United States Geological Survey (USGS) Measurements of Water Quality in San Francisco Bay (SFB) were used to evaluate nutrient concentrations, chlorophyll levels, and temperature at monitoring stations 18 and 19 from 1969 to 2015. Station 19 is situated near the Golden Gate (GG) boundary, albeit with limited observations mainly between 1969 and 1979, while station 18 recorded data until 2015. The comparison of these stations revealed similar nutrient concentrations, although station 19 experienced saltier and warmer conditions due to its proximity to the GG. The differences in model outputs for both stations were found to be consistent across multiple years (2006, 2007, 2011, and 2012). Subsequently, adjustments were made to the CoSiNE model’s average vertical profile at the GG boundary, based on observations from station 18, focusing on depths with observations exceeding 10% of the total data. This correction was uniformly applied through the water column for various variables, such as temperature, salinity, chlorophyll, and nutrients. Additionally, the proportion of diatom chlorophyll to total chlorophyll in the region was estimated at approximately 93%, influencing the specification of boundary concentrations for diatoms and small phytoplankton. Plankton biomass data from the corrected CoSiNE model and converted to carbon units using Redfield ratios. Some tracers like iron (Fe), were not present in CoSiNE and were approximated differently. E.g., for Fe concentration, in-situ observations inside and outside the SF Bay were considered, and a climatological profile of Fe was computed. Boundary conditions for other elements tracked by BEC were generated using ratios calculated from the mean surface field on the 1-km Sothern California ROMS-BEC simulation (Kessouri et al., 2020).

Terrestrial inputs. Similar to the SCB 300 m configuration, we considered inputs from terrestrial sources, including 8 POTW outfalls and 28 rivers. Data on wastewater effluent from these sources encompass monthly time series for various parameters such as water flow, NH4, NO3, NO2, Dissolved Oxygen (DO), temperature, pH, Total Phosphorus (TP), PO4, Organic Phosphorus (OP), total Fe, SiO3, Organic Nitrogen (ON), Alk (alkalinity), salinity, and Total Organic Carbon (TOC) from 2010 to 2015. Data gaps were filled by using adjacent month values, linear interpolation, or correcting erroneous observations with climatological data. River discharges were determined through a combination of model simulations and monitoring observations. Water constituents were assumed to maintain constant values during the wet and dry seasons. For some observations, hourly output from the Soil and Water Assessment Tool (SWAT) model and monitoring stations were used. Salinity was calculated from conductivity and temperature, and missing DO observations were estimated. DOC was calculated using a typical Redfield ratio, and temperature was interpolated spatially to account for gradients. Finally, outfalls and rivers were implemented in ROMS-BEC as deep and surface point sources, respectively. The compilation is being synthesized for publication (Sutula et al. in prep.). 

Dataset for validation. We compiled environmental observations from several monitoring programs along the U.S. West Coast, including California Cooperative Oceanic Fisheries Investigations (CalCOFI), Monterey Bay Aquarium Research Institute (MBARI), West Coast Ocean Acidification Cruises (WCOA), California Current Ecosystem (CCE), among others. For comparison with the model, observations from 1980 to 2020 were organized into a monthly 15x15 km grid with a regular depth spacing (Figure-SFMB-Nutrient-compilation). We note poor data coverage on the continental shelf facing the SFB. The final analysis performed includes Pearson correlations, time series, and vertical profiles of each variable between model and observations in three main regions; San Francisco (SF), Monterey Bay (MB), and Monterey Bay offshore (MBo). The analysis was performed with both the data climatology (all available observations), and with observations only from years 2011-12. To keep it concise, we focused on important variables such as temperature, salinity, oxygen, chlorophyll, primary production, and nitrate concentration. These variables are crucial in understanding and predicting the dynamics of marine ecosystems, as their representation in the model can provide insights on gas exchange, mixing, circulation, biology, and status of marine ecosystems in the region, including HABs. We also included a series of satellite datasets for model validation. Sea Surface Temperature (SST) comes from the Multi-scale Ultra-high Resolution (MUR) SST Analyses, a global, gap-free, gridded, daily 1 km SST dataset created by merging multiple Level-2 satellite SST datasets. Data is available from 2002 to present. For satellite CHL, we used Kahru et al. [2014] merged product at daily 1 km resolution from 2002 to present, which we average monthly. We selected days when we have more than 10% of the observations, and model output was taken on those same days. Monthly Net Primary Production (NPP) comes from the standard product with the Vertically Generalized Production Model (VGPM) from the Ocean Productivity group at Oregon State University based on NASA MODIS data.

Figure-SFMB Nutrient-compilation: Example of in situ observations compiled for the San Francisco - Monterey Bay region, used for model validation. Left: Map of number of nitrate (NO3-) observations in the SFMB region, averaged over 10x10 km horizontal bins. Right: histogram of the number of observations by program in the SFMB region, as a function of year (top) and month (bottom). 

Simulations set up and attribution of anthropogenic effects. We designed three simulations to parse the role of nutrient inputs from terrestrial sources. All simulations are run from 2011 to 2012, with daily outputs. The first simulation (referred to as GGBC) only considered the natural oceanic cycles of nutrients plus the GG boundary condition. The second simulation (referred to as NO-BIO) used the same water flux, temperature, and salinity, but all biogeochemical tracers through the GG were set to zero, except for DIC, N2O, N2, O2, and ALK, which used the GG boundary conditions. The third simulation (referred to as ANTH) included the GGBC model plus inputs of nutrients, inorganic carbon, alkalinity, and dissolved organic matter (N, C, and P) from terrestrial anthropogenic sources. We plan to extend this simulation to the year 2014. The contrast between the ANTH simulation and the NO-BIO and GGBC models allows attribution of biogeochemical effects and eutrophication to different sources of terrestrial nutrients of anthropogenic origin (GG, point sources).

Main publications: Moreno et al., 2022, Harmful Algae.

Production and accumulation of DA, a secondary metabolite synthesized during periods of low primary metabolism, is triggered by environmental stressors such as nutrient limitation. To quantify and estimate the feedback between DA production and environmental conditions, we designed a simple mechanistic model of PN and DA dynamics and validated it against a series of batch and chemostat culture experiments from Kudela’s lab at UCSC. The model provides a framework for predicting PN DA production and accumulation, encompassing physical and biological modules that can be adapted to laboratory cultures or real environmental conditions. 

The biological module, referred to as BEC-HAB, simulates PN dynamics and DA production based on various environmental factors. BEC-HAB is built upon the BEC marine ecosystem model – a common lower trophic level model that is embedded in both global and regional modeling systems, including UCLA’s ROMS-BEC. BEC represents primary production and biogeochemical cycles in the marine ecosystem, accounting for nutrient uptake, photosynthesis, and organic matter cycling. In this specific implementation, the model focuses on dissolved inorganic nutrients (nitrate, ammonium, phosphate, silicate, iron), diatom biomass, chlorophyll, and particulate organic detritus, ignoring other living functional groups and the dynamics of the inorganic carbon system, oxygen, and nitrite. This simplified model is designed to replicate mono-specific PN lab cultures and their interactions with environmental factors. 

To simulate the cycle of DA, we implement a DA component for diatoms, here Pseudo-nitzschia (DiDA), and the two detrital pools (pDA and dDA for particulate and dissolved DA respectively). The equations and parameters for all components except DA are the same as BEC, as documented in previous publications (Moore et al., 2002, 2004). For DA, we represent the evolution of the three DA pools. DiDA production is directly dependent on the photosynthetic rate of PN, which is proportional to PN growth rate and biomass. We introduce a stoichiometric factor alpha (α) (in units of mol DA: mol diatom N) that encapsulates the relationship between production of DA and production of new biomass by diatoms. We assume that this term is a function of the nutrient limitation state of PN. Accordingly, a maximum ratio of DA to biomass production is modulated by a term that encapsulates the increase of DA production under stress by nutrients other than nitrogen (Figure – HAB model schematic).

Figure – HAB model schematic. Simplified DA formulation model schematic. The new DA formulation is part of a version of the BEC model adapted to represent Pseudo-nitzschia mono-specific laboratory cultures. BEC ac- counts for carbon (C), nitrogen (N), phosphorus (P), and iron (Fe). In diatoms (dark green circles), we account for all macronutrients in terms of N primarily, followed by a conversion to domoic acid (DA). Black arrows represent mortality and losses. Purple arrows represent remineralization steps. The green arrow represents primary production (PP) and uptake. The gray arrow represents sinking. After Moreno et al., 2022.

The physical module defines the experimental setup, including nutrient concentrations, water flow, temperature, and light. We configure the model to replicates batch culture experiments, optimize DA production in BEC-HAB against chemostat observations, and explore how environmental drivers affect DA production (Figure – HAB model setup). 

Figure – HAB model setup. Schematic of the 0-D model setup, illustrating different physical and biological components, and the ability to run configurations designed to represent batch and chemostat lab cultures, or an oceanic mixed layer. The model can be also used for optimization against observations. See description in Moreno et al., 2022.

The BEC-HAB model, coupled with a simple nutrient-dependent DA production formulation, effectively represents observed DA production in both batch and chemostat setups (Figure–chemostat results). The main model parameter that controls DA production is the maximum rate model identifies as a key parameter controlling cellular DA synthesis quota, and linearly affects total DA production. The half-saturation constants for nutrient limitation exhibit limited sensitivity, suggesting that once Si limitation is established, DA production quickly saturates to a maximum rate. Accordingly, the model reproduces the positive relationship between Si limitation and DA production and performs well in replicating the partitioning of DA between dissolved and particulate phases. It also suggests a high Si:N requirement in PN cells, potentially explaining the strong impact of Si limitation on DA production. Light and dilution rates also significantly affect DA concentrations. 

Figure – chemostat results. Top row: Model comparison against chemostat (n = 36) experimental observations. Model output against observed chlorophyll (A), macronutrients (B), and domoic acid (C) concentrations from chemostat experiments. Solid filled circles represent Si limited experiments, whereas open circles represent N limited experiments. Solid black lines show the one-to-one match. Bottom row: Model sensitivity to dilution rate and nutrient limitation. (D) Increased DA production with increased Si limitation. (E) Chlorophyll concentration as a function of dilution rate. (F) Seocific DA production rate as a function of Si limitation. Black colors represent chemostat experiment observations. Red colors model output. Filled circles show Si limited chemostats and open circles show N limited. Solid lines show typical model sensitivity. After Moreno et al., 2022.

Overall, the model represents a powerful framework for formulating and evaluating PN growth and DA production models, with the potential for future improvements. The model's simplicity allows it to be integrated into more complex ocean-ecosystem models for simulating DA-driven HABs, although regional adjustments may be needed for accuracy (Section 2.5 and 3.4). The model developed here represents a step toward realistic, regional models for DA and HAB prediction. We also identify additional factors such as trace metal (e.g., Fe) limitation and ocean acidification as future directions that should be investigated for a better understanding of PN HABs in a changing ocean.

Main publications: Sandoval Belmar et al., in preparation-b.

Our analysis of patterns and drivers of DA HABs in the Southern California Bight (Section 1) shows that this region is prone to HABs due to various factors, including natural upwelling and nutrient inputs from terrestrial sources. Furthermore, our modeling work with the 300-m SCB ROMS-BEC model (Sections 2.2 and 3.2) indicated that anthropogenic nutrient inputs increased phytoplankton biomass and eutrophication along the coastal band, and potentially, further offshore. Here, we expand that work to include a direct simulation of DA HAB dynamics by embedding into ROMS-BEC the mechanistic model by Moreno et al., 2022 (Section 2.4), and conducting realistic 3-D simulations in the SCB. This requires adapting the model by Moreno et al., 2022, to the more complex 3-D ROMS-BEC biogeochemical formulation and evaluating the model against the DA observations discussed in Section 1.

Implementation of the Model by Moreno et al., 2022 into ROMS-BEC. Several modifications were made to the DA formulation to enhance accuracy and account for processes not considered by Moreno et al., who focused on simulating chemostat and batch lab culture experiments. In BEC, we assumed that large phytoplankton, representing diatoms, always included potentially toxic PN species, capable of generating DA. This is a reasonable assumption for the US West Coast, where PN species are common. We converted biomass units from nitrogen to carbon, adjusting the stoichiometric factor 𝛽 for DA production based on in situ data (Section XX), reducing the allocation of cell breakdown to dissolved DA. We added a state variable to track DA accumulation in zooplankton (ZDA), introducing a parameter to control DA assimilation by zooplankton. Essentially, this parameter can be used to represent bioaccumulation and biomagnification in zooplankton. We also added a representation of the sinking flux of particulate DA (pDA) and its accumulation in the sediment (sediment DA). These modifications aimed to improve the representation of DA dynamics and were compared to observational data for validation.

Simulation setup for the SCB. We set up two domains for the experiments. The first domain, used for testing purposes, represents the Santa Barbara Channel (SBC), north of the SCB. The second domain is for the entire SCB, from Tijuana to Pismo Beach, covering the same region as Kessouri et al., 2021a,b. To reduce computational costs, both model configurations have a horizontal resolution of ∼1-km. However, experiments at 300-m resolution will be performed as part of future work. The lateral boundary conditions and initial condition for the SBC model are from the 1-km simulation of the Southern US West Coast (Kessouri et al. 2020) while for the Bight domain we use a coarser simulation (4km; Renault et al. 2021; Deutsch et al., 2021). Forcings are identical to Kessouri et al., 2021a,b (Sections 2.1-2.2). The SBC simulation is run for the years 2008 to 2012 and is used as a test for fine-tuning the DA formulation without freshwater inputs. The SCB configuration is run for the same period, with two experiments. The first configuration (CTRL) solely represents the natural oceanic cycles of nutrients, while the second (ANTH) supplements these cycles with nutrients, inorganic carbon, alkalinity, and dissolved organic matter inputs from terrestrial sources (Sutula et al., 2021), including anthropogenic inputs (Section 2.2).

Main publications: Kessouri et al., 2020, Global Biogeochemical Cycles; Kessouri et al., 2022, Journal of Geophysical Research; Damien et al., 2022, Global Biogeochemical Cycles; Damien et al., 2023, Geophysical Research Letters.

The 1-km ROMS-BEC simulations spanning the US West Coast, and the 300-m simulations for the Southern California Bight demonstrated an excellent ability to reproduce natural patterns of primary production nutrient, carbon, and oxygen cycles, and revealed the importance of fine-scale processes (such as sub-mesoscale circulation) in modulating phytoplankton distribution and primary production. Analysis of these simulations was essential to reveal the importance of capturing fine scales for properly simulating the patterns of phytoplankton biomass, productivity, and, consequently, HABs along the California coast. Key findings include:

• Modulation of phytoplankton production by submesoscale circulation (Kessouri et al., 2020, Global Biogeochemical Cycles). With this study, we quantified the impact of submesoscale circulation on nutrient and organic matter cycles in a realistic representation of the CCS, a typical EBUS. We found that: (1) in the coastal region, submesoscale eddies increase nutrient and organic matter subduction and quenching of productivity, further counteracting wind-driven upwelling; (2) in the offshore oligotrophic region, submesoscale eddies enhance the delivery of nutrients to the surface, fueling an increase in new production; (3) these submesoscale effects are modulated by the seasonal cycle, becoming more intense in the coastal band during the upwelling season, and in the offshore band during wintertime. The intensification of vertical transport by submesoscale eddies drives a readjustment of the planktonic ecosystem, with a reduction of phytoplankton and zooplankton biomass, productivity and size in the coastal region following upwelling, and an increase further offshore. In contrast, the export of organic matter from the surface layers, both as sinking and suspended particles and in the dissolved phase increases nearly everywhere. Increase in resolution to scales of 100s of meters or smaller would produce changes in the same direction and greater magnitude. For example, results with the 300-m configurations nearly double the range of vertical eddy fluxes of nutrients and generate much larger pockets of high phytoplankton biomass concentration and primary production. For this reason, we chose to conduct anthropogenic nutrient input assessments at 300-m, to adequately resolve the influence of submesoscale processes on nutrient and material transport, and phytoplankton accumulations, which are critical for HABs.
• Contrasting impacts of mesoscale and submesoscale eddies on primary production (Damien et al., 2023, Geophysical Research Letters). Prior studies, including our results in Kessouri et al., 2020, show that eddies play a crucial role in shaping ocean dynamics by affecting material transport and generating spatio-temporal heterogeneity. However, how eddies at different scales modulate biogeochemical transformation rates remains an open question. In this study, we applied a multi-scale decomposition to our coast-wide 1-km ROMS-BEC simulations (Section 2.1) and investigated the respective impact of mesoscale and submesoscale eddies on nutrient transport and biogeochemical cycling in the California Current System. Using a Reynolds decomposition of model tracers and biogeochemical rates, we showed that the non-linear nature of nutrient uptake by phytoplankton results in a 50% reduction in primary production in the presence of mesoscale and submesoscale eddies. Eddies shape the vertical transport of nutrients with a strong compensation between mesoscale and submesoscale. The eddy effect on primary production is controlled by the covariance of temperature, nutrient and phytoplankton fluctuations caused by eddies. These findings challenge the adequacy of non-eddy resolving global models to accurately represent phytoplankton primary production, and, consequently, HABs. Since the dynamics of pelagic ecosystems is govern by a variety of non-linear processes, from food-web interactions to respiration and microbial dynamics, eddy effects could greatly alter ecosystem dynamics and marine habitats, especially in environments naturally sensitive to multiple stressors such as ocean acidification, warming and deoxygenation and HABs.
• Fine-scale winds drive high productivity, low oxygen, and low Ph conditions in the Santa Barbara Channel (Kessouri et al., 2022, Journal of Geophysical Research). The Santa Barbara Channel is one of the most productive regions of the California Current System, and a HAB hotspot (Section 1). Yet, the physical processes that sustain this high productivity remain unclear. With this study, we used our 300-m ROMS-BEC SCB configuration to show that submesoscale eddies generated by islands are energized by orographic effects on the wind, with significant impacts on nutrient, carbon, and oxygen cycles. These eddies are modulated by two co-occurring air-sea-land interactions: transfer of wind energy to ocean currents that intensifies ocean eddies, and wind-current feedback that tends to dampen them. Our new analysis shows that the dampening is overwhelmed by fine scale wind patterns induced by the presence of surrounding capes and islands. The fine-scale winds cause an additional transfer of momentum from the atmosphere to the ocean that energizes submesoscale eddies. This drives upward doming of isopycnals in the center of the channel, allowing a more efficient injection of nutrients to the surface, and triggering intense phytoplankton blooms that nearly double productivity relative to the case without fine-scale winds. Since the frequency and severity of DA outbreaks are strongly correlated to chlorophyll (Section 1), these processes have clear consequences for HABs in the region. The intensification of the doming effect by the wind-curl and submesoscale eddies also pumps deep low oxygen, acidic waters to the center of the cyclonic eddies. These eddies are then transported away from the Channel into the California Current, where they impact a wider area along the central coast, with potential ecological consequences. Thus, this study highlights the important role of air-sea-land interactions in modulating coastal processes and suggests that submesoscale resolving models are required to correctly represent coastal processes and their ecological manifestations.

• Fine scale shelf dynamics enhances primary production on the US West Coast (Damien et al., 2022 Global Biogeochemical Cycles). The highest primary productivity and phytoplankton accumulations are observed along continental margins, which not only play a major role on carbon and nutrient cycles but are also hot-spots of HAB occurrence and impacts. However, assessments of biogeochemical cycles along continental margins remain uncertain. This uncertainty arises from the large variability over a broad range of temporal and spatial scales of the physical circulation and biogeochemical transformations that characterize these environments. Here, we use our 1-km resolution ROMS-BEC simulations for the entire US West Coast to examine the processes and balances for carbon, oxygen, and nitrogen cycles along the U.S. West Coast. We describe and quantify the biogeochemical cycles on the continental shelf, with a focus on primary production, and the connection of shelf cycles to the broader regional context encompassing the California Current System. On the shelf, coastal and wind stress curl upwelling drive a vigorous overturning circulation that supports primary production, biogeochemical rates, and fluxes that are approximately twice as large as offshore. Exchanges with the proximate sediments, submesoscale shelf currents, bottom boundary layer transport, and intensified cross-shelf export of shelf-produced materials modulate coastal and open-ocean primary production and other biogeochemical rates. The analysis approach laid out in this study is particularly powerful to identify the dominant physical-biogeochemical dynamics along continental shelf settings and asses their large-scale consequences.

Altogether, this body of work advances our understanding of the processes and dynamics controlling natural primary production and phytoplankton accumulations in coastal regions, with consequences for HABs. They also indicate that high-resolution ocean biogeochemical models that resolve the submesoscale are essential to study phytoplankton dynamics and HABs along continental shelves.

Main publications: Sutula et al., 2021 Marine Pollution Bulletin; Kessouri et al., 2021a, Journal of Advances in Models of Earth Systems; Kessouri et al., 2021b, Proceeding of the National Academy of Science, Kessouri et al., 2023 in review; Ho et al., 2023 in review; Frieder et al., 2023, in review; Hoel et al., in preparation, Kessouri et al. in preparation. 

The development and analysis of the 300-m SCB ROMS-BEC configuration with the inclusion of terrestrial nutrient inputs represented a major advance for the study of anthropogenic effects on the coastal ecosystem of the region, and revealed enhancement of anthropogenic stressors such as eutrophication, oxygen loss, acidification and HABs in the SCB. Key findings include:

• Anthropogenic nutrient inputs drive coastal eutrophication, acidification, oxygen loss, and ecosystem change in the SCB (Sutula et al., 2021 Marine Pollution Bulletin; Kessouri et al., 2021a, Journal of Advances in Models of Earth Systems; Kessouri et al., 2021b, Proceeding of the National Academy of Science). In this study, we use simulations with the 300-m resolution ROMS-BEC SCB configuration described in Section 2.2 to quantify the link between terrestrial and atmospheric nutrients, organic matter, and carbon inputs and biogeochemical change in the coastal waters of the Southern California Bight. The model is forced by large-scale climatic drivers and the reconstruction of local inputs via rivers, wastewater outfalls, and atmospheric deposition described in Section 2.2 (Sutula et al., 2021). The model captures the fine scales of ocean circulation along the shelf and it is validated against a large collection of physical and biogeochemical observations (Kessouri et al., 2021a). Based on comparison of simulations with and without terrestrial nutrient inputs, we show that local land-based and atmospheric inputs, enhanced by anthropogenic sources, drive a 79% increase in phytoplankton biomass, a 23% increase in primary production, and a nearly 44% increase in sub-surface respiration rates along a 15km coastal band during the peak upwelling period in summer, reshaping the biogeochemistry of the Southern California Bight. Seasonal reductions in subsurface oxygen, pH, and aragonite saturation state, by up to 50 mmol/m3, 0.09, and 0.47, respectively, are comparable to the global open-ocean oxygen loss and acidification since the preindustrial period. The results indicate consequences for zooplankton, local fisheries, water clarity, and submerged aquatic vegetation. These findings are consistent with the observational analysis discussed in Section 1 (Sandoval Belmar et al., 2023), which demonstrates a positive correlation in the SPC region between the frequency and severity of DA HABs and the volume flows from rivers and POTWs.
• Cross-shore physical transport promotes large-scale offshore response to urban eutrophication (Kessouri et al., 2023 in review). A key control on the magnitude of coastal eutrophication is the degree to which currents quickly transport nitrogen away from the coast to the open ocean before eutrophication symptoms develop. In Kessouri et al., 2021b, we showed that, in the SCB, anthropogenic nitrogen inputs increase algal biomass and productivity, and cause subsurface acidification and deoxygenation along the coast. However, the extent of anthropogenic influence beyond the coastal band, and the physical transport mechanisms responsible for these effects, were not previously documented. With this study, we extended the 300-m ROMS-BEC SCB simulations (Section 2.1) to the recent period (2013-2017) and investigated the transport of anthropogenic nitrogen and its effects on SCB offshore habitats. We found that anthropogenic nutrient inputs promote an increase in productivity and respiration offshore, with recurrent oxygen loss and pH decline in a region located 30 – 90 km from the mainland. Over 2013 to 2017, peak losses up to 14.2 mmol m−3 O2 persisted 4 to 6 months of the year over an area of 278,400 km2 (30% of SCB area). These recurrent features are associated with eddy cross-shore transport of nutrients and plankton biomass, and their accumulation and retention within persistent eddies offshore counteract the dilution and dispersion by mean currents that transport nitrogen, organic matter, and phytoplankton out of the SCB.
• Evaluation of the effect of ocean outfall discharge volume and dissolved inorganic nitrogen load on urban eutrophication outcomes in the Southern California Bight (Ho et al., 2023 in review). Because climate change is increasing drought severity worldwide, ocean discharges of municipal wastewater are a target for potable water recycling. However, potable water recycling would reduce wastewater volume, while maintaining nitrogen loading constant through discharge of reverse osmosis concentrate. Therefore, this practice has the potential to influence spatial patterns in coastal eutrophication and HABs. In this study, we applied the 300-m ROMS-BEC SCB model configuration (Section 2.1) to understand the influence of nitrogen management and potable wastewater recycling on net primary productivity (NPP), pH, and oxygen. We altered the wastewater discharges of 19 outfalls through combinations of dissolved inorganic nitrogen (DIN) reductions from 50 to 85% and recycling from 0 to 90% in the Southern California Bight. Under no recycling, NPP, acidification, and oxygen loss declined proportionally with DIN reductions. This simulated habitat volume expansion for pelagic calcifiers and aerobic taxa. Recycling scenarios under intermediate DIN reduction showed patchier areas of pH and oxygen loss with steeper vertical declines relative to a "no recycling" scenario. These extremes diminished under 85% DIN reduction across all recycling levels, suggesting nitrogen management lowers eutrophication risk even with concentrated discharges. These findings represent a valuable application of our ROMS-BEC model framework to investigate the regional effects of outfall management on eutrophication. Additional work is needed to investigate the non-linear interactions of nitrogen management with water recycling and to contextualize the benefit of these management actions, given accelerating HABs, acidification, and oxygen loss from human influences.
• Urban eutrophication affects pelagic habitat capacity in the Southern California Bight (Frieder et al., 2023, in review). Our results linked land-based nutrient inputs to increased coastal productivity, subsurface acidification and O2 loss in the SCB. However, whether eutrophication alters the capacity to support key taxa has yet to be evaluated for this region. In this study, we assessed the impact of land-based nutrient inputs on the availability of aerobic and calcifying habitat for key pelagic taxa using output from the 300-m ROMS-BEC SCB model. We found that acute, lethal conditions are not commonly induced in epipelagic surface waters, but that sublethal, ecologically relevant changes are pervasive. Land-based nutrient inputs reduce the potential aerobic and calcifier habitat during late summer, when viable habitat is at its seasonal minimum. A region of annually recurring habitat compression is found 30 – 90 km from the mainland, southeast of Santa Catalina Island. Here, both aerobic and calcifier habitat is vertically compressed by, on average, 25%, but can be as much as 60%. This effect can be traced to enhanced remineralization of organic matter that originates from the coast. These findings suggest that effects of land-based nutrients are not restricted to chemistry but extend to habitat capacity for multiple taxa of ecological and economic importance. Considerable uncertainty exists, however, in how this habitat compression translates to population-level effects.
• Natural oceanic and anthropogenic nutrient loading modulate the risk of toxic harmful algal blooms in the Southern California Bight (Kessouri et al. in preparation). In this study, we employed ROMS-BEC to investigate the relative role of natural oceanic versus land-based anthropogenic factors in influencing diatom production in the Southern California Bight and quantify how this modulates the “window of opportunity” for PN HABs to manifest. Specifically, we sought to: (1) describe the spatial and temporal patterns of productivity in diatoms versus other plankton assemblages in the SCB, (2) explore the relative contribution of natural, oceanic versus anthropogenic factors in driving this productivity, and (3) quantify how these natural versus anthropogenic factors alter the risk of PN HABs.

Main publications: Sandoval Belmar et al., in preparation-a.

We conducted an evaluation of the 300-m ROMS-BEC simulation for the SFMB (Section 2.3) including exchanges through the Golden Gate Strait and terrestrial inputs from rivers and POTW, focusing on the year 2011 (with planned extension to 2014). The model validation for 2011 demonstrates reasonable agreement with observations for temperature (R² ~0.9) and salinity, with a slight cold bias of up to 1°C at depth. The representation of nutrients such as nitrate as well as for oxygen is also consistent with observations (Figure SFMB ROMS-BEC validation). The model successfully captures the seasonality of chlorophyll, with moderate correlations with in situ data. Moreover, there is a weaker offshore gradient in chlorophyll, possibly linked to the representation of Fe release from the sediment in the model.

Figure – SFMB ROMS-BEC validation: Results from the 300m resolution ROMS-BEC configuration for the San Francisco-Monterey Bay region for year 2011. The model ANTH configuration includes exchanges from the Golden Gate Strait, rivers, and POTW). Top panels show a comparison of model and satellite chlorophyll. Bottom panels show a comparison of modeled and observed nitrate, as scatterplot, timeseries, and profile, for the Gulf of Farallones. The model captures the main patterns of nutrients and chlorophyll.

By comparing a run with all terrestrial inputs (ANTH) with a run with no terrestrial nutrient inputs (NOBIO) and a run with only inputs from the Golden Gate (GG) we can assess the effects of different terrestrial sources. The influence of inputs through the Golden Gate is substantial, potentially doubling surface nutrient and chlorophyll concentrations in regions adjacent to the Golden Gate plume and resulting in an increase of surface nitrate as large as ~15 mmol m⁻³ (about 30%). This enrichment extends offshore but declines rapidly after few tens of km. Near the estuary mouth, chlorophyll increases by up to ~10 mg m⁻³ (100%) when these inputs are included. 

In the Monterey Bay, terrestrial inputs from rivers and wastewater pipes also enhance the representation of chlorophyll, aligning it more closely to satellite data compared to a model with only freshwater inputs (NOBIO). Dissolved inorganic nitrogen increases by 80% when all terrestrial source are included, increasing coastal chlorophyll concentrations (Figure SFMB ROMS-BEC effects). Subsurface ammonia also shows improvement and appears more realistic. The overall increase in chlorophyll is around 7% in the Monterey Bay area and 30% outside the San Francisco Bay mouth, comparing ANTH and NOBIO. This difference is well seen in all the Greater Farallones National Marine Sanctuary around San Francisco Bay and the Monterey Bay National Marine Sanctuary, highlighting the potential ecological significance of the addition of nutrients from rivers and pipes and consequent increase in chlorophyll. 

Figure – SFMB ROMS-BEC anthropogenic effects: Effects of terrestrial nutrient inputs from the Golden Gate, rivers, and POTWs in the San Francisco-Monterey region. Maps and sections in the Gulf of Farallones (top row) and Monterey Bay (bottom row) show the difference between the model with all terrestrial inputs (ANTH) and a model where only freshwater terrestrial inputs are included (NOBIO). In each row, the first panel shows a map of the change in dissolved inorganic nitrogen (DIN), the second panel a map of the change in chlorophyll, and the third panel a vertical section of the change in O2 as a function of the distance from the coast.

A significant consequence of this increase in chlorophyll is an increase in dissolved oxygen of ~6 mmol m⁻³ from the surface to around 40 m depth due to boosted algal production. However, there is a rapid decrease below that depth by up to 5 mmol O2 m⁻³ all the way to the bottom because of enhanced respiration. Dissolved inorganic carbon mirrors the oxygen picture, showing an increase at depth near the coast about the same magnitude as oxygen. Finally, the comparison of the ANTH and NOBIO simulations reveals that the boost in chlorophyll is due to an increase in all functional groups, with diatoms increasing the probability of chlorophyll concentrations above 4 mg m⁻³, a threshold identified by Sandoval Belmar et al. as conducive to DA accumulations (Section XX). Given the known occurrence of Pseudo-nitzschia HABs in Monterey Bay, the increase in diatom biomass and chlorophyll suggests an increased occurrence of HABs in the region as driven by terrestrial inputs.

Main publications: Sandoval Belmar et al., in preparation-b. 

Initial model adjustments and validation. A series of initial simulations with the fully coupled ROMS-BEC-HAB model for the Santa Barbara Channel (SBC), without anthropogenic inputs, and Southern California Bight (SCB), with anthropogenic inputs, at 1km (Section 2.1-2.2) were conducted and analyzed against our biogeochemical and DA compilations (Section 1). Initial findings indicate that, to enhance the accuracy of DA production in the model, it was essential to improve representation of the Silica-to-Carbon ratio of diatoms. Extensive model runs covering 2008 to 2017 in the Southern California Bight revealed that the diatom component required an adjustment to better capture the increased uptake of silicon by PN during the upwelling season. This adjustment aimed to align the Si* (silica deficit over nitrate) values with observations. The model, without this improvement, tended to overestimate Si* due to higher nitrate (NO3) values compared to observations. Instead of altering the simulation of NO3, which was well-validated in previous studies, but not adjusted for river inputs, this approach altered Si(OH)4 values to improve Si* representation. The Si:C ratio for diatoms was increased from 0.137 to 0.274. Additionally, a systematic analysis identified a bias of 1 mmol/m3 in Si(OH4) throughout the coast-wide “parent” 4-km and 1-km simulations, which provides boundary conditions for these local simulations. To rectify this, the bias was removed at the boundary conditions during each time-step, resulting in improved Si* accuracy and a more precise seasonal cycle representation with reduced root mean square error. Overall, these initial simulations without anthropogenic inputs show a reasonable skill of the model to capture climatological DA observations in the SCB, although interannual variability remains more difficult to capture (Figure SCB ROMS-BEC-HAB results). 

Figure - SCB ROMS-BEC-HAB results: Results of ROMS-BEC-HAB simulation for the Southern California Bight. The model includes an explicit DA cycle (production from phytoplankton and accumulation in the food web and as particulate organic matter). Top panels: average surface DA for two model months (3/2008 and 4/2008) with observations from the same period shown as colored dots. Bottom panels: model-data comparison for the Santa Barbara Channel region. A model simulation (black line) is compared with climatological observations (red line) and monthly-varying observations (red line) averaged over the Santa Barbara Channel. Note that the model does not include anthropogenic nutrient inputs.

We applied the model improvements identified with the SBC configuration to a 1-km simulation for the entire SCB. A comparison of the results for the years 2008 and 2009 between the control simulation without anthropogenic inputs and the simulation with full terrestrial inputs reveals that the inclusion of terrestrial inputs from rivers and POTW outfalls, by significantly increasing the concentrations of NO3 and NH4 in the model, alleviates nitrogen limitation, allowing for greater Si(OH)4 drawdown and an increase in diatom biomass, which intensifies silica limitation. Consequently, pDA concentrations notably rise across the model domain, especially in the vicinity of POTW outfalls in the San Pedro Channel region. At the same time, improved spatial correlations with observational data are observed in both the Santa Barbara Channel and San Pedro Channel. This finding indicates that integrating river and POTW nutrient inputs is essential to capture DA dynamics. Despite these improvements, the model exhibits some level of DA accumulation in months where observational data do not indicate DA presence, e.g., December. Following episodic wintertime mixing, the model can produce high phytoplankton biomass, although in situ evidence for strong wintertime blooms remains limited. If these blooms lead to silica limitation, they can trigger DA production in the model. Assessing the realism of these mechanisms may need further research and observations and could lead to model improvements such as inclusion of temperature dependence for DA production.

Main publications: Smith et al., 2023, submitted.

As a complement to the synthesis and analysis of field DA observations discussed in Section 1 (Sandoval et al., 2023) and the anthropogenic pollution impact described above, we conducted statistical analyses to better link field observations of DA to ecosystem impacts, specifically targeting marine mammal strandings as key indicator of regional HAB impacts. In this analysis, we identify the spatial and temporal linkages between routine pier-based monitoring of DA to the stranding rates of adult and subadult California sea lions (Zalophus californianus). This was a leveraged effort with an additional $30K match of SCCWRP internal funding aimed at better identifying biologically relevant DA concentrations in the southern California region. To conduct this work, we paired 5 years of weekly pier monitoring data from four stations within the Southern California Bight (piers in Santa Barbara Channel, Santa Monica Bay, San Pedro Bay and San Diego) with the case records of animals rescued and treated by the Pacific Marine Mammal Center (PMMC) along the Orange County coastline. Basic patient characteristics including sex, seizure activity and other DA exposure symptomology was collated. From this, we conducted analyses to quantify the significance of observations from each of the monitoring sites across a region, determined if specific time-lags between DA observations and strandings occurred and identified which DA concentrations are associated with an elevated risk of stranding events (Figure – Stranding probability). We found that DA concentrations at Stearns Wharf, nearly 200 km from stranding locations, best predicted increased total strandings, and strandings of sea lions with domoic acid intoxication symptoms. Particulate DA concentrations greater than 0.05 μg/L at Stearns Wharf led to a detectable increase in stranding probability in Orange County, and concentrations over 0.25 μg/L resulted in a nearly 1.6-fold increase in stranding probabilities for a given week.

Figure – Stranding Probability: Observations and Poisson regression models of the relationship between domoic acid concentration and stranding scenarios, (A) total strandings and (B) sea lions exhibiting behavioral symptoms of domoic acid intoxication. Points indicate weekly observations. The x-axis corresponds to DA concentrations and the Y axis and point color reflect numbers of strandings. Points at or near the y-value of 0 correspond to weeks in which there were no strandings. Black points with y value at or near 1 correspond to weeks with one stranding. Blue points indicate weeks with two or more strandings. Black and blue bands indicate the probability, predicted by a Poisson regression model, of one or more (black) and to or more strandings (blue) at different domoic acid concentrations. Lines indicate the maximum likelihood probability and bands indicate two standard errors of that mean value

Main publications: Yamaguchi et al. 2022, Nature Climate Change; Moscoso et al., 2022, Ecosystem Modeling.

In parallel with the main activities of the grant, we have engaged in a series of collaborative efforts across a range of projects that have significantly enhanced our understanding of ocean ecosystem models and their sensitivity to climate change, both at regional and global scales. 

With the study by Yamaguchi et al., 2022, we focused on the global controls of phytoplankton blooms and their phenology under a changing climate. Using a 30-member Large Ensemble of climate change projections, we showed that earlier bloom initiation is expected by the end of the century in most ocean regions, with wide regional variations in bloom peak timing changes. Shifts in both initiation and peak timing are induced by decoupling between altered phytoplankton growth and zooplankton grazing, with increased zooplankton predation (a top-down control) playing an important role in altered bloom peak timing over much of the ocean. In the mid-to-high latitudes, phenological shifts will exceed background natural variability by the end of the twenty-first century, with implications for HABs and marine food webs.

With the study by Moscoso et al., 2022, we investigated how allometric (i.e., size-dependent) dependences in fundamental processes, such as growth and grazing rates, influence the size structure of planktonic communities. The work is based on an idealized size structured ecosystem model (SSEM) of the type used to represent the behavior of pelagic ecosystems. Specifically, we investigated the phenomenon of “quantization” in size structure, whereby localized biomass “peaks” emerge along the size spectrum (i.e., as dominant sizes in plankton communities). We found that that size selectivity of zooplankton predation controls the placement of these peaks, while variations in nutrient supply control the total magnitude of biomass. Leveraging these localized biomass peaks could lead to a more efficient discretization of continuous planktonic size spectra in models with a reduced number of size classes.

With the study by Moscoso et al., 2023, we investigate the role of different physical drivers in establishing observed cross-shore patterns of phytoplankton biomass and primary productivity in eastern boundary upwelling systems (EBUS). Previous research highlighted the significance of factors like winds, mesoscale eddies, and offshore nutrient distributions. In this study, we used a computationally cheap two-dimensional model to examines the impact of these factors on the California Current ecosystem. We find strengthening the wind stress maximum increases the average planktonic size in the coastal upwelling area, while the ecosystem is less responsive to wind stress curl changes. Deepening the nutricline shifts phytoplankton blooms toward the shore but also supports larger phytoplankton sizes. Eddy stirring of nutrients suppresses coastal primary productivity, while eddy restratification has minor effects. These results highlight the importance of wind stress maxima, isopycnal eddy diffusion, and nutricline depth in influencing coastal ecosystems in EBUSs, offering insights into differences among EBUSs and their responses to climate shifts.