Variable organic matter stoichiometry enhances the biological drawdown of CO₂ in the northwest European shelf seas

Variations in the elemental ratios of carbon, nitrogen, and phosphorus in marine organic matter (OM) and their influence on the marine carbon cycle remain poorly understood for both the open and coastal oceans. Observations consistently show an enrichment of carbon and a depletion of phosphorus relative to elemental Redfield ratios. However, many biogeochemical models are constrained to Redfield stoichiometry, neglecting the effects of variable stoichiometry on carbon cycling and typically underestimating biological carbon fixation. This impedes the accurate representation of OM cycling and the resulting carbon fluxes, especially in productive temperate shelf seas such as the northwest European shelf seas (NWES). Here, the efficiency of oceanic CO₂ uptake strongly depends on the biological uptake of inorganic carbon and its export to the North Atlantic, both of which are influenced by OM stoichiometry. In this study, we provide a first comprehensive and quantitative assessment of the effects of variable OM stoichiometry on carbon cycling in the NWES. For this purpose, we integrate two pathways for variable OM stoichiometry, motivated by observational and experimental results, into the regional high-resolution coupled 3D physical–biogeochemical modeling system SCHISM-ECOSMO-CO₂ (Semi-implicit Cross-scale Hydroscience Integrated System Model – ECOSystem MOdel): first, the release of carbon-enriched dissolved OM under nutrient limitation and, second, the preferential remineralization of organic nitrogen and phosphorus. With these extensions we reproduce the observed OM stoichiometry and evaluate its impact on marine carbon cycling, with a focus on OM cycling and the resulting air–sea CO₂ exchange. Compared to the reference simulation with fixed Redfield stoichiometry, the variable stoichiometry configurations show an increase in the annual net CO₂ uptake of 10 %–33 % in the North Sea and 9 %–31 % in the entire NWES, depending on the relative contribution of the two new implementations. As the main driver of the additional CO₂ uptake, we identify a corresponding intensification of annual and seasonal OM cycling, resulting in higher net autotrophy in surface waters and higher net heterotrophy in sub-surface layers. This enhanced gradient in net community production leads to an increased biological drawdown of inorganic carbon, most pronounced in the Norwegian Trench. By increasing the biological control on the surface partial pressure of CO₂, this leads to higher summer and lower winter uptake. Our results highlight the importance of variable stoichiometry for an accurate representation of the shelf carbon pump mechanism in the NWES, as it significantly influences the efficiency of carbon sequestration. Since the response depends largely on regional physical conditions and pre-existing carbon export mechanisms, regional assessments are essential to understand the sensitivity of the carbon cycle to OM stoichiometry, which should be included in global models to accurately represent the coastal carbon cycle.

1 Introduction

Marine organic matter (OM) is a major global reservoir of reduced carbon (Friedlingstein et al., 2023; Hansell et al., 2009) and plays a key role in sequestering atmospheric CO₂ in the ocean (Carlson and Hansell, 2015; Falkowski et al., 1998; Hansell et al., 2009). The balance between OM production and decomposition controls the transformation and transports of both organic and inorganic carbon (del Giorgio and Duarte, 2002; Smith and Hollibaugh, 1993), especially in the shallow and highly productive temperate shelf seas (Barrón and Duarte, 2015; Bauer et al., 2013; Canuel et al., 2012; Dai et al., 2022; Gattuso et al., 1998). Through its chemical composition, OM further regulates the coupling between the elemental cycles of carbon, nitrogen, and phosphorus, with implications for the relative magnitudes of elemental fluxes, nutrient availability, and OM production (Hopkinson and Vallino, 2005). In this way, variations in OM composition may affect the biological carbon drawdown and thereby the air–sea CO₂ exchange through both changes in carbon fixation and variations in the composition of the OM that is available for lateral and vertical transports (Hopkinson and Vallino, 2005; Loh and Bauer, 2000; Lønborg and Álvarez-Salgado, 2012; Tanioka et al., 2021; Williams, 1995; Zakem and Levine, 2019).

Despite observational evidence for large global and regional variations in elemental $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratios in dissolved (DOM) and particulate (POM) organic matter (Hopkinson and Vallino, 2005; Liang et al., 2023; Loh and Bauer, 2000; Lønborg et al., 2024; Martiny et al., 2013; Tanioka et al., 2022b), many biogeochemical models assume a constant Redfield stoichiometry of $\mathrm{C}:\mathrm{N}:\mathrm{P}=\mathrm{106}:\mathrm{16}:\mathrm{1}$ or similar, originally derived as an average ratio for the entire global ocean (Redfield, 1963). As a consequence, these models tend to underestimate carbon fixation, especially during nutrient limitation in summer, which is indicated by their deficiency in representing the seasonality of the dissolved inorganic carbon (DIC) concentration and partial pressure of CO₂ (pCO₂) in surface waters (Bozec et al., 2006; Kähler and Koeve, 2001; Prowe et al., 2009). Hence, variable stoichiometry in OM production and decomposition is essential for the adequate representation of biogeochemical fluxes and especially their seasonality (Anderson and Pondaven, 2003; Bozec et al., 2006). Previous experiments using Earth system models suggest that the globally observed carbon enrichment in both DOM and POM significantly enhances global carbon and nitrogen fixation, drawdown of CO₂, and export production when compared to Redfield-stoichiometry-based estimates (Letscher et al., 2015; Letscher and Moore, 2015; Tanioka et al., 2021). However, the regional impact on coastal seas and its spatiotemporal variability remain uncertain due to the scarcity of measurements and the resulting missing consideration in regional studies (Aricò et al., 2021).

The global coastal ocean disproportionately contributes to the sequestration of atmospheric CO₂ in relation to its areal extent through carbon transport to the deep ocean and burial in sediments (Bauer et al., 2013; Dai et al., 2022; Laruelle et al., 2014). Its role is suggested to have changed from a weak pre-industrial net source to a significant contemporary net sink of CO₂ due to both rising atmospheric CO₂ concentrations and enhanced productivity driven by higher terrestrial nutrient inputs (Bauer et al., 2013; Lacroix et al., 2021; Laruelle et al., 2018; Mathis et al., 2024). With a regionally significant contribution (Frankignoulle and Borges, 2001), the northwest European shelf seas (NWES) efficiently export atmospheric CO₂ to the North Atlantic Ocean through the shelf carbon pump mechanism (Tsunogai et al., 1999; Thomas et al., 2004; Bozec et al., 2005; Legge et al., 2020).

In the seasonally stratified parts of the central and northern North Sea, the spatial separation of production in the surface mixed layer and respiration in the sub-surface layer lead to a seasonally occurring drawdown of CO₂, resulting in a vertical gradient of DIC. Subsequently, the DIC-enriched sub-surface waters are transported across the shelf edge, leading to a net export of DIC (Thomas et al., 2004). This DIC export is further driven by riverine and Baltic Sea discharge of DIC and net heterotrophy, where the net import of organic carbon, largely from cross-shelf transports, increases net respiration and thereby adds DIC to the system (Thomas et al., 2005). With a smaller contribution, particulate organic carbon (POC) is exported by net deposition and burial in sediments, mainly in the deep regions of the Norwegian Trench, the Skagerrak, and along the shelf edge (Legge et al., 2020; Thomas et al., 2005). The sedimentation of POC is associated with high uncertainties and may play a larger role on long timescales (Legge et al., 2020; Thomas et al., 2005).

This regional carbon export mechanism in the NWES may strongly depend on the observed variations in OM stoichiometry, which makes them essential for regional model-based carbon budget estimates (Chaichana et al., 2019; Davis et al., 2019). Reproducing the observed OM stoichiometry may account for the typically underestimated carbon fixation (Humphreys et al., 2019; Prowe et al., 2009) and consequently affect the estimated net air–sea CO₂ exchange, cross-shelf transport of DIC, and vertical transport of organic carbon in the NWES through changes in the biological carbon pump and the resulting drawdown of DIC.

In this work, we assess the effects of observed variations in OM stoichiometry on carbon fixation, respiration, biological drawdown of inorganic carbon, and the implications for the air–sea CO₂ exchange in the NWES. This serves as a first step in identifying the regional impact on marine carbon cycling, with a focus on the continental shelf carbon pump mechanism. For this, we integrate two pathways of variable OM stoichiometry, motivated by observational and experimental results, into a regional high-resolution coupled 3D physical–biogeochemical ocean model.

Firstly, we introduce an extracellular release of carbon-enriched DOM under nutrient limitation, modified from Neumann et al. (2022). While primary production at Redfield stoichiometry dominates biological carbon fixation under nutrient-replete conditions, the extracellular release allows for additional carbon fixation beyond nutrient limitation and, with variable OM stoichiometry, depending on nutrient availability. This release of carbon-enriched DOM has been observed as an extracellular release of carbohydrates under nutrient stress by phytoplankton cells in mesocosm experiments (Børsheim et al., 2005; Fajon et al., 1999; Søndergaard et al., 2000), as well as field measurements during bloom periods (Myklestad, 1995; Wear et al., 2015b; Williams, 1995). Estimates suggest that the extracellular release of DOC on average contributes 10 %–20 % but temporarily up to 80 % of primary productivity (Wear et al., 2015a). Accordingly, global observations suggest that stoichiometric variations depend on macronutrient availability, specifically nitrogen-versus-phosphorous stress, in addition to temperature (Tanioka et al., 2022b).

Secondly, we include the preferential remineralization of organic nitrogen and phosphorus with the sequence P>N>C for both DOM and POM. Both observations (Clark et al., 1998; Hopkinson et al., 1997, 2002; Loh and Bauer, 2000; Thomas et al., 1999; Williams et al., 1980) and experimental evidence (Hach et al., 2020; Lønborg et al., 2009; Lønborg and Álvarez-Salgado, 2012) suggest the preferential utilization of organic nitrogen and phosphorus over carbon due to higher bioavailability as a mechanism for variable OM stoichiometry. Assuming remineralization at Redfield ratios thus implies underestimating the recycling of nutrients, new production, and carbon export (Hach et al., 2020; Lønborg and Álvarez-Salgado, 2012). This process is thought to contribute to the slow formation of strongly nutrient-depleted and carbon-enriched refractory DOM from labile DOM (Cai and Jiao, 2023; Hach et al., 2020), with refractory DOM having an estimated average lifetime of approximately 6000 years (Bauer et al., 1992).

While previous studies have simulated OM cycling with variable OM stoichiometry for the North Sea (Butenschön et al., 2016; Lorkowski et al., 2012) and the Baltic Sea (Neumann et al., 2022), this work represents the first comprehensive and quantitative assessment of the effects of observed variations in OM stoichiometry on marine carbon cycling in the NWES. This study quantifies the regional impacts of two pathways for variable OM stoichiometry on OM cycling and the resulting air–sea CO₂ exchange in the NWES. In doing so, it will improve our understanding of these implications for marine carbon cycling and the limitations of Redfield-stoichiometry-based estimates. To achieve an accurate representation of the magnitude of both mechanisms, we account for the observed range of POM and labile DOM stoichiometry. We assess the effects of the individual and combined mechanisms by comparing four hindcast simulations: a reference configuration using Redfield stoichiometry, separate configurations for each mechanism, and a final one that integrates both with a reduced contribution. These simulations are conducted using the SCHISM-ECOSMO-CO₂ modeling system from Kossack et al. (2023, 2024) and span the period from 1995 to 2010, covering the NWES, the Baltic Sea, and parts of the northeastern Atlantic Ocean.

2 Materials and methods

2.1 SCHISM-ECOSMO-CO₂ modeling system

The SCHISM-ECOSMO-CO₂ modeling system couples 3D hydrodynamic, biogeochemical, and carbonate system models through the Framework for Aquatic Biogeochemical Models (FABM; Bruggeman and Bolding, 2014). This one-way coupling allows us to consider the control of dynamic forcing conditions on ecosystem and carbon dynamics, as well as the effects of biological production, respiration, and the physical state on carbon cycling (Fig. 1). The Semi-implicit Cross-scale Hydroscience Integrated System Model (SCHISM; Zhang et al., 2016b) simulates the physical state and the resulting tracer transports over the ocean domain. Under these physical forcing conditions, the here-introduced variable stoichiometry version of the ecosystem model ECOSMO II (ECOSystem MOdel; Daewel and Schrum, 2013; Schrum et al., 2006) prognostically calculates tracer concentrations through local sources and sinks from biogeochemical processes. Based on simulated temperature, pressure, salinity, and two prognostic carbonate system variables provided by ECOSMO II, the carbonate system model by Blackford and Gilbert (2007) diagnostically computes the missing carbonate system variables and the air–sea CO₂ exchange under equilibrium conditions. The carbonate system calculations follow the HALTAFALL algorithm (Ingri et al., 1967) for the composition of equilibrium mixtures.

Figure 1 Schematic diagram of the coupled physical–biogeochemical modeling system SCHISM-ECOSMO-CO₂ based on the FABM framework. The modeling system couples the 3D baroclinic circulation model SCHISM, an extended version of the lower-trophic-level ecosystem model ECOSMO II, and a carbonate system model based on the HALTAFALL algorithm to calculate the composition of equilibrium mixtures.

2.1.1 Hydrodynamic model SCHISM

The hydrodynamical model component SCHISM is capable of efficiently resolving 3D baroclinic ocean circulation in the NWES and adjacent areas across a range of depths from the shallow southern North Sea to the deep waters of the northeastern Atlantic (Zhang et al., 2016a). This is achieved through the use of an unstructured triangular horizontal grid (Zhang et al., 2016b) in combination with localized sigma coordinates with shaved cells (LSC²) for the vertical grid structure (Zhang et al., 2015). Previous studies have demonstrated that SCHISM is able to accurately resolve transports across shelf slopes and canyons (Wang et al., 2022). This makes it particularly suitable for simulating cross-shelf transports and resolving the Norwegian Trench, the Skagerrak, and shelf slope areas in the NWES. The adaptive vertical resolution further enhances the representation of bottom and slope currents, with relevance for the resulting benthic pelagic coupling. By employing a third-order finite-volume transport scheme built upon the weighted essentially non-oscillatory (WENO) formalism, SCHISM is capable of simulating mesoscale eddies in the open ocean (Ye et al., 2019), thereby facilitating cross-scale applications.

The model simulations for this study employ the NWES-LR grid, initialization, and forcing configuration from Kossack et al. (2023, 2024) and Porz et al. (2024), which successfully reproduced temperature, nutrient, and carbonate system variability in the NWES. Our simulations cover the ocean domain from 40–66° N and 20° W–30° E, including the NWES, the Baltic Sea, and parts of the northeastern Atlantic Ocean for the years 1995 to 2010, with a time step of 200 s. For the purposes of this analysis, the initial 5 years are excluded as a spin-up period. With 140 152 triangular grid elements and a maximum of 53 vertical layers, the horizontal resolution ranges from about 4.5 to 10 km, with increasing depth on the continental shelf, and reaches up to 15 km in the Atlantic Ocean. While the Baltic Sea and parts of the northeast Atlantic are included for the explicit resolution of exchange fluxes, the focus of this study is on the North Sea and the entire NWES. In accordance with Kossack et al. (2023, 2024), we subdivide the NWES into sub-areas based on a combination of bathymetry and both physical and biogeochemical variability (Fig. 2). The boundary of the NWES along the shelf edge is defined by the 200 m isobath, with the exception of the Norwegian Trench, which is included because of its important role in cross-shelf exchange fluxes.

Figure 2 Sub-regions and bathymetry of the northwest European shelf seas (NWES) based on a combination of bathymetric features and both physical and biogeochemical variability. The division of the NWES into sub-regions with relatively homogeneous physical and biogeochemical characteristics allows for the identification of consistent carbon cycle responses.

2.1.2 Ecosystem model ECOSMO II

The lower-trophic-level ecosystem model ECOSMO II describes the interactions of nutrient cycles and functional groups of phytoplankton, zooplankton, detritus, and surface sediments with 16 state variables (Daewel and Schrum, 2013). Here, nutrient cycles include phosphate (PO₄), oxygen (O₂), silicate (SiO₂), and nitrogen as both ammonium (NH₄) and nitrate (NO₃). ECOSMO II represents three functional groups of phytoplankton: flagellates, limited by nitrogen and phosphate; diatoms, additionally limited by silicate; and cyanobacteria with the ability to fix nitrogen from the atmosphere at the water surface. Light limits primary production for all functional groups, with the consideration of self-shading by plankton and organic matter. Cyanobacteria production is additionally temperature-dependent and only possible under sufficient light and under low-salinity conditions in the Baltic Sea. The two zooplankton groups differ in their herbivorous and omnivorous feeding behavior.

Detritus is present in the form of POM, DOM, and biogenic opal (SiO₂ × 2H₂O), with different rates of temperature-dependent remineralization. POM and DOM contain carbon, nitrogen, and phosphorus at Redfield ratios of $\mathrm{C}:\mathrm{N}:\mathrm{P}=\mathrm{106}:\mathrm{16}:\mathrm{1}$ (Redfield, 1963). New detritus from assimilation losses and mortality separates into 60 % POM and 40 % DOM, with diatom losses further enriching biogenic opal. While POM and biogenic opal sink at constant velocities, DOM is only transported through advection and diffusion. Below a critical bottom-shear stress, both POM and opal accumulate by deposition in two corresponding sediment pools. Above this critical shear stress, they are resuspended back into the water column. Burial permanently removes both POM and silicate from the sediment pools at a constant percentage rate. Benthic remineralization releases DIC and silicate to the water column as a function of temperature only, whereas nitrogen and phosphorous release from sedimentary POM also depends on oxygen concentrations. Under oxic conditions, remineralization of particulate organic nitrogen (PON) from sediments releases ammonium. Under anoxic conditions, the release of ammonium doubles, and denitrification additionally converts nitrate as an oxidation agent into gaseous nitrogen, which becomes unavailable for photosynthesis. Remineralization converts particulate organic phosphorus (POP) into iron-bound particulate inorganic phosphorous (PIP), the third sediment pool. Dissolution of this PIP then releases phosphate back to the water column. This process increases with increasing temperature and decreasing oxygen concentrations, reaching its maximum under anoxic conditions.

2.1.3 Variable stoichiometry in ECOSMO II

We here extend the ECOSMO II formulation of pelagic and benthic OM, formerly constrained to elemental Redfield ratios, to represent variable OM stoichiometry and its consequences for carbon and nutrient cycling. For this, we introduce independent state variables for pelagic dissolved organic carbon (DOC), nitrogen (DON), and phosphorus (DOP), as well as pelagic and benthic particulate organic carbon (POC), nitrogen (PON), and phosphorus (POP). However, the internal stoichiometry of phyto- and zooplankton biomass remains at constant Redfield ratios. In accordance with Neumann et al. (2022), we assume that this provides a reasonable simplification, as healthy phytoplankton cells only marginally deviate from Redfield ratios compared to DOM and POM (Ho et al., 2003), while also limiting computational cost with reduced model complexity. However, regional variations and future trends in ecosystem stoichiometry may motivate further model developments (Finkel et al., 2010; Geider and La Roche, 2002; Sardans et al., 2021). Concentrations and rates continue to use units of grams of carbon equivalent according to Redfield ratios in conformity with the model implementation. Rates depending on OM concentrations now scale with the respective C, N, or P content. Based on these changes, we implement two pathways for variable OM stoichiometry in both carbon fixation and respiration, motivated by observational evidence.

Firstly, we include an extracellular release of carbon-enriched DOM adapted from Neumann et al. (2022). While primary and secondary production are still constrained to Redfield stoichiometry, the extracellular release allows for DOM production and hence carbon fixation beyond nutrient limitation, with variable stoichiometry depending on the nutrient conditions. The release of DOM results in a corresponding uptake of DIC and nutrients independent of phytoplankton biomass, which is described as part of the full set of ecosystem model equations in the Supplement (Sect. S2). The extracellular release base rate E defined in Eq. (1) is a function of phytoplankton biomass $C_{P_j}$, the respective maximum uptake rates σ_j, an optional temperature dependence ${\mathit{\beta }}_{T_j}$(here only relevant for cyanobacteria), and a constant scaling factor B_ER, which controls the range of stoichiometric variations and organic carbon concentrations to ensure an accurate representation. It is summed over all phytoplankton groups P_j and ensures that extracellular release only persists as long as phytoplankton biomass is present and remains on a comparable scale to primary production. As the extracellular release does not generate additional phytoplankton biomass, it will eventually decline in absence of primary production.

In contrast to primary production, the extracellular release is not co-limited by the nutrient availability of nitrogen and phosphorus. Instead, nutrient and light conditions control the balance between the release of DOC, DON, and DOP. The release rates E_DON and E_DOP defined in Eqs. (2) and (3) are limited by the respective nutrient limitations β_N and β_P of nitrogen and phosphorus but also by the inverse nutrient limitations 1−β_P and 1.1−β_N, such that the release is highest in the absence of the other nutrient source. Consequently, when primary production is limited only by one nutrient source, extracellular release increases for the less limiting nutrient source. In the absence of both nutrient sources, the extracellular release shifts to release only DOC at a rate E_DOC, limited by both inverse nutrient limitations, as defined in Eq. (4). Any release of DON and DOP also releases DOC, such that the release of DOM is always carbon-rich relative to elemental Redfield ratios. As it is inherent to the implementation, all limiting terms and maximum uptake rates have been adapted to align with the existing terms for primary production in ECOSMO II.

$$\begin{array}{}\text{(1)}& {\displaystyle }E=B_{\mathrm{ER}}\times {\mathrm{\Sigma }}_{j=\mathrm{1}}^{\mathrm{3}}\left({\mathit{\sigma }}_j{C}_{P_j}{\mathit{\beta }}_{Tj}\right)\text{(2)}& {\displaystyle }{E}_{\mathrm{DON}}=\left\{\begin{array}{l}E\times min\left(\mathrm{1}-{\mathit{\beta }}_{\mathrm{P}},{\mathit{\beta }}_{\mathrm{N}},\mathit{\alpha }\left(I\right)\right),\\ {\mathit{\beta }}_N>\mathrm{0.1}\mathrm{and}{\mathit{\beta }}_{\mathrm{P}}<\mathrm{1}\\ \mathrm{0},{\mathit{\beta }}_N\le \mathrm{0.1}\mathrm{or}{\mathit{\beta }}_{\mathrm{P}}\ge \mathrm{1}\end{array}\right.\text{(3)}& {\displaystyle }{E}_{\mathrm{DOP}}=\left\{\begin{array}{l}E\times min\left({\mathit{\beta }}_{\mathrm{P}},\mathrm{1.1}-{\mathit{\beta }}_{\mathrm{N}},\mathit{\alpha }\left(I\right)\right),\\ {\mathit{\beta }}_{\mathrm{P}}>\mathrm{0.1}\mathrm{and}{\mathit{\beta }}_{\mathrm{N}}<\mathrm{1.1}\\ \mathrm{0},{\mathit{\beta }}_{\mathrm{P}}\le \mathrm{0.1}\mathrm{or}{\mathit{\beta }}_{\mathrm{N}}\ge \mathrm{1.1}\end{array}\right.\text{(4)}& {\displaystyle }\begin{array}{rl}{E}_{\mathrm{DOC}}& =E_{\mathrm{DON}}+E_{\mathrm{DOP}}\\ & +\left\{\begin{array}{l}E\times min\left(\max (\mathrm{1}-{\mathit{\beta }}_{\mathrm{P}},\mathrm{1.1}-{\mathit{\beta }}_{\mathrm{N}}),\mathit{\alpha }\left(I\right)\right),\\ {\mathit{\beta }}_{\mathrm{P}}<\mathrm{1}\mathrm{and}{\mathit{\beta }}_{\mathrm{N}}<\mathrm{1.1}\\ \mathrm{0},{\mathit{\beta }}_{\mathrm{P}}\ge \mathrm{1}\mathrm{or}{\mathit{\beta }}_{\mathrm{N}}\ge \mathrm{1.1}\end{array}\right.\end{array}\end{array}$$

To account for the formation of transparent exopolymer particles (Engel, 2002; Mari et al., 2017) from the released DOM, as proposed by Neumann et al. (2022), we include a flocculation rate from DOM to POM. The particle formation is thought to linearly depend on the resulting DIC uptake (Engel, 2002; Neumann et al., 2022), for which we adapt the formulation to be a function of the same limitations as the extracellular release itself. The base rate of flocculation F in Eq. (5) is defined by a constant rate F_DOM2POM and the scaling factor B_ER, both of which are specified in the configuration. The resulting fluxes F_DON2PON, F_DOP2POP, and F_DOC2POC as defined in Eqs. (6)–(8) are depended upon the base rate F and the respective tracer concentrations C_DON, C_DOP, and C_DOC, as well as the same limiting factors as those that apply to the extracellular release itself. To account for the increased carbon contents in observed POM concentrations, we add the flocculation rate for nitrogen and phosphorus to the rate F_DOC2POC for carbon. Consequently, the POC formation is the largest under any conditions, as is the release of DOC itself.

$$\begin{array}{}\text{(5)}& {\displaystyle }F=F_{\mathrm{DOM}\mathrm{2}\mathrm{POM}}\times B_{\mathrm{ER}}\text{(6)}& {\displaystyle }{F}_{\mathrm{DON}\mathrm{2}\mathrm{PON}}=\left\{\begin{array}{l}F\times {\mathrm{C}}_{\mathrm{DON}}\times min\left(\mathrm{1}-{\mathit{\beta }}_{\mathrm{P}},{\mathit{\beta }}_{\mathrm{N}},\mathit{\alpha }\left(I\right)\right),\\ {\mathit{\beta }}_{\mathrm{N}}>\mathrm{0.1}\mathrm{and}{\mathit{\beta }}_{\mathrm{P}}<\mathrm{1}\\ \mathrm{0},{\mathit{\beta }}_{\mathrm{N}}\le \mathrm{0.1}\mathrm{or}{\mathit{\beta }}_{\mathrm{P}}\ge \mathrm{1}\end{array}\right.\text{(7)}& {\displaystyle }{F}_{\mathrm{DOP}\mathrm{2}\mathrm{POP}}=\left\{\begin{array}{l}F\times {\mathrm{C}}_{\mathrm{DOP}}\times min\left({\mathit{\beta }}_{\mathrm{P}},\mathrm{1.1}-{\mathit{\beta }}_{\mathrm{N}},\mathit{\alpha }\left(I\right)\right),\\ {\mathit{\beta }}_{\mathrm{P}}>\mathrm{0.1}\mathrm{and}{\mathit{\beta }}_{\mathrm{N}}<\mathrm{1.1}\\ \mathrm{0},{\mathit{\beta }}_{\mathrm{P}}\le \mathrm{0.1}\mathrm{or}{\mathit{\beta }}_{\mathrm{N}}\ge \mathrm{1.1}\end{array}\right.\text{(8)}& {\displaystyle }\begin{array}{rl}{F}_{\mathrm{DOC}\mathrm{2}\mathrm{POC}}& =F_{\mathrm{DON}\mathrm{2}\mathrm{PON}}+F_{\mathrm{DOP}\mathrm{2}\mathrm{POP}}\\ & +\left\{\begin{array}{l}F\times {\mathrm{C}}_{\mathrm{DOC}}\\ \times min\left(\max (\mathrm{1}-{\mathit{\beta }}_{\mathrm{P}},\mathrm{1.1}-{\mathit{\beta }}_{\mathrm{N}}),\mathit{\alpha }\left(I\right)\right),\\ {\mathit{\beta }}_{\mathrm{P}}<\mathrm{1}\mathrm{and}{\mathit{\beta }}_{\mathrm{N}}<\mathrm{1.1}\\ \mathrm{0},{\mathit{\beta }}_{\mathrm{P}}\ge \mathrm{1}\mathrm{or}{\mathit{\beta }}_{\mathrm{N}}\ge \mathrm{1.1}\end{array}\right.\end{array}\end{array}$$

Under strongly nutrient-limited conditions, β_N≤0.1 or β_P≤0.1, the extracellular release and flocculation for DON or DOP are reduced to zero to prevent an excessive depletion of nutrients. Should 1−β_P or 1.1−β_N fall below zero, the respective terms are also set to zero, as the formulation would otherwise permit negative production. In practice, this state is, however, not reached since the limitations stay within these bounds.

As a second mechanism for variable OM stoichiometry, we incorporate the preferential remineralization of organic nitrogen and phosphorus expressed through higher relative remineralization rates. The remineralization rates ϵ for all of the OM variables are defined in Eqs. (9)–(14). All rates depend on the same temperature function, here expressed in Eq. (9) for POC with the constant reference temperature T_ref and the simulated temperature T. The remineralization rate for DOC is defined relative to POC at a constant ratio, ϵ_DOC:POC. Similarly, the ratios between the percentage rates of PON, POP, DON, and DOP compared to POC and DOC are defined by the constant ratios ϵ_POC:PON, ϵ_POC:POP, ϵ_DOC:DON, and ϵ_DOC:DOP, all of which are specified in the configuration. In the model configurations for this study, we selected values between zero and 1 to reflect the greater bioavailability of nitrogen and phosphorus relative to carbon. The lower the ratio, the higher the respective bioavailability. This allows for the approximation of the observed mean DOM and POM stoichiometry through a faster recycling of nutrients compared to carbon.

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathit{ϵ}}_{\mathrm{POC}}\left(T\right)=\mathrm{0.006}\times \left(\mathrm{1}+\mathrm{20}\times \left({\displaystyle \frac{T^{\mathrm{2}}}{T_{\mathrm{ref}}^{\mathrm{2}}+T^{\mathrm{2}}}}\right)\right)\text{(10)}& {\displaystyle }{\mathit{ϵ}}_{\mathrm{PON}}\left(T\right)={\displaystyle \frac{{\mathit{ϵ}}_{\mathrm{POC}}\left(T\right)}{{\mathit{ϵ}}_{\mathrm{POC}:\mathrm{PON}}}}\text{(11)}& {\displaystyle }{\mathit{ϵ}}_{\mathrm{POP}}\left(T\right)={\displaystyle \frac{{\mathit{ϵ}}_{\mathrm{POC}}\left(T\right)}{{\mathit{ϵ}}_{\mathrm{POC}:\mathrm{POP}}}}\text{(12)}& {\displaystyle }{\mathit{ϵ}}_{\mathrm{DOC}}\left(T\right)={\mathit{ϵ}}_{\mathrm{DOC}:\mathrm{POC}}\times {\mathit{ϵ}}_{\mathrm{POC}}\left(T\right)\text{(13)}& {\displaystyle }{\mathit{ϵ}}_{\mathrm{DON}}\left(T\right)={\displaystyle \frac{{\mathit{ϵ}}_{\mathrm{DOC}}\left(T\right)}{{\mathit{ϵ}}_{\mathrm{DOC}:\mathrm{DON}}}}\text{(14)}& {\displaystyle }{\mathit{ϵ}}_{\mathrm{DOP}}\left(T\right)={\displaystyle \frac{{\mathit{ϵ}}_{\mathrm{DOC}}\left(T\right)}{{\mathit{ϵ}}_{\mathrm{DOC}:\mathrm{DOP}}}}\end{array}$$

A complementary schematic of the OM cycling in the variable stoichiometry version of ECOSMO II is shown in Fig. 3, which illustrates how these two additional parameterizations are incorporated into the existing model structure. The full model description is provided in the Supplement (Sect. S2). This model version includes the prognostic equations for dissolved inorganic carbon (DIC) and total alkalinity (TA) for the coupling to the carbonate system model described below. These were previously introduced in Kossack et al. (2024).

Figure 3 Schematic diagram showing the organic matter (OM) cycling in the variable stoichiometry ECOSMO II including the preferential remineralization of organic nitrogen and phosphorous, as well as the extracellular release of carbon-enriched DOM with the associated particle formation.

2.1.4 Carbonate system model

The carbonate system model by Blackford and Gilbert (2007) describes marine carbonate chemistry based on four main state variables: DIC, TA, pH, and pCO₂. Additionally, it provides the speciation of DIC into carbonic acid (H₂CO₃), bicarbonate (HCO${}_{\mathrm{3}}^{-}$), and carbonate ions (CO${}_{\mathrm{3}}^{\mathrm{2}-}$). Given two of the main state variables, the model estimates the remaining variables iteratively using the HALTAFALL algorithm (Ingri et al., 1967), a method for calculating the composition of equilibrium mixtures. Here, ECOSMO II prognostically simulates DIC and TA, thereby enabling the carbonate system model to diagnostically calculate pH, pCO₂, and the carbonate speciation based on temperature, pressure, and salinity, as in Kossack et al. (2024). For the carbonate speciation, the necessary dissociation constants are based on Millero et al. (2006).

In ECOSMO II, DIC is released into the water column by two processes: pelagic and benthic remineralization of organic carbon and zooplankton excretion. Conversely, carbon fixation by net primary production and the extracellular release of DOC consume DIC. At the ocean surface, the exchange of CO₂ with the atmosphere results in the addition or removal of DIC from the water column. The air–sea CO₂ exchange depends on the difference between the simulated surface pCO₂ and monthly prescribed atmospheric CO₂ concentrations from the Mace Head station (Lan et al., 2023) and on wind speed as described by the gas transfer parametrization in Wanninkhof (2014). Changes in TA are caused by biological uptake and release of nitrate, ammonium, and phosphate, as described in the prognostic TA formulation by Wolf-Gladrow et al. (2007) and Gustafsson (2013). A complete description of the DIC and TA equations can also be found in the full model description in the Supplement (Sect. S2). A number of previous studies have already examined carbon dynamics in the North Sea and the NWES using versions of this carbonate system model in combination with different physical and ecosystem models (Artioli et al., 2012, 2014; Blackford et al., 2017; Blackford and Gilbert, 2007; Wakelin et al., 2012).

2.1.5 Model configuration and experiments

To assess the individual and combined effects of the release of carbon-enriched DOM and the preferential remineralization of organic nitrogen and phosphorus, we compare four model configurations with identical initial conditions and forcing. The first configuration (RS), which serves as a reference simulation, adheres to Redfield stoichiometry by considering equal remineralization rates for carbon, nitrogen, and phosphorus and by neglecting the extracellular release of DOM. The second configuration (known as ER) considers the DOM release with a scaling factor of B_ER=0.4 and equal remineralization rates. The third configuration (known as PR) includes the preferential remineralization of nitrogen and phosphorus with +60 % and +100 % higher bioavailability for POM and DOM, respectively. Lastly, the fourth configuration (ER&PR) combines the extracellular release and PR, both with a reduced contribution. This is expressed by a lower extracellular release scaling factor, B_ER=0.2, and a reduced increase in bioavailability of +30 % and +50 % for nitrogen and phosphorus, respectively. These configurations have been developed to represent the mean stoichiometry of labile DOM and POM within observational bounds (Sect. 2.2.2).

All relevant parameters for the four configurations are listed in the Supplement (Table S1). The initial conditions and river discharge, open-boundary, and atmospheric forcing conditions are identical to the NWES-LR configuration in Kossack et al. (2023, 2024), except for the surface deposition of ammonium and nitrate. The spatially resolved surface deposition is now prescribed monthly rather than by a climatological mean. The monthly means are derived from EMEP MSC-W data (Simpson et al., 2012). Additional variables for DON, DOP, PON, and POP are calculated from the DOM and POM concentrations assuming Redfield stoichiometry in the initial and forcing conditions. Accordingly, any discrepancies from elemental Redfield ratios can be attributed to the aforementioned parameterizations.

Initial fields and monthly climatological means for the open-boundary conditions of temperature, salinity, nutrients, and oxygen are interpolated from the World Ocean Atlas 2018 (WOA18) data products (Garcia et al., 2018a, b; Locarnini et al., 2018; Zweng et al., 2019). For DIC and TA, the initial and boundary conditions are based on interpolated climatological fields from NNGv2LDEO (Broullón et al., 2019, 2020), a neural network (NN) approach combining the Global Ocean Data Analysis Project version 2 (Gv2) and the Lamont–Doherty Earth Observatory (LDEO) datasets. For the Baltic Sea, the TA initial conditions are instead derived from the salinity relation in Hjalmarsson et al. (2008). To account for globally increasing atmospheric CO₂ concentrations, an interannual trend from global ICON-Coast simulations (Mathis et al., 2022) is added to the DIC open-boundary conditions. Measurements from the Mace Head station (Lan et al., 2023) provide monthly mean atmospheric CO₂ forcing for the air–sea CO₂ exchange calculation. All initial sediment fields are based on ECOSMO simulations conducted by Fabian Werner (personal communication, 2023) and Samuelsen et al. (2022).

River loads of nutrients, DIC, and TA are prescribed along with river discharge for the 172 largest rivers in the model domain, while the remaining loads are set to zero. A regional dataset introduced by Daewel and Schrum (2013) and extended by Zhao et al. (2019) provides daily nutrient loads from a 21-year mean over 1995–2015. Annual mean concentrations from Pätsch and Lenhart (2004) and Neal and Davies (2003) additionally supply DIC, DOC, and TA river loads for the Scheldt, Meuse, Rhine, Ems, Elbe, Humber, Wear, Tweed, Great Ouse, and Thames. Where possible, missing TA loads are derived from the salinity relation in Hjalmarsson et al. (2008) and Artioli et al. (2012). For the remaining rivers, we use average DIC loads of 2700 µM from Pätsch and Lenhart (2004) and an average TA in the NWES of 2050 µM. Missing DOC loads are estimated to be 1/10 of the corresponding DIC concentration. Riverine DOC is assumed to be only 10 % bioavailable, as suggested by Kühn et al. (2010). This fraction is considered to be remineralized and hence is directly added to DIC.

2.2 Organic matter stoichiometry: representation and validation

2.2.1 Observed organic matter stoichiometry

The spatial and temporal variability of $\mathrm{C}:\mathrm{N}:\mathrm{P}$ stoichiometry in DOM and POM remains poorly constrained. However, case studies and global observations permit the estimation of their range of variations. Due to scarcity of measurements within the NWES, we rely on a mixture of global and regional observations. For reference, a literature compilation of both concentrations and the stoichiometry of OM in the global open ocean and various marginal seas is provided in the Supplement (Tables S3–S10). These observations indicate globally consistent compositional differences between DOM and POM, as well as between different lability fractions of DOM (Hopkinson and Vallino, 2005; Liang et al., 2023). Furthermore, some latitudinal and vertical patterns and predictors have previously been identified, including temperature and nutrient availability (Liang et al., 2023; Martiny et al., 2013; Tanioka et al., 2022b). Ecosystem stoichiometry may also change under shifting environmental conditions such as temperature, stratification, and CO₂ concentrations (Finkel et al., 2010; Moreno and Martiny, 2018; Sardans et al., 2012, 2021), which could contribute to future trends in DOM and POM composition.

In the surface layer of the global open ocean, the observed average stoichiometry of bulk DOM exhibits ranges of $\mathrm{DOC}:\mathrm{DON}:\mathrm{DOP}=\mathrm{387}$–1164:26–123:1 (Hopkinson and Vallino, 2005; Letscher and Moore, 2015; Liang et al., 2023; Lønborg and Álvarez-Salgado, 2012), which is more carbon-rich than POM, with an estimated average composition of $\mathrm{DOC}:\mathrm{DON}:\mathrm{DOP}=\mathrm{146}$–294:21–36:1 (Liang et al., 2023; Martiny et al., 2014; Singh et al., 2015; Tanioka et al., 2021, 2022a). For both DOM and POM, there is a persistent vertical trend that indicates increasing C:N, C:P, and N:P ratios with depth. This relation is thought to arise as organic nitrogen and phosphorus is preferentially respired, leaving older and deeper material more carbon-rich (Hopkinson and Vallino, 2005; Letscher and Moore, 2015; Painter et al., 2017; Tanioka et al., 2021).

In the case of DOM, this vertical gradient is more pronounced, which is linked to old and strongly nutrient-depleted refractory material. The continuous spectrum of lability in DOM is typically simplified into two or three pools of (semi-)labile (LDOM) and refractory DOM (RDOM) (Hansell et al., 2009). While the labile and semi-labile fractions have a life time of days to years, the refractory pool has an average age of approximately 6000 years (Bauer et al., 1992; Hansell et al., 2009). Over time, the bioavailable compounds are removed, leaving the remaining refractory material highly carbon-rich, with an observed average stoichiometry of $\mathrm{RDOC}:\mathrm{RDON}:\mathrm{RDOP}=\mathrm{1373}$–3511:10–202:1 (Aminot and Kérouel, 2004; Hopkinson and Vallino, 2005; Liang et al., 2023; Lønborg and Álvarez-Salgado, 2012).

Despite its high carbon content, RDOM only acts as a minor source of inorganic carbon and nutrients on sub-millennial timescales due to its slow formation and decomposition (Hopkinson and Vallino, 2005). In contrast, LDOM, with an observed average stoichiometry of $\mathrm{LDOC}:\mathrm{LDON}:\mathrm{LDOP}=\mathrm{179}$–$\mathrm{199}:\mathrm{20}:\mathrm{1}$ (Hopkinson and Vallino, 2005; Liang et al., 2023) in the global surface ocean, acts as the major control on the recycling of carbon and nutrients on shorter timescales. Consequently, the majority of OM cycling on the annual to decadal timescales considered here occurs at an intermediate stoichiometry between the composition of bulk DOM and elemental Redfield ratios (Hopkinson and Vallino, 2005). The balance between concentrations of RDOM and LDOM to a large extent determines the vertical gradient in bulk DOM stoichiometry (Aminot and Kérouel, 2004). High LDOM concentrations are only observed in surface waters, while RDOM concentrations are instead vertically homogeneous (Aminot and Kérouel, 2004; Hansell and Carlson, 1998). Hence, the fraction of RDOM, and with that the relative carbon content in bulk DOM, increase with depth (Aminot and Kérouel, 2004; Hopkinson and Vallino, 2005; Hung et al., 2003; Liang et al., 2023).

In marginal seas such as the North Sea and Baltic Sea, river discharge and water residence times may further control bulk DOM stoichiometry (Voss et al., 2021). Regions of freshwater influence exhibit higher carbon contents, with $\mathrm{DOC}:\mathrm{DON}:\mathrm{DOP}=\mathrm{1199}$–3714:52–221:1 (Stepanauskas et al., 2002; Voss et al., 2021). This is attributed to the elevated C:N ratios of 30–60:1 observed in riverine OM (Bauer et al., 2013). Long residence times in the Baltic Sea facilitate the accumulation of riverine organic matter, in contrast to the relatively short residence times in the North Sea, which reduce the influence of riverine inputs on bulk OM stoichiometry. These factors contribute to the large regional variations, for which we provide an overview in the literature compilation.

2.2.2 Simulated organic matter stoichiometry

To ascertain a realistic magnitude of the preferential remineralization and the release of carbon-enriched DOM, we identified parameter settings for which the simulated stoichiometry adheres to observational bounds for LDOM and POM. Accordingly, we reduced the contribution of each mechanism in the ER&PR configuration to ensure consistency across all configurations. We found that using the parameter settings from the individual configurations in the combined ER&PR configuration resulted in C:N and C:P ratios significantly exceeding those that were observed and in the other configurations, making their comparison less meaningful.

It should be noted that the DOM simulated in ECOSMO is fully biodegradable, with the majority of DOM produced being respired within 1 year. Consequently, the simulated DOM represents a spectrum of only labile and semi-labile DOM but not RDOM. The latter is not currently considered, as its impact on the carbon and nutrient cycles is assumed to be negligible on the timescales in question. Therefore, we assess the simulated DOM concentrations and stoichiometry with respect to observations of LDOM. Given the large discrepancies in $\mathrm{C}:\mathrm{N}:\mathrm{P}$ stoichiometry of bulk DOM, LDOM, and RDOM, we assume this to be the most reasonable approach.

As a representation of the simulated OM stoichiometry, we present the time-averaged horizontal distribution of column-integrated C:N, C:P, and N:P ratios, as well as vertical profiles of horizontally averaged ratios over the northern, central, and southern North Sea. Figures 4 and 5 show these distributions for the ER, PR, and ER&PR configurations for DOM and POM, respectively. Observations do not allow us to resolve regional vertical profiles. For this reason, we compare the simulated stoichiometry to observed ranges and, where available, the observed mean composition, independent of depth. In the NWES, observations of DOM stoichiometry are limited to bulk DOM, without a distinction between LDOM and RDOM (Chaichana et al., 2019; Davis et al., 2014, 2019; Painter et al., 2018; Suratman et al., 2009). Thus, we compare the vertical distribution of DOM stoichiometry to globally estimated ranges for LDOM from Hopkinson and Vallino (2005) and Liang et al. (2023). In contrast, for POM stoichiometry, we provide both global estimates from Liang et al. (2023), and regional estimates from the Hebrides Shelf (Painter et al., 2017) and the Celtic Sea (Davis et al., 2014, 2019). This is because direct observations for the North Sea are limited to C:N ratios and do not provide estimated mean compositions (Chaichana et al., 2019; Suratman et al., 2009).

Overall, the simulated DOM stoichiometry reproduces the observed range within reasonable limits (Fig. 4). The mean simulated DOC:DON ratios of 10.2–11.2 compare well with the observed means and ranges of 8.9 (5.4–12.0) and 10.7 (8.7–14.1). Similarly, simulated mean DOC:DOP ratios of 198–232 agree with the observational estimates of 179 (83–414) and 199 (154–245). In both cases, the ER&PR configuration is on the upper end of the range, with the ER and PR configurations being more central. In terms of the DON:DOP ratios, all configurations reproduce the observed means of 20 (15–49) and 20 (16–25), with mean simulated DON:DOP ratios of 19.5–20.7. For all three elemental ratios, the PR configuration shows the observed increase with depth as N and P are preferentially removed. This is not the case for the ER configuration, where the additional DOM is only released in the euphotic zone. In terms of the horizontal distribution, the carbon enrichment is highest for the ER configuration in the central and northern parts of the North Sea, as well as in the deep parts of the Atlantic Ocean. This shows a higher DOC accumulation in nutrient depleted areas, as expected from the implementation. For the PR configuration, instead, the carbon contents are highest in the Norwegian Trench and along the shelf edge, where OM respiration is generally high, which amplifies the difference in remineralization rates. The average composition is close to the surface DOM composition, as DOM concentrations rapidly decrease with depth below the euphotic zone.

Figure 4 Mean DOM stoichiometry for the ER, PR, and ER&PR configurations over the simulation period of 2000–2010. The first three columns show the vertically integrated horizontal distribution. For each elemental ratio and configuration, the titles display the mean values and the respective standard deviations over the NWES. The last column shows the vertical distribution of simulated horizontal mean stoichiometry and the corresponding standard deviation over the NNS, CNS, and SNS sub-areas. Observed depth-independent ranges in DOM stoichiometry and the estimated mean composition from Hopkinson and Vallino (2005) and Liang et al. (2023) are represented by the error bars. The corresponding Redfield ratios of C:N = 6.625, C:P = 106, and N:P = 16 are indicated by the dashed blue lines in the vertical profiles. The horizontal means and standard deviations over the entire NWES are listed above the horizontal distributions. For the vertical profiles, the horizontal standard deviation is represented by the shaded area around the solid lines.

The simulated POM stoichiometry is also consistent with the observed ranges but is more challenging for the observed mean ratios, as only one global study explicitly provides an average composition (Fig. 5). Simulated mean POC:PON ratios of 7.4–8.1 are consistent with the global mean of 7.7 (6.1–9.2) and the ranges of 8.3–13.7 and 4.3–11.6 for the Hebrides Shelf and the Celtic Sea. Likewise, the simulated POC:POP ratios of 139–158 are within the center of the observed mean and ranges of 160 (73–295), (142–260), and (63–283), with the PR configuration being closest to the estimated mean. The increased PON:POP ratios compared to Redfield stoichiometry are reproduced in all configurations. Simulated PON:POP ratios of 19.0–19.6 are comparable to the observational estimates, with a mean and ranges of 21 (10–37), (16–26), and (10–47). For all three elemental ratios, the PR configuration shows the highest and the ER configuration the lowest horizontal averages across all vertical levels. As for DOM, the increasing ratios with depth indicated by DOM and POM observations are reproduced by the PR but not the ER configuration. The horizontal distribution for the ER configuration shows the highest carbon enrichment in the central and northern North Sea, the same as for DOM. In contrast, the PR configuration instead shows homogeneous elemental ratios across the NWES but high carbon enrichment in the deep, open Atlantic. Altogether, simulated DOM and POM stoichiometry both reasonably represent the available observations.

Figure 5 Mean POM stoichiometry for the ER, PR, and ER&PR configurations over the simulation period of 2000–2010. The first three columns show the vertically integrated horizontal distribution. For each elemental ratio and configuration, the titles display the mean values and the respective standard deviation over the NWES. The last column shows the vertical distribution of simulated horizontal mean stoichiometry and the corresponding standard deviation over the NNS, CNS, and SNS sub-areas. Observed depth-independent ranges in POM stoichiometry and the estimated mean composition from Liang et al. (2023), Davis et al. (2014, 2019), and Painter et al. (2017) are represented by the error bars. The corresponding Redfield ratios of C:N = 6.625, C:P = 106, and N:P = 16 are indicated by the dashed blue lines in the vertical profiles. The horizontal means and standard deviations over the entire NWES are listed above the horizontal distributions. For the vertical profiles, the horizontal standard deviation is represented by the shaded area around the solid lines.

2.3 Impact on model performance

The ability of SCHISM-ECOSMO-CO₂ to reproduce the observed variability in temperature, nutrients, and carbonate system variables in the NWES has already been established in Kossack et al. (2023, 2024). As the extended OM representation in the variable stoichiometry ECOSMO II version may have consequences for the model validation, we here describe the effects of the model extension on the model validation in the ER, PR, and ER&PR configurations compared to the RS reference configuration. With the primary objective being to investigate the consequences of the proposed modifications on marine carbon cycling, particularly with regard to the biological uptake of inorganic carbon, we provide a more comprehensive account of the carbonate system validation. Differences in the model representation compared to Kossack et al. (2023, 2024) may also be linked to the different data availability in the simulated time period, to the monthly resolved surface deposition, and to a reduced remineralization rate for DOM. The latter was adjusted for a better representation of DOC concentrations, which are essential for assessing the effects of variable stoichiometry. This is associated with a decrease in nutrient availability and hence productivity.

Based on surface, bottle, and pump data from the ocean hydrochemistry data collection of the International Council for the Exploration of the Sea (ICES), we find that there is no significant improvement or deterioration in the model representation of phosphate and nitrate for any of the three variable stoichiometry configurations with respect to the reference configuration. For these variables we provide the horizontal distribution of seasonal model biases at co-located data points and corresponding Taylor diagrams in the Supplement (Figs. S1–S6). The spatial distribution of regional and seasonal biases is persistent across all configurations, which results in overall negligible differences in the resulting model performance. We include temperature and salinity for the RS configuration, as it may have a relevant influence on the carbonate system calculations. All other configurations are excluded for the physical variables, as the physical simulations are equivalent.

In the validation for the carbonate system, we assess DIC, TA, and pCO₂. While we do not discuss TA in the results of this study, it is used to derive the simulated pCO₂ and, as such, is relevant to the model performance. We compare the simulated DIC and TA concentrations to a total of 5685 and 5110 observations from the Global Ocean Data Analysis Project (GLODAP). Additionally, we validate the simulated CO₂ concentrations against 772 202 CO₂ fugacity (fCO₂) measurements from the Surface Ocean CO₂ Atlas (SOCAT). For this, the observed fCO₂ values are converted to pCO₂ using the Python library PyCO2SYS (Humphreys et al., 2022, 2024). Simulated surface CO₂ concentrations are converted to pCO₂ according to Eqs. (15)–(16) from the guide to best practices by Dickson et al. (2007) and an updated water vapor pressure equation as a function of temperature from Huang (2018). Here, xCO₂ is the simulated dry mole fraction, p_air and $p_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ are the air and water vapor pressure, and T is the simulated temperature.

$$\begin{array}{}\text{(15)}& {\displaystyle }p{\mathrm{CO}}_{\mathrm{2}}=x{\mathrm{CO}}_{\mathrm{2}}\times (p_{\mathrm{air}}-p_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}})\text{(16)}& {\displaystyle }\begin{array}{rl}{p}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}& =\mathrm{9.869}\times \exp \left(\mathrm{34.494}-{\displaystyle \frac{\mathrm{4924.99}}{T+\mathrm{237.1}}}\right)\\ & \times {\left(T+\mathrm{105}\right)}^{-\mathrm{1.57}}\end{array}\end{array}$$

In stark contrast to the nutrient validation, there is a significant discrepancy between the model configurations with respect to DIC and pCO₂ related to changes in the inorganic carbon uptake. TA validation, however, does not differ significantly between configurations. This is demonstrated by the changes in the corresponding Taylor diagrams in Fig. 6. The Taylor diagrams represent the model performance in terms of the correlation, normalized standard deviation (NSD) with respect to the observations, and their normalized and centered root-mean-squared difference (RMSD). As complementary indicators of the model performance, we provide the seasonal co-located model biases and the mean percentage biases in the Supplement (Figs. S7–S9 and Table 2).

With regard to DIC, the RS configuration underestimates the observed variability in the NNS, CNS, SNS, and NT, with an NSD of 0.48–0.80 (Fig. 6). Compared to this, the configurations allowing for variable OM stoichiometry show a consistent improvement in the represented variability, with an NSD of 0.70–0.98, 0.55–0.85, and 0.63–0.90 for the ER, PR, and ER&PR configurations, respectively (Fig. 6). With an average increase in NSD of 0.07–0.21, these more closely reproduce the range of observed variability across all sub-regions. The improvement is consistently largest for the ER and smallest for the PR configuration. This range of improvement is proportional to the additional DIC uptake from biological carbon fixation, which is discussed below. Only in the SK sub-region does the NSD exceed 1 in the RS configuration, showing an increase from 1.17 to 1.21–1.32, which adds to the overestimation of the observed variability. Differences in the correlation between the configurations are negligible, with correlation coefficients of 0.81–0.88 for the NNS, CNS, and NT sub-regions (Fig. 6). For the SK and SNS, the correlation is weaker, with correlation coefficients of 0.71–0.72 and 0.60–0.66, respectively. Improvements in the correlation may be limited due to an expected mismatch between the small-scale variability captured in measurements and simulated average concentrations on much larger scales. The model's percentage bias in DIC is generally between 0.6 % and 1.3 % and does not show any conclusive differences between configurations (Table S2).

Differences in the model performance for TA are negligible and inconsistent across sub-regions, ranging from 0.002 to 0.013 in the NSD and 0.001 to 0.003 in the correlation coefficients (Fig. 6). The correlation is generally weaker for TA than for DIC, which can be attributed to both missing constituents in ECOSMO II that contribute to changes in TA (Wolf-Gladrow et al., 2007) and to uncertainties in terrestrial TA sources. In the SNS and NT sub-regions, the correlation coefficients of 0.60 and 0.82 are on the higher end. In comparison, the NNS and CNS have correlation coefficients of 0.29 and 0.40, which are more influenced by exchange with the Atlantic and less by river discharge. This may suggest that not only terrestrial sources contribute to the lower correlation but local alkalinity generation also contributes.

As pCO₂ is largely dependent on DIC, improvements in the representation of DIC are expected to be reflected, to some extent, in the representation of pCO₂. However, given the considerably greater data availability for pCO₂, we provide a validation here explicitly. For the simulated surface pCO₂, as for DIC, there is an improvement in the representation of the observed variability in most of the sub-regions. Additionally, several sub-regions show a significant increase in correlation in combination with a small but consistent decrease of approximately 0.2 %–3.3 % (Table S2) in the percentage bias across all sub-regions except from the SK and SWC regions. The correlation coefficients for the NNS, CNS, SNS, NT, SK, EC, and AS regions increase significantly from 0.66, 0.60, 0.36, 0.21, 0.34, 0.43, and −0.06 in the RS configuration to 0.78, 0.76, 0.53, 0.42, 0.44, 0.59, and 0.39 in the ER configuration (Fig. 6). This increase in correlation is similarly dependent on the increase in inorganic carbon uptake, the same as for the changes in standard deviation for DIC. The PR and ER&PR configurations yield intermediate results, with the lowest increase in correlation for the PR configuration. Only for the SWC does the correlation remain rather consistent. Likewise, the observed standard deviation is more closely represented in the NNS, CNS, SNS, EC, SWC, and AS sub-regions by variable OM stoichiometry. The corresponding NSD increases from 0.67, 0.81, 0.70, 0.65, 0.75, and 0.37 in the RS configuration up to 1.10, 1.08, 0.87, 0.74, 1.14, and 0.48 in the ER configuration (Fig. 6). Variability is generally overestimated in the NT and SK sub-regions, likely due to the representation of rivers as point sources rather than explicitly resolved estuaries.

In summary, the variable stoichiometry configurations consistently improve model performance in terms of the represented range of variability for DIC and pCO₂, represented by the NSD, as well as the correlation with observations for pCO₂. The magnitude of the improvement is lowest in the PR and highest in the ER configuration, which is related to the corresponding increase in biological uptake of inorganic carbon, as discussed below. This supports the need to account for carbon overconsumption compared to estimates based on Redfield stoichiometry. While variable OM stoichiometry here improves the representation in terms of the range of variability, some challenges remain in representing the timing of DIC release in autumn. This is further discussed below with respect to surface DIC concentrations. Overall, the changes in TA are negligible, indicating that the implementation has a minimal impact on TA.

Figure 6 Taylor diagrams for simulated dissolved inorganic carbon (DIC), total alkalinity (TA), and partial pressure of CO₂ (pCO₂) in comparison to the GLODAPv2022 and SOCATv2023 datasets, including all four model configurations. The model performance is represented in terms of correlation, normalized standard deviation, and normalized and centered root mean square difference (RMSD), with respect to the observations. Regions were selected based on sufficient data availability. The results are differentiated by configuration, as indicated by color, and sub-region, as indicated by shape.

2.4 Net community production

Net community production (NCP), defined as the difference between carbon fixation and respiration, is a metric used to assess whether a system is net autotrophic (NCP >0), with a net production of organic carbon, or net heterotrophic (NCP <0), with a net respiration of organic carbon (Smith and Hollibaugh, 1993). Net autotrophic regions typically act as a sink for inorganic carbon, while net heterotrophic regions increase the inorganic carbon content respired from additional imported organic carbon (Thomas et al., 2005). For this study we define NCP as carbon fixation through primary production and the extracellular release of DOC minus carbon respiration, which consists of pelagic and benthic remineralization as well as zooplankton excretion contributing to DIC release (Fig. 3).

3 Results

3.1 Annual carbon fluxes

As a first step towards quantifying the large-scale effects of both pathways for variable OM stoichiometry on OM cycling and the resulting changes in the air–sea CO₂ exchange, we assess their influence on key annual carbon fluxes and their seasonal variability. In particular, we examine the effects on vertically integrated carbon fixation, carbon respiration, NCP, and the air–sea CO₂ exchange. We begin by assessing changes in the balance between organic carbon production and respiration to understand how the associated uptake and release of inorganic carbon contributes to resulting changes in the oceanic CO₂ uptake.

3.1.1 Biological carbon fixation

Biological carbon fixation exhibits a consistent percentage increase between the North Sea and the entire NWES for each of the variable stoichiometry configurations (Fig. 7). However, the magnitude of this increase varies considerably, with the PR configuration showing the smallest and the ER configuration showing the largest increase. In comparison to the vertically integrated carbon fixation in the North Sea of 4.5 Tmol C yr⁻¹ for the RS configuration, the PR, ER&PR, and ER configurations exhibit a notable increase, reaching 4.9 (+10 %), 6.0 (+34 %), and 6.9 Tmol C yr⁻¹ (+55 %), respectively. For the entire NWES, the corresponding reference carbon fixation of 9.7 Tmol C yr⁻¹ increases to 10.6 (+9 %), 12.8 (+32 %), and 14.8 Tmol C yr⁻¹ (+52 %).

Figure 7 Simulated annual vertically integrated carbon fixation, pelagic and benthic remineralization, net community production (NCP), and net air–sea CO₂ exchange for all four model configurations, averaged over the simulated period of 2000–2010. The air–sea CO₂ flux is defined downward such that positive values indicate an oceanic uptake and negative values an outgassing of CO₂. The error bars show the corresponding standard deviation. White bars show the literature estimates and ranges for NCP and the air–sea CO₂ exchange in the North Sea from Thomas et al. (2005) and the estimated range of the air–sea CO₂ exchange in the NWES from Legge et al. (2020), with reference to both modeling and observational studies. The annual time series for selected sub-regions and the average fluxes for all sub-regions are provided in the Supplement (Figs. S10–S14 and Tables S11–S15).

In the RS reference configuration, carbon fixation is highest in the shallow and well-mixed shelf regions (Fig. 8a), including the SNS, EC, AS, CS, IS, ISS, and NWA sub-areas (Fig. 2). The highest annual sub-region mean is reached in the Celtic Sea, with a mean value of 9.1 mol C m⁻² yr⁻¹. In contrast, the lowest annual mean of carbon fixation is found in the deeper central North Sea, with a value of 7.2 mol C m⁻² yr⁻¹. This is due to the seasonal stratification, which promotes nutrient limitation. The spatial distribution of simulated carbon fixation is consistent with the established physical controls on primary production and its characteristic spatial pattern in the NWES (Holt et al., 2012).

Across the North Sea area, the carbon fixation response differs considerably between the PR and ER configurations. For the ER configuration, the largest increase in carbon fixation here is observed in the CNS, NNS, and NT sub-regions, with values ranging from 3.9 to 4.7 mol C m⁻² yr⁻¹ (Fig. 8a). In contrast, the SNS sub-region exhibits a smaller increase, with a value of 3.2 mol C m⁻² yr⁻¹. Conversely, the PR configuration exhibits the highest increase of 0.9 mol C m⁻² yr⁻¹ in the SNS and smaller increases of 0.4–0.7 mol C m⁻² yr⁻¹ in the CNS, NNS, and NT sub-regions. Consequently, the north–south gradient in productivity is weakened in the ER configuration but strengthened in the PR configuration. There is an evenly distributed positive response in the outer shelf regions, while some inner shelf regions, including the EC, IS, and ISS, show a much weaker response. All configurations further show a large increase in carbon fixation in the AS region, from 1.2 to 5.6 mol C m⁻² yr⁻¹. Consistent with the mean annual changes, the amplitude of these changes is largest in the ER and smallest in the PR configuration.

Figure 8 Spatial distribution of annual mean vertically integrated (a) carbon fixation and (b) respiration, including pelagic and benthic remineralization, for the RS configuration and differences for the ER, PR, and ER&PR configurations over the simulation period of 2000–2010. The respective seasonal means for carbon fixation and respiration are provided in the Supplement (Figs. S15 and S16).

3.1.2 Carbon respiration

The rate of carbon respiration exceeds the rate of carbon fixation in the North Sea and the entire NWES in all configurations (Fig. 7). In the RS configuration, the magnitude of pelagic remineralization is comparable to that of carbon fixation, with annual means of 4.5 Tmol C yr⁻¹ in the North Sea and 9.6 Tmol C yr⁻¹ across the entire NWES. The additional benthic remineralization of 0.37 and 0.70 Tmol C yr⁻¹, however, results in a total respiration that surpasses carbon fixation. In comparison, the variable stoichiometry configurations exhibit higher pelagic remineralization, reaching values of 5.0 (+10 %), 6.1 (+35 %), and 7.1 Tmol C yr⁻¹ (+56 %) for the North Sea in the PR, ER&PR, and ER configurations, respectively. In the entire the NWES, the corresponding values are 10.5 (+9 %), 12.8 (+33 %), and 14.9 Tmol C yr⁻¹ (+55 %). In contrast, the increases in benthic remineralization are equal across configurations and thus do not match this percentage increase. The annual means of benthic remineralization are 0.44 (+20 %), 0.43 (+17 %), and 0.45 Tmol C yr⁻¹ (+21 %) in the North Sea and 0.83 (+18 %), 0.81 (+15 %), and 0.85 Tmol C yr⁻¹ (+20 %) in the NWES. This is likely linked to a comparable increase in POC concentrations, which is further discussed below.

The pattern of the RS reference carbon respiration (Fig. 8b) resembles the distribution of carbon fixation, with the exception of the NT region. Similar to carbon fixation, respiration is particularly low in the CNS, 7.5 mol C m⁻² yr⁻¹, and particularly high in the western outer shelf regions, up to 9.3 mol C m⁻² yr⁻¹ in the SWC. However, respiration reaches its maximum in the NT, with an average rate of 11.1 mol C m⁻² yr⁻¹ despite only moderate carbon fixation. The north–south gradient in the impact of carbon respiration across the North Sea is positive in the PR configuration and negative in the ER and ER&PR configurations. The preferential remineralization thus increases the north–south gradient in production and respiration. Conversely, the extracellular release of DOC promotes production and respiration in the otherwise nutrient-depleted and less biologically productive areas. The strong respiration in the NT is further enhanced by 16.1 mol C m⁻² yr⁻¹ in the variable stoichiometry configurations. Otherwise, the differences are consistent with the increases in carbon fixation.

Shifts in the balance between carbon fixation and respiration rather than either process alone play a crucial role in air–sea CO₂ exchange, influencing whether specific locations act as sources or sinks of DIC. Therefore, we will further assess the annual NCP, focusing on both vertically integrated values and the vertical gradient between the surface and sub-surface layers.

3.1.3 NCP

Overall, the North Sea and the entire NWES are net heterotrophic in all configurations, with a corresponding NCP of −0.43 and −0.63 Tmol C yr⁻¹ in the RS configuration (Fig. 7). For the North Sea, this value is situated at the lower end of the uncertainty range from observations, with an estimated NCP of $-\mathrm{0.34}\pm \mathrm{0.11}$ Tmol C yr⁻¹ reported in Thomas et al. (2005). The PR, ER&PR, and ER configurations exhibit a stronger net heterotrophy (i.e., lower NCP), with NCP values of −0.49 (−15 %), −0.54 (−27 %), and −0.59 Tmol C yr⁻¹ (−40 %) in the North Sea. In the NWES, the percentage decrease in NCP is even more pronounced, with an annual NCP of −0.75 (−18 %), −0.86 (−36 %), and −0.96 Tmol C yr⁻¹ (−52 %).

In the RS configuration, net autotrophic regions are found in the shallow coastal areas of the SNS, EC, CS, IS, and NWA sub-regions; the shelf edge; and the Dogger Bank (Fig. 9a). The most pronounced net autotrophic region is a broad band in the North Sea along the Norwegian Trench. Conversely, beyond the 200 m isobath, the NT and the deep parts of the SK region are the most net heterotrophic regions across the entire shelf. This results in a large NCP gradient across the slope. Similarly, portions of the shelf slope and the deeper central IS, ISS, and NWA sub-regions exhibit comparatively elevated net heterotrophy, with net autotrophy observed in adjacent, shallower areas. The remaining shelf area is largely weakly net heterotrophic or nearly neutral.

The most notable NCP response is the intensification of the aforementioned gradient between net autotrophy along the edge of the NT region and the net heterotrophy inside the NT region. This feature indicates an increase in OM that is produced in the CNS, transported into the NT, and subsequently remineralized. The variable stoichiometry configurations demonstrate a reinforcement of the pre-existing pattern of NCP across the entire shelf area. This is expressed by an increase in NCP in net autotrophic and a decrease in NCP in net heterotrophic regions. The RS reference NCP is lowest in the NT and SK regions, with annual means of −3.2 and −1.0 mol C m⁻² yr⁻¹, respectively. These regions also exhibit the most pronounced NCP response, with decreases of −0.3 to −1.6 mol C m⁻² yr⁻¹ in the NT and −0.3 to −1.0 mol C m⁻² yr⁻¹ in the SK sub-region. In accordance with the response in carbon fixation and respiration, the PR configuration overall shows the weakest decrease, while the ER configuration shows the strongest decrease. Nevertheless, even in the PR configuration, the largest decrease in NCP is shown in the NT and SK regions, thereby amplifying the pre-existing NCP gradient. In the SNS, IS, and ISS, the impact is intermediate, with values ranging from −0.1 to −0.2 mol C m⁻² yr⁻¹, and it is less significant in the remaining regions. The ER configuration shows an increase in the CNS and AS, with a range of 3–0.05 mol C m⁻² yr⁻¹. A significant decrease can be seen in the NNS, SNS, SWC, IS, ISS, HS, and MS, with a range of −0.2 to −0.4 mol C m⁻² yr⁻¹.

Figure 9 (a) Spatial distribution of annual vertically integrated NCP for the RS configuration and differences for the ER, PR, and ER&PR configurations over the simulation period of 2000–2010. (b) Respective monthly 11-year mean of carbon fixation (CF, dashed); carbon respiration (CR, dotted), including benthic and pelagic remineralization; and NCP (dash–dot). The respective seasonal means of NCP are provided in the Supplement (Fig. S17).

While the increase in carbon fixation is limited to the euphotic zone, where sufficient light and nutrients are available, respiration increases in both the surface and sub-surface layers (Fig. 10). In the upper 50 m, pelagic remineralization is approximately 10 %–61 % higher in the variable stoichiometry configurations. Below 50 m, the remaining water column shows a smaller increase in pelagic remineralization of approximately 5 %–41 %. The separation into carbon fixation at the surface only and respiration across various depths leads to opposite changes in the pelagic NCP (i.e., excluding benthic remineralization) between the surface and sub-surface layers. The net autotrophic surface layer above 50 m becomes increasingly autotrophic, by 7 %–34 %, while the remaining water column becomes more net heterotrophic, by 6 %–41 %. This enhanced vertical gradient in NCP suggests a stronger depletion of DIC in surface waters and a stronger release of DIC at depth, implying a stronger biological control on surface pCO₂ and an increase in the resulting oceanic CO₂ uptake.

Figure 10 Simulated annual vertically integrated pelagic remineralization and pelagic NCP, separated into above and below 50 m depth and averaged over the simulated period of 2000–2010 for all four model configurations. The error bars show the corresponding standard deviation. In contrast to the NCP shown with all vertically integrated fluxes in Fig. 9, pelagic NCP here does not include benthic remineralization.

3.1.4 Enhanced seasonal organic matter cycling

The impacts of variable stoichiometry on the annual means result from marked changes in the seasonality of carbon fixation and respiration. Across all configurations, carbon fixation in the North Sea and the NWES exceeds respiration between mid-January and the end of May, reaching its maximum at the beginning of May (Fig. 9b). Carbon respiration accordingly surpasses carbon fixation at the end of May and remains higher until mid-January, with the maximum reached in July. This results in a range of approximately 4.5 months of net autotrophy, with the maximum NCP in April, and 7.5 months of net heterotrophy, with minimum NCP in October.

Considering the observed OM stoichiometry in our simulations leads to an intensification of seasonal carbon fixation, respiration, and NCP. As with the annual mean fluxes, the intensification is lowest for the PR configuration and highest for the ER configuration. The maximum reference carbon fixation of 16.3 mol C m⁻² yr⁻¹ reached in May increases to a range of 18.3 (+12 %) to 26.2 mol C m⁻² yr⁻¹ (+60 %) in the NWES. Similarly, the maximum carbon respiration of 15.1 mol C m⁻² yr⁻¹ in July increases to values of 16.9 (+12 %) to 23.3 mol C m⁻² yr⁻¹ (+54 %). For NCP, both the maximum of 7.1 mol C m⁻² yr⁻¹ in April and the minimum of −4.3 mol C m⁻² yr⁻¹ in October are enhanced. The seasonal maximum in NCP increases to a range of 7.7 (+8 %) to 12.3 mol C m⁻² yr⁻¹ (+72 %) in the variable stoichiometry configurations. At the same time, the seasonal minimum decreases to an amplitude of −4.8 (−13 %) to −7.3 mol C m⁻² yr⁻¹ (−72 %). The response is nearly identical in the North Sea. In conclusion, these results indicate that variable OM stoichiometry intensifies both the seasonal biological production of organic carbon and its subsequent degradation. In combination with the lateral and vertical gradient in NCP, this finding suggests an increased biological control of the seasonal pCO₂ and the resulting air–sea CO₂ exchange.

3.1.5 Net air–sea CO₂ exchange

In response to the enhanced OM cycling, all variable stoichiometry configurations result in an increased net oceanic CO₂ uptake from the atmosphere. The magnitude of this increase in the CO₂ uptake, compared across configurations, is proportional to the increase in carbon fixation and respiration (Figs. 7–9) and the resulting increase in the vertical gradient of NCP (Fig. 10). The reference simulation yields an average net CO₂ uptake by the ocean of 0.67 Tmol C yr⁻¹ for the North Sea and 1.11 Tmol C yr⁻¹ for the entire NWES. These values are slightly below observational estimates of 0.79 Tmol C yr⁻¹ for the North Sea (Thomas et al., 2005) and 1.3 to 3.3 Tmol C yr⁻¹ for the NWES (Legge et al., 2020). In the variable stoichiometry configurations, the enhanced CO₂ uptake is in good agreement with the literature estimates in the North Sea, with an annual CO₂ uptake of 0.74 (+10 %) to 0.89 Tmol C yr⁻¹ (+33 %). In the NWES, the annual CO₂-uptake is on the lower end of the literature estimates, 1.21 (+9 %) to 1.44 Tmol C yr⁻¹ (+31 %). Thus, considering the observed stoichiometry improves the representation of the annual net air–sea CO₂ exchange both for the North Sea and the entire NWES.

In terms of the spatial distribution, the increase in CO₂ uptake is not uniform across the shelf but instead amplifies existing regional differences (Fig. 11). The reference configuration reproduces the observed north–south gradient across the North Sea, with the annual CO₂ uptake decreasing from the NNS and NT sub-regions towards the SNS sub-region (Kitidis et al., 2019; Schiettecatte et al., 2007; Thomas et al., 2004). All variable stoichiometry configurations consistently strengthen this latitudinal gradient in the air–sea CO₂ exchange. However, the amplitude of this effect varies consistently with the magnitude of the increased annual uptake. The largest increase in CO₂ uptake is consistently shown in the deeper central shelf areas and parts of the outer shelf regions including the NNS, CNS, NT, SWC, MS, AS, and WI. In the ER configuration, the mean CO₂ uptake increases by 0.3 to 0.5 mol C m⁻² yr⁻¹ in these regions. Here, the CNS shows the greatest average increase, ranging from 0.1 mol C m⁻² yr⁻¹ in the PR to 0.5 mol C m⁻² yr⁻¹ in the ER configuration. Smaller but consistent increases are displayed for the southern and near-coastal regions including the SNS, SK, and SWC regions. A persistent increase in CO₂ outgassing or decrease in CO₂ uptake is only shown in a narrow band along the Norwegian coast, near the Elbe estuary, and in inner shelf regions including the EC and IS. In the Atlantic Ocean, the CO₂ uptake consistently increases. Overall, the variable OM stoichiometry configurations demonstrate an increase in the annual net CO₂ uptake, with an intensification of the north–south gradient across the North Sea. This is consistent with the enhanced NCP gradient across the slope of the Norwegian Trench (Fig. 9) and between surface and sub-surface NCP (Fig. 10). The latter suggest additional drawdown of DIC with potential for a stronger CO₂ uptake, especially in but not limited to this region.

Figure 11 Annual mean air–sea CO₂ exchange for the RS configuration and differences for the ER, PR, and ER&PR configurations over the simulation period of 2000–2010. The direction of the flux is defined as downward, meaning a positive air–sea CO₂ exchange indicates an oceanic uptake and a negative one an outgassing of CO₂. Similarly, positive differences represent an increase in uptake or a decrease in outgassing and vice versa.

3.1.6 Seasonal air–sea CO₂ exchange

In addition to the net annual changes, the seasonal intensification of OM cycling also results in a shifted seasonality of the air–sea CO₂ exchange. This can be attributed to an intensification of the biological contribution to changes in pCO₂ relative to the influence of temperature, which is further discussed below. The RS reference configuration shows the strongest CO₂ uptake during winter and spring (Fig. 12), with averages of 1.7 and 2.4 mol C m⁻² yr⁻¹, respectively. In winter, the highest uptake occurs in the NNS and SNS, with the CNS exhibiting the lowest uptake. For the remaining year, there is a consistent north–south gradient, with the highest uptake rates shown in the northern half and weaker uptake or outgassing in the southern half. In summer and autumn, the southern parts and inner shelf areas of the NWES including the SNS, SK, EC, IS, and parts of the CNS display an outgassing of CO₂. The extent of this outgassing is greatest in summer, while the intensity is greatest in autumn. On average, there is a net outgassing of −0.2 mol C m⁻² yr⁻¹ in summer and autumn. This qualitative pattern with a pronounced north–south gradient is overall in good agreement with previous monthly estimates from observations (Kitidis et al., 2019). However, the seasonality indicates that in the reference configuration, the effect of biological carbon fixation on the air–sea CO₂ exchange is underestimated relative to the influence of temperature. This is demonstrated by the majority of uptake occurring during winter and spring, although observations also suggest that significant uptake continues during summer (Kitidis et al., 2019; Thomas et al., 2004).

The seasonal impact of variable stoichiometry is consistent across configurations, with differences in amplitude proportional to the annual changes. For each season, the smallest differences in uptake and outgassing are shown in the PR configuration, while the largest are seen in the ER configuration. In winter, the CO₂ uptake is reduced by 0.2 to 0.8 mol C m⁻² yr⁻¹ across the entire NWES. The strongest reduction is shown in the central and outer shelf areas, as well as the Norwegian Trench. During the spring and summer months, the CO₂ uptake is significantly stronger, with the most pronounced increases in the central North Sea and the Norwegian Trench. The average increase is 0.1 to 0.7 mol C m⁻² yr⁻¹ in spring and 0.4 to 1.3 mol C m⁻² yr⁻¹ in summer, resulting in a shift from net outgassing to net uptake in summer. In autumn, the gradient between uptake in the northern and outgassing in the southern parts is enhanced, with an average decrease in uptake or increase in outgassing of up to 0.2 mol C m⁻² yr⁻¹. This is consistent with the increased autotrophy in spring and summer and increased heterotrophy during autumn and winter. These shifts in the seasonal air–sea CO₂ exchange indicate that the additional drawdown of DIC through increased carbon fixation enhances CO₂ uptake and decreases outgassing during the productive season, while additional respiration reduces uptake and promotes outgassing in autumn and winter. This implies that the air–sea CO₂ flux becomes less temperature controlled and more biologically controlled. This indication is supported by the intensified seasonal carbon fixation and respiration (Figs. 9 and S15–S17), as well as the increased vertical gradient in NCP (Fig. 10). In the absence of biological carbon fixation, seasonal temperature variations would lead to a high CO₂ uptake due to the low temperatures in winter and a lower uptake or higher outgassing due to high temperatures in summer. By depleting the surface waters of DIC during the productive season, the biological pump drives the seasonality in the opposite direction. Thus, with the increased carbon fixation in the variable stoichiometry configurations, the CO₂ uptake shifts towards the productive season.

Figure 12 Seasonal mean air–sea CO₂ exchange for the RS configuration and differences for the ER, PR, and ER&PR configurations over the simulation period of 2000–2010. The CO₂ exchange for each configuration is shown in the Supplement (Fig. S18).

3.2 Organic and inorganic carbon concentrations

The seasonal variations in organic and inorganic carbon concentrations are regulated by the balance between biological carbon fixation and respiration. Consequently, the intensification of the annual and seasonal OM cycling as described above has implications for the transformation and transport of both organic and inorganic carbon. In turn, these carbon pools determine carbon fluxes, which are of relevance for the regional carbon budget. Therefore, we assess the changes in organic and inorganic carbon contents in the context of the enhanced OM cycling and identify further potential consequences for marine carbon cycling in the NWES. Our analysis focuses on changes in the surface DIC and pCO₂, which are relevant for the air–sea CO₂ exchange, as well as seasonal vertical gradients in DIC, which influence cross-shelf exchanges of inorganic carbon. With regard to organic carbon contents, we assess the seasonal vertical distribution of both DOC and POC, which represent how long inorganic carbon uptake is seasonally retained in OM and control cross-shelf transports and sedimentation of organic carbon.

3.2.1 Surface inorganic carbon seasonality

Biogeochemical models constrained to elemental Redfield ratios tend to underestimate the seasonal uptake of DIC by biological carbon fixation in surface waters (Bozec et al., 2006; Kähler and Koeve, 2001; Prowe et al., 2009). Here, the RS reference configuration demonstrates this underestimation of seasonal variability in surface DIC concentrations and ΔpCO₂ (Fig. 13). While ΔpCO₂ in winter agrees well with observational estimates, the seasonal minimum between April and October is strongly underestimated in the RS configuration for both DIC and ΔpCO₂.

The variable stoichiometry configurations demonstrate a higher reduction in surface DIC concentrations during summer in proportion to the increased seasonal carbon fixation described above. The decreased summer ΔpCO₂ relative to the RS configuration increases the potential for CO₂ uptake in summer, which is consistent with the abovementioned shift from the maximum CO₂ uptake in winter and spring towards spring and summer, connected to a more pronounced biological contribution. The additional DIC uptake leads to a reduction in the RMSD of ΔpCO₂ between simulated and observed monthly means, dropping from 57, 21, 41, and 66 µatm in the NT, NNS, CNS, and SNS regions to 37–50, 16–19, 21–32, and 47–58 µatm, respectively. This improvement is most significant in the ER configuration and least significant in the PR configuration. In contrast, the RMSD for DIC shows smaller, less consistent differences. This is likely due to a small general overestimation of surface DIC concentrations and an overestimation of DIC during the second half of the year.

Although reproducing the observed OM stoichiometry improves the representation of summer minima in ΔpCO₂ and DIC, even with variable stoichiometry, simulated DIC increases more rapidly than observed levels do in the second half of the year. This DIC overestimation can be attributed to a lack of organic carbon retention, such that respiration exceeds production too early, thereby releasing DIC back into the water column. One factor limiting the retention of organic carbon may be the missing representation of higher trophic levels including pelagic fish stocks and macrobenthos. Incorporating these and explicitly representing the microbial loop instead of solely relying on bulk remineralization rates for DOM and POM may further improve the model's ability to capture the seasonal dynamics of inorganic carbon, particularly in late summer and autumn.

Figure 13 Simulated daily 11-year mean of horizontal means and standard deviations of DIC concentrations and differences between simulated surface pCO₂ and prescribed atmospheric concentrations from the Mace Head measurement station (Lan et al., 2023) for all four model configurations. The error bars indicate monthly means and standard deviations of corresponding GLODAP and SOCAT observations for DIC and pCO₂. The conversion from simulated CO₂ concentrations and observed fCO₂ to pCO₂ was conducted in accordance with the methodology described above for the carbonate system validation. For the derivation of the ΔpCO₂ from observed pCO₂, the same monthly atmospheric concentrations were subtracted as for the simulated ΔpCO₂. Both variables are shown for the Norwegian Trench (NT) as well as the northern (NNS), central (CNS), and southern North Sea (SNS).

3.2.2 Seasonal vertical DIC gradient

In the deep and seasonally stratified central and northern parts of the North Sea, the separation of OM production in the euphotic zone and respiration in deeper layers leads to a seasonal vertical DIC gradient (Thomas et al., 2004). The high productivity in spring and summer reduce DIC in the surface waters, exporting OM to the sub-surface layer, where subsequent respiration increases DIC concentrations at depths. Vertical mixing in autumn and winter then reduces the gradient and finally restores the nearly homogeneous vertical distribution in winter. This qualitative pattern is reproduced in the RS reference configuration (Fig. 14) but with an underestimated seasonal inorganic carbon drawdown. Because observations are limited and represent only parts of each season, often covering a single month, direct comparisons of seasonally averaged gradients are of limited value. Therefore, we focus here on the representation of the maximum DIC drawdown in summer and the implications of variable stoichiometry for DIC drawdown to sub-surface layers.

In both sub-regions, the vertical DIC gradient in summer is more closely reproduced in the variable stoichiometry configurations, with the closest representation in the ER configuration. Compared to the RS configuration, the additional carbon fixation leads to a reduction in the vertically integrated DIC content in the upper 50 m during summer by 122 to 564 mmol m⁻² in the NT and 169 to 820 mmol m⁻² in the combined NNS and CNS sub-regions. Conversely, in autumn and winter, bottom DIC concentrations are significantly higher compared to those in the RS configuration due to additional respiration. This increase in bottom DIC concentrations is more pronounced in the NT region. The vertically integrated autumn DIC concentrations below 50 m in the NT show a large increase of 1367 to 3855 mmol m⁻² compared to the reference configuration. In the NNS and CNS regions, this increase only ranges from 176 to 384 mmol m⁻². The much larger increase in bottom DIC concentrations in the NT indicates additional external sources of DIC to this region from additional downwelling of organic carbon. An intensification of the vertical DIC gradient overall suggests a higher biological drawdown of DIC, with the potential for higher DIC concentrations in waters exported to the North Atlantic.

Figure 14 Simulated and observed horizontally averaged seasonal DIC gradients in the Norwegian Trench (NT) and the combined northern (NNS) and central North Sea (CNS) sub-regions for all four model configurations. Solid lines represent the simulated concentrations, while dashed lines and the shaded area indicate observed concentrations and their corresponding standard deviation. The observed vertical profiles were generated by firstly calculating the mean concentrations and standard deviation within every meter of depth and, secondly, generating a fifth-order polynomial fit for the vertical profiles of means and standard deviations using numpy.polyfit for each season. The polynomial coefficients are listed in the Supplement (Table S16).

3.2.3 Seasonal organic carbon production

The enhanced seasonal carbon fixation in the variable stoichiometry configurations also results in an increased seasonal accumulation of DOC in both the well-mixed and seasonally stratified regions, with the greatest accumulation occurring in the euphotic zone (Fig. 15). For the purpose of comparison, we present the means and ranges of derived (semi-)labile DOC (LDOC) concentrations in surface and bottom waters in August from 2 consecutive years. These observational ranges demonstrate considerable interannual variability. The ranges from both years demonstrate that the simulated DOC concentrations are a reasonable representation of concentrations during the summer months. The RS stoichiometry is situated at the lower end of the estimated range, while the ER configuration is situated at the upper end.

In the combined NNS and CNS region, the water-column-integrated DOC content during summer exhibits a 68 % to 123 % increase in the ER&PR and ER configurations relative to the reference simulation. In contrast, the PR configuration exhibits a small decrease of 1.2 % in the total summer DOC content, with concentrations increasing in the surface layer but decreasing in the sub-surface layer. A consistent increase in the entire water column is observed in the variable stoichiometry configurations, resulting in a 4 %, 48 %, and 82 % increase in the vertically integrated DOC content in the SNS during summer for the PR, ER&PR, and ER configurations, respectively. The surface DOC content approximately doubles in the SNS and triples in the NNS and CNS during all seasons for the ER configuration, with a maximum increase of approximately 25 % in the PR and an intermediate increase in the combined configuration. The increased DOC concentrations are likely to contribute to carbon export from surface waters through downwelling transport. By retaining more organic carbon in the water column that can be redistributed and respired elsewhere, this likely contributes to the greater depletion of surface DIC concentrations and hence higher CO₂ uptake.

Figure 15 Vertical distribution of the horizontally averaged seasonal mean (solid) and standard deviation (shaded) of DOC concentrations for all four model configurations over the simulation period from 2000 to 2010 in the southern North Sea (SNS) and over the combined northern (NNS) and central North Sea (CNS) sub-regions. The error bars indicate derived estimates of mean LDOC concentrations and observed ranges. These estimates were derived from observed seasonal bulk DOC concentrations in the seasonally stratified and well-mixed regions of the North Sea from Chaichana et al. (2019) by subtracting an average surface ocean refractory or recalcitrant DOC (RDOC) concentration of 44.3 µM from three observational estimates (Aminot and Kérouel, 2004; Hopkinson and Vallino, 2005; Liang et al., 2023).

The enhanced seasonal carbon fixation in the variable stoichiometry configurations further increases the POC content during spring and summer (Fig. 16). Simulated summer POC concentrations are consistent with observations across all configurations in the CNS. In the SNS, the reference configuration underestimates the mean POC concentrations, which are more closely represented in all variable stoichiometry configurations. The vertically integrated POC content is 22 % to 31 % higher in the CNS, with the largest increase in the PR and the smallest increase in the ER configuration. In the SNS, increases are larger and more consistent across configurations, with a range of 55 % to 66 %. The increases are smaller in spring and autumn and negligible in winter. The overall larger POC concentrations in the variable stoichiometry configurations suggest a higher availability of organic carbon for lateral and vertical transports, as well as sedimentation of POC. An increased export from surface waters and subsequent respiration of POC likely contribute to the abovementioned higher drawdown of DIC and the resulting increase in CO₂ uptake. The comparable increase in POC concentrations for the variable stoichiometry configurations with a varying increase in CO₂ uptake indicates that both POC and DOC production contribute to the drawdown of DIC, which is relevant for the enhanced CO₂ uptake. This also corroborates the high additional accumulation of DIC from downwelling into and subsequent respiration of organic carbon from the Norwegian Trench.

Figure 16 Vertical distribution of horizontally averaged seasonal mean and standard deviation of POC concentrations in the central (CNS) and southern North Sea (SNS) for all four model configurations. The error bars indicate estimated mean POC concentrations and observed ranges from Chaichana et al. (2019).

4 Discussion

For the first time, this study provides a comprehensive modeling assessment of how observed variations in the $\mathrm{C}:\mathrm{N}:\mathrm{P}$ stoichiometry of OM impact model-based carbon cycle estimates in the NWES compared to OM cycling at fixed Redfield ratios. We focus on changes in OM cycling in terms of carbon fixation, respiration, NCP, the resulting export of DIC to sub-surface layers, and its consequences for the air–sea CO₂ exchange.

By implementing two pathways for variable $\mathrm{C}:\mathrm{N}:\mathrm{P}$ stoichiometry in the production and respiration of OM (Fig. 3), we demonstrate that variable stoichiometry consistently results in an increased oceanic CO₂ uptake, with a range of 10 %–33 % in the North Sea and 9 %–31 % in the entire NWES. This additional CO₂ uptake is driven by an enhanced biological drawdown of CO₂ through intensified surface production and sub-surface respiration of organic carbon (Figs. 7–10). With the increase in OM cycling, both the North Sea and the entire NWES become more net heterotrophic (Figs. 7–9), with an increased NCP in net autotrophic regions and a decreased NCP in net heterotrophic regions (Fig. 9). This amplifies the pre-existing gradient in NCP, especially along the border between the central North Sea and the Norwegian Trench. Further, the vertical gradient between the net autotrophic surface layer and net heterotrophic sub-surface layer increases, leading to increased DIC uptake at the surface and increased DIC release at depth (Fig. 10). Seasonally, there is an increase in net autotrophy during spring and early summer, while net heterotrophy increases during the rest of the year. This increase in the seasonal amplitude in NCP appears without a notable shift in the timing of the transition from net autotrophy to net heterotrophy. The resulting changes in both organic and inorganic carbon concentrations (Figs. 13–16) have implications for marine carbon cycling in the NWES.

Firstly, the reduction in surface DIC concentrations and hence in ΔpCO₂ enhances the annual net CO₂ uptake (Figs. 7 and 11) and additionally shifts the seasonality to be more biologically controlled and therefore less dominated by temperature. This is expressed by an increased CO₂ uptake in spring and summer and a reduced CO₂ uptake in winter. Furthermore, there is a stronger north–south gradient, with the maximum uptake in the stratified northern regions where export of organic carbon to deeper layers is most efficient (Fig. 12). Secondly, a higher seasonal accumulation of both DOC and POC indicates greater availability of OM contents for lateral transport, the export from surface waters through sinking and downwelling, and the sedimentation of POC (Figs. 13–16). Since comparable simulated POC concentrations with different DOC concentrations show large differences in air–sea CO₂ exchange, these results suggest that both additional POC and DOC production and the resulting carbon export from surface waters contribute significantly to seasonal DIC reduction in surface waters.

In particular, an increase in sub-surface DIC in seasonally stratified areas (Fig. 14) in combination with higher availability of organic carbon (Figs. 15–16) and a stronger vertical (Fig. 10) and lateral gradient in NCP across the slope of the Norwegian Trench (Fig. 9) indicate an intensification of the shelf carbon pump mechanism. This is likely due to additional degradation of OM in the Norwegian Trench imported from the North Atlantic and other parts of the North Sea, resulting in DIC enrichment of deeper shelf waters, which are subsequently exported to the North Atlantic. The increased export of DIC to the open ocean and potentially to sediments is likely responsible for the net uptake of atmospheric CO₂.

The reason for the increased productivity varies between the two mechanisms. The preferential remineralization of organic nitrogen and phosphorus increases production through higher nutrient availability, while the extracellular release of DOM allows for additional carbon fixation beyond nutrient limitation, which additionally reduces surface nutrient concentrations (Supplement Fig. S19). For a comparable range in OM stoichiometry, which reflects the observed variations in the composition of LDOM and POM (Figs. 4 and 5), the impact described is of a higher amplitude for the extracellular release of DOM than for the preferential remineralization of organic nitrogen and phosphorus across all of these carbon fluxes (Fig. 7). As expected, the combined configuration, with a reduced contribution of each mechanism, consistently reproduces intermediate changes between the two individual configurations. While the extracellular release of DOM more effectively accounts for the underestimated reduction in DIC in surface waters, the preferential remineralization better reflects the vertical increase in C:N and C:P ratios with depth and, as such, age of the material. Since both mechanisms have been suggested by observational and experimental studies and reproduce relevant features of OM cycling, the consideration of variable stoichiometry in both production and decomposition of OM is likely necessary for an adequate model representation of marine carbon cycling.

While previous studies have considered carbon fixation beyond elemental Redfield ratios in the North Sea and Baltic Sea (Neumann et al., 2022; Prowe et al., 2009), here we explicitly quantify the impact on regional carbon cycling across the NWES and assess the individual and combined contributions of two contributing mechanisms. Our findings regarding the regional impact on marine carbon cycling are consistent with previous global Earth system model experiments considering the preferential remineralization of organic nitrogen and phosphorus, as well as with proposed consequences of higher C:N and C:P ratios in OM from observational studies. In particular, our results regionally support the suggestion that increased carbon content relative to elemental Redfield ratios would imply higher export production, an associated drawdown of inorganic carbon, and increased availability of organic carbon for lateral carbon transport (Hopkinson and Vallino, 2005; Letscher and Moore, 2015; Tanioka et al., 2021; Zakem and Levine, 2019). Beyond the large-scale effects, we highlight how the regional carbon export mechanism determines the response to variable OM stoichiometry in terms of spatial distribution and seasonality and, as such, needs to be considered regionally as well as globally.

Given the significant impact on major carbon fluxes, we show that estimates of regional carbon budgets are sensitive to a reasonable representation of OM stoichiometry. These findings not only underscore the importance of the elemental composition of OM in regional carbon cycling but also suggest that global carbon budget models may need to be re-evaluated to account for its accurate representation. The pronounced seasonal shifts in CO₂ uptake driven by variable OM stoichiometry highlight the need for its consideration in carbon modeling in temperate shelf seas globally.

Future research efforts should focus on identifying the contribution of different processes to the observed OM stoichiometry and their spatiotemporal variability. Other effects on OM composition that are beyond the scope of this study include the quality of imported OM from land (Painter et al., 2018), the effect of variable stoichiometry in prey on higher trophic levels (Schindler and Eby, 1997), and pathways for variable stoichiometry in phyto- and zooplankton (Moreno and Martiny, 2018). The latter pathways include the direct uptake of DON and DOP by phytoplankton (Fitzsimons et al., 2020), as well as variable $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratios in zooplankton grazing and excretion (Anderson et al., 2005; Elser and Urabe, 1999). Another important factor to consider is the role of community composition in shaping ecosystem stoichiometry (Kwiatkowski et al., 2018). All of these contributions may be subject to regional and interannual variability, as well as long-term trends imposed by climate change (Sardans et al., 2012). Additionally, spatiotemporal variations in ecosystem processes, such as the partitioning of new detritus into POM and DOM, grazing rates, and temperature-dependent remineralization rates, could introduce regional and interannual variability not explicitly accounted for in this study. While these processes contribute to organic matter cycling and could influence our findings, the model's accurate representation of OM concentrations and stoichiometry ensures that the broader conclusions regarding regional carbon fluxes remain robust.

While this study offers valuable insights into the drawdown of inorganic carbon and explores changes in both organic and inorganic carbon concentrations, it does not fully quantify the impact on the carbon budget of the NWES. Specifically, the study has not yet established key fluxes explicitly, such as lateral cross-shelf transports of dissolved organic carbon (DOC) and particulate organic carbon (POC), vertical POC sedimentation rates, and the export of dissolved inorganic carbon (DIC) to the open ocean. Future research should aim to quantify these fluxes to provide a more comprehensive understanding of the impact on the full carbon budget of the NWES. Setting these impacts in the context of multiple anthropogenic stressors including climate change impacts will be crucial for refining carbon budget estimates and predicting long-term carbon sequestration in this region and their significance for global carbon cycling.

5 Conclusions

This study highlights the significant impact of the observed variations in $\mathrm{C}:\mathrm{N}:\mathrm{P}$ stoichiometry of OM on carbon cycling in the NWES. By incorporating two pathways for variable OM stoichiometry into the regional 3D physical–biogeochemical modeling system SCHISM-ECOSMO-CO₂, we reproduce the observed stoichiometry of LDOM and POM. In comparison to fixed Redfield stoichiometry, our results show an enhanced biological drawdown of DIC and a resulting increase in the oceanic CO₂ uptake. This additional CO₂ uptake is driven by an intensification of OM cycling and the resulting increase in seasonal, lateral, and vertical NCP gradients, which enhances the seasonal DIC gradient. With an increased biological contribution to the seasonality of pCO₂, CO₂ uptake shifts from winter towards spring and summer. This is particularly evident in the deep and seasonally stratified regions of the North Sea. Despite the difference in the pathways for variable stoichiometry in the production and respiration of OM, the model response shows a consistent impact on OM cycling and the resulting air–sea CO₂ exchange, with noticeable differences only in their magnitude. To determine the amplitude of these changes, a deeper understanding of the individual process contributions to variable stoichiometry is required.

Our findings underscore the limitations of using fixed Redfield stoichiometry for simulating OM cycling, revealing that the observed variations in OM stoichiometry have a profound influence on marine carbon cycling in the NWES. Incorporating variable OM stoichiometry into biogeochemical models will be essential for an accurate representation of carbon cycling, regional carbon budgets, and their spatiotemporal variations. This will be a crucial aspect when evaluating the impact of multiple anthropogenic stressors on the marine ecosystem and the role of OM stoichiometry in future carbon budgets under changing environmental conditions. As the results are largely determined by the pre-existing regional physical controls and the resulting structure of biogeochemical cycles, we emphasize that regional assessments are key to understanding the role of OM stoichiometry in regional carbon cycling and should be incorporated into the representation of the coastal ocean in global modeling studies.