# Biogeochemical variability in the California Current System

Curtis Deutsch†<fnref n="1" />, Hartmut Frenzel†, James C. McWilliams*, Lionel Renault*, Faycal Kessouri*#,
Evan Howard†, Jun-Hong Liang%, Daniele Bianchi*, and Simon Yang*

* Department of Atmospheric and Oceanic Sciences,
University of California, Los Angeles, CA 90095-1565

† University of Washington, School of Oceanography,
Box 357940 Seattle, WA 98195-7940

# Southern California Coastal Water Research Project,
3535 Harbor Blvd, Suite 110, Costa Mesa, CA 92626

% Dept. of Oceanography & Coastal Sciences, College of the
Coast & Environment, Louisiana State University, LA 70803

February 13, 2021

## Abstract

The biological productivity and diversity of the California Current System (CCS) is at the leading edge of major emerging climate trends, including hypoxia and acidification. We present results from a hindcast simulation (reanalysis) of an eddy-resolving oceanic physical-biogeochemical model of the CCS, to characterize its mean state and its patterns and drivers of variability in marine biogeochemical and ecosystem processes from 1995-2010. This is a companion paper to a physical analysis in Renault et al. (2021). The model reproduces long-term mean distributions of key ecosystem metrics, including surface nutrients and productivity and subsurface O₂ and carbonate undersaturation. The spatial patterns of Net Primary Productivity (NPP) are broadly consistent with measured and remotely sensed rates, and they reflect a predominant limitation by nitrogen, with seasonal and episodic limitation by *Fe* nearshore in the central CCS, and in the open ocean northern CCS. The vertical distribution of NPP is governed by the trade-off between nutrient and light limitation, a balance that reproduces and explains the observed spatial variations in the depth of the deep *Chl* maximum. The seasonal to interannual variability of biogeochemical properties and rates is also well captured by model simulations. Because of the prevailing nutrient limitation, fluctuations in the depth of the pycnocline and associated nutricline are the leading single factor explaining interannual variability in the interior biogeochemical state, and the relationships between density and biogeochemical rates and tracers are consistent between model and observations. The magnitude and relationship between density structure and biogeochemical processes is illustrated by the 1997-98 El Niño event, which faithfully reproduces the single largest deviation from the mean state in the simulated period. A slower decadal shoaling of the pycnocline also accounts for the concomitant trends in hypoxic and corrosive conditions on the shelf. The resulting variability is key to understanding the vulnerability of marine species to oceanic change, and to the detection of such changes, soon projected to exceed the range of conditions in the past century.

<footnote n="1">¹Corresponding author: cdeutsch@uw.edu</footnote>
38

# 1 Introduction

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

Coastal upwelling systems along the eastern boundaries of the subtropical oceanic basins are some
of the most climatically and biologically dynamic regions of the world's oceans (Carr and Kearns,
2003; Kudela et al., 2008). In the California Current System (CCS) of the Eastern North Pacific,
the seasonal cycle of alongshore winds and offshore surface currents yields a direct conduit for
nutrient-rich water to rise from the deep ocean to the sunlit surface (Chavez and Messie, 2009).
The upwelling of cold water near the coast is evidenced in satellite remote sensing of sea sur-
face temperature, expansive offshore marine stratocumulus clouds, and a nearshore ribbon of high
chlorophyll from phytoplankton. The resulting cascade of phytoplankton biomass up the food
chain supports high biodiversity and productive fisheries (FAO, 2009).

The physical state of the Northeastern Pacific varies on time scales of days to decades. This variability includes mesoscale eddies, seasonal cycle, interannual El Niño Southern Oscillation (ENSO) variability, and lower frequency climate fluctuations characterized by the Pacific Decadal Oscillation (PDO, Mantua et al. (1997)) and the North Pacific Gyre Oscillation (NPGO, Di Lorenzo et al. (2008)). Ecosystems throughout the Eastern Pacific respond strongly to physical forcing at each of these timescales, through the physical influence of winds, light, and heat, and their effects on the supply of nutrients to phytoplankton and oxygen to marine animals. Interannual climate cycles associated with ENSO in particular are major perturbations to these parameters and thus to plankton productivity (Chavez et al., 2002; Bograd and Lynn, 2001; Turi et al., 2018). Because thermocline waters entering the CCS are, like other subtropical Eastern Boundary Upwelling Systems, far from the sites of atmospheric ventilation in the west, the rising waters are “old” and bear the signature of decades of biogeochemical process that yield low O₂, high nutrients, and low pH. The upwelling in eastern boundary systems also generates energetic mesoscale eddy fields (Capet et al., 2008). These eddies can transport the extreme properties of CCS thermocline waters far offshore and down into the subtropical interior, connecting the biogeochemistry of the coastal CCS with the adjacent oceanic gyres (Nagai et al., 2015; Gruber et al., 2011; Renault et al., 2016a; Frenger et al., 2018).

Oceanic acidification and deoxygenation are also emerging trends in the CCS ecosystem (Gru-
ber et al., 2012; Chan et al., 2008). Anthropogenic CO2 has been detected in coastal subsurface
waters off Northern California (Feely et al., 2008). Decadal trends in oxygen have also been ob-
served in the California Current (McClatchie et al., 2010; Pierce et al., 2012; Dussin et al., 2019)
and have altered the proportions of biologically important nutrients (Deutsch et al., 2011). In the
Northern California Current, massive benthic die-offs have been attributed to episodes of extreme
hypoxia along the Oregon shelf (Chan et al., 2008). Fluctuating abundance of species in upper
trophic levels observed over decadal and longer time scales arise from climate variability, but the
specific mechanisms remain obscure (Chavez et al., 2003; Rykaczewski and Checkley, 2008). In
addition to climate forcing, the coastal ocean is subjected to anthropogenic pollution that could
locally exacerbate hypoxia and acidification (Doney et al., 2007), and may contribute to the in-
creasing frequency and toxicity of harmful algal blooms along California's coast (Andersson et al.,
2008).

The continental margin of western North America has a rich endowment of historical and on-
going observational programs aiding evaluation of climate-ecosystem interactions, including one
of the most extensive time series programs (the California Cooperative Ocean Fisheries Investi-
gations, CalCOFI) anywhere in the world's oceans (Ducklow et al., 2009). The age of shallow

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81
thermocline waters also make them prone to large amplitudes of low-frequency variability (Ito and Deutsch, 2010) that confound the detection of climate and other anthropogenic changes impinging through the oceanic surface. Thus, even in this relatively well-studied region of the ocean, it is difficult to distinguish long-term trends associated with anthropogenic climate change from the low-frequency variability that pervades oceanic properties in regions of close contact with old, deep waters. For example, it remains unclear whether the changes in hypoxia in the California Current are driven by internal climate variability, or are an early sign of long-term climate warming (McClatchie et al., 2010; Pierce et al., 2012; Long et al., 2016).

Models of this dynamic region often exhibit substantial biases in the mean state and unknown fidelity in representing historical variability and its causal mechanisms. The purpose of this paper is to report a system-wide validation of eddy-resolving, regional model fields and property relationships through comparison to a variety of hydrographic, experimental, and remote sensing observations. Here we focus on biogeochemical aspects of the model solution; the physical dynamics and its model validation are discussed in a companion paper (Renault et al., 2021). In this work, we demonstrate the fidelity of modeled spatial patterns and seasonal to interannual variability to observational datasets. We further probe the mechanisms underlying this variability–better understanding these drivers is critical to attributing and projecting the biogeochemical responses of the CCS to natural climate fluctuations and anthropogenic change. Section 2 provides a description of the model, its boundary conditions, the simulations performed, and the datasets used for model evaluation. These constitute the core results described in Section 3 and summarized in the final Section 4.

# 2 Methods

## 2.1 Model description

The ecosystem and biogeochemical cycles are simulated in the Regional Ocean Modeling System (ROMS, Shchepetkin and McWilliams (2005)). As in Renault et al. (2016b), the primary domain extends from 144.7°W to 112.5°W and from 22.7°N to 51.1°N. Its grid is 437 x 662 points with a horizontal resolution of $dx = 4$ km and 60 vertical levels. Initial and horizontal boundary data for temperature, salinity, surface elevation, and horizontal velocity are taken from the quarter-degree, daily-averaged Mercator Glorys2V3 data-assimilating ocean reanalysis (accessible via http://www.myocean.eu; described further at https://www.mercator-ocean.fr/en/science-publications/glorys/). In order to maintain a realistic representation of the variability in water mass properties and transport into the model domain over time, monthly anomalies from the Mercator data are added to the mean monthly climatology from the World Ocean Atlas (WOA, Locarnini et al. (2013); Zweng et al. (2013)) over the period 1995-2004. The freshwater, turbulent heat, and momentum fluxes are estimated using bulk formulae (Large, 2006) and the atmospheric fields derived from an uncoupled simulation with the Weather Research and Forecasting model (WRF). Heat and momentum fluxes are computed from bulk formulae, as detailed in Renault et al. (2021). The freshwater flux from river runoff is included as surface precipitation and is spread using a Gaussian distribution over the grid cells that fall within the range from the coast to 150 km offshore; this excludes a detailed representation of river plumes. The river-runoff forcing dataset we use is a monthly climatology from Dai et al. (2009). The river inputs are assumed to carry negli-
gible chemical concentrations, except the outflow from the Juan de Fuca Strait (see below). This statistically equilibrated solution, named USW4, is integrated over the period 1995-2010 after a spin up of 1 year (with initial conditions derived from the World Ocean Atlas).

Further details are described in the companion paper (Renault et al., 2021).

The coastal biogeochemical dynamics are simulated using an ecosystem model (the Biogeochemical Elemental Cycling (BEC) model, Fig. 1 and Appendix). This model includes both phytoplankton and zooplankton, and dissolved, suspended, and sinking particulates (Moore et al., 2004). The model includes four phytoplankton functional groups (picoplankton, diatoms, coccolithophores, and diazotrophs) characterized by distinct biogeochemical functions (nutrient recycling, silicification, calcification, and $N_2$ fixation, respectively). Four nutrient cycles (nitrogen, silicic acid, phosphate, and iron) are simulated and are coupled through a fixed phytoplankton stoichiometry, except for iron (Fe), which varies in proportion to the other nutrients (see Equations A100-109; Moore et al. (2002, 2004)). The ecosystem is linked to an oceanic biogeochemistry module that includes total inorganic carbon (DIC), alkalinity, iron, and dissolved $O_2$. Remineralization of sinking organic matter is parameterized according to the mineral ballast model of Armstrong et al. (2001). Gas exchange fluxes for $O_2$ and $CO_2$ are based on Wanninkhof (1992). The BEC equations are listed in the Appendix, and the model code, including parameter settings, are available through the GitHub repository (see the remark at the end of the paper).

<figure><img src="image_1.png" /><figcaption>Figure 1: Schematic structure and physical configuration of ROMS-BEC biogeochemical model. (a) The main ecosystem state variables and fluxes. (b) Geographic scale of simulation, and sources of surface forcing, open boundary condition data and representations of benthic nutrient fluxes.</figcaption></figure>

The iron Fe cycle includes dissolved iron, scavenged iron, and iron associated with organic matter pools and dust particles, but only dissolved iron and organically bound iron are explicitly modeled as state variables. For dissolved iron, four processes are considered: atmospheric deposition, biological uptake and remineralization, scavenging by sinking particles, and release by sediments. Atmospheric iron deposition is based on the dust climatology of Mahowald et al. (2006). We implemented a sedimentary iron source based on benthic flux chamber measurements in the California margin. An equation relating sediment Fe release as a function of bottom water $O_2$ ($\log_{10}[F_{fe}] = 2.5 - 0.0165 \cdot O_2$, where $O_2$ is in mmol m$^{-3}$ and the efflux units are $\mu$ mol m$^{-2}$
$d^{-1}$) is derived from data compiled by Severmann et al. (2010). The resulting rates of Fe supply from sediments (Fig. 2) exceed those from atmospheric dust deposition throughout the model domain.

We also added an anoxic nitrogen cycle, with losses to the sediments and water column. Bottom water nitrate is removed using a statistical description of sediment denitrification proposed by Middelburg et al. (1996), based on a vertically resolved diagenetic model that predicts the primary dependence of benthic denitrification to be on organic carbon sedimentation rate, with a secondary sensitivity to bottom water oxygen concentration. This statistical description of the complete diagenetic model reproduces basic controls on observed sediment fluxes, without the considerable computational cost of a sedimentary submodel. The predicted rates of $NO_3^-$ loss from this sedimentary sink (Fig. 2) amount to a small loss of $\approx 3 \times 10^{12}$ gN yr$^{-1}$. Denitrification in the water column is also modeled, but its integrated removal rate is an order of magnitude smaller than sedimentary losses, and has negligible impact on the results because $O_2$ in the model domain rarely falls below the threshold ($5$ mmol m$^{-3}$) assumed for this process. The higher $O_2$ thresholds associated with anaerobic particle micro-environments could increase the importance of anaerobic processes in the CCS, but they are not represented in this model (Bianchi et al., 2018). The removal of $NO_3^-$ by denitrification also acts as a sink for alkalinity.

<figure><img src="image_2.png" /><figcaption>Figure 2: Parameterized fluxes of iron and nitrate between the water column and sediments. (a) The Fe efflux (in mmol m$^{-2}$ yr$^{-1}$) to the water column from the sediments. (b) The $NO_3^-$ flux from the water column to the sediments due to net denitrification in sedimentary pore waters. Both fluxes are parameterized as a function of bottom water $O_2$, and denitrification is additionally parameterized as a function of organic matter flux to the seafloor. These maps are therefore part of the model solutions, and not prescribed forcings.</figcaption></figure>
166

## 2.2 Biogeochemical forcing and validation data

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

The biogeochemical model components and the physical (Renault et al., 2021) and biogeochemical forcings are schematically represented in Fig. 1. The biogeochemical boundary conditions for nutrients ($NO_3^-$, $PO_4^{3-}$, $Si(OH)_4$) and $O_2$ are taken from monthly climatological observations in the 2013 World Ocean Atlas (WOA) (Garcia-Reyes et al., 2014). Additionally, lateral $NO_3^-$ fluxes derived from prior model simulations (Davis et al., 2014) of nutrient exchange between the Strait of Juan de Fuca and the coastal ocean were imposed as landward boundary conditions in the Northern CCS. Non-nitrogenous nutrients were not available for inclusion in fluxes from the Strait of Juan de Fuca. However, we imposed an Fe concentration at this boundary by scaling it to nitrate. The scaling factor ($Fe=NO_3^-*3e-5$) was chosen to be similar to that of the surrounding coastal water. This ensures that Juan de Fuca nutrient inputs will not alter the locally limiting nutrient, absent data to support such an alteration. Boundary condition data for Fe is taken from global simulations with the Community Earth System Model (CESM) that used an earlier version of the same BEC ecosystem model. The $NH_4^+$ boundary concentrations, being small in nature, are set to zero, but adjust rapidly to the ecosystem processes in the interior of the domain. Time-dependent carbon cycle parameters, DIC and Alk, are taken from GLODAP (Key et al., 2004), with a reference year of 1995. An imposed trend at the boundary scales the anthropogenic component of DIC in proportion to the rise of atmospheric $CO_2$ since 1995. Time-dependent atmospheric $pCO_2$ is also used as a surface boundary condition for air-sea gas exchange. Aside from the boundary carbonate system parameters, the only non-stationary forcing of the model solution comes from the physical boundary conditions and surface forcing (Renault et al., 2021).

In order to ensure the integrity of tracer relationships along isopycnal surfaces, we map the biogeochemical boundary conditions from source data to the model grid using density rather than depth as the vertical coordinate, while retaining the mean-seasonal values of *T* and *S* (hence density) along the boundary as specified in the physical conditions (Renault et al., 2021). This prevents any errors in the depth of isopycnal surfaces inherited from the physical boundary data (Merca tor) from biasing the biogeochemical properties along that surface. Concurrently, this results in biogeochemical boundary conditions responding to interannual variability in isopycnal depth, despite being climatologically fixed along isopycnal surfaces at the boundary. Thus interannual bio- geochemical variability is propagated into the model domain because of the covariance between property isopleths and isopycnal displacements rather than changing biogeochemical water mass properties on density surfaces. However, biogeochemical variability along isopycnal surfaces in the interior domain can still arise from varying the proportions (mixture) of water masses enter- ing from different boundaries, or from time variable rates of biogeochemical transformation on isopycnals in the model interior.

The CCS is among the best-sampled regions of the world's oceans. Hydrographic sampling and biological rate measurements have been conducted repeatedly if not routinely along several sections off the West coast, most notably in the CalCOFI sampling area in the Southern Califor- nia Bight, off Monterey Bay, California, off Newport, Oregon, and Line P off Victoria, British Columbia at the northern edge of the 4 km model domain. Despite the abundant datasets from this region, data density is still sparse for much of the central California coast and for many biogeo- chemical properties of interest (e.g., *Fe*). The total number of profiles in the 2013 World Ocean Database (WOD; downloaded from https://www.nodc.noaa.gov/OC5/WOD13/ and including the hydrographic line data) are plotted for $NO_3^-$ and $O_2$ over the entire historical data period (1955-
2013; Fig. 3).

To facilitate comparison of model outputs to data, and particularly seasonal cycles, we defined 6 regions (Fig. 3), dividing the CCS by distance from the coast into nearshore (0-100 km) and offshore (100-500 km) regions, and by latitude into the Northern, Central and Southern CCS. These designations are somewhat arbitrary, but are based on a combination of topographic delinquations and to ensure adequate data coverage in each region. We focus our validation efforts on broad-scale measures that can be evaluated from climatological databases, namely WOD and its objectively-mapped climatological representation, the World Ocean Atlas (WOA; Garcia-Reyes et al. (2014)). We further evaluate the vertical and cross-shelf structure of biogeochemical variables at greater resolution at the hydrographic line locations. Higher frequency biogeochemical measurements from moorings are generally available only for more recent periods, and primarily from nearshore environments. Model comparison to mooring data is left for planned downscaling of these simulations better suited to examining high-frequency variability.

<figure><img src="image_3.png" /></figure>

Figure 3: Hydrographic data density used for model validation. Observations in the World Ocean Database are binned in a regular 1° latitude/longitude grid for each month over the entire historical data period (1955-2013). Total number of months with a profile are shown for $NO_3^-$ (left) and $O_2$ (right). Nominal station locations for major repeat hydrographic lines used for validation (see Fig. 16) are shown (red circles), along with the boundaries (black lines) used for regional time series comparisons in Figs. 6, 7, 14, and 21.
## 3 Results

We describe the spatial patterns and temporal variability of model biogeochemical solutions and their fidelity to observational datasets. Of the numerous properties and rates of the biogeochemical system that are predicted by the model, we focus on those that are most important to ecosystem primary productivity and the overall elemental cycling of carbon, nitrogen, and oxygen, and are best observed over scales captured by the model. This analysis begins with the photic zone, with particular emphasis on the factors driving NPP variability at multiple scales of space and time. We also evaluate the export of this productivity to depth. Second, we present results from the ocean thermocline, where the respiration of exported surface productivity contributes to hypoxic and corrosive conditions. Aspects of the model solution that are not presented include nutrients that are not limiting, rates that lack large-scale and climatological datasets, and variability that is poorly resolved by a $dx = 4$ km model (e.g., submesoscale and nearshore phenomena).

### 3.1 Photic Zone

#### 3.1.1 Chlorophyll and Net Primary Production

We begin with an evaluation of model distributions of Chlorophyll-a (hereafter *Chl*), and Net Primary Productivity (NPP), both of which can be estimated from remote sensing of ocean color. While NPP is of greater biological significance, its estimation is less direct than for *Chl*.

The model *Chl* concentrations are governed by the product of biomass and the *C : Chl* ratios. Biomass is subject to advection and to ecosystem transformations (see Appendix). The *C : Chl* ratio is determined by photoacclimation, or the amounts of light-harvesting pigments and photo-protective compounds produced by phytoplankton in response to their growth environment. This process is included in the ROMS-BEC representation of phytoplankton physiology following the model of Geider et al. (1998), which relates changes in chlorophyll synthesis and nutrient uptake in response to changing PAR. The dominant patterns of *Chl* are also found in biomass (see below), indicating that photoacclimation is not the leading factor in *Chl* variability. While biomass may be a more ecologically meaningful comparison, we validated model solutions using *Chl* because it is more directly estimated from ocean color sensing. The frequency distribution of *Chl* in both ROMS and SeaWIFS remote sensing data is approximately log-normal and is mapped after logarithmic transformation.

The annual mean concentrations of *Chl* vary most strongly in the cross-shore direction, with relatively weak alongshore gradients, a well-known pattern in observations (Banas and Hickey, 2008) that is well represented by the model (Fig. 4a,b). The offshore drop-off in *Chl* is somewhat weaker in model simulations, resulting in a wider band of high coastal *Chl*, a tendency that is not reflected in NPP (discussed immediately below). Kessouri et al. (2020) shows that there is some sensitivity in these distributions to model resolution, with higher resolution increasing the nearshore biomass and productivity. The leading pattern of variability in climatological *Chl* is characterized by a seasonal cycle that is also largely synchronous along the coast.
<figure><img src="image_4.png" /></figure>

Figure 4: Mean annual Chlorophyll (Chl) and its seasonal cycle. Mean annual concentrations (upper panels) are shown in mg m⁻³ averaged over the simulation period in both model output and SeaWiFS remote sensing level 3 product. The seasonal cycle is shown as the spatial loading patterns (middle panels) and time series (bottom panel) of the first EOF of climatological values of log 10[Chl]. In both the mean annual and seasonal variations of Chl, the dominant variations are cross-shore and at Point Conception, which separates the central CCS from the SCB. High coastal Chl extends further offshore in model solutions than in observations, a bias that is not found in productivity (Fig. 5).
The leading Empirical Orthogonal Function (EOF) (Fig. 4c,d,e), has *Chl* reaching peak values in late summer, and it accounts for the large majority of the climatological variance in both observations (EOF1=64% variance) and model solutions (EOF1 = 69% variance). The second EOF (not shown) also has a similar loading pattern, but with more meridional structure offshore and a minimum in spring–this accounts for only ≈13% and 19% of the variance in observations and the model, respectively. The seasonality of *Chl* reveals anti-phased cycles between nearshore high-*Chl* band and the lower *Chl* offshore. While near-shore *Chl* peaks in late summer, the offshore surface *Chl* has a minimum. This pattern is not found in the EOFs of depth-integrated NPP, indicating that the dipole structure of the *Chl* pattern results from a vertical redistribution of *Chl* to greater depths in offshore waters as they become more oligotrophic due to nutrient uptake during summer months. This interpretation is confirmed in the analysis of the vertical *Chl* maximum (Sec. 3.1.2).

In BEC, NPP depends on the sum of *j* model phytoplankton biomasses ($B_j$), their maximum growth rates ($\mu_{max}(T) = \mu_0 T^{1.06}$), and the limitation of those rates by light ($0 \le \gamma_j(I) \le 1$) and the minimum Michaelis-Menten function $\lambda_{i,j}$ ($0 \le \lambda_{i,j}(N_i) \le 1$) among the *i* nutrients with half-saturation $K_{i,j}$ ($0 \le \lambda_{i,j}(N_i) \le 1$), written as:

$$NPP = \sum_{i,j} \mu_{max,j}(T) \gamma_j(I) \min[\lambda_{i,j}(N_i)] B_j. \quad (1)$$

The spatial patterns of modeled NPP fall within the range of satellite-derived estimates (Fig. 5). The two commonly used satellite-based algorithms are the vertically generalized production model (VGPM) (Behrenfeld and Falkowski, 1997) and the carbon-based productivity model (CbPM) (Westberry et al., 2008). The VGPM estimates productivity on the basis of light and chlorophyll concentrations, calibrated to a predominantly coastal radiocarbon incubation dataset. The CbPM additionally incorporates phytoplankton backscattering and growth rate relationships in order to estimate productivity as a function of phytoplankton biomass, and it is calibrated to subtropical gyre radiocarbon incubations. The two algorithms exhibit a relatively wide range for the CCS region, reflecting the considerable uncertainty in “empirical” NPP estimates as well as differences in the measurements underlying each. The VGPM algorithm has a larger offshore gradient, with higher coastal values, and lower values in the open ocean, compared to the carbon-based CbPM. NPP rates from ROMS-BEC fall between the two remote sensing products, but are generally closer to the values of the VGPM algorithm, supporting higher near-shore rates, lower offshore rates, and increased seasonality relative to the CbPM. The VGPM has been explicitly calibrated against radiocarbon bottle incubations from the CalCOFI program, and it is therefore likely to be more accurate in this region (Kahru et al., 2009). Indeed, we find that ROMS-BEC rates and distributions of productivity are also consistent with direct estimates from ship-based data both from CalCOFI (Fig. 5; Munro et al. (2013)) and the broader subtropical Northeast Pacific (Palevsky et al., 2016).

Productivity in the northern CCS has been consistently biased in regional models (Banas and Hickey, 2008), including in our initial simulations. We conducted simulations with and without lateral nitrogen fluxes at the Strait of Juan de Fuca imposed as boundary condition from model simulations by Davis et al. (2014). Consistent with that study, without the nutrient inputs from the Salish Sea, NPP was biased low by > 50%. The inclusion of the effects of nitrogen inputs at the Strait of Juan de Fuca brought the model much closer to satellite-based empirical models. The inclusion of these inputs is consistent with the study by Davis et al. (2014), and it is used in all results reported here.
<figure><img src="image_5.png" /></figure>

Figure 5: Spatial distribution of annual Net Primary Productivity (mol C m⁻² yr⁻¹) integrated over the depth of the photic zone, from (top left) ROMS, (top right) CalCOFI, (bottom left) VGPM, and (bottom right) CbPM.

The seasonal cycle of NPP is also well captured by the model (Fig. 6). In all 6 regions of the CCS, the climatological NPP, integrated over the depth of the photic zone and averaged over the regional mask, exhibits an amplitude and phasing that is within the range of satellite-based empirical models. The most notable exception is the offshore Northern domain, where a spring bloom is predicted to be of stronger magnitude than estimated by either satellite product. This model result is consistent with measured geochemical tracers, which also indicate the spring bloom in the offshore Northeast Pacific is greater than estimated from the satellite algorithms (Palevsky et al., 2016). A smaller discrepancy occurs in the southern nearshore region, where ROMS-BEC generates greater summer production than either of the satellite algorithms. Model NPP in the oligotrophic part of the domain is lower than satellite estimates, but is more consistent with the most offshore values in the depth-integrated rates based on radiocarbon bottle incubations from CalCOFI. Overall, where and when ROMS-BEC and satellite algorithms for NPP disagree, ROMS-BEC output is generally closer to the available observational data.
<figure><img src="image_6.png" /></figure>

Figure 6: Seasonal cycle of annual Net Primary Productivity. The NPP rate (mol C m⁻² yr⁻¹) from ROMS-BEC (black), and two satellite algorithms (VGPM, blue; CbPM, red) are integrated over the depth of the photic zone, and averaged over 6 regions (see Fig. 3) from northern (top row), central (middle row), and southern (bottom row) of the CCS, and divided by distance from the coast into a nearshore (0-100 km; right column) and offshore (100-500 km; left column) region.

### 3.1.2 Seasonal limitation of productivity by light and nutrients

To evaluate the role of environmental factors shaping the seasonal cycle and regional differences in rates of productivity, we computed monthly mean limitation factors for each of the environmental variables that modulate the maximum growth rates, including macronutrients ($NO_3^-$, $PO_4^{3-}$, $Si(OH)_4$), Fe, and light; see (1). By construction, growth rates are limited by only one nutrient at a time (Liebig's Law of the Minimum), such that only the lowest value has an influence on rates. Light operates as a multiplicative factor on nutrient limitation, reducing growth relative to the light-saturated photosynthesis rate (Geider et al., 1998). Temperature influences maximum growth rate, but is not considered a limiting factor in the upper ocean, so is not analyzed here.

Over a climatological seasonal cycle, small plankton growth rates are almost always reduced by light more so than by nutrients, regardless of season or location (Fig. 7). The small plankton are assumed to have a lower half-saturation constant for nutrients, and the resulting higher affinity makes them less prone to nutrient limitation than large plankton are. Thus, in the inshore regions where nutrients are high, light is always limiting. Offshore, nutrients can limit small phytoplank-
ton to a similar degree as light in the summer. However, large phytoplankton make up ≈ 90% of modeled NPP on average, and large phytoplankton (“diatoms”) are primarily limited by nutrient availability. In the model, this results from a nutrient supply by combined wind-induced upwelling and eddy-induced subduction (Gruber et al., 2011; Nagai et al., 2015; Renault et al., 2016a) that is unable to saturate the potential net uptake of nitrogen by phytoplankton at prevailing light levels. Significant eddy fluxes also occur due to submesoscale eddies and fronts in the CCS (Kessouri et al., 2020), but are not resolved in this model. Large phytoplankton only experience light limitation in the northern CCS, where the seasonal cycle alternates between winter light limitation and nutrient limitation for the rest of the year.

<figure><img src="image_7.png" /><figcaption>Figure 7: Seasonal cycle of growth limitation factors for light (γ term in (1), blue) and for nutrients (λ term in (1): $NO_3^-$ (green) and $Fe$ (red) for diatoms (solid) and small phytoplankton (dashed)). Factors are NPP-weighted and averaged over the photic zone for each region shown in Fig. 3. Limitation factors close to 1 mean no limitation; values close to 0 mean complete limitation.</figcaption></figure>

While on the regional scales used for this analysis nitrogen limitation appears more stringent than $Fe$ limitation throughout the CCS, significant $Fe$ limitation occurs on smaller scales and shorter durations (see below). The seasonal amplitude of plankton growth rates is relatively small (≈ 10%), indicating that the amplitude of seasonal production (≈ 100%) is governed by seasonal controls on biomass rather than growth rates. The trade-off between light and nutrient limitation spatially and seasonally is a ubiquitous feature of phytoplankton distributions and phenology. Moving deeper in the water column, light becomes more limiting as photosynthetically active radiation
(PAR, light of 400 - 700 nm wavelengths) is absorbed and scattered, while nutrient concentrations are greatest below the surface mixed layer and rapidly decrease upward across the seasonally variable pycnocline; i.e., there is a well-defined nutricline.

The competing influences of nutrient and light limitation on depth of optimal plankton growth (i.e., highest growth rates) are reflected in the depth of the deep chlorophyll maximum (DCM). In the CCS, the depth of the DCM deepens from the coast to the open ocean, suggesting that the growth-maximizing combination of light and nutrients is found deeper offshore, consistent with a deepening nutricline that intensifies nutrient limitation at the surface and light limitation where nutrients are abundant, for both small and large phytoplankton. We use the observed pattern of the DCM depth as an indicator of whether the model achieves a realistic trade-off between these two countervailing growth condition gradients (Fig. 8).

However, for DCM depth to be a reliable indicator of phytoplankton growth rate trade-offs, we must rule out two alternative interpretations relating to isopycnal advection and photoacclimation. First, the offshore deepening of the DCM can have a physical origin. Because it closely follows the plunging of isopycnal surfaces offshore, the vertical peak in biomass and associated chlorophyll could be caused by eddy subduction carrying high surface chlorophyll away from the coast along deepening isopycnals (Gruber et al., 2011; Nagai et al., 2015; Renault et al., 2016a). To evaluate this possibility, we compared the depths of maximum chlorophyll concentration and phytoplankton biomass to the depth at which the product of light and nutrient limitation factors are maximized (Fig. 8). These maps are virtually indistinguishable, suggesting the DCM follows growth rates rather than advection. As a more stringent test, we performed a short (5-year) simulation in which surface PAR was reduced by 10%. The results revealed a significant shoaling of both the biomass and *Chl* peaks, but no detectable change in isopycnal depths, confirming that these depths do in fact reflect a nutrient-light trade-off rather than advection along density surfaces. Second, the peak depth of *Chl* may also be decoupled from that of biomass and growth rates due to photoacclimation, or a shift in the amounts of light-harvesting pigments and photoprotective compounds produced by phytoplankton in response to the light-environment. This process is included in the ROMS-BEC representation of phytoplankton physiology following the model of Geider et al. (1998), which relates changes in chlorophyll synthesis and nutrient uptake in response to changing PAR. Indeed, we find that the depth of the DCM is slightly deeper than that of the maximum plankton biomass. However, the offshore and latitude gradients of the depth of peak biomass and chlorophyll are very similar. DCM deepening offshore is consistent with optimized growth conditions in the model, and reproduces the pattern observed in the available CalCOFI data (Fig. 8).

In summary, the *Chl* maximum does not primarily reflect photoacclimation or isopycnal transport, although it is affected by those processes. Instead, it is found at approximately the same depth as that of peak biomass (panel c), which in turn is found at the depth that maximizes growth rate through the combined impacts of nutrients and light (panel d). This is the same depth at which NPP, the product of biomass and growth rate, is also maximum. Thus, these correspondences indicate that the model DCM is a reflection of the essential trade-off between light and nutrient limitation, and the fidelity to the observed DCM implies this trade-off is adequately captured.

Growth rates are modulated by a complex and evolving pattern of nutrient limitation by reactive nitrogen ($NO_3^- + NH_4^+$) and soluble dissolved *Fe*, with no appreciable limitation by $Si(OH)_4$ and $PO_4^{3-}$ in the CCS. The limitation factors are mapped as a biomass-weighted fraction of time that each of the nutrients is most limiting (Fig. 9). The spatial pattern among nutrients largely reflects
<figure><img src="image_8.png" /></figure>

Figure 8: Depth of the vertical maximum of Chlorophyll in CalCOFI (top left) and ROMS (top right), and of model diatom biomass (bottom left) and nutrient limitation factor (bottom right). The correspondence between model fields demonstrates that the model DCM tracks nutrient limitation, while the fidelity to observations implies the model captures a realistic trade-off between nutrient and light limitation.

the areas where $NO_3^-$ supply routinely exceeds maximum potential uptake seasonally. Thus, the waters entering the CCS from the subarctic High Nutrient – Low Chlorophyl (HNLC) region are most frequently $Fe$ limited. Along the coast, the seasonal upwelling of excess nitrogen leads to significant periods of $Fe$ limitation as well. In the coastal zone off Monterey Bay, $Fe$ limitation has been diagnosed via incubation experiments, in a band of water slightly offshore, with nitrogen limitation both in more shoreward and open coastal zones (Firme et al., 2003). This pattern is consistent with that predicted by the model (Fig. 9b inset). Most of the rest of the model domain is perpetually nitrogen limited.
<figure><img src="image_9.png" /></figure>

Figure 9: Frequency of limitation by $NO_3^-$ (left) and $Fe$ (right) for the model’s dominant primary producer, diatoms. The limitation factors are weighted by biomass (as in Fig. 7), using 5-day output, and plotted as the fraction of time that each nutrient limitation factor is the lowest among nutrients. The inset shows offshore band of relatively frequent $Fe$ limitation along the central CCS, similar to that observed by (Firme et al., 2003).

### 3.1.3 Nutrient concentrations

Surface nutrient concentrations provide an important measurable test of system behavior. For nutrients that limit phytoplankton growth, accurately simulating their distributions is a necessary condition for a mechanistic prediction of NPP. Moreover, they provide an integrative measure of net community production (equal to NPP less community heterotrophic respiration) and export of organic matter to the thermocline. We therefore compared model predicted distributions of the two primary limiting nutrients, $NO_3^-$ and $Fe$, to available observations.

Coastal measurements of dissolved $Fe$ reveal a spatially patchy distribution, reflecting its short residence time with respect to removal by plankton uptake and particle-active scavenging. Existing data are too sparse to yield a clear climatological pattern for model validation. However, the primary coastal region where the model predicts most frequent $Fe$ limitation, in the central CCS, has been relatively well sampled, including on two cruises off Monterey Bay that also tested for $Fe$ limitation (Firme et al. (2003); see Fig. 9, inset). Given the lack of a clear large-scale pattern of surface $Fe$ levels, we used a more statistically-based validation metric, focusing on the relative frequency of $Fe$ measurements versus concentration and distance from the shore (Fig. 10). On average, the data and model both show a decline in the mean and median $Fe$ levels with offshore distance. This reduction is driven largely by the decreasing frequency of high concentrations, while the most commonly observed $Fe$ levels remain consistent at $\approx 0.5 \cdot 10^{-6} \text{ mol m}^{-3}$ regardless of distance from shore. The thinning tail of high concentrations in $Fe$ distribution with cross-shelf distance occurs in both modeled and observed fields, but is more pronounced in the measurements.

In the CCS, the nutrient most often limiting NPP is reactive nitrogen, of which by far the largest and most commonly measured pool is $NO_3^-$.
<figure><img src="image_10.png" /></figure>

Figure 10: Histogram of surface *Fe* concentrations (10<sup>-6</sup> mol m<sup>-3</sup>) from observations (left) and ROMS (right). Because *Fe* concentrations are patchy in nature and sparsely observed, values are binned by distance from shore (0-20 km, upper row; 20-40 km middle row; 40-60 km lower row) to reveal a cross-shelf gradient, in the same latitude band off Monterey Bay for which observations were made on summer cruises (Firme et al., 2003).

of NO₃⁻ to climatological values from the World Ocean Database (Fig. 11). A depth of 50 m is chosen because it approximates the average depth of maximum biomass and NPP (see Fig. 8) and is generally near the base of the photic zone on the continental shelf. We included all historical measurements for this analysis because the data density in the model period (1994-2010) was much sparser, and no significant differences were found between the average NO₃⁻ in this period relative to 1955-2013.

ROMS-BEC captures regional patterns well for NO₃⁻ (Fig. 11). Annual mean concentrations of ≈ 15 mmol m⁻³ along the coast decline to values below the half-saturation level for model diatom growth (NO₃⁻ < 1 mmol m⁻³) within 500 km from shore. The offshore gradient is similar throughout most of the CCS, except in the Southern California Bight (SCB), where coastal surface values are much lower. Similar model fidelity was found for other macronutrients (Si(OH)₄, PO₄³⁻), but not shown because they do not reach limiting concentrations. The coastal zone exhibits strong variability in NO₃⁻ at 50 m, with standard deviations of 5-10 mmol m⁻³ throughout most
of the coastal zone, but with a slight northward increase in variance (Fig. 11c). The variability of $NO_3^-$ in the climatology (WOA) exhibits a similar spatial pattern, but with ≈ 50% of the magnitude. Thus, approximately half of the variation in surface nutrients most commonly limiting NPP is associated not with the seasonal cycle, but with interannual variability. The model also reproduces observed magnitudes and patterns of $NO_3^-$ variability (Fig. 11d). We use the interannual anomaly fields in the model and in observations to test the importance of nutrient supply as a mechanism driving changes in NPP over time.

<figure><img src="image_11.png" /></figure>

Figure 11: Long-term mean (upper row) and historical variability (bottom row) of nitrate (mmol m⁻³) near the base of the photic zone (50 m) in ROMS (left column), WOD (center column), and their correlations (right column). Variability is mapped as the standard deviation and consists of roughly equal contributions from seasonal and interannual variability (see text). The full period of the World Ocean Database (1955-2013) is used to yield the most robust estimate of variance. The correlations between ROMS and WOD are highly significant (p<0.01) for both mean and variability, with squared Pearson correlation coefficients (R²) of 0.84 and 0.80, respectively.
3.1.4 Interannual variability in surface nutrients and NPP

Factors that limit NPP during the mean seasonal cycle may drive interannual and longer term productivity changes. We examined correlations between NPP and both light and nutrient concentrations in model simulations and observations, where available. Interannual anomalies in $NO_3^-$ in ROMS-BEC are found to be well correlated ($R^2 \approx 0.5$) to the density of water at 100 m (Fig. 12). This reflects the role of pycnocline heave and vertical mixing in supplying macronutrients, and alleviating local nutrient limitation. Observations show a similar magnitude of correlation in the southern and central CCS, but a weaker correlation to the north. This may reflect the role of nutrient supply processes that are either missing or represented only climatologically in our model, and not connected to pycnocline heave. Because the weaker correlations are in the northern domain where nutrients can enter from subarctic surface waters, the climatological $NO_3^-$ used for the boundary conditions is a likely culprit. However, the presence of river $NO_3^-$ sources or a variable Juan de Fuca flux could also weaken the correlation in the observations relative to the model.

<figure><img src="image_12.png" /><figcaption>Figure 12: Correlation ($R^2$) between NPP and density, from ROMS-BEC (left) and CalCOFI measurements. The NPP rate is integrated over depth from $^{14}$C bottle incubation data, and density is interpolated to 100 m as an index of nutrient supply (see Fig. 13). In both variables, the mean seasonal cycle is removed leaving interannual variations. Relationships have similar strength in data and model, and indicate that ≈ 50% of interannual variability in NPP can be attributed to anomalous nutrient supply due to pycnocline heaving.</figcaption></figure>

Predicted correlations between NPP and density can be tested directly by combining bottle incubations and hydrographic observations in the southern CCS (Fig. 13). Relationships between nutrient and density anomalies (subtracting the mean seasonal cycle), are of similar strength, ac-
counting for ≈ 50% of the variability in both CalCOFI observations and model simulations. Inter-annual anomalies in NPP in ROMS-BEC are also significantly correlated with surface PAR (not shown) due to variable cloudiness, though it accounts for a smaller fraction of the variance (≈ 20%). The role of light is confined closer to the coastal upwelling where surface $NO_3^-$ is high, and light availability thus limits phytoplankton growth.

<figure><img src="image_13.png" /><figcaption>Figure 13: Correlation ($R^2$) between $NO_3^-$ and density at 100 m depth in ROMS-BEC (left) and WOD (right). For both variables, the mean seasonal cycle is removed leaving interannual variations. Interannual variations in subsurface (100 m) $NO_3^-$ highly correlated with density ($R^2 \ge 0.8$) in most of CCS. Correlation is weaker in northern CCS in data than in model.</figcaption></figure>

### 3.1.5 Carbon fluxes from the photic zone

Of the net production of organic matter by phytoplankton, a substantial fraction can be respired by zooplankton and higher trophic levels. The remainder is available to be transported away, by particle sinking and transport of dissolved organic matter, e.g., via eddy subduction (Sec. 3.1.2). The fraction of NPP that is regenerated within the surface ocean depends on food web processes, such as grazing rates. Although data is not available to evaluate large-scale patterns of grazing rates, an indirect comparison can be made through the export flux and the fraction of NPP that is exported in sinking particles rather than recycled (the so-called pe-ratio (Dunne et al., 2005; Murray et al., 1996)).

The fraction of NPP that is exported varies from 5-25%, consistent with the range of values inferred in field studies (Fig. 14). The model predicts highest pe-ratios in the coastal zone, where
productivity is high and sea surface temperature relatively low. These dependencies are also con-
sistent with those inferred from a global analysis of rate measurements for NPP and net community
production (NCP, assumed equal to total export of both particulate and dissolved organic carbon)
(Dunne et al., 2005). For a more quantitative comparison, we compared modeled pe-ratios to those
predicted by a statistical model fit to global observations by Dunne et al. (2005), over the sea-
sonal cycle in each of our 6 standard CCS regions (Fig. 14). In both the empirical model and in
ROMS-BEC, the mean value, phasing, and seasonal amplitude of changes in pe-ratio are similar.
Empirically based estimates generate a similar result that remains consistent with model simula-
tions, even when sea surface temperature is held constant. This suggests that temperature and its
impacts on the relative growth rates of phytoplankton and their grazers are not the essential cause
of variable pe-ratios.

<figure><img src="image_14.png" /><figcaption>Figure 14: Annual mean and seasonal cycle of fraction of NPP that is exported in sinking particles (i.e., pe-ratio; colored lines), and the export flux (mol C m⁻² yr⁻¹; black line). Values diagnosed from model simulations (red line) are compared to an empirically derived algorithm (Dunne et al., 2005) based on *Chl* and SST (blue line), and from the same algorithm applied with constant SST (green line). Observed annual net community production, which should approximate export on an annual basis are measured to 5-10 mol C m⁻² yr⁻¹ (Munro et al., 2013).</figcaption></figure>

Annual mean export flux represents the transfer of biogenic material from the surface to depth,
and thus the influence of CCS productivity on air-sea CO₂ flux and thermocline properties. Model simulated export production ranges from ≈ 1 – 10 mol C m⁻² y⁻¹ (Fig. 14b). Few measurements of export or net community production are available to evaluate the overall pattern of this flux. In the SCB region, measurements of O₂ have been used to estimate NCP rates of 3 - 17 (mean 6.5) mol C m⁻² y⁻¹ (Munro et al., 2013), in line with modeled rates and with similar spatial patterns of export (greatest along the Northern coast of the SCB). Munro et al. (2013) also combine information from ¹⁵N uptake experiments and nitrate based new production ratios (Dugdale et al., 1992; Eppley et al., 1992) with radiocarbon incubations to generate very similar estimates of the magnitude and spatial variability of production in the SCB. Particle-based estimates of export (neglecting the role of dissolved carbon) are lower than observed NCP as well as ROMS-BEC export, 2 - 4 mol C m⁻² y⁻¹ (Collins et al., 2011; Stukel et al., 2011), though enhanced particle export associated with mesoscale fronts (2 - 3 times greater rates over smaller spatial and temporal scales (Stukel et al., 2017)) highlights the potentially important role of subduction by eddies and fronts in explaining mismatches between observations. Importantly, ROMS-BEC generates such mesoscale features which contribute to model export estimates. While the model assumes particulate matter is redistributed vertically without being transported by the lateral circulation, the magnitudes and pattern of export are not substantially different from models that include explicit 3-dimensional particle transport (Frischknecht et al., 2018).

Model export production is also similar to regional nutrient budget analyses which suggest NCP averaging 7 - 9 mol C m⁻² y⁻¹ over the broader CalCOFI region (Roemmich and McCallister, 1989; Bograd and Lynn, 2001). Similar nutrient budget analysis indicates that annual NCP should be ≈ 17 mol C m⁻² y⁻¹ off Monterrey Bay, again in line with ROMS-BEC estimates of export for that region (Fig. 14b). As noted above, the fidelity of modeled nitrate distributions across the model domain provides a critical broad-scale measure of net biological drawdown, and thus of net community and export production.

Surface ocean CO₂ partial pressure and sea-to-air flux (Fig. 15) is reduced by net community production, increased by surface heat flux, reduced by freshwater fluxes, and modulated by upwelling and lateral circulation. It thus represents another important metric of overall system function. The role of the CCS and its sub-regions in the atmosphere-ocean balance of CO₂ has previously been investigated in several studies, both empirically (e.g., Hales et al. (2012)) and in models (e.g., Fiechter et al. (2014); Turi et al. (2014) and references therein). We evaluated the patterns of annual mean and summertime surface pCO₂ in the model hindcast simulation against observations in the SOCATv6 database from 1995-2010 (Bakker and coauthors, 2016), and the associated air-sea fluxes over the modeled coastal region. Similar to observations over the simulated period, strong CO₂ supersaturations are simulated in a narrow band of coastal water within about 100 km of the shore (Fig. 15). The values are highest along the central coast (35N-43N), lowest in the northern CCS, and intermediate in the southern domain. The highest values are associated with major topographic features, as previously noted by Fiechter et al. (2014).

The lack of apparent regional or seasonal bias in the surface pCO₂, together with the good model-data agreement in surface buoyancy fluxes (Renault et al., 2021) and biological rates in these simulations, suggests that the balance of processes regulating the model's surface CO₂ fluxes is reliably captured. Consistent with previous studies, the net integrated CO₂ flux to the atmosphere (an uptake of 1.4 Tg C yr⁻¹) is found to be a relatively small residual of larger compensating outgassing and ingassing fluxes. A detailed accounting of factors driving pCO₂ and air-sea flux variability has been described (Fiechter et al., 2014; Turi et al., 2014). Our model predicts a
somewhat larger net uptake than previous studies, likely because our model domain extends farther
529

into the northern CCS where the combination of heat loss and fresh water forcing suppresses
530

surface pCO2. An extension of similarly detailed analysis of air-sea flux variability to the northern
531

CCS domain is left to future work.
532
<figure><img src="image_15.png" /></figure>

Figure 15: Surface ocean CO₂ partial pressure (ppm) and air-sea flux (mol m⁻² yr⁻¹). (a) Annual surface pCO₂ from ROMS. (b) Surface pCO₂ from SOCATv6 gridded coastal dataset, averaged over all months from 1995-2010. (c) Annual air-sea CO₂ flux from ROMS. Negative values indicate outgassing to the atmosphere.
### 3.2 Thermocline

Here we test the model's representation of biogeochemical properties below the mixed layer and photic zone, i.e., in the thermocline. We focus on distributions and variability of $O_2$ and aragonite saturation state ($\Omega_A$) as these properties influence habitability for calcification and aerobic respiration by marine animals. The values of $\Omega_A$ are calculated from model dissolved inorganic carbon (DIC) and total alkalinity (Alk) using carbonate system equilibrium equations (CO2SYS) (van Heuven et al., 2011). The variability of $O_2$ has been analyzed in greater detail using previous simulations of this model with different atmospheric and physical boundary conditions (Durski et al., 2017).

To evaluate the vertical structure of biogeochemical tracers, we turn to repeat hydrographic sections (Fig. 16). Transects through three cross-sections spanning the Southern California Bight (CalCOFI line 80), the central California coast (MBARI line 67), and the central Oregon coast (Newport line) show the typical vertical and cross-shelf gradients of $NO_3^-$ and $O_2$. The downward enhancement of $NO_3^-$ and depletion of $O_2$ is a signal of the broad-scale vertical redistribution of these elements by the formation and degradation of organic matter within the CCS, as well as the gradients imported from the Pacific basin through the boundary conditions. The shoaling of the isopleths of both quantities follows that of the isopycnal surfaces by upwelling along the coast. The distributions of $NO_3^-$ and $O_2$ are generally well reproduced by ROMS-BEC. The model somewhat underestimates the slope of these isopleths very nearshore. This tendency is also reflected in, and likely derived from, the same underestimate in the zonal tilt of isopycnal surfaces (Renault et al., 2021).

Along isopycnal surfaces, $O_2$ generally increases with latitude and with distance from shore, reflecting the contrasting properties carried by the broad offshore California current from the $O_2$-rich subarctic, and the narrow near-shore California Undercurrent that transports low-$O_2$ waters of tropical origin northward along the slope. Both northern and southern end-member water types can be seen on the isopycnal surface 26.5 (Fig. 17), which also comprises the source of water upwelling onto the continental shelf along much of the US west coast (Pierce et al., 2012). On this and other density surfaces, the distribution of $O_2$ in ROMS-BEC is consistent with climatological observations, suggesting that the balance of distinct water masses and the respiratory modifications they experience in the interior of the domain are relatively well represented in the model. Thermocline nutrient distributions exhibit a similar model skill (not shown).

The $O_2$ in the thermocline of the CCS is highly variable, with standard deviations of $\approx 20$ mmol m$^{-3}$ that are on average 15-20% of the mean $O_2$ across the historical measurements (Fig. 17). The magnitudes and patterns of variance are well captured by the ROMS hindcast. In both the model hindcast and in observations, the standard deviation of monthly $O_2$ is $\approx 5$ times larger than that of the climatological seasonal cycle. Thus, a large majority of the $O_2$ variability is explained by non-seasonal time-scales, including large eddy-driven fluctuations (Frenger et al., 2018) and low-frequency climate variability (Buil and Lorenzo, 2017; McClatchie et al., 2010). Both model and observations indicate that $O_2$ variation peaks slightly offshore, in a pattern resembling that of eddy kinetic energy (Renault et al., 2021), and reflecting the role of eddies in transporting hypoxic waters offshore (Frenger et al., 2018). The high interannual to decadal $O_2$ variability observed in the central North Pacific, which reaches its maximum on the isopycnal 26.5, is not included in the climatological boundary conditions, and thus likely accounts for the model bias toward low $O_2$ variance in the north of the domain.
<figure><img src="image_16.png" /></figure>

Figure 16: Vertical sections of annual mean $NO_3^-$ (left) and $O_2$ (right) (mmol m$^{-3}$), from ROMS (left) and WOD (right) at the latitudes with regular observations by repeat hydrographic surveys. The lines span the northern CCS (upper row, ≈ 44.5°N, nearest the Newport OR), the central CCS (middle row, MBARI line 67), and the southern CCS (bottom row, CalCOFI line 80). Locations of observations are shown in Fig. 3.

Variations in $O_2$ on depth surfaces in the thermocline are also dominated by interannual anomalies. The fluctuating $O_2$ at depths of 100-200 m are highly correlated with density ($R^2 \ge 0.5$) throughout the CCS, reflecting the importance of isopycnal heaving of the background $O_2$ gradient (Ito et al., 2019). Similar correlations are observed in ROMS and the World Ocean Database (Fig. 18). In the model, the largest such anomaly is associated with the ENSO event in 1997-98, in which deepened isopycnals yield high $O_2$ conditions that last for ≈ 1 year. The signal is recorded in the central CCS as well, although the magnitude of the anomalies is reduced by ≈ 50% relative to the better-sampled event in the SCB.

The variability of $O_2$ within an isopycnal surface can be used to account for this portion of variance, leaving only lateral circulation and respiration. We find that along $\sigma_\theta = 26.5$ kg m$^{-3}$, a large fraction of $O_2$ variability is correlated with salinity (Fig. 19), commonly used as a proxy for tropical low-oxygen and high-salinity water transported poleward in the coastal undercurrent (Meinvielle and Johnson, 2013). Interannual variability in respiration rates on this surface also accounts for ≈ 20% of isopycnal $O_2$ variance, and is in turn correlated with the depth of the density layer (Deutsch et al., 2011). Historical observations show a declining strength of correlation between S and $O_2$ with latitude, suggesting that variability from dynamics other than the CUC become an increasingly important source of $O_2$ variability to the north. Indeed, the variability of $O_2$ in central mode water from the open North Pacific is most pronounced on this density surface (?Buil
<figure><img src="image_17.png" /><figcaption>Figure 17: Thermocline $O_2$ in model simulations and observations of the climatological mean (upper) and standard deviation (lower). All maps are interpolated to a potential density surface (26.5) surface chosen as the density class of waters upwelling onto the shelf in summer (a.k.a. “source water”). The mean $O_2$ maps are from summer months (JJA), however other seasons reveal similar model fidelity. Observed mean summer $O_2$ is from objectively mapped climatology (WOA). Variability is mapped as the standard deviation of monthly values from the World Ocean Database covering 1955-2013 (lower right), and it is predominantly due to interannual variability rather than the seasonal cycle (see text). The correlations between ROMS and WOD are highly significant ($p \ll 0.01$) for both mean $O_2$ and its variability, with squared Pearson correlation coefficients (R$^2$) of 0.95 for the mean state and 0.35 for the spatial pattern of temporal variability.</figcaption></figure>

and Lorenzo, 2017), and the North Pacific Current is thought to transport isopycnal $O_2$ anomalies into the northern CCS and contribute to the significant interdecadal fluctuations observed in the CCS [insert citation in comment here]. However, this longer term variance along isopycnals is not represented in the climatological biogeochemical boundary conditions in the model. Thus the relatively constant correlation between S and $O_2$ along isopycnal surfaces across latitude in the model may stem from the lack of $O_2$ variability on isopycnal surfaces as those waters enter the domain from the open North Pacific. A complete attribution of the observed magnitudes of $O_2$ variance
<figure><img src="image_18.png" /></figure>

Figure 18: Correlation between $O_2$ and density anomalies at 100 m, in ROMS-BEC (left) and WOD (right). Interannual variations in subsurface $O_2$ are highly correlated with density ($R^2 \ge 0.5$) throughout the CCS, and the strength of the relationship is similar in model and observations.

in source waters to the CCS will require inclusion of anomalies entering from the broader North Pacific (Deutsch et al., 2006; Kwon et al., 2016) and is left for future work.

In addition to low coastal $O_2$, the CCS is characterized by shallow depth horizons for carbonate saturation (Fig. 20). Below ≈ 100 m depth, carbonate concentrations are commonly undersaturated with respect to aragonite mineral formation, and thus inhibit shell formation by calcifying organisms. The aragonite saturation state, $\Omega_A$, is predicted to fall below saturation ($\Omega_A < 1$) along most of the coast in summer, consistent with observations in NOAA coastal surveys (e.g., Feely et al. (2008)). Coastal hydrographic surveys reveal a strong mesoscale patchiness to the carbonate saturation state, likely reflecting mesoscale eddies and submesoscale features. In the multi-annual mean distribution of $\Omega_A$, the model hindcast captures the scale and intensity of undersaturated conditions well.

Both low $O_2$ and low $\Omega_A$ have been implicated as primary factors mediating the influence of climate on organism fitness and species habitat in the CCS (Howard et al., 2020; Busch and McElhany, 2016). We compare decadal trends in both these properties from the hindcast simulations to the observed changes over time. For $\Omega_A$, the measurements are too sparse and the distribution too patchy to define a robust trend, even over the short model period. For each property, time series are shown for the regions with the most data coverage (Fig. 21): northern CCS for aragonite saturation, and southern CCS for $O_2$. In both cases, the trend in the data is within the uncertainty in the measurements.

As a metric of variability in these habitat constraints, we computed the volume of water subject
<figure><img src="image_19.png" /></figure>

Figure 19: Correlation coefficient between $O_2$ and salinity anomalies in (left) ROMS, and (right) WOD along the isopycnal surface $\sigma_{\theta} = 26.5 \text{ kg m}^{-3}$. Prevailing negative values indicate that high S occurs when $O_2$ is low. To maximize data availability, the observational correlation is based on monthly anomalies in WOD from the period 1955-2013. The correlations are of similar magnitude when confined to the simulated period, 1995-2010, but are only available in the Southern California Bight.

to hypoxic or corrosive conditions. We use a constant $O_2$ level of 100 mmol m$^{-3}$ as a simple indicator of hypoxic constraints, recognizing that this value varies among species, and depends on other factors, including temperature (Deutsch et al., 2015). Corrosive conditions are defined by simple thermodynamic undersaturation ($\Omega_A < 1$), though biological sensitivities may begin at higher thresholds. Water volumes are computed as the sum of grid cell volumes with $O_2 < 100$ mmol m$^{-3}$ or $\Omega_A < 1$ that are on the continental shelf ($-z < 200$ m).

The volume of hypoxic and corrosive water in the CCS varies strongly over latitude and time (Fig. 22). For both properties, restriction of putative habitat volume is stronger to the north of Pt. Conception, opposite the latitudinal gradient of $O_2$ and $\Omega_A$. The corrosive volumes are much larger than hypoxic volumes, exceeding 90% of water volume in northern latitudes during the summer upwelling season, consistent with NOAA survey data. On an annual basis, waters with a more stringent criterion for calcification ($\Omega_A < 2$) are about twice as voluminous still, primarily because the length of the season with low carbonate is broadened. Hypoxic conditions occupy a smaller fraction of shelf waters, but reach 30% of shelf water volume over a broad latitude range. The fractional coverage by hypoxia peaks around 45°N, on the Oregon coast, and it is quite small in the Southern California Bight, where $O_2$ declines most sharply below 200 m rather than on the shelf.

An analogous figure is shown in Renault et al. (2021) for the along-coast and temporal vari-
<figure><img src="image_20.png" /></figure>

Figure 20: Aragonite saturation state ($\Omega_A$) at 100 m from observations (left) and ROMS-BEC (right). Observations are from large-scale objectively analyzed fields (color field, GLODAP2 (Lau-vset and coauthors, 2016)), and from NOAA coastal surveys in the summer of 2007 (circles, Feely et al. (2008)). Model distribution is averaged over summer (JJA) from 2004-2010 climatology.

ability of the sea-surface height and depth of the $\sigma_\theta = 26.5 \text{ kg m}^{-3}$ depth anomalies. Relative to the quantities shown in Fig. 22, they exhibit more along-coast coherence and a more dominant seasonal cycle, with less evident interannual variability than shown here, apart from the 1997-98 ENSO event. This indicates somewhat smoother physical fields than biogeochemical ones, reflective of non-conservative biogeochemical processes acting on top of the broader patterns of physical circulation influence.

Variability in both habitat constraints is largely synchronous (Fig. 22c,d), reflecting the strong control on both $O_2$ and $\Omega_A$ by the effects of cumulative organic matter respiration. For both volumes, the fractional variation is similar and substantial, reaching ≈ 50% of the mean across much of the latitude range. Variability of corrosive volumes is greatly attenuated south of Pt. Conception. Simulated hypoxic volumes increased sharply during the early 2000s off Oregon, when major ecosystem die-offs were attributed to the onset of extreme hypoxia there (Chan et al., 2008).
<figure><img src="image_21.png" /><figcaption>Figure 21: Trends in thermocline $O_2$ and $\Omega_A$ over the simulated period. Maps of the linear trend are shown in upper panels. Time series (lower panels) are shown for the regions with the most data coverage for each tracer: northern CCS for carbonate, and southern CCS for $O_2$. For $\Omega_A$ all available profiles are shown. For $O_2$, the mean value and standard deviation are plotted as box-whisker for each month in the WOD.</figcaption></figure>
<figure><img src="image_22.png" /></figure>

Mean volume of low Ω_A

<figure><img src="image_23.png" /></figure>

Mean volume of low O₂

<figure><img src="image_24.png" /></figure>

Volume-averaged temperature

<figure><img src="image_25.png" /></figure>

Volume fraction anomaly of aragonite undersaturation

<figure><img src="image_26.png" /></figure>

Volume fraction anomaly of low oxygen

<figure><img src="image_27.png" /></figure>

Trend of temperature

<figure><img src="image_28.png" /></figure>

Trend of low-Ω_A vol. fraction

<figure><img src="image_29.png" /></figure>

Trend of low-O₂ vol. fraction

<figure><img src="image_30.png" /><figcaption>Figure 22: Top row: mean temperature and temporal trend (°C year⁻¹) of water on the shelf (−z < 200 m). Middle and bottom rows: volume of corrosive and hypoxic water over time and latitude. Water volumes are computed as sum of grid cell volumes with O₂ < 100 mmol m⁻³ (bottom) and Ω_A < 1 (middle) that are on the shelf (−z < 200 m). Mean values for each month are shown in the left column, and anomalies, computed as a fractional deviation from the climatological mean monthly volume at each latitude, in the middle column. Trends over time are shown in the right column.</figcaption></figure>
652
# 4 Conclusions

653
We present model simulations of ecosystem and biogeochemical cycles in the CCS that reproduce
654
the broad patterns of processes and states observed in this region over the past couple of decades.
655
Our results demonstrate that productivity of the CCS reflects a complex interplay of factors. The
656
limitation by the physical supply and removal of macronutrients (nitrate) provides the dominant
657
seasonal and spatial pattern of NPP, but with significant constraints from light and Fe on a seasonal
658
basis, especially in the northern CCS. Interannual variations in NPP are reasonably well predicted
659
by fluctuations in pycnocline depth that modulate the rates of surface nutrient supply. Expanded
660
datasets on near-surface Fe concentrations are needed to better establish its role as a limiting
661
factor for growth in the CCS. A significant correlation between model NPP and surface irradiance
662
suggests that changes in light are also influential. Together, these results suggest that an index of
663
NPP that accounts for both regional pycnocline structure and cloud cover would be more skillful
664
than one based only on coastal winds (e.g., Bakun (1990); Jacox et al. (2018)). Our results highlight
665
the value of continued measurements of the depth of the chlorophyll maximum.

666
Biogeochemical properties of subsurface waters in the CCS are also well reproduced by model
667
simulations. The amplitude of interannual variability in NO3- at the base of the photic zone and
668
of O2 in the thermocline are also both strongly correlated to undulations of the pycnocline. The
669
largest such anomalies in our simulation period were associated with the 1997-98 ENSO event,
670
whose amplitude of density and O2 anomalies remains coherent over a wide latitude band, albeit
671
with declining magnitude. For NO3-, the overall variance is somewhat lower, and the strength
672
of density correlations is somewhat higher in model output than in observations in the northern
673
domain. This suggests an important role for anomalies entering the CCS from the subarctic North
674
Pacific, an HNLC region. Basin-scale changes in biogeochemical properties are known to be ex-
675
ceptionally high at the gyre boundary ≈ 45°N (Mecking et al., 2008), and these remote anomalies
676
are likely to play an important but uncertain role in the variability observed in the CCS. Similarly,
677
inputs from terrestrial and riverine sources of nutrients and organic carbon could contribute to in-
678
terannual variability in the system that is currently difficult to constrain from observational data.
679
The impact of nearshore influences on the mean state suggests that its contribution to variability
680
could also be substantial. Evaluating these remote influences from both the open ocean and from
681
boundary inputs, using empirically-based time-dependent biogeochemical boundary conditions is
682
an important avenue for future research.

683
The variability of biogeochemical properties leads to significant changes in the volume of wa-
684
ters characterized by biologically stressful conditions of hypoxia and carbonate undersaturation. In
685
the volume anomalies for both habitat constraints, there is strong coherence across the CCS. Years
686
with unusually large volumes of hypoxic or corrosive water offer few obvious latitudinal refuges.
687
The onset of these conditions tends to propagate from the central CCS (≈ 40°N), arriving in the
688
northern CCS with a 2-3 month delay. Thus, monitoring hypoxia and CO2 system parameters in
689
the central CCS may offer some seasonal predictability for northern ecosystem impacts.
**Codes and Simulation Data** The physical and biogeochemical codes used for our simulations are at https://github.com/UCLA-ROMS/Code. Simulation model output archive data can be made available by email requests to the Corresponding Authors.

**Acknowledgments** This work was funded by grants from the National Science Foundation (OCE-1419323), the National Oceanic and Atmospheric Administration (DOC-NOAA NA15NOS4780186), and the Gordon and Betty Moore Foundation (GBMF #3775). Model simulations were carried out using the Extreme Science and Engineering Discovery Environment (XSEDE) and Yellowstone computers supported by the National Science Foundation at San Diego Supercomputer Center and NCAR, respectively. We gratefully acknowledge the PIs that supply high quality data to national centers, and to the programs that compile, manage and distribute these data for large-scale analyses, including the National Ocean Data Center (NOAA), the remote sensing products from NASA, and the California Cooperative Oceanic Fisheries Investigations (CalCOFI).
## 5 Appendix: Biogeochemical Model

Here, for completeness, we summarize the equations of the Biogeochemical Elemental Cycling (BEC) model in the implementation used for this work. This formulation is based on the original version presented in Moore et al. (2004).

### 5.1 Variables and parameters

#### 5.1.1 Prognostic variables

<table><thead><tr><th>Name</th><th>Description</th><th>Units</th></tr></thead><tbody><tr><td>N<sub>NH4</sub></td><td>ammonium</td><td>mmol N/m<sup>3</sup></td></tr><tr><td>N<sub>NO2</sub></td><td>nitrite</td><td>mmol N/m<sup>3</sup></td></tr><tr><td>N<sub>NO3</sub></td><td>nitrate</td><td>mmol N/m<sup>3</sup></td></tr><tr><td>N<sub>N2O</sub></td><td>nitrous oxide</td><td>mmol N/m<sup>3</sup></td></tr><tr><td>N<sub>N2</sub></td><td>nitrogen</td><td>mmol N/m<sup>3</sup></td></tr><tr><td>N<sub>don</sub></td><td>dissolved organic nitrogen</td><td>mmol N/m<sup>3</sup></td></tr><tr><td>P<sub>PO4</sub></td><td>Phosphorus</td><td>mmol P/m<sup>3</sup></td></tr><tr><td>P<sub>dop</sub></td><td>dissolved organic phosphorus</td><td>mmol P/m<sup>3</sup></td></tr><tr><td>O<sub>2</sub></td><td>dissolved oxygen</td><td>mmol O<sub>2</sub>/m<sup>3</sup></td></tr><tr><td>Fe</td><td>iron</td><td>nmol Fe/m<sup>3</sup></td></tr><tr><td>Fe<sub>sp</sub></td><td>small phytoplankton iron</td><td>mmol Fe/m<sup>3</sup></td></tr><tr><td>Fe<sub>diat</sub></td><td>diatom iron</td><td>mmol Fe/m<sup>3</sup></td></tr><tr><td>Fe<sub>diaz</sub></td><td>diazotroph iron</td><td>mmol Fe/m<sup>3</sup></td></tr><tr><td>Fe<sub>dofe</sub></td><td>dissolved organic iron</td><td>mmol Fe/m<sup>3</sup></td></tr><tr><td>Ca<sub>sp</sub></td><td>small phytoplankton calcium carbonate</td><td>mmol CaCO<sub>3</sub>/m<sup>3</sup></td></tr><tr><td>Si<sub>diat</sub></td><td>diatom silicate</td><td>mmol SiO<sub>2</sub>/m<sup>3</sup></td></tr><tr><td>Si<sub>SiO2</sub></td><td>Silicate</td><td>mmol SiO<sub>2</sub>/m<sup>3</sup></td></tr><tr><td>Alk</td><td>total alkalinity</td><td>mmol/m<sup>3</sup></td></tr><tr><td>C<sub>DIC</sub></td><td>dissolved inorganic carbon</td><td>mmol C/m<sup>3</sup></td></tr><tr><td>C<sub>sp</sub></td><td>small phytoplankton carbon</td><td>mmol C/m<sup>3</sup></td></tr><tr><td>C<sub>diat</sub></td><td>diatom carbon</td><td>mmol C/m<sup>3</sup></td></tr><tr><td>C<sub>diaz</sub></td><td>diazotroph carbon</td><td>mmol C/m<sup>3</sup></td></tr><tr><td>C<sub>doc</sub></td><td>dissolved organic carbon</td><td>mmol C/m<sup>3</sup></td></tr><tr><td>Chl<sub>sp</sub></td><td>small phytoplankton chlorophyll</td><td>mmol Chl/m<sup>3</sup></td></tr><tr><td>Chl<sub>diat</sub></td><td>diatom chlorophyll</td><td>mmol Chl/m<sup>3</sup></td></tr><tr><td>Chl<sub>diaz</sub></td><td>diazotroph chlorophyll</td><td>mmol Chl/m<sup>3</sup></td></tr><tr><td>C<sub>zoo</sub></td><td>zooplankton carbon</td><td>mmol C/m<sup>3</sup></td></tr><tr><td>POC<sub>sed</sub></td><td>organic carbon in sediment</td><td>mmol C/m<sup>2</sup></td></tr><tr><td>CaCO<sub>3, sed</sub></td><td>inorganic carbon in sediment</td><td>mmol C/m<sup>2</sup></td></tr><tr><td>Si<sub>sed</sub></td><td>Silicate in sediment</td><td>mmol Si/m<sup>2</sup></td></tr><tr><td>I<sub>SW</sub></td><td>Penetrative solar heat flux</td><td>W/m<sup>2</sup></td></tr></tbody></table>
709 5.1.2 Local variables

<table><thead><tr><th>Name</th><th>Description</th><th>Units</th></tr></thead><tbody><tr><td>J<sup>u</sup><sub>sp,no3</sub></td><td>nitrate uptake by small phytoplankton</td><td>mmol N/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>u</sup><sub>sp,nh4</sub></td><td>ammonium uptake by small phytoplankton</td><td>mmol N/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>u</sup><sub>diat,no3</sub></td><td>nitrate uptake by diatoms</td><td>mmol N/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>u</sup><sub>diat,nh4</sub></td><td>ammonium uptake by diatoms</td><td>mmol N/m<sup>3</sup>/sec</td></tr><tr><td>J<sub>ammox</sub></td><td>rate of NH<sub>4</sub> oxidation to NO<sub>2</sub></td><td>mmol N/m<sup>3</sup>/sec</td></tr><tr><td>J<sub>nit</sub></td><td>rate of NO<sub>2</sub> oxidation to NO<sub>3</sub></td><td>mmol N/m<sup>3</sup>/sec</td></tr><tr><td>J<sub>denit</sub><sub>no3,no2</sub></td><td>amount of NO<sub>3</sub> converted to NO<sub>2</sub> by water column denitrif.</td><td>mmol N/m<sup>3</sup>/sec</td></tr><tr><td>J<sub>denit</sub><sub>sed</sub></td><td>amount of NO<sub>3</sub> converted to N<sub>2</sub> by benthic denitrif.</td><td>mmol N/m<sup>3</sup>/sec</td></tr><tr><td>J<sub>denit</sub><sub>no2,n2o</sub></td><td>amount of NO<sub>2</sub> converted to N<sub>2</sub>O by water column denitrif.</td><td>mmol N/m<sup>3</sup>/sec</td></tr><tr><td>J<sub>denit</sub><sub>no2,n2</sub></td><td>amount of NO<sub>2</sub> converted to N<sub>2</sub> by water column denitrif.</td><td>mmol N/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>Nfix</sup><sub>diaz</sub></td><td>total N<sub>2</sub>-fixation by diazotrophs</td><td>mmol N/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>Nexcrete</sup><sub>diaz</sub></td><td>N excreted by diazotrophs</td><td>mmol N/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>photoN</sup><sub>diaz</sub></td><td>N fixed by diazotroph (non-excreted)</td><td>mmol N/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>u</sup><sub>sp,po4</sub></td><td>phosphate uptake by small phytoplankton</td><td>mmol P/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>u</sup><sub>diat,po4</sub></td><td>phosphate uptake by diatoms</td><td>mmol P/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>u</sup><sub>diaz,po4</sub></td><td>phosphate uptake by diazotrophs</td><td>mmol P/m<sup>3</sup>/sec</td></tr><tr><td>P<sup>remain</sup><sub>diaz</sub></td><td>remaining diazotrophs phosphate</td><td>mmol P/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>loss</sup><sub>diaz,dop</sub></td><td>non-grazing mortality of diazotrophs routed to POC</td><td>mmol P/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>loss</sup><sub>diaz,dip</sub></td><td>non-grazing mortality of diazotrophs routed to PIC</td><td>mmol P/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>prod</sup><sub>sp,caco3</sub></td><td>CaCO<sub>3</sub> production by small phytoplankton</td><td>mmol CaCO<sub>3</sub>/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>photo</sup><sub>sp,C</sub></td><td>carbon uptake by photosynthesis in small phytoplankton</td><td>mmol C/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>photo</sup><sub>diat,C</sub></td><td>carbon uptake by photosynthesis in diatoms</td><td>mmol C/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>photo</sup><sub>diaz,C</sub></td><td>carbon uptake by photosynthesis in diazotrophs</td><td>mmol C/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>grz</sup><sub>sp</sub></td><td>grazing loss for small phytoplankton</td><td>mmol C/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>grz</sup><sub>sp,dic</sub></td><td>grazed small phytoplankton routed to DIC</td><td>mmol C/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>grz</sup><sub>sp,doc</sub></td><td>grazed small phytoplankton routed to DOC</td><td>mmol C/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>grz</sup><sub>sp,poc</sub></td><td>grazed small phytoplankton routed to POC</td><td>mmol C/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>grz</sup><sub>sp,zoo</sub></td><td>grazed small phytoplankton routed to new zooplankton biomass</td><td>mmol C/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>grz</sup><sub>diat</sub></td><td>grazing loss for diatoms</td><td>mmol C/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>grz</sup><sub>diat,dic</sub></td><td>grazed diatoms routed to DIC</td><td>mmol C/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>grz</sup><sub>diat,doc</sub></td><td>grazed diatoms routed to DOC</td><td>mmol C/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>grz</sup><sub>diat,poc</sub></td><td>grazed diatoms routed to POC</td><td>mmol C/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>grz</sup><sub>diat,zoo</sub></td><td>grazed diatoms routed to new zooplankton biomass</td><td>mmol C/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>grz</sup><sub>diaz</sub></td><td>grazing loss for diazotrophs</td><td>mmol C/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>grz</sup><sub>diaz,dic</sub></td><td>grazed diazotrophs routed to DIC</td><td>mmol C/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>grz</sup><sub>diaz,doc</sub></td><td>grazed diazotrophs routed to DOC</td><td>mmol C/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>grz</sup><sub>diaz,poc</sub></td><td>grazed diazotrophs routed to POC</td><td>mmol C/m<sup>3</sup>/sec</td></tr><tr><td>J<sup>grz</sup><sub>diaz,zoo</sub></td><td>grazed diazotrophs routed to new zoo biomass</td><td>mmol C/m<sup>3</sup>/sec</td></tr></tbody></table>

710
<table>
  <tr>
    <td>J<sub>sp</sub><sup>agg</sup></td>
    <td>aggregation loss of small phytoplankton</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp</sub><sup>agg</sup></td>
    <td>aggregation loss of diatoms</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp</sub><sup>l</sup></td>
    <td>non-grazing mortality of small phytoplankton</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,dic</sub><sup>l</sup></td>
    <td>non-grazing mortality of small phytoplankton routed to DIC</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,dic</sub><sup>l</sup></td>
    <td>non-grazing mortality of small phytoplankton routed to DOC</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,poc</sub><sup>l</sup></td>
    <td>non-grazing mortality of small phytoplankton routed to POC</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,poc</sub><sup>l</sup></td>
    <td>non-grazing mortality of diatoms</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,dic</sub><sup>l</sup></td>
    <td>non-grazing mortality of diatoms routed to DIC</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,dic</sub><sup>l</sup></td>
    <td>non-grazing mortality of diatoms routed to DOC</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,poc</sub><sup>l</sup></td>
    <td>non-grazing mortality of diatoms routed to POC</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,dic</sub><sup>l</sup></td>
    <td>non-grazing mortality of diazotrophs</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,dic</sub><sup>l</sup></td>
    <td>non-grazing mortality of diazotrophs routed to DIC</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,dic</sub><sup>l</sup></td>
    <td>non-grazing mortality of diazotrophs routed to DOC</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,poc</sub><sup>l</sup></td>
    <td>non-grazing mortality of diazotrophs routed to POC</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,zoo</sub><sup>l</sup></td>
    <td>zooplankton mortality</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,zoo,dic</sub><sup>l</sup></td>
    <td>zooplankton mortality routed to DIC</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,zoo,dic</sub><sup>l</sup></td>
    <td>zooplankton mortality routed to DOC</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,zoo,poc</sub><sup>l</sup></td>
    <td>zooplankton mortality routed to POC</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,zoo,poc</sub><sup>l</sup></td>
    <td>fractional factor for routing of zoo losses</td>
    <td>no units</td>
  </tr>
  <tr>
    <td>J<sub>sp</sub><sup>agg</sup></td>
    <td>aggregation of small phytoplankton</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp</sub><sup>agg</sup></td>
    <td>aggregation of diatoms</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp</sub><sup>agg</sup></td>
    <td>aggregation of diazotrophs</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,dic</sub><sup>don,nh4</sup></td>
    <td>amount of dissolved organic N remineralized</td>
    <td>mmol N/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,dic</sub><sup>don,nh4</sup></td>
    <td>amount of dissolved organic C remineralized</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,dic</sub><sup>don,nh4</sup></td>
    <td>amount of dissolved organic P remineralized</td>
    <td>mmol P/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,dic</sub><sup>other</sup></td>
    <td>amount of dissolved Fe remineralized</td>
    <td>mmol Fe/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>sp,remin</sub></td>
    <td>remineralization in the sediments by processes other than oxic remin. and denitrif.</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>poc</sub><sup>prod</sup></td>
    <td>amount of particulate organic C produced</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>poc</sub><sup>prod</sup></td>
    <td>amount of particulate inorganic C produced</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>pic</sub><sup>prod</sup></td>
    <td>amount of particulate organic Si produced</td>
    <td>mmol Si/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>pic</sub><sup>prod</sup></td>
    <td>amount of particulate organic Fe produced</td>
    <td>mmol Fe/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>prod</sub><sup>prod</sup></td>
    <td>remineralized particulate organic C</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>prod</sub><sup>prod</sup></td>
    <td>remineralized particulate inorganic C</td>
    <td>mmol C/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>prod</sub><sup>prod</sup></td>
    <td>remineralized particulate organic Si</td>
    <td>mmol SiO<sub>2</sub>/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>prod</sub><sup>prod</sup></td>
    <td>remineralized particulate organic Fe</td>
    <td>mmol Fe/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>prod</sub><sup>scav</sup></td>
    <td>scavenging of inorganic Fe</td>
    <td>mmol Fe/m<sup>3</sup>/sec</td>
  </tr>
  <tr>
    <td>J<sub>prod</sub><sup>scav</sup></td>
    <td>remineralized Fe from sediment</td>
    <td>mmol Fe/m<sup>3</sup>/sec</td>
  </tr>
</table>
<table><tr><td>PAR</td><td>photosynthetically active radiation</td><td>W/m<sup>2</sup></td></tr><tr><td>Φ<sup>mineral</sup><sub>poc</sub></td><td>incoming mineral associated particulate C flux</td><td>mmol C cm/m<sup>3</sup>/sec</td></tr><tr><td>Φ<sup>out</sup><sub>poc</sub></td><td>particulate C-flux buried in the sediments</td><td>mmol C /m<sup>2</sup>/sec</td></tr><tr><td>Φ<sup>free</sup><sub>poc</sub></td><td>free associated particulate C flux</td><td>mmol C /m<sup>2</sup>/sec</td></tr><tr><td>Φ<sup>sol</sup><sub>pic</sub></td><td>soluble and pic associated particulate C flux</td><td>mmol C /m<sup>2</sup>/sec</td></tr><tr><td>Φ<sup>nonsol</sup><sub>pic</sub></td><td>non-soluble and pic associated particulate C flux</td><td>mmol C /m<sup>2</sup>/sec</td></tr><tr><td>Φ<sup>sol</sup><sub>psi</sub></td><td>soluble and SiO2 associated particulate C flux</td><td>mmol SiO<sub>2</sub> /m<sup>2</sup>/sec</td></tr><tr><td>Φ<sup>nonsol</sup><sub>psi</sub></td><td>non-soluble and SiO2 associated particulate C flux</td><td>mmol SiO<sub>2</sub> /m<sup>2</sup>/sec</td></tr><tr><td>Φ<sup>sol</sup><sub>dust</sub></td><td>soluble and dust associated particulate C flux</td><td>mmol C /m<sup>2</sup>/sec</td></tr><tr><td>Φ<sup>nonsol</sup><sub>dust</sub></td><td>non-soluble and dust associated particulate C flux</td><td>mmol C /m<sup>2</sup>/sec</td></tr><tr><td>Φ<sup>sol</sup><sub>pfe</sub></td><td>soluble and dust associated particulate C flux</td><td>mmol Fe /m<sup>2</sup>/sec</td></tr><tr><td>Φ<sup>nonsol</sup><sub>pfe</sub></td><td>non-soluble and dust associated particulate C flux</td><td>mmol Fe /m<sup>2</sup>/sec</td></tr><tr><td>Φ<sup>remin,sed</sup><sub>poc</sub></td><td>remineralized organic C flux from sediment</td><td>mmol C/m<sup>2</sup>/sec</td></tr><tr><td>Φ<sup>remin,sed</sup><sub>pic</sub></td><td>remineralized inorganic C flux from sediment</td><td>mmol C/m<sup>2</sup>/sec</td></tr><tr><td>Φ<sup>remin,sed</sup><sub>Si</sub></td><td>remineralized Si flux from sediment</td><td>mmol C/m<sup>2</sup>/sec</td></tr><tr><td>Φ<sup>air</sup><sub>O2</sub></td><td>O<sub>2</sub> air-sea flux</td><td>mmol O<sub>2</sub>/m<sup>2</sup>/sec</td></tr><tr><td>Φ<sup>air</sup><sub>N2O</sub></td><td>N<sub>2</sub>O air-sea flux</td><td>mmol N/m<sup>2</sup>/sec</td></tr><tr><td>Φ<sup>air</sup><sub>N2</sub></td><td>N<sub>2</sub> air-sea flux</td><td>mmol N/m<sup>2</sup>/sec</td></tr><tr><td>Φ<sup>air</sup><sub>CO2</sub></td><td>CO<sub>2</sub> air-sea flux</td><td>mmol C/m<sup>2</sup>/sec</td></tr><tr><td>Sc</td><td>Schmidt Number</td><td>-</td></tr><tr><td>Ws</td><td>Wind speed at 10 m</td><td>m/sec</td></tr><tr><td>PV</td><td>Piston Velocity</td><td>m/sec</td></tr></table>

## 713 5.2 Ecosystem parameters

<table><thead><tr><th>Parameters</th><th>Description</th><th>Values</th><th>Units</th></tr></thead><tbody><tr><td colspan="4"><strong>Model grid</strong></td></tr><tr><td>Δz</td><td>Model layer thickness</td><td>variable</td><td>meters</td></tr><tr><td>k<sub>ρ</sub></td><td>Index of model vertical level at tracer points</td><td>1-60</td><td>no units</td></tr><tr><td colspan="4"><strong>Carbon</strong></td></tr><tr><td>PC<sup>sp</sup><sub>ref</sub></td><td>max phyto C-specific growth rate at Tref (GD98) for small phytoplankton</td><td>3.0</td><td>1/d</td></tr><tr><td>PC<sup>diat</sup><sub>ref</sub></td><td>max phyto C-specific growth rate at Tref (GD98) for diatoms</td><td>3.0</td><td>1/d</td></tr></tbody></table>
<table><tr><td>PC<sup>diaz</sup><sub>ref</sub></td><td>max phyto C-specific growth rate at Tref (GD98) for diazotrophs</td><td>0.4</td><td>1/d</td></tr><tr><td>α<sub>chl</sub></td><td>chlorophyll-specific initial slope of P vs. I curve for diatoms and small phytoplankton</td><td>0.3</td><td>mmol C m<sup>2</sup>/(mg Chl W day)</td></tr><tr><td>α<sup>diaz</sup><sub>chl</sub></td><td>chlorophyll-specific initial slope of P vs. I curve for diazotrophs</td><td>0.036</td><td>mmol C m<sup>2</sup>/(mg Chl W day)</td></tr></table>

Nutrient limitation

<table><tr><td>k<sup>sp</sup><sub>NO<sub>3</sub></sub></td><td>Small phyto. half saturation constant for NO<sub>3</sub> uptake</td><td>0.5</td><td>mmol N/m<sup>3</sup></td></tr><tr><td>k<sup>diat</sup><sub>NO<sub>3</sub></sub></td><td>Diatom half saturation constant for NO<sub>3</sub> uptake</td><td>2.5</td><td>mmol N/m<sup>3</sup></td></tr><tr><td>k<sup>sp</sup><sub>NH<sub>4</sub></sub></td><td>Small phyto. half saturation constant for NH<sub>4</sub> uptake</td><td>0.01</td><td>mmol N/m<sup>3</sup></td></tr><tr><td>k<sup>diat</sup><sub>NH<sub>4</sub></sub></td><td>Diatom half saturation constant for NH<sub>4</sub> uptake</td><td>0.1</td><td>mmol N/m<sup>3</sup></td></tr><tr><td>k<sup>sp</sup><sub>Fe</sub></td><td>Small phyto. half saturation constant for Fe uptake</td><td>0.035e-3</td><td>mmol Fe/m<sup>3</sup></td></tr><tr><td>k<sup>diat</sup><sub>Fe</sub></td><td>Diat half saturation constant for Fe uptake</td><td>0.08e-3</td><td>mmol Fe/m<sup>3</sup></td></tr><tr><td>k<sup>diaz</sup><sub>Fe</sub></td><td>Diazotroph half saturation constant for Fe uptake</td><td>0.1e-3</td><td>mmol Fe/m<sup>3</sup></td></tr><tr><td>k<sup>sp</sup><sub>PO<sub>4</sub></sub></td><td>Small phyto. half saturation constant for PO<sub>4</sub> uptake</td><td>0.01</td><td>mmol P/m<sup>3</sup></td></tr><tr><td>k<sup>diat</sup><sub>PO<sub>4</sub></sub></td><td>Diatom half saturation constant for PO<sub>4</sub> uptake</td><td>0.1</td><td>mmol P/m<sup>3</sup></td></tr><tr><td>k<sup>diaz</sup><sub>PO<sub>4</sub></sub></td><td>Diazotroph half saturation constant for PO<sub>4</sub> uptake</td><td>0.005</td><td>mmol P/m<sup>3</sup></td></tr></table>
<table><tr><td>k<sup>diaz</sup><sub>SiO<sub>2</sub></sub></td><td>Diatom half saturation constant for SiO<sub>2</sub> uptake</td><td>1</td><td>mmol SiO<sub>2</sub>/m<sup>3</sup></td></tr><tr><td colspan="4"><strong>Fixed Stoichiometry</strong></td></tr><tr><td>Q<sub>N:C</sub></td><td>Nitrogen to Carbon ratio</td><td>0.137</td><td>mmol N/mmol C</td></tr><tr><td>Q<sub>P:C</sub></td><td>Small phyto. and diatom P:C ratio</td><td>0.00855</td><td>mmol P/mmol C</td></tr><tr><td>Q<sup>diaz</sup><sub>P:C</sub></td><td>Diazotroph P:C ratio</td><td>0.002735</td><td>mmol P/mmol C</td></tr><tr><td colspan="4"><strong>Dissolved organic matter</strong></td></tr><tr><td>τ<sub>dom</sub></td><td>Dissolved organic matter remineralization inverse timescale</td><td>0.01</td><td>1/d</td></tr><tr><td>τ<sub>remin,sed</sub></td><td>remineralization inverse timescale in sediment</td><td>0.003</td><td>1/d</td></tr><tr><td colspan="4"><strong>N<sub>2</sub> fixation</strong></td></tr><tr><td>α<sup>diaz</sup><sub>ex</sub></td><td>Ratio of N excreted by diazotrophs to total N fixed</td><td>0.3</td><td>no units</td></tr><tr><td colspan="4"><strong>Chl stoichiometry and production</strong></td></tr><tr><td>Q<sup>sp,growth</sup><sub>N:chl,max</sub></td><td>Max ratio of Chl produced to N photosynthesized for small phyto.</td><td>2.5</td><td>mg Chl/mmol N</td></tr><tr><td>Q<sup>diaz,growth</sup><sub>N:chl,max</sub></td><td>Max ratio of Chl produced to N photosynthesized for diatom</td><td>4.0</td><td>mg Chl/mmol N</td></tr><tr><td>Q<sup>diaz,growth</sup><sub>N:chl,max</sub></td><td>Max ratio of Chl produced to N photosynthesized for diaz.</td><td>3.4</td><td>mg Chl/mmol N</td></tr></table>
Light

<table><tr><td>f<sub>Isw</sub></td><td>Fraction of incoming radiation used for photosynthesis</td><td>0.45</td><td>no units</td></tr><tr><td>μ<sub>w</sub></td><td>Radiation attenuation coefficient per unit chlorophyll</td><td>0.03</td><td>1/m/(mg Chl/m<sup>3</sup>)</td></tr><tr><td>μ<sub>Chl</sub></td><td>Radiation attenuation coefficient for water</td><td>0.04</td><td>1/m</td></tr></table>

Fe stoichiometry and cycling

<table><tr><td>Q<sup>sp,g0</sup><sub>Fe:C</sub></td><td>Ratio used in the calculation of the Fe to C ratio of uptake for small phyto.</td><td>6e-6</td><td>mmol Fe/mmol C</td></tr><tr><td>Q<sup>sp,g1</sup><sub>Fe:C</sub></td><td>Maximum Fe:C ratio of uptake by small phyto. if Fe &lt; 2 &middot; k<sup>sp</sup><sub>Fe</sub>.</td><td>2.5e-6</td><td>mmol Fe/mmol C</td></tr><tr><td>Q<sup>diat,g0</sup><sub>Fe:C</sub></td><td>Ratio used in the calculation of the Fe to C ratio of uptake for diatoms</td><td>6e-6</td><td>mmol Fe/mmol C</td></tr><tr><td>Q<sup>diat,g1</sup><sub>Fe:C</sub></td><td>Maximum Fe:C ratio of uptake by diatoms if Fe &lt; 2 &middot; k<sup>diat</sup><sub>Fe</sub></td><td>2.5e-6</td><td>mmol Fe/mmol C</td></tr><tr><td>Q<sup>diaz,g0</sup><sub>Fe:C</sub></td><td>Ratio used in the calculation of the Fe to C ratio of uptake for diazotrophs</td><td>42e-6</td><td>mmol Fe/mmol C</td></tr><tr><td>Q<sup>diaz,g1</sup><sub>Fe:C</sub></td><td>Maximum Fe:C ratio of uptake by diazotrophs if Fe &lt; 2 &middot; k<sup>diaz</sup><sub>Fe</sub></td><td>14e-6</td><td>mmol Fe/mmol C</td></tr><tr><td>β<sup>max,scale1</sup><sub>fe</sub></td><td>Scaling parameter used in the calculation of Fe scavenging</td><td>3.0</td><td>no units</td></tr><tr><td>β<sup>thresh,1</sup><sub>fe</sub></td><td>Fe threshold parameter used in the calculation of Fe scavenging</td><td>0.6e-3</td><td>mmol Fe/m<sup>3</sup></td></tr></table>
<table><tr><td>β<sub>fe</sub><sup>thresh,2</sup></td><td>Fe threshold parameter used in the calculation of Fe scavenging</td><td>0.5e-3</td><td>mmol Fe/m<sup>3</sup></td></tr><tr><td>J<sub>Fe</sub><sup>scaav,0</sup></td><td>Fe scavenging reference rate</td><td>0.12</td><td>mmolFe/m<sup>3</sup>/s</td></tr></table>

## CaCO<sub>3</sub> stoichiometry and production

<table><tr><td>Q<sub>CaCO<sub>3</sub>:C</sub><sup>sp,max</sup></td><td>Maximum calcification to C photosynthesis ratio</td><td>0.4</td><td>mmol CaCO<sub>3</sub>/mmol C</td></tr><tr><td>f<sub>CaCO<sub>3</sub></sub><sup>prod</sup></td><td>Initial calcification to C photosynthesis ratio</td><td>0.026</td><td>mmol CaCO<sub>3</sub>/mmol C</td></tr><tr><td>T<sub>1</sub><sup>CaCO<sub>3</sub></sup></td><td>Temperature parameter used to modify the initial calcification rate</td><td>1.</td><td>°C</td></tr><tr><td>T<sub>2</sub><sup>CaCO<sub>3</sub></sup></td><td>Temperature parameter used to modify the initial calcification rate</td><td>-2.</td><td>°C</td></tr></table>

## Si stoichiometry

<table><tr><td>Q<sub>Si:C</sub><sup>diat,max</sup></td><td>Maximum Q<sub>Si:C</sub> ratio for diatoms</td><td>0.685</td><td>mmol Si/mmol C</td></tr><tr><td>Q<sub>Si:C</sub><sup>diat,g0</sup></td><td>Default Si:C ratio of growth for diatoms</td><td>0.137</td><td>mmol Si/mmol C</td></tr><tr><td>Q<sub>Si:C</sub><sup>diat,g1</sup></td><td>Ratio used to calculate the Si to C multiplicative ratio of growth for diatoms</td><td>2.5</td><td>no units</td></tr></table>

## Grazing

<table><tr><td>J<sub>sp</sub><sup>g,max</sup></td><td>maximum grazing loss for small phytoplankton</td><td>2.5</td><td>1/d</td></tr><tr><td>J<sub>diat</sub><sup>g,max</sup></td><td>maximum grazing loss for diatoms</td><td>1.95</td><td>1/d</td></tr><tr><td>J<sub>diaz</sub><sup>g,max</sup></td><td>maximum grazing loss for diazotrophs</td><td>1.2</td><td>1/d</td></tr></table>
<table>
  <tr>
    <td>β<sub>z</sub><sup>grz</sup></td>
    <td>grazing coefficient, used in density dependent grazing modification</td>
    <td>1.05</td>
    <td>mmol C/m<sup>3</sup></td>
  </tr>
  <tr>
    <td>β<sub>thres</sub><sup>grz,sp,0</sup></td>
    <td>Small phytoplankton threshold concentration for grazing</td>
    <td>0.001</td>
    <td>mmol C/m<sup>3</sup></td>
  </tr>
  <tr>
    <td>β<sub>thres</sub><sup>grz,diat,0</sup></td>
    <td>Diatom threshold concentration for grazing</td>
    <td>0.02</td>
    <td>mmol C/m<sup>3</sup></td>
  </tr>
  <tr>
    <td>β<sub>thres</sub><sup>grz,diaz,0</sup></td>
    <td>Diazotroph threshold concentration for grazing</td>
    <td>0.01</td>
    <td>mmol C/m<sup>3</sup></td>
  </tr>
  <tr>
    <td>α<sub>sp</sub><sup>grz,zoo</sup></td>
    <td>Fraction of small phyto. grazing going to zooplankton</td>
    <td>0.3</td>
    <td>no units</td>
  </tr>
  <tr>
    <td>α<sub>diat</sub><sup>grz,zoo</sup></td>
    <td>Fraction of diatom grazing going to zooplankton</td>
    <td>0.3</td>
    <td>no units</td>
  </tr>
  <tr>
    <td>α<sub>diaz</sub><sup>grz,zoo</sup></td>
    <td>Fraction of diazotroph grazing going to zooplankton</td>
    <td>0.21</td>
    <td>no units</td>
  </tr>
  <tr>
    <td>α<sub>sp</sub><sup>grz,poc</sup></td>
    <td>Default fraction of small phytop. grazing going to POC</td>
    <td>0.22</td>
    <td>no units</td>
  </tr>
  <tr>
    <td>α<sub>diat</sub><sup>grz,poc</sup></td>
    <td>Fraction of diatom grazing going to POC</td>
    <td>0.26</td>
    <td>no units</td>
  </tr>
  <tr>
    <td>α<sub>diaz</sub><sup>grz,poc</sup></td>
    <td>Fraction of diazotroph grazing going to POC</td>
    <td>0.0</td>
    <td>no units</td>
  </tr>
  <tr>
    <td>α<sub>sp</sub><sup>grz,doc</sup></td>
    <td>Fraction of small phyto. grazing going to DOC</td>
    <td>0.34</td>
    <td>no units</td>
  </tr>
  <tr>
    <td>α<sub>diat</sub><sup>grz,doc</sup></td>
    <td>Fraction of diatom grazing going to DOC</td>
    <td>0.13</td>
    <td>no units</td>
  </tr>
  <tr>
    <td>α<sub>diaz</sub><sup>grz,doc</sup></td>
    <td>Fraction of diazotroph grazing going to DOC</td>
    <td>0.24</td>
    <td>no units</td>
  </tr>
  <tr>
    <td>α<sub>sp</sub><sup>grz,dic</sup></td>
    <td>Fraction of small phytop. grazing going to DIC</td>
    <td>0.36</td>
    <td>no units</td>
  </tr>
  <tr>
    <td>α<sub>diat</sub><sup>grz,dic</sup></td>
    <td>Fraction of diatom grazing going to DIC</td>
    <td>0.31</td>
    <td>no units</td>
  </tr>
  <tr>
    <td>α<sub>diaz</sub><sup>grz,dic</sup></td>
    <td>Fraction of diazotroph grazing going to DIC</td>
    <td>0.55</td>
    <td>no units</td>
  </tr>
</table>

Losses
<table><tr><td>λ<sup>mort</sup><sub>sp</sub></td><td>Small phyto. mortality 1</td><td>0.15</td><td>1/d</td></tr><tr><td>λ<sup>mort</sup><sub>diat</sub></td><td>Diatom mortality 1</td><td>0.15</td><td>1/d</td></tr><tr><td>λ<sup>mort</sup><sub>diaz</sub></td><td>Diazotroph mortality 1</td><td>0.16</td><td>1/d</td></tr><tr><td>α<sup>l,poc</sup><sub>diaz</sub></td><td>Fraction of diatom loss going to POC</td><td>0.05</td><td>no units</td></tr><tr><td>α<sup>l,poc</sup><sub>diaz</sub></td><td>Fraction of diazotroph loss going to POC</td><td>0</td><td>no units</td></tr><tr><td>f<sub>labile</sub></td><td>Fraction of Labile dissolved organic matter for loss calculations</td><td>0.70</td><td>no units</td></tr><tr><td>λ<sup>mort</sup><sub>zoo</sub></td><td>Zooplankton linear mortality</td><td>0.08</td><td>1/d</td></tr><tr><td>λ<sup>mort2</sup><sub>zoo</sub></td><td>Zooplankton quadratic mortality</td><td>0.42</td><td>1/(mmol C m<sup>3</sup> d)</td></tr><tr><td>β<sup>thres0,l</sup><sub>zoo</sub></td><td>Zooplankton threshold concentrations for mortality</td><td>0.03</td><td>mmol C/m<sup>3</sup></td></tr></table>

## Aggregation

<table><tr><td>τ<sup>agg,min</sup><sub>diat</sub></td><td>Minimum aggregation rate for diatoms</td><td>0.01</td><td>1/d</td></tr><tr><td>τ<sup>agg,min</sup><sub>diaz</sub></td><td>Minimum aggregation rate for diazotrophs</td><td>0.01</td><td>1/d</td></tr><tr><td>τ<sup>agg,max</sup><sub>sp</sub></td><td>Maximum aggregation rate for small phyto.</td><td>0.75</td><td>1/d</td></tr><tr><td>τ<sup>agg,max</sup><sub>diat</sub></td><td>Maximum aggregation rate for diatoms</td><td>0.75</td><td>1/d</td></tr><tr><td>τ<sup>agg,max</sup><sub>diaz</sub></td><td>Maximum aggregation rate for diazotrophs</td><td>0.75</td><td>1/d</td></tr><tr><td>λ<sup>mort2</sup><sub>sp</sub></td><td>Small phyto. quadratic mortality</td><td>0.0035</td><td>1/(mmol C m<sup>3</sup> d)</td></tr><tr><td>λ<sup>mort2</sup><sub>diat</sub></td><td>diatom quadratic mortality</td><td>0.0035</td><td>1/(mmol C m<sup>3</sup> d)</td></tr><tr><td>λ<sup>mort2</sup><sub>diaz</sub></td><td>diazotroph quadratic mortality</td><td>0.16</td><td>1/(mmol C m<sup>3</sup> d)</td></tr></table>

## N cycle rates

<table><tr><td>τ<sup>ammox</sup></td><td>NH<sub>4</sub> oxidation inverse timescale</td><td>0.06</td><td>1/d</td></tr><tr><td>τ<sup>nit</sup></td><td>NO<sub>2</sub> oxidation inverse timescale</td><td>0.33</td><td>1/d</td></tr></table>
<table>
  <tr>
    <td>PAR<sub>max</sub><sup>nitrif</sup></td>
    <td>Light threshold for onset of nitrification</td>
    <td>4.0</td>
    <td>W / m<sup>2</sup></td>
  </tr>
  <tr>
    <td>Q<sub>N:C</sub><sup>denit</sup></td>
    <td>Ratio of inorganic N consumed to POC mineralized during denitrification</td>
    <td>104/106</td>
    <td>mmol N / mmol C</td>
  </tr>
  <tr>
    <td>f<sub>n2o</sub><sup>denit</sup></td>
    <td>Fraction of denitrification that goes to N<sub>2</sub>O</td>
    <td>0.95</td>
    <td>no units</td>
  </tr>
  <tr>
    <td>τ<sub>n2o</sub><sup>cons</sup></td>
    <td>N<sub>2</sub>O reduction inverse timescale</td>
    <td>0.33333</td>
    <td>1/d</td>
  </tr>
  <tr>
    <td>a<sub>0</sub><sup>bohl</sup></td>
    <td>Parameter used to calculate sedimentary denitrification (Bohlen et al., 2012)</td>
    <td>0.06</td>
    <td>no units</td>
  </tr>
  <tr>
    <td>a<sub>1</sub><sup>bohl</sup></td>
    <td>Parameter used to calculate sedimentary denitrification (Bohlen et al., 2012)</td>
    <td>0.19</td>
    <td>no units</td>
  </tr>
  <tr>
    <td>a<sub>2</sub><sup>bohl</sup></td>
    <td>Parameter used to calculate sedimentary denitrification (Bohlen et al., 2012)</td>
    <td>0.99</td>
    <td>no units</td>
  </tr>
</table>

Particle cycling

<table>
  <tr>
    <td>ρ<sub>pic</sub></td>
    <td>Organic carbon to inorganic carbon mass ratio in PIC</td>
    <td>0.07 · <sup>100.09</sup>/<sub>12.01</sub></td>
    <td>no units</td>
  </tr>
  <tr>
    <td>ρ<sub>psi</sub></td>
    <td>Organic carbon to inorganic Si mass ratio in PSiO<sub>2</sub></td>
    <td>0.035 · <sup>60.08</sup>/<sub>12.01</sub></td>
    <td>mmol SiO<sub>2</sub> / mmol C</td>
  </tr>
  <tr>
    <td>ρ<sub>dust</sub></td>
    <td>Organic carbon to dust mass ratio in dust</td>
    <td>0.07 · <sup>10<sup>6</sup></sup>/<sub>12.01</sub></td>
    <td>no units</td>
  </tr>
  <tr>
    <td>λ<sub>poc</sub></td>
    <td>remin. length scale for poc, modified by T</td>
    <td>130</td>
    <td>meters</td>
  </tr>
  <tr>
    <td>λ<sub>pic</sub></td>
    <td>remin. length scale for pic</td>
    <td>600</td>
    <td>meters</td>
  </tr>
  <tr>
    <td>λ<sub>psi</sub></td>
    <td>remin. length scale for psi, modified by T</td>
    <td>210</td>
    <td>meters</td>
  </tr>
  <tr>
    <td>λ<sub>dust</sub></td>
    <td>remin. length scale for dust</td>
    <td>600</td>
    <td>meters</td>
  </tr>
  <tr>
    <td>λ<sub>nonsol</sub></td>
    <td>remin. length scale for non soluble material</td>
    <td>40000</td>
    <td>meters</td>
  </tr>
  <tr>
    <td>Q<sub>10,poc</sub></td>
    <td>Temperature dependency factor for poc, non soluble material</td>
    <td>1.12</td>
    <td>no units</td>
  </tr>
</table>
<table>
  <tr>
    <td>Q<sub>10,psi</sub></td>
    <td>Temperature dependency factor for psi<br>remin. length scale</td>
    <td>4.0</td>
    <td>no units</td>
  </tr>
  <tr>
    <td>γ<sub>pic</sub></td>
    <td>Production fraction<br>routed to the nonsoluble subclass</td>
    <td>0.55</td>
    <td>no units</td>
  </tr>
  <tr>
    <td>γ<sub>psi</sub></td>
    <td>Production fraction<br>routed to the nonsoluble subclass</td>
    <td>0.37</td>
    <td>no units</td>
  </tr>
  <tr>
    <td>Φ<sup>ref</sup><sub>poc</sub></td>
    <td>particulate C flux of<br>reference</td>
    <td>2.0e-3</td>
    <td>mmolC/m<sup>2</sup>/sec</td>
  </tr>
</table>

Oxygen

<table>
  <tr>
    <td>O<sub>2</sub><sup>min</sup></td>
    <td>Oxygen minimum</td>
    <td>1.0</td>
    <td>mmol O<sub>2</sub>/m<sup>3</sup></td>
  </tr>
  <tr>
    <td></td>
    <td>threshold</td>
    <td></td>
    <td></td>
  </tr>
</table>

714

**5.3 Model equations**

715

**5.3.1 Tracer equations**

716

Here the symbol *d/dt* denotes the sum of the local time derivative and the physical transport.

$$
\frac{d}{dt}(N_{no3}) = J^{nit} - J_{no3,no2}^{denit} - J_{sed}^{denit} - (J_{diat,no3}^{u} + J_{sp,no3}^{u}) \quad (A1)
$$

717

718

$$
\begin{equation}
\begin{split}
\frac{d}{dt}(N_{nh4}) = {}& Q_{N:C}J_{poc}^{remin} + J_{don,nh4} + Q_{N:C}(J_{diat,dic}^{l} + J_{sp,dic}^{l} \\
& \qquad + J_{diaz,dic}^{l}) + Q_{N:C}J_{zoo,dic}^{l} + Q_{N:C}(J_{diat,dic}^{grz} + J_{sp,dic}^{grz} + J_{diaz,dic}^{grz}) \\
& \qquad - J^{ammox} - (J_{diat,nh4}^{u} + J_{sp,nh4}^{u}) + Q_{N:C}\Phi_{poc}^{remin, sed}/\Delta z
\end{split}
\tag{A2}
\end{equation}
$$

719

720

$$
\frac{d}{dt}(N_{no2}) = J^{ammox} - J^{nit} - J_{no2,n2o}^{denit} + J_{no3,no2}^{denit} \quad (A3)
$$

721

722

$$
\frac{d}{dt}(N_{n2o}) = 0.5 J_{no2,n2o}^{denit} - J_{n2o,n2}^{denit} + \Phi_{N_2O}^{air} / \Delta z \quad (A4)
$$

723

724

$$
\frac{d}{dt}(N_{n2}) = J_{n2o,n2}^{denit} + 0.5 J_{sed}^{denit} + \Phi_{N_2}^{air} / \Delta z \quad (A5)
$$
$$
\begin{align}
\frac{d}{dt}(Fe) ={}& J_{dofe,fe} - J_{Fe}^{scau} + (Q_{Fe:C}^{zoo} J_{zoo,dic}^{l} + J_{pFe}^{remin} - J_{diaz,Fe}^{u} \nonumber \\
& + Q_{Fe:C}^{sp}(J_{sp,dic}^{l} + J_{sp,dic}^{grz}) + Q_{Fe:C}^{diat}(J_{diat,dic}^{l} + J_{diat,dic}^{grz}) \nonumber \\
& + Q_{Fe:C}^{diaz}(J_{diaz,dic}^{l} + J_{diaz,dic}^{grz}) - J_{sp,Fe}^{u} - J_{diat,Fe}^{u} \nonumber \\
& + J_{sp,zoo}^{grz}(Q_{Fe:C}^{sp} - Q_{Fe:C}^{zoo}) + J_{diat,zoo}^{grz}(Q_{Fe:C}^{diat} - Q_{Fe:C}^{zoo}) \nonumber \\
& + J_{diaz,zoo}^{grz}(Q_{Fe:C}^{diaz} - Q_{Fe:C}^{zoo}) + J_{Fe}^{Sed} \tag{A6}
\end{align}
$$

$$
\begin{align}
\frac{d}{dt}(P_{po4}) ={}& Q_{P:C}(J_{poc}^{remin} + J_{zoo,dic}^{l} + J_{sp,dic}^{l} + J_{diat,dic}^{l} + J_{sp,dic}^{grz} + J_{diat,dic}^{grz}) \nonumber \\
& + J_{diaz,dip}^{loss} - (J_{sp,P}^{u} + J_{diat,P}^{u} + J_{diaz,P}^{u}) + J_{dop,po4} + Q_{P:C} \frac{\Phi_{poc}^{remin,sed}}{\Delta z} \tag{A7}
\end{align}
$$

$$
\begin{align}
\frac{d}{dt}(Si_{sio2}) ={}& Q_{Si:C}(0.5 J_{diat}^{grz} + 0.95 J_{diat}^{l}) \nonumber \\
& - J_{diat,Si}^{u} + J_{pSi}^{remin} + \Phi_{Si}^{remin,sed} / \Delta z \tag{A8}
\end{align}
$$

If $(O_2 > O_{2min})$

$$
\begin{equation}
\begin{aligned}
\frac{d}{dt}(O_2) = & \left( J_{sp,C}^{\text{photo}} + J_{diat,C}^{\text{photo}} + J_{diaz,C}^{\text{photo}} \right) / R_{C:O}^d + \left( \begin{aligned}[t]
& -(-J_{poc}^{\text{remin}} - J_{doc,dic} - J_{zoo,dic} - J_{sp,dic}^{l} \\
& - J_{sp,dic}^{grz} - J_{diat,dic}^{l} - J_{diat,dic}^{grz} - J_{diaz,dic}^{l} - J_{diaz,dic}^{grz} - J_{diaz,dic}^{grz}) / R_{C:O}^d \\
& - 1.5 J^{\text{ammox}} - 0.5 J^{\text{nit}} - \Phi_{poc}^{\text{remin,sed}} / (R_{C:O}^d \Delta z) + \Phi_{O_2}^{\text{air}} / \Delta z
\end{aligned} \right)
\end{aligned}
\tag{A9}
\end{equation}
$$

If $(O_2 \le O_{2min})$

$$
\frac{d}{dt}(O_2) = \frac{(J_{sp,C}^{\text{photo}} + J_{diat,C}^{\text{photo}} + J_{diaz,C}^{\text{photo}})}{R_{C:O}^{d} + \Phi_{O_2}^{\text{air}} / \Delta z} \quad (A10)
$$

$$
\begin{align}
\frac{d}{dt}(DIC) ={}& J_{doc,dic} + 0.33 \quad J_{sp}^{grz} \quad Q_{C:caco3} \nonumber \\
& + J_{zoo,dic}^{l} + J_{sp,dic}^{l} + J_{sp,dic}^{grz} + J_{diat,dic}^{l} + J_{diat,dic}^{grz} - J_{sp,C}^{photo} - J_{diat,C}^{photo} \nonumber \\
& - J_{sp,caco3}^{prod} + J_{diaz,dic}^{grz} + J_{diaz,dic}^{l} - J_{diaz,C}^{photo} + (J_{poc}^{remin} + J_{pic}^{remin}) \nonumber \\
& + (\Phi_{poc}^{remin,sed} + \Phi_{pic}^{remin,sed}) / \Delta z + \Phi_{CO_2}^{air} / \Delta z \tag{A11}
\end{align}
$$
$$ \frac{d}{dt}(Alk) = -\frac{dN_{no3}}{dt} + \frac{dN_{nh4}}{dt} + 2J_{pic}^{remin} + 0.33 \quad J_{sp}^{grz} \quad Q_{C:caco3} - J_{sp,caco3}^{prod} \\ + 2\Phi_{pic}^{remin, sed} / \Delta z \quad (A12) $$

$$ \frac{d}{dt}(C_{sp}) = J_{sp,C}^{photo} - (J_{sp}^{grz} + J_{sp}^{l} + J_{sp}^{agg}) \quad (A13) $$

$$ \frac{d}{dt}(Chl_{sp}) = J_{sp,chl}^{pa} - Q_{Chl:C}^{sp}(J_{sp}^{grz} + J_{sp}^{l} + J_{sp}^{agg}) \quad (A14) $$

$$ \frac{d}{dt}(Fe_{sp}) = gQ_{Fe:C}^{sp}J_{sp,C}^{photo} - Q_{Fe:C}^{sp}(J_{sp}^{grz} + J_{sp}^{l} + J_{sp}^{agg}) \quad (A15) $$

$$ \frac{d}{dt}(Ca_{sp}) = J_{sp,caco3}^{prod} - Q_{caco3:C}^{sp}(J_{sp}^{grz} + J_{sp}^{l} + J_{sp}^{agg}) \quad (A16) $$

$$ \frac{d}{dt}(C_{diat}) = J_{diat,C}^{photo} - (J_{diat}^{grz} + J_{diat}^{l} + J_{diat}^{agg}) \quad (A17) $$

$$ \frac{d}{dt}(Fe_{diat}) = gQ_{Fe:C}^{diat}J_{diat,C}^{photo} - Q_{Fe:C}^{diat}(J_{diat}^{grz} + J_{diat}^{l} + J_{diat}^{agg}) \quad (A18) $$

$$ \frac{d}{dt}(Chl_{diat}) = J_{diat,chl}^{pa} - Q_{Chl:C}^{diat}(J_{diat}^{grz} + J_{diat}^{l} + J_{diat}^{agg}) \quad (A19) $$

$$ \frac{d}{dt}(Si_{diat}) = gQ_{Si:C}J_{diat,C}^{photo} - Q_{Si:C}^{diat}(J_{diat}^{grz} + J_{diat}^{l} + J_{diat}^{agg}) \quad (A20) $$
$$ \frac{d}{dt}(C_{diaz}) = J_{diaz,C}^{photo} - (J_{diaz}^{grz} + J_{diaz}^{l} + J_{diaz}^{agg}) \quad (A21) $$

$$ \frac{d}{dt}(Fe_{diaz}) = gQ_{Fe:C}^{diaz}J_{diaz,C}^{photo} - Q_{Fe:C}^{diaz}(J_{diaz}^{grz} + J_{diaz}^{l} + J_{diaz}^{agg}) \quad (A22) $$

$$ \frac{d}{dt}(Chl_{diaz}) = J_{diaz,chl}^{pa} - Q_{Chl:C}^{diaz}(J_{diaz}^{grz} + J_{diaz}^{l} + J_{diaz}^{agg}) \quad (A23) $$

$$ \frac{d}{dt}(C_{zoo}) = (J_{sp,zoo}^{grz} + J_{diat,zoo}^{grz} + J_{diaz,zoo}^{grz}) - J_{zoo}^{l} \quad (A24) $$

$$ \frac{d}{dt}(C_{doc}) = J_{sp,doc}^{l} + J_{diat,doc}^{l} + J_{diaz,doc}^{l} + J_{zoo,doc}^{l} + (J_{sp,doc}^{grz} + J_{diat,doc}^{grz} + J_{diaz,doc}^{grz}) - J_{doc,dic} \quad (A25) $$

$$ \frac{d}{dt}(N_{don}) = Q_{N:C} (J_{sp,doc}^{l} + J_{diat,doc}^{l} + J_{diaz,doc}^{l} + J_{zoo,doc}^{l} + J_{sp,doc}^{grz} + J_{diat,doc}^{grz} + J_{diaz,doc}^{grz}) - J_{don,nh4} + J_{diaz,N}^{excrete} \quad (A26) $$

$$ \frac{d}{dt}(P_{dop}) = Q_{P:C} (J_{sp,doc}^{l} + J_{diat,doc}^{l} + J_{sp,doc}^{grz} + J_{diat,doc}^{grz} + J_{zoo,doc}^{l}) + J_{diaz,dop}^{l} - J_{dop,po4} \quad (A27) $$

$$ \frac{d}{dt}(Fe_{dofe}) = Q_{Fe:C}^{sp}(J_{sp,doc}^{l} + J_{sp,doc}^{grz}) + Q_{Fe:C}^{diat}(J_{diat,doc}^{l} + J_{diat,doc}^{grz}) + Q_{Fe:C}^{diaz}(J_{diaz,doc}^{l} + J_{diaz,doc}^{grz}) + Q_{Fe:C}^{zoo}J_{zoo,doc}^{l} - J_{dofe,fe} \quad (A28) $$
### 5.3.2 Treatment of particulate organic matter

#### General model

Particulate organic matter is produced and instantaneously distributed over the depth of the water column following the exponential solution to the steady-state 1-dimensional production-remineralization equation:

$$ \frac{\partial \Phi(z)}{\partial z} = J^{\text{remin}}(z) - J^{\text{prod}}(z) = -\frac{\Phi(z)}{\lambda} + J^{\text{prod}}(z), \qquad (\text{A29}) $$

where $\Phi$ is a flux, $\lambda$ is the remineralization length-scale, $J^{\text{prod}}$ is the known production rate within the layer, and $J^{\text{remin}}$ the remineralization rate within the layer, which needs to be determined. For a single layer, assuming the flux at the top of the layer $\Phi(k)$ is known, and the production $J^{\text{prod}}(k)$ is constant within the layer, the solution to equation A29, can be cast to determine the flux out of the layer, $\Phi(k - 1)$, for each element $i$:

$$ \Phi_i(k-1) = \Phi_i(k) e^{-\frac{\Delta z}{\lambda_i}} + J^{\text{prod}}(k) (1 - e^{-\frac{\Delta z}{\lambda_i}})\lambda, \qquad (\text{A30}) $$

Particulate organic carbon (POC) is partitioned between a free and mineral associated component:

$$ \Phi_{\text{poc}}(z) = \Phi_{\text{poc}}^{\text{free}}(z) + \Phi_{\text{poc}}^{\text{mineral}}(z) \qquad (\text{A31}) $$

$\Phi_{\text{poc}}^{\text{mineral}}$ can be associated with CaCO$_3$, SiO$_2$ or dust. Each mineral-associated POC flux is further partitioned into a “soluble” component, which remineralizes with the length-scale of the associated mineral, and a “non-soluble” component which remineralizes with a length-scale of 40,000 m.

$$ \Phi_{\text{poc}}^{\text{mineral}}(z) = \rho_{\text{pic}} (\Phi_{\text{pic}}^{\text{sol}}(z) + \Phi_{\text{pic}}^{\text{nonsol}}(z)) \\ + \rho_{\text{psi}} (\Phi_{\text{psi}}^{\text{sol}}(z) + \Phi_{\text{psi}}^{\text{nonsol}}(z)) \\ + \rho_{\text{dust}} (\Phi_{\text{dust}}^{\text{sol}}(z) + \Phi_{\text{dust}}^{\text{nonsol}}(z)) \qquad (\text{A32}) $$

For all components of the fluxes except Fe, the flux out of the layer, $\Phi(k-1)$, is computed first, with knowledge of the source within the layer, $J^{\text{prod}}(k)$, and of the remineralization length-scale. Remineralization in the layer is then calculated from conservation, i.e., from A29. Below, for each component, we list the equations used to determine the production terms $J^{\text{prod}}(k)$, followed by the fluxes out, $\Phi(k-1)$, and finally the remineralization terms, $J^{\text{remin}}(z)$, which enter the tracer conservation equations.

#### Production

$$ J_{\text{poc}}^{\text{prod}} = (J_{\text{sp,poc}}^{\text{grz}} + J_{\text{diat,poc}}^{\text{grz}} + J_{\text{diaz,poc}}^{\text{grz}}) + (J_{\text{sp}}^{\text{agg}} + J_{\text{diat}}^{\text{agg}}) \\ + (J_{\text{sp,poc}}^{l} + J_{\text{diat,poc}}^{l} + J_{\text{diaz,poc}}^{l}) + J_{\text{zoo,poc}}^{l} \qquad (\text{A33}) $$
$$J_{pic}^{prod} = (0.67 J_{sp}^{grz} + J_{sp}^{agg} + J_{sp}^{l}) Q_{caco3:C}^{sp} \quad (A34)$$

$$J_{psi}^{prod} = (0.5 J_{diat}^{grz} + 0.05 J_{diat}^{l} + J_{diat}^{agg}) Q_{Si:C}^{diat} \quad (A35)$$

$$J_{dust}^{prod} = 0 \quad (A36)$$

Available production for free POC is then:

$$J_{poc,avail}^{prod} = J_{poc}^{prod} - \rho_{pic}J_{pic}^{prod} - \rho_{psi}J_{psi}^{prod} \quad (A37)$$

## Fluxes out

Temperature dependency is used to modify the remineralization length scales of particulate organic carbon, POC, and opal, SiO₂:

$$T_{poc} = (Q_{10,poc})^{\frac{T-30}{10}} \quad (A38)$$

$$T_{Psi} = (Q_{10,psi})^{\frac{T-30}{10}} \quad (A39)$$

Free POC flux equation:

$$\Phi_{poc}^{free}(k-1) = \Phi(k)_{poc}^{free} e^{-\frac{\Delta z T_{poc}}{\lambda_{poc}}} + J_{poc,avail}^{prod}(k) (1 - e^{-\frac{\Delta z T_{poc}}{\lambda_{poc}}}) \frac{\lambda_{poc}}{T_{poc}} \quad (A40)$$

Soluble mineral-associated POC flux equation:

$$\Phi_{pic}^{sol}(k-1) = \Phi(k)_{pic}^{sol} e^{-\frac{\Delta z}{\lambda_{pic}}} + J_{pic}^{prod}(k) (1 - \gamma_{pic}) (1 - e^{-\frac{\Delta z}{\lambda_{pic}}}) \lambda_{pic} \quad (A41)$$

$$\Phi_{psi}^{sol}(k-1) = \Phi(k)_{psi}^{sol} e^{-\frac{\Delta z T_{psi}}{\lambda_{psi}}} + J_{psi}^{prod}(k) (1 - \gamma_{psi}) (1 - e^{-\frac{\Delta z T_{psi}}{\lambda_{psi}}}) \frac{\lambda_{psi}}{T_{psi}} \quad (A42)$$

$$\Phi_{dust}^{sol}(k-1) = \Phi(k)_{dust}^{sol} e^{-\frac{\Delta z}{\lambda_{dust}}} \quad (A43)$$

Non-soluble mineral-associated POC flux equation:

$$\Phi_{pic}^{nonsol}(k-1) = \Phi(k)_{pic}^{nonsol} e^{-\frac{\Delta z}{\lambda_{nonsol}}} + J_{pic}^{prod}(k) \gamma_{pic} \Delta z \quad (A44)$$

$$\Phi_{psi}^{nonsol}(k-1) = \Phi(k)_{psi}^{nonsol} e^{-\frac{\Delta z}{\lambda_{nonsol}}} + J_{psi}^{prod}(k) \gamma_{psi} \Delta z \quad (A45)$$

$$\Phi_{dust}^{nonsol}(k-1) = \Phi(k)_{dust}^{nonsol} e^{-\frac{\Delta z}{\lambda_{nonsol}}} \quad (A46)$$
806 **Remineralization**

807 Remineralization is computed from conservation, i.e, A29:

$$
J^{\text{remin}} = J^{\text{prod}} + \frac{\partial \Phi(z)}{\partial z} \qquad (\text{A47})
$$

808 Numerically, for each individual layer, we have:

$$
J_{poc}^{\text{remin}}(k) = J_{poc}^{\text{prod}}(k) + \frac{\Phi_{poc}^{\text{free}}(k) - \Phi_{poc}^{\text{free}}(k-1) + \Phi_{poc}^{\text{mineral}}(k) - \Phi_{poc}^{\text{mineral}}(k-1)}{\Delta z} \quad (\text{A48})
$$

809

$$
J_{pic}^{\text{remin}}(k) = J_{pic}^{\text{prod}}(k) + \frac{\Phi_{pic}^{\text{sol}}(k) - \Phi_{pic}^{\text{sol}}(k-1) + \Phi_{pic}^{\text{nonsol}}(k) - \Phi_{pic}^{\text{nonsol}}(k-1)}{\Delta z} \quad (\text{A49})
$$

810

$$
J_{\text{psi}}^{\text{remin}}(k) = J_{\text{psi}}^{\text{prod}}(k) + \frac{\Phi_{\text{psi}}^{\text{sol}}(k) - \Phi_{\text{psi}}^{\text{sol}}(k-1) + \Phi_{\text{psi}}^{\text{nonsol}}(k) - \Phi_{\text{psi}}^{\text{nonsol}}(k-1)}{\Delta z} \quad (\text{A50})
$$

811

$$
J_{\text{dust}}^{\text{remin}}(k) = \frac{\Phi_{\text{dust}}^{\text{sol}}(k) - \Phi_{\text{dust}}^{\text{sol}}(k-1) + \Phi_{\text{dust}}^{\text{nonsol}}(k) - \Phi_{\text{dust}}^{\text{nonsol}}(k-1)}{\Delta z} \quad (\text{A51})
$$

812 Particulate Fe

813 Production is as follows:

$$
\begin{align}
J_{pfe}^{\text{prod}} ={}& (J_{sp}^{\text{agg}} + J_{sp,poc}^{\text{grz}} + J_{sp,poc}^{l}) Q_{Fe:C}^{sp} \nonumber \\
& + (J_{diat}^{\text{agg}} + J_{diat,poc}^{\text{grz}} + J_{diat,poc}^{l}) Q_{Fe:C}^{diat} \nonumber \\
& + (J_{diaz}^{\text{agg}} + J_{diaz,poc}^{\text{grz}}) Q_{Fe:C}^{diaz} + 0.1 J_{Fe}^{\text{scav}} \nonumber \\
& + J_{zoo}^{l} f_{zoo}^{d} Q_{Fe:C}^{zoo} \tag{A52}
\end{align}
$$

814

815

Particulate Fe remineralization is assumed to be proportional to POC remineralization plus a
release from dust :

$$
J_{pfe}^{\text{remin}}(k) = J_{poc}^{\text{remin}}(k) \frac{\Phi_{pfe}^{\text{sol}}(k) + \Phi_{pfe}^{\text{nonsol}}(k)}{\Phi_{poc}^{\text{free}}(k) + \Phi_{poc}^{\text{mineral}}(k)} + 626.712 J_{\text{dust}}^{\text{remin}}(k) \quad (\text{A53})
$$

816

The flux out can then be computed from conservation:

$$
\Phi_{pfe}^{\text{sol}}(k-1) = \Phi_{pfe}^{\text{sol}}(k) + \Delta z \left( J_{pfe}^{\text{prod}}(k) - J_{poc}^{\text{remin}}(k) \frac{\Phi_{pfe}^{\text{sol}}(k) + \Phi_{pfe}^{\text{nonsol}}(k)}{\Phi_{poc}^{\text{free}}(k) + \Phi_{poc}^{\text{mineral}}(k)} \right) \quad (\text{A54})
$$

$$
\Phi_{pfe}^{\text{nonsol}}(k-1) = \Phi_{pfe}^{\text{nonsol}}(k) \qquad (\text{A55})
$$

817 Bottom flux to sediment

$$
\Phi_{\text{poc}}^{\text{out}} = \Phi_{\text{poc}}(k_{\text{bottom}}) \qquad (\text{A56})
$$
### 818 5.3.3 Biogeochemical rates

#### 819 Carbon

$$T_{func} = 2^{0.1*T-3} \qquad (A57)$$

$$J_{sp,C}^{photo} = PC_{ref}^{sp} f_{nut}^{sp} T_{func} \left( 1 - e^{-\frac{\alpha_{chl} Q_{chl:C}^{sp} PAR}{PC_{ref}^{sp} f_{nut}^{sp} T_{func}}} \right) C_{sp} \qquad (A58)$$

$$J_{diat,C}^{photo} = PC_{ref}^{diat} f_{nut}^{diat} T_{func} \left( 1 - e^{-\frac{\alpha_{chl} Q_{chl:C}^{diat} PAR}{PC_{ref}^{diat} f_{nut}^{diat} T_{func}}} \right) C_{diat} \qquad (A59)$$

$$J_{diaz,C}^{photo} = PC_{ref}^{diaz} f_{nut}^{diaz} T_{func} \left( 1 - e^{-\frac{\alpha_{chl} Q_{chl:C}^{diaz} PAR}{PC_{ref}^{diaz} f_{nut}^{diaz} T_{func}}} \right) C_{diaz} \qquad (A60)$$

#### 823 Remineralization

$$J_{don,nh4} = \tau_{dom} N_{don} \qquad (A61)$$

$$J_{doc,dic} = \tau_{dom} C_{doc} \qquad (A62)$$

$$J_{dop,po4} = \tau_{dom} P_{dop} \qquad (A63)$$

$$J_{dofe,fe} = \tau_{dom} F_{e} d_{ofe} \qquad (A64)$$

$$\Phi_{poc}^{remin,sed} = \begin{cases} \tau_{remin,sed} POC_{sed}, & \text{if } k = k_{bottom} \\ 0, & \text{elsewhere} \end{cases} \qquad (A65)$$

$$\Phi_{pic}^{remin,sed} = \begin{cases} \tau_{remin,sed} CaCO_3_{sed}, & \text{if } k = k_{bottom} \\ 0, & \text{elsewhere} \end{cases} \qquad (A66)$$

$$\Phi_{Si}^{remin,sed} = \begin{cases} \tau_{remin,sed} Si_{sed}, & \text{if } k = k_{bottom} \\ 0, & \text{elsewhere} \end{cases} \qquad (A67)$$

#### 830 Nutrient limitation

$$f_{nut}^{sp} = \min(V_{NO_3}^{sp} + V_{NH_4}^{sp}, V_{Fe}^{sp}, V_{PO_4}^{sp}) \qquad (A68)$$

$$f_{nut}^{diat} = \min(V_{NO_3}^{diat} + V_{NH_4}^{diat}, V_{Fe}^{diat}, V_{SiO_3}^{diat}, V_{PO_4}^{diat}) \qquad (A69)$$

$$f_{nut}^{diaz} = \min(V_{Fe}^{diaz}, V_{PO_4}^{diaz}) \qquad (A70)$$
833

$$V_{NO_3}^{sp} = \frac{N_{NO_3}/k_{no3}^{sp}}{1 + N_{NO_3}/k_{no3}^{sp} + N_{NH_4}/k_{nh4}^{sp}} \quad (A71)$$

834

$$V_{NO_3}^{diat} = \frac{N_{NO_3}/k_{no3}^{diat}}{1 + N_{NO_3}/k_{no3}^{diat} + N_{NH_4}/k_{nh4}^{diat}} \quad (A72)$$

835

$$V_{NH_4}^{sp} = \frac{N_{NH_4}/k_{nh4}^{sp}}{1 + N_{NO_3}/k_{no3}^{sp} + N_{NH_4}/k_{nh4}^{sp}} \quad (A73)$$

836

$$V_{NH_4}^{diat} = \frac{N_{NH_4}/k_{nh4}^{diat}}{1 + N_{NO_3}/k_{no3}^{diat} + N_{NH_4}/k_{nh4}^{diat}} \quad (A74)$$

837

$$V_{Fe}^{sp} = \frac{Fe}{Fe + k_{fe}^{sp}} \quad (A75)$$

838

$$V_{Fe}^{diat} = \frac{Fe}{Fe + k_{fe}^{diat}} \quad (A76)$$

839

$$V_{Fe}^{diaz} = \frac{Fe}{Fe + k_{fe}^{diaz}} \quad (A77)$$

840

$$V_{PO_4}^{sp} = \frac{PO_4}{PO_4 + k_{po4}^{sp}} \quad (A78)$$

841

$$V_{PO_4}^{diat} = \frac{PO_4}{PO_4 + k_{po4}^{diat}} \quad (A79)$$

842

$$V_{PO_4}^{diaz} = \frac{PO_4}{PO_4 + k_{po4}^{diaz}} \quad (A80)$$

843

$$V_{SiO_2}^{diat} = \frac{Si_{sio3}}{Si_{sio2} + k_{sio2}^{diat}} \quad (A81)$$

844

NO₃ and NH₄ uptake

$$J_{sp,no3}^{u} = Q_{N:C} \frac{V_{NO_3}^{sp}}{V_{NO_3}^{sp} + V_{NH_4}^{sp}} J_{sp,C}^{photo} \quad (A82)$$

845

$$J_{sp,nh4}^{u} = Q_{N:C} \frac{V_{NH_4}^{sp}}{V_{NO_3}^{sp} + V_{NH_4}^{sp}} J_{sp,C}^{photo} \quad (A83)$$

846

$$J_{diat,no3}^{u} = Q_{N:C} \frac{V_{NO_3}^{diat}}{V_{NO_3}^{diat} + V_{NH_4}^{diat}} J_{diat,C}^{photo} \quad (A84)$$

847

$$J_{diat,nh4}^{u} = Q_{N:C} \frac{V_{NH_4}^{diat}}{V_{NO_3}^{diat} + V_{NH_4}^{diat}} J_{diat,C}^{photo} \quad (A85)$$
848 N₂ fixation

$$J_{\text{diaz,N}}^{\text{photo}} = Q_{N:C} J_{\text{diaz,C}}^{\text{photo}} \qquad (\text{A86})$$

849

$$J_{\text{diaz}}^{\text{nfix}} = J_{\text{diaz,N}}^{\text{photo}}/(1 - \alpha_{\text{ex}}^{\text{diaz}}) \qquad (\text{A87})$$

850

$$J_{\text{diaz,N}}^{\text{excrete}} = \alpha_{\text{ex}}^{\text{diaz}} J_{\text{diaz}}^{\text{nfix}} \qquad (\text{A88})$$

851 PO₄ uptake

$$J_{\text{sp,po4}}^{u} = Q_{P:C} J_{\text{sp,C}}^{\text{photo}} \qquad (\text{A89})$$

852

$$J_{\text{diat,po4}}^{u} = Q_{P:C} J_{\text{diat,C}}^{\text{photo}} \qquad (\text{A90})$$

853

$$J_{\text{diaz,po4}}^{u} = Q_{P:C} J_{\text{diaz,C}}^{\text{photo}} \qquad (\text{A91})$$

854 Chl stoichiometry and production

$$Q_{\text{Chl:C}}^{\text{sp}} = \text{Chl}_{\text{sp}}/C_{\text{sp}} \qquad (\text{A92})$$

855

$$Q_{\text{Chl:C}}^{\text{diat}} = \text{Chl}_{\text{diat}}/C_{\text{diat}} \qquad (\text{A93})$$

856

$$Q_{\text{Chl:C}}^{\text{diaz}} = \text{Chl}_{\text{diaz}}/C_{\text{diaz}} \qquad (\text{A94})$$

$$J_{\text{sp,chl}}^{\text{pa}} = Q_{\text{N:chl,max}}^{\text{sp,growth}} \frac{J_{\text{sp,C}}^{\text{photo}}/C_{\text{sp}}}{\alpha_{\text{chl}} Q_{\text{Chl:C}}^{\text{sp}} \text{PAR}} Q_{\text{N:C}} J_{\text{sp,C}}^{\text{photo}} \qquad (\text{A95})$$

$$J_{\text{diat,chl}}^{\text{pa}} = Q_{\text{N:chl,max}}^{\text{diat,growth}} \frac{J_{\text{diat,C}}^{\text{photo}}/C_{\text{diat}}}{\alpha_{\text{chl}} Q_{\text{Chl:C}}^{\text{diat}} \text{PAR}} Q_{\text{N:c}} J_{\text{diat,C}}^{\text{photo}} \qquad (\text{A96})$$

$$J_{\text{diaz,chl}}^{\text{pa}} = Q_{\text{N:chl,max}}^{\text{diaz,growth}} \frac{J_{\text{diaz,C}}^{\text{photo}}/C_{\text{diaz}}}{\alpha_{\text{chl}} Q_{\text{Chl:C}}^{\text{diaz}} \text{PAR}} Q_{\text{N:c}} J_{\text{diaz,C}}^{\text{photo}} \qquad (\text{A97})$$

857 Light

$$\text{PAR}_{\text{B}}(z) = \max(0, f_{\text{Isw}} I_{\text{SW}}) e^{\int_{z}^{0} (\mu_{\text{chl}}(\text{Chl}_{\text{sp}} + \text{Chl}_{\text{diat}} + \text{Chl}_{\text{diaz}}) + \mu_{w}) dz} \qquad (\text{A98})$$

858

$$\text{PAR} = \frac{1}{\mu_{\text{chl}}(\text{Chl}_{\text{sp}} + \text{Chl}_{\text{diat}} + \text{Chl}_{\text{diaz}}) + \mu_{w})} \frac{\partial \text{PAR}_{\text{B}}(z)}{\partial z} \qquad (\text{A99})$$
859 Fe stoichiometry and cycling

$$Q_{Fe:C}^{sp} = Fe_{sp}/C_{sp} \qquad (A100)$$

860

$$Q_{Fe:C}^{diat} = Fe_{diat}/C_{diat} \qquad (A101)$$

861

$$Q_{Fe:C}^{diaz} = Fe_{diaz}/C_{diaz} \qquad (A102)$$

862

$$Q_{Fe:C}^{Zoo} = 2.5e^{-6} \qquad (A103)$$

863

$$gQ_{Fe:C}^{sp} = \begin{cases} Q_{Fe:C}^{sp,g0}, & \text{if } Fe \ge 2 k_{fe}^{sp} \\ \max(Q_{Fe:C}^{sp,g1}, \frac{Q_{Fe:C}^{sp,g0} Fe}{2 k_{fe}^{sp}}), & \text{otherwise} \end{cases} \qquad (A104)$$

864

$$gQ_{Fe:C}^{diat} = \begin{cases} Q_{Fe:C}^{diat,g0}, & \text{if } Fe \ge 2 k_{fe}^{diat} \\ \max(Q_{Fe:C}^{diat,g1}, \frac{Q_{Fe:C}^{sp,g0} Fe}{2 k_{fe}^{diat}}), & \text{otherwise} \end{cases} \qquad (A105)$$

865

$$gQ_{Fe:C}^{diaz} = \begin{cases} Q_{Fe:C}^{diaz,g0}, & \text{if } Fe \ge 2 k_{fe}^{diaz} \\ \max(Q_{Fe:C}^{diaz,g1}, \frac{Q_{Fe:C}^{diaz,g0} Fe}{2 k_{fe}^{diaz}}), & \text{otherwise} \end{cases} \qquad (A106)$$

$$J_{sp,Fe}^{u} = gQ_{Fe:C}^{sp} J_{sp,C}^{photo} \qquad (A107)$$

866

$$J_{diat,Fe}^{u} = gQ_{Fe:C}^{diat} J_{diat,C}^{photo} \qquad (A108)$$

867

$$J_{diaz,Fe}^{u} = gQ_{Fe:C}^{diaz} J_{diaz,C}^{photo} \qquad (A109)$$

$$J_{Fe}^{scav} = \begin{cases} (J_{Fe}^{scav,0} \min((\Phi_{poc} + 8.33 \cdot 10^4 \Phi_{dust})/\Phi_{poc}^{ref}, \beta_{fe}^{max,scale1}) + (Fe - \beta_{fe}^{thresh,1}) 6.0/1.4e^{-3}) 3.1709792e^{-8} Fe, & \text{if } Fe > \beta_{fe}^{thresh,1} \\ (J_{Fe}^{scav,0} \min((\Phi_{poc} + 8.33 \cdot 10^4 \Phi_{dust})/\Phi_{poc}^{ref}, \beta_{fe}^{max,scale1}) + \frac{Fe}{\beta_{fe}^{thresh,2}}) 3.1709792e^{-8} Fe, & \text{if } Fe < \beta_{fe}^{thresh,2} \end{cases} \qquad (A110)$$

$$J_{Fe}^{Sed} = \begin{cases} 0.75 \cdot 10^{(2.5-0.0165 O_2)} (0.001/86400)/\Delta z, & \text{if } k = k_{bottom} \text{ and } O_2 > O_2^{min} \\ 0.25 \cdot 10^{(2.5-0.0165 O_2)} (0.001/86400)/\Delta z, & \text{if } k = k_{bottom-1} \text{ and } O_2 > O_2^{min} \\ 10^{(2.5-0.0165 O_2)} (0.001/86400)/\Delta z & \text{elsewhere} \end{cases} \qquad (A111)$$

868 75% of $J_{Fe}^{Sed}$ is released in the bottom layer and 25% in the layer above.
869 CacO₃ stoichiometry production

$$Q_{caco3:C}^{sp} = \begin{cases} CaCO_3_{sp}/C_{sp} & \text{if } CaCO_3_{sp}/C_{sp} \le Q_{caco3:C}^{sp,max} \\ Q_{caco3:C}^{sp,max} & \text{otherwise} \end{cases} \qquad (A112)$$

$$J_{sp,caco3}^{prod,max} = f_{caco3}^{prod} J_{sp,C}^{photo} f_{nut}^{sp} \qquad (A113)$$

$$J_{sp,caco3}^{prodT} = \begin{cases} J_{sp,caco3}^{prodmax} & \text{if } T \ge T_1^{caco3} \\ J_{sp,caco3}^{prodmax} \frac{\max[T-T_2^{caco3}, 0]}{T_1^{caco3}-T_2^{caco3}} & \text{if } T < T_1^{caco3} \end{cases} \qquad (A114)$$

$$J_{sp,caco3}^{prod} = \begin{cases} \min[(J_{sp,caco3}^{prodT}C_{sp}/3, 0.4J_{sp,C}^{photo}]) & \text{where } C_{sp} > 3.0 \\ J_{sp,caco3}^{prod,T} & \text{elsewhere} \end{cases} \qquad (A115)$$

870 Si stoichiometry

$$Q_{Si:C}^{diat} = \max(Si_{diat}/C_{diat}, Q_{Si:C}^{diat,max}) \qquad (A116)$$

$$gQ_{Si:C}^{diat} = \begin{cases} \min(Q_{Si:C}^{diat,g0}(2Q_{Si:C}^{diat,g1}k_{fe}^{diat}/Fe - Q_{Si:C}^{diat,g1} + 1), Q_{Si:C}^{diat,max}), & \text{if } 0 < Fe < 2k_{fe}^{diat} \text{ and } Si_{sio2} > 2k_{sio2}^{diat} \\ Q_{Si:C}^{diat,max} & \text{if } Fe = 0 \\ Q_{Si:C}^{diat,g0} & \text{otherwise} \end{cases} \qquad (A117)$$

$$J_{diat,Si}^{u} = gQ_{Si:C}^{diat} J_{diat,C}^{photo} \qquad (A118)$$

871 Grazing

$$\beta_{thres}^{grz,sp} = \begin{cases} \beta_{thres,0}^{grz,sp} & z \ge -100 \text{ m} \\ \beta_{thres,0}^{grz,sp} \cdot (200 + z)/100 & -100 \text{ m} > z > -200 \text{ m} \\ 0 & z \le -200 \text{ m} \end{cases} \qquad (A119)$$

$$P'_{sp} = \max(C_{sp} - \beta_{thres}^{grz,sp}, 0) \qquad (A120)$$

872

$$\beta_{thres}^{grz,diat} = \begin{cases} \beta_{thres,0}^{grz,diat} & z \ge -100 \text{ m} \\ \beta_{thres,0}^{grz,diat} \cdot (200 + z)/100 & -100 \text{ m} > z > -200 \text{ m} \\ 0 & z \le -200 \text{ m} \end{cases} \qquad (A121)$$

$$P'_{diat} = \max(C_{diat} - \beta_{thres}^{grz,diat}, 0) \qquad (A122)$$

$$\beta_{thres}^{grz,diaz} = \begin{cases} \beta_{thres,0}^{grz,diaz} & z \ge -100 \text{ m} \\ \beta_{thres,0}^{grz,diaz} \cdot (200 + z)/100 & -100 \text{ m} > z > -200 \text{ m} \\ 0 & z \le -200 \text{ m} \end{cases} \qquad (A123)$$
$$P'_{\text{diaz}} = \max(C_{\text{diaz}} - \beta_{\text{thres}}^{\text{grz, diaz}}, 0) \qquad (\text{A124})$$

$$J_{sp}^{\text{grz}} = J_{sp}^{g,max} T_{\text{func}} \frac{P'_{sp}}{P'_{sp} + (\beta_{z}^{\text{grz}})^2} C_{\text{zoo}} \qquad (\text{A125})$$

$$J_{\text{diat}}^{\text{grz}} = J_{\text{diat}}^{g,max} T_{\text{func}} \frac{P'_{\text{diat}}}{P'_{\text{diat}} + 0.81(\beta_{z}^{\text{grz}})^2} C_{\text{zoo}} \qquad (\text{A126})$$

$$J_{\text{diaz}}^{\text{grz}} = J_{\text{diaz}}^{g,max} T_{\text{func}} \frac{P'_{\text{diaz}}}{P'_{\text{diaz}} + (\beta_{z}^{\text{grz}})^2} C_{\text{zoo}} \qquad (\text{A127})$$

$$J_{sp}^{g,zoo} = \alpha_{sp}^{\text{grz,zoo}} J_{sp}^{\text{grz}} \qquad (\text{A128})$$

$$J_{\text{diat}}^{g,zoo} = \alpha_{\text{diat}}^{\text{grz,zoo}} J_{\text{diat}}^{\text{grz}} \qquad (\text{A129})$$

$$J_{\text{diaz}}^{g,zoo} = \alpha_{\text{diaz}}^{\text{grz,zoo}} J_{\text{diaz}}^{\text{grz}} \qquad (\text{A130})$$

$$J_{sp}^{g,poc} = \max(Q_{caco3:C}^{sp,max}, Q_{caco3:C}^{sp}, \min(0.18 P'_{sp}, \alpha_{sp}^{\text{grz,poc}})) J_{sp}^{\text{grz}} \qquad (\text{A131})$$

$$J_{\text{diat}}^{g,poc} = \alpha_{\text{diat}}^{\text{grz,poc}} J_{\text{diat}}^{\text{grz}} \qquad (\text{A132})$$

$$J_{\text{diaz}}^{g,poc} = \alpha_{\text{diaz}}^{\text{grz,poc}} J_{\text{diaz}}^{\text{grz}} \qquad (\text{A133})$$

$$J_{sp}^{g,doc} = \alpha_{sp}^{\text{grz,doc}} J_{sp}^{\text{grz}} - J_{sp}^{\text{grz,poc}} \qquad (\text{A134})$$

$$J_{\text{diat}}^{g,doc} = \alpha_{\text{diat}}^{\text{grz,doc}} J_{\text{diat}}^{\text{grz}} \qquad (\text{A135})$$

$$J_{\text{diaz}}^{g,doc} = \alpha_{\text{diaz}}^{\text{grz,doc}} J_{\text{diaz}}^{\text{grz}} \qquad (\text{A136})$$

$$J_{sp}^{g,dic} = \alpha_{sp}^{\text{grz,dic}} J_{sp}^{\text{grz}} \qquad (\text{A137})$$

$$J_{\text{diat}}^{g,dic} = \alpha_{\text{diat}}^{\text{grz,dic}} J_{\text{diat}}^{\text{grz}} \qquad (\text{A138})$$

$$J_{\text{diaz}}^{g,dic} = \alpha_{\text{diaz}}^{\text{grz,dic}} J_{\text{diaz}}^{\text{grz}} \qquad (\text{A139})$$

873

874

875

876

877

878

879

880
$$J_{sp}^{l} = \lambda_{sp}^{mort} P'_{sp} \qquad (A140)$$

882

$$J_{diat}^{l} = \lambda_{diat}^{mort} P'_{diat} \qquad (A141)$$

883

$$J_{diaz}^{l} = \lambda_{diaz}^{mort} P'_{diaz} \qquad (A142)$$

$$J_{sp}^{l,poc} = Q_{caco3:C}^{sp} J_{sp}^{l} \qquad (A143)$$

884

$$J_{diat}^{l,poc} = \alpha_{diat}^{l,poc} J_{diat}^{l} \qquad (A144)$$

885

$$J_{diaz}^{l,poc} = \alpha_{diaz}^{l,poc} J_{diaz}^{l} \qquad (A145)$$

$$J_{sp}^{l,doc} = (1 - f_{labile})(J_{sp}^{l} - J_{sp}^{l,poc}) \qquad (A146)$$

886

$$J_{diat}^{l,doc} = (1 - f_{labile})(J_{diat}^{l} - J_{diat}^{l,poc}) \qquad (A147)$$

887

$$J_{diaz}^{l,doc} = (1 - f_{labile})(J_{diaz}^{l} - J_{diaz}^{l,poc}) \qquad (A148)$$

$$J_{sp}^{l,dic} = f_{labile}(J_{sp}^{l} - J_{sp}^{l,poc}) \qquad (A149)$$

888

$$J_{diat}^{l,dic} = f_{labile}(J_{diat}^{l} - J_{diat}^{l,poc}) \qquad (A150)$$

889

$$J_{diaz}^{l,dic} = f_{labile}(J_{diaz}^{l} - J_{diaz}^{l,poc}) \qquad (A151)$$

$$P_{diaz}^{remain} = Q_{P:C}^{diaz} (J_{diaz}^{grz} + J_{diaz}^{l} + J_{diaz}^{agg}) - Q_{P:C} (J_{diaz,poc}^{grz} + J_{diaz,zoo}^{grz}) \qquad (A152)$$

890

$$J_{diaz,dop}^{l} = (1 - f_{labile}) P_{diaz}^{remain} \qquad (A153)$$

891

$$J_{diaz,dip}^{l} = f_{labile} P_{diaz}^{remain} \qquad (A154)$$

$$\beta_{zoo}^{thres,l} = \begin{cases} \beta_{zoo}^{thres0,l} & z \geq -100m \\ \beta_{zoo}^{thres0,l} \cdot (200 + z)/100 & -100 m > z > -200 m \\ 0 & z \leq -200 m \end{cases} \qquad (A155)$$

$$Z'_{zoo} = \max(C_{zoo} - \beta_{zoo}^{thres,l}, 0) \qquad (A156)$$
$$J_{zoo}^{l} = \lambda_{zoo}^{\text{mort2}} T_{\text{func}} Z_{zoo}^{\prime 2} + \lambda_{zoo}^{\text{mort}} T_{\text{func}} Z_{zoo}^{\prime} \quad (\text{A157})$$

$$f_{zoo}^{d} = \frac{0.1333 J_{diat}^{grz} + 0.0333 J_{sp}^{grz}}{J_{diat}^{grz} + J_{sp}^{grz} + J_{diaz}^{grz}} \quad (\text{A158})$$

$$J_{zoo,doc}^{l} = (1 - f_{\text{labile}}) (1 - f_{zoo}^{d}) J_{zoo}^{l} \quad (\text{A159})$$

$$J_{zoo,dic}^{l} = f_{\text{labile}} (1 - f_{zoo}^{d}) J_{zoo}^{l} \quad (\text{A160})$$

$$J_{zoo,poc}^{l} = f_{zoo}^{d} J_{zoo}^{l} \quad (\text{A161})$$

892

893

## 894 Aggregation

$$J_{sp}^{agg} = \min(\tau_{sp}^{agg,max} P'_{sp}, \lambda_{sp}^{\text{mort2}} P'_{sp}^{\prime 2}) \quad (\text{A162})$$

895

$$J_{diat}^{agg} = \max(\tau_{diat}^{agg,min} P'_{diat}, \min(\tau_{diat}^{agg,max} P'_{diat}, \lambda_{diat}^{\text{mort2}} P'_{diat}^{\prime 2})) \quad (\text{A163})$$

896

$$J_{diaz}^{agg} = 0.0 \quad (\text{A164})$$

## 897 N-cycle rates

$$J^{\text{ammox}} = \begin{cases} \tau^{\text{ammox}} N_{nh4}, & \text{where } \text{PAR}_B(k) < \text{PAR}_{\text{max}}^{\text{nitrif}} \\ 0, & \text{elsewhere} \end{cases} \quad (\text{A165})$$

$$J^{\text{nit}} = \begin{cases} \tau^{\text{nit}} N_{no2}, & \text{where } \text{PAR}_B(k) < \text{PAR}_{\text{max}}^{\text{nitrif}} \text{ and } O_2 > O_2^{\text{min}} \\ 0, & \text{elsewhere} \end{cases} \quad (\text{A166})$$

$$\begin{aligned} J_{o2}^{\text{cons}} = & (J_{poc}^{\text{remin}} + J_{doc,dic} + J_{zoo,dic}^{l} + J_{sp,dic}^{l} \\ & + J_{sp,dic}^{grz} + J_{diat,dic}^{l} + J_{diaz,dic}^{grz} \\ & + J_{diaz,dic}^{l} + J_{diaz,dic}^{grz} + J_{remin}^{poc})/Q_{C:O} \\ & + 1.5 J^{\text{ammox}} + 0.5 J^{\text{nit}} \end{aligned} \quad (\text{A167})$$

$$J_{no3,no2}^{\text{denit}} = \begin{cases} Q_{N:C}^{\text{denit}} J_{poc}^{\text{remin}}, & \text{where } O_2 \le O_2^{\text{min}} \\ 0, & \text{elsewhere} \end{cases} \quad (\text{A168})$$

898

$$J_{no2,n2o}^{\text{denit}} = \begin{cases} f_{n2o}^{\text{denit}} J_{no3,no2}^{\text{denit}}, & \text{where } O_2 \le O_2^{\text{min}} \\ 0, & \text{elsewhere} \end{cases} \quad (\text{A169})$$

$$J_{n2o,n2}^{\text{denit}} = \begin{cases} \tau_{n2o}^{\text{cons}} N_{n2o}, & \text{where } O_2 \le O_2^{\text{min}} \\ 0, & \text{elsewhere} \end{cases} \quad (\text{A170})$$
$$J_{sed}^{denit} = \begin{cases} (a_{0}^{bohl} + a_{1}^{bohl} (a_{2}^{bohl})^{O_{2}-N_{no3}}) \frac{\Phi_{poc}^{out}}{\Delta z}, & \text{where } z = z_{bottom} \\ 0, & \text{elsewhere} \end{cases} \qquad (A171)$$

900 Air-sea fluxes

$$SC_{O_2} = 1638.0 + T_{surf}(-81.83 + T_{surf}(-0.008004(1.483 + T_{surf}))) \qquad (A172)$$

901

$$SC_{CO_2} = 2073.1 + T_{surf}(-125.62 + T_{surf}(-0.043219(3.6276 + T_{surf}))) \qquad (A173)$$

902

$$PV_{O_2} = 8.6 \cdot 10^{-7} W s^2 \sqrt{\frac{660}{Sc_{O_2}}} \qquad (A174)$$

903

$$PV_{CO_2} = 8.6 \cdot 10^{-7} W s^2 \sqrt{\frac{660}{Sc_{CO_2}}} \qquad (A175)$$

904

$$\Phi_{O_2}^{air} = PV_{O_2}(O_{2sat} - O_2) \qquad (A176)$$

905

906

$O_{2sat}$ is the oxygen saturation concentration computed from Garcia and Gordon (1992), page 1310, eq. 8.

$$\Phi_{N_2O}^{air} = -PV_{O_2}N_2O \qquad (A177)$$

907

$$\Phi_{N_2}^{air} = -PV_{O_2}N_2 \qquad (A178)$$

908

$$\Phi_{CO_2}^{air} = PV_{CO_2}(CO_{2air}^{*} - CO_2^{*}) \qquad (A179)$$

909

910

with $CO_{2air}^*$ and $CO_2^*$ the concentrations corresponding to the partial pressure in air and water computed from Dickson and Goyet (1994) (SOP No. 3, p25-26).

911 Concentration in sediment

$$POC_{sed} = \Phi_{poc}^{out} - J_{poc}^{remin,sed} \qquad (A180)$$

912

$$CaCO_3_{sed} = \Phi_{pic}^{out} - J_{pic}^{remin,sed} \qquad (A181)$$

913

$$Si_{sed} = \Phi_{Si}^{out} - J_{Si}^{remin,sed} \qquad (A182)$$
References
914

915 Andersson, A. J., F. T. Mackenzie, and N. R. Bates, 2008: Life on the margin: implications of
916 ocean acidification on Mg-calcite, high latitude and cold-water marine calcifiers. *Mar. Ecol.
917 Progr. Ser.*, **373**, 265–273.

918 Armstrong, R. A., C. Lee, J. I. Hedges, S. Honjo, and S. G. Wakeham, 2001: A new, mechanistic
919 model for organic carbon fluxes in the ocean based on the quantitative association of POC with
920 ballast minerals. *Deep-Sea Res. II*, **49**, 219–236.

921 Bakker, D., and coauthors, 2016: A multi-decade record of high-quality fCO2 data in ver-
922 sion 3 of the Surface Ocean CO2 Atlas (SOCAT). *Earth System Science Data*, **8**, 383–413,
923 doi:https://doi.org/10.5194/essd-8-383-2016.

924 Bakun, A., 1990: Global climate change and intensification of coastal ocean upwelling. *Science*, **247**, 198–201.

926 Banas, N., and B. Hickey, 2008: Why is the northern end of the California Current System so
927 productive? *Oceanography*, **21**, 90–107.

928 Behrenfeld, M. J., and P. G. Falkowski, 1997: Photosynthetic rates derived from satellite-based
929 chlorophyll concentration. *Limnol. Oceanogr.*, **42**, 1–20.

930 Bianchi, D., T. S. Weber, R. Kiko, and C. Deutsch, 2018: Global niche of marine anaer-
931 obic metabolisms expanded by particle microenvironments. *Nature Geosci.*, **11**, 263–268,
932 doi.org/10.1038/s41561-018-0081-0.

933 Bograd, S. J., and R. J. Lynn, 2001: Physical-biological coupling in the California Current during
934 the 1997–99 El Niño-La Niña cycle. *Geophys. Res. Lett.*, **28**, 275–278.

935 Bohlen, L., A. W. Dale, and K. Wallmann, 2012: Simple transfer functions for calculating benthic
936 fixed nitrogen losses and C:N:P regeneration ratios in global biogeochemical models. *Global
937 Biogeochem. Cycles*, **26**, GB3029, 1–16.

938 Buil, M. P., and E. D. Lorenzo, 2017: Decadal dynamics and predictability of oxygen and
939 subsurface tracers in the California Current System. *Geophys. Res. Lett.*, **44**, 4204–4213,
940 doi:10.1002/2017GL072931.

941 Busch, D. S., and P. McElhany, 2016: Estimates of the direct effect of seawater pH
942 on the survival rate of species groups in the California Current ecosystem. *PLOS One*,
943 doi:10.1371/journal.pone.0160669.

944 Capet, X., J. McWilliams, M. Molemaker, and A. Shchepetkin, 2008: Mesoscale to submesoscale
945 transition in the California Current System. Part I: Flow structure, eddy flux, and observational
946 tests. *Journal of Physical Oceanography*, **38** (1), 29–43.

947 Carr, M.-E., and E. J. Kearns, 2003: Production regimes in four Eastern Boundary Current Sys-
948 tems. *Deep Sea Res. II*, **50**, 3199–3221.
949 Chan, F., J. Barth, J. Lubchenco, A. Kirincich, H. Weeks, W. T. Peterson, and B. Menge, 2008: Emergence of anoxia in the California Current large marine ecosystem. *Science*, **319**, 920–920.

951 Chavez, F., and Coauthors, 2002: Biological and chemical consequences of the 1997–1998 El Niño in central California waters. *Progr. Oceanogr.*, **54**, 205–232.

953 Chavez, F. P., and M. Messie, 2009: A comparison of eastern boundary upwelling ecosystems. *Progr. Oceanogr.*, **83**, 80–96.

955 Chavez, F. P., J. Ryan, S. E. Lluch-Cota, and M. Ñiquen, 2003: From anchovies to sardines and back: multidecadal change in the Pacific Ocean. *Science*, **299**, 217–221.

957 Collins, W., and Coauthors, 2011: Development and evaluation of an Earth-System model–HadGEM2. *Geosci. Model Devel.*, **4**, 1051–1075.

959 Davis, K. A., N. S. Banas, S. N. Giddings, S. A. Siedlecki, P. MacCready, E. J. Lessard, R. M. Kudela, and B. M. Hickey, 2014: Estuary-enhanced upwelling of marine nutrients fuels coastal productivity in the U.S. Pacific Northwest. *J. Geophys. Res. Oceans*, **119**, 8778–8799.

962 Deutsch, C., H. Brix, T. Ito, H. Frenzel, and L. Thompson, 2011: Climate-forced variability of ocean hypoxia. *Science*, **333**, 336–339.

964 Deutsch, C., S. Emerson, and L. Thompson, 2006: Fingerprints of climate change in North Pacific oxygen. *Geophys. Res. Lett.*, **32**, 1–4.

966 Deutsch, C., A. Ferrel, B. Seibel, H.-O. Pörtner, and R. B. Huey, 2015: Climate change tightens a metabolic constraint on marine habitats. *Science*, **348**, 1132–1135.

968 Di Lorenzo, E., and Coauthors, 2008: North Pacific Gyre Oscillation links ocean climate and ecosystem change. *Geophys. Res. Lett.*, **35**, doi:10.1029/2007GL032838.

970 Dickson, A. G., and C. Goyet, 1994: Handbook of methods for the analysis of the various parameters of the carbon dioxide system in sea water. version 2.

972 Doney, S. C., N. Mahowald, I. Lima, R. A. Feely, F. T. Mackenzie, J.-F. Lamarque, and P. J. Rasch, 2007: Impact of anthropogenic atmospheric nitrogen and sulfur deposition on ocean acidification and the inorganic carbon system. *Proc. Nat. Acad. Sci.*, **104**, 14 580–14 585.

975 Ducklow, H. W., S. C. Doney, and D. K. Steinberg, 2009: Contributions of long-term research and time-series observations to marine ecology and biogeochemistry. *Annu. Rev. Marine. Sci.*, **1**, 279–302.

978 Dugdale, R. C., F. P. Wilkerson, R. T. Barber, and F. P. Chavez, 1992: Estimating new production in the equatorial Pacific Ocean at 150 W. *J. Geophys. Res. Oceans*, **97**, 681–686.

980 Dunne, J. P., R. A. Armstrong, A. Gnanadesikan, and J. L. Sarmiento, 2005: Empirical and mechanistic models for the particle export ratio. *Global Biogeochem. Cycles*, **19**, doi:10.1029/2004GB002390.
Durski, S. M., J. A. Barth, J. C. McWilliams, H. Frenzel, and C. Deutsch, 2017: The influence of variable slope-water characteristics on dissolved oxygen levels in the northern California Current System. *J. Geophys. Res. Oceans*, **122**, 7674–7697, doi:10.1002/2017JC013089.

Dussin, R., E. N. Curchitser, C. A. Stock, and N. V. Oostende, 2019: Biogeochemical drivers of changing hypoxia in the California Current Ecosystem. *Deep-Sea Res. II*, **169-170**, doi:10.1016/j.dsr2.2019.05.013.

Eppley, R. W., F. P. Chavez, and R. T. Barber, 1992: Standing stocks of particulate carbon and nitrogen in the equatorial Pacific at 150° W. *J. Geophys. Res. Oceans*, **97**, 655–661.

FAO, 2009: State of the World's Fisheries and Aquaculture 2008. *Food and Agriculture Organization of the United Nations*.

Feely, R. A., C. L. Sabine, J. M. Hernandez-Ayon, D. Ianson, and B. Hales, 2008: Evidence for upwelling of corrosive “acidified” water onto the continental shelf. *Science*, **320**, 1490–1492.

Fiechter, J., E. N. Curchitser, C. A. Edwards, F. Chai, N. L. Goebel, and F. P. Chavez, 2014: Air-sea CO₂ fluxes in the California Current: Impacts of model resolution and coastal topography. *Global Biogeochem. Cycles*, **28**, doi:10.1002/2013GB004683.

Firme, G. F., E. L. Rue, D. A. Weeks, K. W. Bruland, and D. A. Hutchins, 2003: Spatial and temporal variability in phytoplankton iron limitation along the California coast and consequences for Si, N, and C biogeochemistry. *Global Biogeochem. Cycles*, **17**, doi:10.1029/2001GB001824.

Frenger, I., D. Bianchi, C. Stuhrenberg, A. Oschlies, J. Dunne, C. Deutsch, E. Galbraith, and F. Schutte, 2018: Biogeochemical role of subsurface coherent eddies in the ocean: Tracer cannonballs, hypoxic storms, and microbial stewpots? *Global Biogeochem. Cycles*, **32**, doi:10.1002/2017GB005743.

Frischknecht, M., M. Münnich, and N. Gruber, 2018: Origin, transformation, and fate: The three-dimensional biological pump in the California Current System. *J. Geophys. Res. Oceans*, **123**, doi:10.1029/2018JC013934.

Garcia, H. E., and L. I. Gordon, 1992: Oxygen solubility in seawater: Better fitting equations. *Limnol. Oceanogr.*, **37 (6)**, 1307–1312.

Garcia-Reyes, M., J. L. Largier, and W. J. Sydeman, 2014: Synoptic-scale upwelling indices and predictions of phyto- and zooplankton populations. *Progr. Oceanogr.*, **120**, 177–188.

Geider, R. J., H. L. Maclntyre, and T. M. Kana, 1998: A dynamic regulatory model of phytoplanktonic acclimation to light, nutrients, and temperature. *Limnol. Oceanogr.*, **43**, 679–694.

Gruber, N., C. Hauri, Z. Lachkar, D. Loher, T. L. Frölicher, and G.-K. Plattner, 2012: Rapid progression of ocean acidification in the California Current System. *Science*, **337**, 220–223.

Gruber, N., Z. Lachkar, H. Frenzel, P. Marchesiello, M. Münnich, J. C. McWilliams, T. Nagai, and G.-K. Plattner, 2011: Eddy-induced reduction of biological production in eastern boundary upwelling systems. *Nature Geosci.*, **4**, 787–792.
1019 Hales, B., P. G. Strutton, M. Saraceno, R. Letelier, T. Takahashi, R. Feely, C. Sabine, and F. Chavez, 2012: Satellite-based prediction of pCO2 in coastal waters of the eastern North Pacific. *Progr. Oceanogr.*, doi:10.1016/j.pocean.2012.03.001.

1022 Howard, E. M., and Coauthors, 2020: Climate-driven aerobic habitat loss in the California Current System. *Science Adv.*, **6**, doi:eaay3188.

1024 Ito, T., and C. Deutsch, 2010: A conceptual model for the temporal spectrum of oceanic oxygen variability. *Geophys. Res. Lett.*, **37**, doi:10.1029/2009GL041595.

1026 Ito, T., M. C. Long, C. Deutsch, S. Minobe, and D. Sun, 2019: Mechanisms of low-frequency oxygen variability in the North Pacific. *Global Biogeochem. Cycles*, **33**, doi:10.1029/2018GB005987.

1029 Jacox, M. G., C. A. Edwards, E. L. Hazen, and S. J. Bograd, 2018: Coastal Upwelling Revisited: Ekman, Bakun, and Improved Upwelling Indices for the U.S. West Coast. *J. Geophys. Res. Oceans*, **123**, 7332–7350, doi:10.1029/2018JC014187.

1032 Kahru, M., R. Kudela, M. Manzano-Sarabia, and B. G. Mitchell, 2009: Trends in primary production in the California Current detected with satellite data. *J. Geophys. Res. Oceans*, **114**, doi:10.1029/2008JC004979.

1035 Kessouri, F., D. Bianchi, L. Renault, J. C. McWilliams, H. Frenzel, and C. Deutsch, 2020: Sub-mesoscale currents modulate the seasonal cycle of nutrients and productivity in the California Current System. *Global Biogeochem. Cycles*, **34**, doi:10.1029/2020GB006578.

1038 Key, R. M., and Coauthors, 2004: A global ocean carbon climatology: Results from Global Data Analysis Project (GLODAP). *Global Biogeochem. Cycles*, **18**, doi:10.1029/GB002247.

1040 Kudela, R. M., and Coauthors, 2008: New insights into the controls and mechanisms of plankton productivity in coastal upwelling waters of the Northern California Current System. *Oceanography*, **21**, 46–59.

1043 Kwon, E. Y., C. Deutsch, S.-P. Xie, S. Schmidtko, and Y.-K. Cho, 2016: The North Pacific oxygen uptake rates over the past half century. *J. Climate*, **29**, doi:10.1175/JCLI-D-14-00157.1.

1045 Large, W. B., 2006: Surface fluxes for practitioners of global ocean data assimilation. *Ocean Weather Forecasting*, Springer, 229–270.

1047 Lauvset, S. K., and coauthors, 2016: A new global interior ocean mapped climatology: the 1° × 1° GLODAP version 2. *Earth System Science Data*, **8**, 325–340, doi:10.5194/essd-8-325-2016.

1049 Locarnini, R., and Coauthors, 2013: World Ocean Atlas 2013: Temperature. *NOAA NESDIS*, **1**.

1050 Long, M. C., C. Deutsch, and T. Ito, 2016: Finding forced trends in oceanic oxygen. *Global Biogeochem. Cycles*, **30**, 381–397, doi:10.1002/2015gb005310.
1052
1053
1054
1055

Mahowald, N. M., M. Yoshioka, W. D. Collins, A. J. Conley, D. W. Fillmore, and D. B. Coleman, 2006: Climate response and radiative forcing from mineral aerosols during the last glacial maximum, pre-industrial, current and doubled-carbon dioxide climates. *Geophys. Res. Lett.*, **33**, doi:10.1029/2006GL026126.

1056
1057

Mantua, N. J., S. R. Hare, Y. Zhang, J. M. Wallace, and R. C. Francis, 1997: A Pacific interdecadal climate oscillation with impacts on salmon production. *Bull. Am. Meteor. Soc.*, **78**, 1069–1080.

1058
1059
1060

McClatchie, S., R. Goericke, R. Cosgrove, G. Auad, and R. Vetter, 2010: Oxygen in the Southern California Bight: Multidecadal trends and implications for demersal fisheries. *Geophys. Res. Lett.*, **37**, doi:10.1029/2010GL044497.

1061
1062
1063
1064

Mecking, S., C. Langdon, R. A. Feely, C. L. Sabine, C. A. Deutsch, and D.-H. Min, 2008: Climate variability in the North Pacific thermocline diagnosed from oxygen measurements: An update based on the US CLIVAR/CO2 Repeat Hydrography cruises. *Global Biogeochem. Cycles*, **22**, doi:10.1029/2007GB003101.

1065
1066

Meinvielle, M., and G. C. Johnson, 2013: Decadal water-property trends in the California Under-current, with implications for ocean acidification. *J. Geophys. Res. Oceans*, **118**, 6687–6703.

1067
1068

Middelburg, J. J., K. Soetaert, P. M. Herman, and C. H. Heip, 1996: Denitrification in marine sediments: A model study. *Global Biogeochem. Cycles*, **10**, 661–673.

1069
1070

Moore, J. K., S. C. Doney, J. A. Kleypas, D. M. Glover, and I. Y. Fung, 2002: An intermediate complexity marine ecosystem model for the global domain. *Deep-Sea Res. II*, **49**, 403–462.

1071
1072
1073

Moore, J. K., S. C. Doney, and K. Lindsay, 2004: Upper ocean ecosystem dynamics and iron cycling in a global three-dimensional model. *Global Biogeochem. Cycles*, **18**, GB4028, doi:10.1029/2004GB002220.

1074
1075
1076

Munro, D. R., P. D. Quay, L. W. Juranek, and R. Goericke, 2013: Biological production rates off the Southern California coast estimated from triple O₂ isotopes and O₂: Ar gas ratios. *Limnol. Oceanogr.*, **58**, 1312–1328.

1077
1078
1079

Murray, J. W., J. Young, J. Newton, J. Dunne, T. Chapin, B. Paul, and J. J. McCarthy, 1996: Export flux of particulate organic carbon from the central equatorial Pacific determined using a combined drifting trap-²³⁴Th approach. *Deep-Sea Res. II*, **43**, 1095–1132.

1080
1081
1082

Nagai, T., N. Gruber, H. Frenzel, Z. Lachkar, J. C. McWilliams, and G. K. Plattner, 2015: Dominant role of eddies and filaments in the offshore transport of carbon and nutrients in the California Current System. *J. Geophys. Res. Oceans*, **120**, 5318–5341.

1083
1084
1085

Palevsky, H. I., P. D. Quay, and D. P. Nicholson, 2016: Discrepant estimates of primary and export production from satellite algorithms, a biogeochemical model, and geochemical tracer measurements in the North Pacific Ocean. *Geophys. Res. Lett.*, **43**, 8645–8653.

1086
1087
1088

Pierce, D. W., P. J. Gleckler, T. P. Barnett, B. D. Santer, and P. J. Durack, 2012: The fingerprint of human-induced changes in the ocean's salinity and temperature fields. *Geophys. Res. Lett.*, **39**, L21704, doi:10.1029/2012GL053389.
1089
1090
1091

Renault, L., C. Deutsch, J. C. McWilliams, F. Colas, H. Frenzel, and J. H. Liang, 2016a: Partial decoupling of primary productivity from upwelling in the California Current System. *Nature Geosci.*, **9**, 505–508.

1092
1093
1094

Renault, L., J. C. McWilliams, A. Jousse, C. Deutsch, H. Frenzel, F. Kessouri, and R. Chen, 2021: The physical structure and behavior of the California Current System. *Progr. Oceanogr.*, submitted.

1095
1096
1097

Renault, L., M. J. Molemaker, J. C. McWilliams, A. F. Shchepetkin, F. Lemarié, D. Chelton, S. Illig, and A. Hall, 2016b: Modulation of wind work by oceanic current interaction with the atmosphere. *J. Phys. Ocean.*, **46**, 1685–1704.

1098
1099

Roemmich, D., and T. McCallister, 1989: Large scale circulation of the North Pacific Ocean. *Progr. Oceanogr.*, **22**, 171–204.

1100
1101

Rykaczewski, R. R., and D. M. Checkley, 2008: Influence of ocean winds on the pelagic ecosystem in upwelling regions. *Proc. Nat. Acad. Sci.*, **105**, 1965–1970.

1102
1103
1104

Severmann, S., J. McManus, W. M. Berelson, and D. E. Hammond, 2010: The continental shelf benthic iron flux and its isotope composition. *Geochimica et Cosmochimica Acta*, **74**, 3984–4004.

1105
1106
1107

Shchepetkin, A. F., and J. C. McWilliams, 2005: The Regional Oceanic Modeling System (ROMS): A split-explicit, free-surface, topography-following-coordinate oceanic model. *Ocean Modelling*, **9**, 347–404.

1108
1109
1110

Stukel, M. R., M. R. Landry, C. R. Benitez-Nelson, and R. Goericke, 2011: Trophic cycling and carbon export relationships in the California Current Ecosystem. *Limnol. Oceanogr.*, **56**, 1866–1878.

1111
1112

Stukel, M. R., and Coauthors, 2017: Mesoscale ocean fronts enhance carbon export due to gravitational sinking and subduction. *Proc. Nat. Acad. Sci.*, **114**, 1252–1257.

1113
1114
1115

Turi, G., Z. Lachkar, and N. Gruber, 2014: Spatiotemporal variability and drivers of pCO₂ and air–sea CO₂ fluxes in the California Current System: an eddy-resolving modeling study. *Bio-geosci.*, **11**, 671–690, doi:https://doi.org/10.5194/bg-11-671-2014.

1116
1117
1118

Turi, G., and Coauthors, 2018: Response of O₂ and pH to ENSO in the California Current System in a high-resolution global climate model. *Ocean Science*, **14**, 69–86, doi:10.5194/os-14-69-2018.

1119
1120
1121
1122

van Heuven, S., D. Pierrot, J. Rae, E. Lewis, and D. Wallace, 2011: *MATLAB Program Developed for CO₂ System Calculations*. ORNL/CDIAC-105b. Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy, Oak Ridge, Tennessee, doi:10.3334/CDIAC/otg.CO2SYS_MATLAB_v1.1.

1123
1124

Wanninkhof, R., 1992: Relationship between wind speed and gas exchange over the ocean. *J. Geophys. Res. Oceans*, **97**, 7373–7382, doi:10.1029/92JC00188.
1125 Westberry, T., M. Behrenfeld, D. Siegel, and E. Boss, 2008: Carbon-based primary productivity modeling with vertically resolved photoacclimation. *Global Biogeochem. Cycles*, **22**, GB2024, doi:10.1029/2007GB003078.

1128 Zweng, M., and Coauthors, 2013: World Ocean Atlas 2013,: Salinity. *NOAA NESDIS*, **2**.