Mesoscale Eddy Variability Enhances Fixed Nitrogen Loss and Suppresses Nitrous Oxide Production in Oxygen Minimum Zones

Within oxygen minimum zones, anaerobic processes transform bioavailable nitrogen (N) into the gases dinitrogen (N2) and nitrous oxide (N2O), a potent greenhouse gas. Mesoscale eddies in these regions create heterogeneity in dissolved N tracers and O2 concentrations, influencing nonlinear N cycle reactions that depend on them. Here, we use an eddy‐resolving model of the Eastern Tropical South Pacific to show that eddies enhance N2 production by between 43% and 64% at the expense of reducing N2O production by between 94% and 104% due to both the steep increase of progressive denitrification steps at vanishing oxygen, and the more effective inhibition of N2O consumption relative to production. Our findings reveal the critical role of eddies in shaping the N cycle of oxygen minimum zones, which is not currently represented by coarse models used for climate studies.

As O 2 approaches anoxia (here defined as O 2 < 5 mmol m −3 ), aerobic remineralization ceases and is followed by stepwise heterotrophic denitrification ( , which generates both N 2 O (as an intermediary) and N 2 (Babbin et al., 2015;Kalvelage et al., 2013;Ward et al., 2009).The release of  NH + 4 from remineralization further stimulates fixed N loss.In oxygenated waters, small amounts of N 2 O are generated as a byproduct of  NH + 4 oxidation to  NO − 2 (the first step of nitrification), with higher yields at lower O 2 concentrations Abstract Within oxygen minimum zones, anaerobic processes transform bioavailable nitrogen (N) into the gases dinitrogen (N 2 ) and nitrous oxide (N 2 O), a potent greenhouse gas.Mesoscale eddies in these regions create heterogeneity in dissolved N tracers and O 2 concentrations, influencing nonlinear N cycle reactions that depend on them.Here, we use an eddy-resolving model of the Eastern Tropical South Pacific to show that eddies enhance N 2 production by between 43% and 64% at the expense of reducing N 2 O production by between 94% and 104% due to both the steep increase of progressive denitrification steps at vanishing oxygen, and the more effective inhibition of N 2 O consumption relative to production.Our findings reveal the critical role of eddies in shaping the N cycle of oxygen minimum zones, which is not currently represented by coarse models used for climate studies.
Plain Language Summary Nitrogen limits photosynthesis over large swaths of the ocean.Oxygen minimum zones are important oceanic environments where the presence of low O 2 concentrations allow anaerobic microbial communities to transform bioavailable forms of nitrogen into the gases dinitrogen (N 2 ) and nitrous oxide (N 2 O, a powerful greenhouse gas), which together lead to nitrogen loss from the oceans.The large scale ocean circulation acts to maintain these oxygen minimum zones, and the gradients in nitrogen and oxygen concentrations surrounding them.Turbulence also generates eddies and filaments which stir these gradients, promoting chemical heterogeneity.However, it is uncertain how eddies impact the nitrogen transformation rates dependent on nitrogen and oxygen tracer concentrations.Using a high resolution ocean model of the Eastern Tropical South Pacific oxygen minimum zone, we explore the contributions from eddies to total N 2 and N 2 O production.We highlight that roughly half of total N 2 production is caused by the presence of eddies, while N 2 O production is drastically reduced or even consumed.These findings expand our understanding of the ways in which eddies alter nitrogen loss from the oceans, and challenge the ability of non-eddy resolving models to accurately represent N 2 and N 2 O production in a changing ocean.
The ETSP hosts the second largest OMZ by area in the global ocean (Bianchi et al., 2012;Fuenzalida et al., 2009).Mesoscale eddies of   10-100 km radii are regularly generated in this region, with no obvious preference for polarity, particularly poleward of 15°S where eddy frequency reaches between 30% and 50% (Chaigneau et al., 2011;Chelton et al., 2011).Cyclones are typically surface intensified, with cores located between 100 and 200 m due to near-surface generation via instability of the Peru Coastal Current.Anticyclones, in contrast, are shed from the Peru-Chile Undercurrent and thus have deeper cores found between 200 and 400 m (Chaigneau et al., 2008).As their vertical influence overlaps with the depth of the OMZ (Chaigneau et al., 2008), mesoscale variability shapes the transport of biogeochemical tracers due to eddy-scale correlations between currents and tracer concentrations, for example, by intermittently supplying O 2 to the OMZ (Bettencourt et al., 2015;Czeschel et al., 2011;Gnanadesikan et al., 2013).
A recognition of the importance of eddy-induced transport has generated a variety of techniques to represent their unresolved influence in global models used for long-term climate and ocean biogeochemical studies (e.g., Gent and McWilliams (1990) and Fox-Kemper et al. (2008)).However, since many biogeochemical transformations depend nonlinearly on tracer concentrations, and occur over timescales comparable to those of turbulent fluctuations (Mahadevan & Campbell, 2002;Pasquero, 2005;Prend et al., 2021), chemical heterogeneity can also lead to amplification or reduction of transformation rates at the scale of eddies (herein referred to as "eddy reactions") (Goodman & Robinson, 2008;Lévy et al., 2012;Martin et al., 2015;Prend et al., 2021;Wallhead et al., 2013).As a consequence, biogeochemical rates estimated from a "mean field approximation," that is, from properties averaged over temporal or spatial scales greater than those of eddies, can fail to capture the magnitude of transformation rates in a fluctuating environment (Brentnall, 2003;Lévy et al., 2013;Martin et al., 2015;Rovinsky et al., 1997).For example, Lévy et al. (2013) used an eddy-resolving biogeochemical model of the North Atlantic to suggest that eddy reactions reduce primary production and increase grazing rates by between 5% and 30% relative to the mean state, owing to correlations between nutrients, phytoplankton and zooplankton caused by eddies.
N cycle reactions responsible for N 2 and N 2 O production are highly nonlinear in their response to environmental variability, in particular oxygen concentrations, suggesting a significant role for eddy reactions.In this study, we examine the role of eddy heterogeneity on fixed N loss and N 2 O production, which is usually ignored by global ocean biogeochemistry and Earth system models used to study OMZ processes (J.J. Busecke et al., 2022;Cabré et al., 2015;Martinez-Rey et al., 2015;Kwiatkowski et al., 2020).To this end, we apply an eddy-mean decomposition formalism (Lévy et al., 2013(Lévy et al., , 2024) ) to a realistic eddy-resolving simulation of the ETSP OMZ (McCoy, Damien, Clements, et al., 2023) to quantify the importance of eddy-scale reactions on N 2 and N 2 O production.

Model Configuration
We employ an identical model configuration for the ETSP as in McCoy, Damien, Clements, et al. (2023), but with a finer horizontal resolution of 5 km (from 10 km) to better resolve mesoscale eddies.The domain extends from −111.38°W to −66.62°W and from 42.52°S to 3.41°N, encompassing the ETSP OMZ, the equatorial current system, and the wind-driven subtropical gyre circulation (Figure 1a).
The physical component consists of the Regional Ocean Modeling System (ROMS), a primitive-equation, hydrostatic, topography-following ocean general circulation model (Shchepetkin, 2015;Shchepetkin & McWilliams, 2005).ROMS is coupled online to the Biogeochemical Elemental Cycling (BEC) model (Deutsch et al., 2021;Moore et al., 2004) with an expanded representation of OMZ biogeochemistry (NitrOMZ) (Bianchi et al., 2022).We analyze daily output from a 10-year simulation initialized with the 10 km configuration (McCoy, Damien, Clements, et al., 2023), following a one year adjustment period.Model setup and forcings, and a summary of NitrOMZ's formulation and equations, are presented in Sections S1 and S2 in Supporting Information S1, and described in detail in previous publications (Bianchi et al., 2022;McCoy, Damien, Clements, et al., 2023).

N 2 and N 2 O Production
In the model, the net N 2 O and N 2 production rates (  N 2 O and  N 2 , in units of mmol N m −3 s −1 ) are given by: (2) Here,  (Goreau et al., 1980;Ji et al., 2018).

Eddy-Mean Decomposition of Tracer Equations
The conservation equation for a biogeochemical tracer (C) can be written as: (3) Here, u and κ are the velocity vector (u, v, w) and the vertical eddy diffusivity, resulting in the advective transport and turbulent diffusion terms, and J C encapsulates biogeochemical sources and sinks.
The eddy-mean decomposition approach requires filtering Equation 3 to separate contributions from large scale, low frequency "mean" fields from higher frequency "eddy" fluctuations.For example, the decomposition of C follows: Here, the mean filter consists of monthly climatological averages from 10 years of model output.Eddy fluctuations (C′) therefore represent tracer heterogeneity at sub-monthly time scales, which is mostly driven by mesoscale variability (Le Traon, 1991;Morrow & Le Traon, 2012).
Following the approach of Lévy et al. (2013), substituting C for  C + C ′ and u for   +  ′ in Equation 3, and applying the mean filter leads to: (5)

Sign and Amplitude of Eddy Reactions
To attribute eddy reactions to correlations between specific tracer pairs, we follow the approach of Lévy et al. (2013), and approximate monthly-averaged reaction rates (herein referred to as "total" reactions,   ) with a Taylor series expansion: .Note that we ignore third order and higher terms, assuming small fluctuations; this approximation is supported by our results (Section 3.4).Equation 6allows to quantify the contribution of different tracer correlations to the eddy reaction rate.

Mean State and Eddy Variability
The model produces a realistic climatological representation of the physical circulation and biogeochemistry of the ETSP; a detailed validation is presented in McCoy, Damien, Clements, et al. (2023).Briefly, it simulates a persistent OMZ centered at roughly 8°S that extends offshore to nearly 100°W (Figure 1a).Oxygen concentrations shape N tracer gradients in the OMZ (Figures 1d-1g).The oxycline is observed at roughly 100 m depth, coinciding with a maximum in  NH + 4 , and overlies an OMZ ranging in thickness from roughly 200 m offshore to nearly 400 m near the coast (Figures 1d and 1e).In the steep suboxic gradients surrounding the OMZ, high concentrations of N 2 O accumulate (Figure 1g).As O 2 approaches anoxia, characteristic  NO − 2 maxima and N 2 O minima appear (Figures 1f and 1g).
The model generates a vigorous mesoscale eddy field, particularly around 15°S where eddy kinetic energy (EKE) reaches a maximum, in agreement with satellite-based estimates (Figure S1 in Supporting Information S1).Long-lived anticyclonic subsurface eddies, also known as submesoscale coherent vortices (McCoy et al., 2020;McWilliams, 1985), are frequently generated at the coast and propagate eastward into the subtropical gyre.Larger, more ephemeral surface-intensified cyclonic eddies are typically spawned between 5 and 15°S.The monthly filter successfully separates eddy-driven tracer heterogeneity from mean distributions set by the large-scale circulation.The resulting tracer fluctuations (e.g., Figure 1b and Figures S2e-S2h in Supporting Information S1) reveal the variability caused by mesoscale turbulence within the OMZ and along the fringe of the subtropical gyre.Because each tracer is characterized by different mean gradients, eddy fluctuations lead to distinctive patterns of tracer variance (Figures 1h-1k).Comparing the magnitude of these fluctuations to the respective tracer means (Figure 10.1029/2023GL106179 6 of 13 S3 in Supporting Information S1) highlights the vertical and horizontal impact of eddies, which supply O 2 to the anoxic OMZ core, and export N tracers from the OMZ to the adjacent ocean.

Mean and Eddy Reactions
Despite a higher N 2 O yield at decreasing O 2 (Equation S9 in Supporting Information S1), N 2 O production from  NH + 4 oxidation (Figure 2a) reaches a maximum in the upper oxycline and is only partially enhanced within the OMZ.The mean reaction (Figure 2e), which reflects the effect of large-scale tracer distributions (Figures 1d-1g), nearly captures the total reaction apart from an overestimation at the boundary of the anoxic OMZ core, as revealed by the negative sign of the eddy reaction (Figure 2i).Anammox (Figure 2h) is also well captured by the mean tracer distribution (Figure 2d), with only a weak positive contribution from eddies along the upper oxycline, and a weak negative contribution below it (Figure 2l).Integrated over the OMZ domain (outlined in Figure 1a), mesoscale heterogeneity reduces N 2 O production from  NH + 4 oxidation by −7 to −28% (here and in the following, the range encapsulates minimum and maximum monthly mean values) (Figures 2i and 2m), whereas N 2 production from anammox is altered by −17 to +2% (Figures 2l and 2p).
The two denitrification steps are tightly coupled and similarly peak in the upper oxycline, weakening offshore of 85°W at depth (Figures 2b and 2c).In contrast to anammox and  NH + 4 oxidation to N 2 O, their eddy components are positive (Figures 2j and 2k) and greater in magnitude.Thus, eddies significantly amplify both denitrification steps along the fringe of the anoxic core and along the upper oxycline.Integrated over the OMZ, eddy  NO − 2 reduction contributes between 48% and 53% to the total rate (Figure 2n), while eddy N 2 O reduction shows a slightly larger contribution (61%-65%, Figure 2o).
The eddy-mean decomposition reveals that mean N 2 O production overestimates total N 2 O production by an order of magnitude, and is closely compensated by N 2 O consumption at the scale of eddies (Figures 3c and 3e).This result emerges from the much larger difference between mean  NO − 2 and N 2 O reduction rates (Figures 2f and 2g) in the suboxic gradients surrounding the OMZ.Thus, while the mean field approximation underestimates the expression of individual denitrification rates due to the lack of eddies, it overestimates the magnitude of N 2 O production from incomplete denitrification.As a consequence, eddies cause a dramatic decrease (94%-104%) in net N 2 O production compared to the mean field approximation, leading to net N 2 O consumption between June and August (Figure 3g).
Because N 2 O reduction is both a sink of N 2 O and a source of N 2 (Equations 1 and 2), these results also highlight the role of eddy heterogeneity in regulating N 2 production.Despite the negative effect of eddies on anammox (Figure 2l), eddy stimulation of N 2 O reduction enhances N 2 production by 43%-64% (Figure 3h).

Attributing Eddy Effects to Oxygen Variability
To attribute the sign and magnitude of eddy N 2 and N 2 O production to specific aspects of eddy-driven heterogeneity, we apply the Taylor series approach outlined in Section 2.4.Analysis of individual tracer-pair contributions in Equation 6 reveals the dominant role of O 2 fluctuations, such that the eddy N 2 and N 2 O production patterns (Figures 3e and 3f) are nearly recovered using the O 2 variance term alone (Figures 4a and 4b; see also Figures S6-S9 in Supporting Information S1 for contributions from other tracer pairs to the eddy reactions in Figures 2i-2l, and Figure S10 in Supporting Information S1 for full reconstructions).
The role of O 2 fluctuations is explained schematically in Figure 4c.At low, but not completely anoxic O 2 concentrations, N 2 O reduction (blue square) is nearly suppressed, whereas  NO − 2 reduction (red square) is significantly less inhibited.This establishes a "window," that is, an excess of N 2 O production over consumption, where net N 2 O production emerges from incomplete denitrification.This picture is altered in the presence of eddy fluctuations.Because of a higher sensitivity to O 2 (Table S1 in Supporting Information S1), N 2 O reduction (blue curve) increases faster than  NO − 2 reduction (red curve) as O 2 declines below the mean value.Thus, after averaging over the eddy fluctuations, both reactions are enhanced relative to the mean.However the amplification of N 2 O reduction (blue circle, from 2% to 19% inhibited) is greater than that for  NO − 2 reduction (red circle, from 42% to 51% inhibited), and the "window" for N 2 O production by incomplete denitrification is reduced, in agreement with the results in Figure 3e.Effectively, in the presence of eddies, the last two steps of denitrification become more tightly coupled across a broader swath of the OMZ, leading to greater N 2 O consumption.Figure 4c also explains the weaker impact of eddies on anammox (green curve): because anammox is less sensitive to O 2 (Table S1 in Supporting Information S1), fluctuations only slightly amplify the reaction.
While these results rely on the nitrogen cycle model formulation employed herein (Bianchi et al., 2022), they are grounded in the observed progressive intolerance of anaerobic metabolisms to O 2 (Babbin et al., 2015;Bristow et al., 2016;Dalsgaard et al., 2014;Sun et al., 2017), and are robust across a range of O 2 inhibition parameter values (see Section S3 and Figures S11-S14 in Supporting Information S1).Thus, eddy-driven O 2 heterogeneity can be expected to non-uniformly enhance anaerobic processes, leading to an increase in N 2 production in OMZs at the expense of reduced N 2 O production.

Discussion and Conclusions
Numerical models are essential tools for extrapolating biogeochemical processes to global scales.However, their limited spatial resolution hinders accurate representation of the central role played by mesoscale processes (Mahadevan et al., 2012;McGillicuddy et al., 2007).Notably, eddies enhance chemical variability and tracer transport (Damien et al., 2023;Gruber et al., 2011), causing both intrusions of oxygenated waters into the core of the OMZ, and anoxic fluctuations along the OMZ margins .This in turn broadens the overlap between aerobic (e.g., Figure 2a) and anaerobic metabolisms (Figures 2b-2d; see also Buchanan et al. (2023)), and favors metabolite exchange.
Using an eddy-resolving model, we show that eddy heterogeneity slightly suppresses aerobic and anaerobic  NH + 4 oxidation rates, while greatly stimulating denitrification rates.As a result, denitrification rates occurring at eddy scales are of comparable magnitude to the rates dictated by mean tracer fields, and, in the case of N 2 O reduction, even greater (Figures 2j and 2k).This decoupling between denitrification steps is governed by their progressive inhibition at increasing O 2 , while the weaker stimulation of anammox by eddies reflects a greater tolerance to O 2 .
Stimulation of N loss by mesoscale processes may have important consequences in future oceans.As EKE decreases sharply with depth (Ni et al., 2023), changes in the position of the oxycline via increased stratification or wind variations (Bakun, 1990;Deutsch et al., 2014; J. J. M. Busecke et al., 2019) could alter the impact of eddy fluctuations on N 2 and N 2 O production (Yang et al., 2020), and the balance between N fixation and fixed-N losses in the ocean (Wang et al., 2019).
A primary goal of ocean modeling is to upscale physical and biogeochemical processes at unresolved scales based on resolved fields, enabling long-term global simulations that are still prohibitively expensive with eddy-resolving models.Physical scale-dependent parameterizations that capture unresolved eddy transports are widely adopted (Fox-Kemper et al., 2008;Gent & McWilliams, 1990;Jansen et al., 2019).However, analogous parameterization for eddy reactions are missing.Our results add to a growing literature that indicates that neglecting subgrid-scale eddy covariances in tracers provides a biased picture of biogeochemical transformations (Damien et al., 2023;Lévy et al., 2013;Wallhead et al., 2013).We suggest that these issues may be more critical in OMZs, because of the nonlinear dependence of the N cycle to O 2 concentrations.
Note that the choice of horizontal resolution affects both the magnitude and expression of eddy reactions induced by the mesoscale circulation (Lévy et al., 2013), in part because eddy activity strengthens at higher resolution (Marchesiello et al., 2003).Our simulations capture mesoscale eddies, and therefore the bulk of eddy kinetic energy in the ETSP (Chelton et al., 2011).However, unresolved submesoscale processes can be locally important (Hauschildt et al., 2021), in particular near the coast, where the oxycline shoals closer to the mixed layer (Thomsen, Kanzow, Colas, et al., 2016).There, intense vertical velocities induced by submesoscale cross-frontal circulations (Lévy et al., 2022;McWilliams, 2016) likely enhance O 2 intrusions into the OMZ.To what extent these currents can effectively ventilate anoxic waters remains unclear (Thomsen, Kanzow, Krahmann, et al., 2016).It is likely that increased heterogeneity at submesoscale enhances the eddy effects discussed here.A dedicated characterization of submesoscales will require simulations at sub-kilometer resolution (Damien et al., 2023).
Biogeochemical biases are typically corrected by adjusting model parameters, enabling non-eddy resolving models to attain realistic biogeochemical transformation rates.However, this may not always robustly capture the underlying dynamics (Wallhead et al., 2013), such that global models tuned to large-scale observations may fail in predicting responses to perturbations in a changing ocean.These include expansion of OMZs (Schmidtko et al., 2017;Stramma et al., 2008), and shifts between anoxic and suboxic volumes (J.J. Busecke et al., 2022), which continue to be poorly resolved in global Earth system models (Cabré et al., 2015;Séférian et al., 2020).Despite these challenges, a number of methods to upscale biogeochemistry show promise.Higher resolution simulations can provide statistical information such as tracer variances and pairwise covariances driven by eddies (e.g., Figures 1h-1k).These could be validated by high-resolution biogeochemical observations (Martin et al., 2015), or by large-scale observational syntheses (Olsen et al., 2016).With knowledge of the functional forms that represent biogeochemical transformations, this information can be used to reconstruct eddy reactions via a Taylor series approximation (e.g., Section 3.4; Figures 4a and 4b and Figures S8e-S8h in Supporting Infprmation S1), correcting biases in coarse simulations (Wallhead et al., 2013).Our study is a step forward toward developing scale-dependent parameterizations of the effects of eddies on biogeochemistry.This work was supported by NSF Grant OCE-1847687.D.B. acknowledges support from the Alfred P. Sloan Foundation.This work used the Expanse system at the San Diego Supercomputer Center through allocation TG-OCE170017 from the Advanced Cyber infrastructure Coordination Ecosystem: Services and Support (ACCESS) program, which is supported by National Science Foundation Grants 2138259, 2138286, 2138307, 2137603, and 2138296.Satellite-based eddy kinetic energy estimates were processed by SSALTO/DUACS and distributed by AVISO+ (https://www.aviso.altimetry.fr) with support from CNES.

Figure 2 .
Figure 2. (a,e,i) Model transects of total, mean, and eddy  NH + 4 oxidation to N 2 O from model years 2-11.Panels (b,f,j) show the same results, but for  NO − 2 reduction to N 2 O, panels (c,g,k) for N 2 O reduction to N 2 , and panels (d,h,l) for the anammox reaction to N 2 .(m-p) Time-series of integrated total (solid) and mean (dashed) rates from the OMZ domain (see dashed box in Figure1a, extending to 1,000 m).Blue/red shading highlights when eddy reactions contribute negatively/positively to the total reaction rates, respectively.

Figure 3 .
Figure 3. Same as in Figure 2, but for (a, c, e, g) net N 2 O production and (b, d, f, h) net N 2 production.Note that panel (a) is multiplied by a factor of 10 due to the strong compensation between mean (c) and eddy (e) terms.

Figure 4 .
Figure 4. (a), (b) Reconstructed eddy net N 2 O and N 2 production (Figures 3e and 3f) following the Taylor series expansion in Equation 6 using only the contribution from O 2 variance (Figure 1h).(c) Schematic of eddy effects on the exponential inhibition by O 2 for anammox (   2  , green curve),  NO − 2 reduction (   2 2 , red curve), and N 2 O reduction (   2 3 , blue curve); see Equations S4 and S8 in Supporting Information S1.An example mean field approximation, here assuming  O2 of 2 mmol m −3 , is shown with square markers, while the total inhibition (circle markers) is estimated assuming O 2 fluctuations (  O ′ 2 ) of ±1.5 mmol m −3 .N 2 O production "windows" (mean and total) are defined as the difference in O 2 -dependent inhibition between  NO − 2 reduction (which produces N 2 O) and N 2 O reduction (which consumes it).
x ) ⏟⏞⏞⏟⏞⏞⏟  = 0 , the linear terms also reduce to zero, and the monthly-averaged eddy reaction rates can be estimated with knowledge of eddy tracer variances and pairwise (x1,...,x ), and total rates would therefore depend on mean reactions only.For nonlinear reactions, since  x ′