Biogeochemical feedbacks may amplify ongoing and future ocean deoxygenation: a case study from the Peruvian oxygen minimum zone

A new box model is employed to simulate the oxygen-dependent cycling of nutrients in the Peruvian oxygen minimum zone (OMZ). Model results and data for the present state of the OMZ indicate that dissolved iron is the limiting nutrient for primary production and is provided by the release of dissolved ferrous iron from shelf and slope sediments. Most of the removal of reactive nitrogen occurs by anaerobic oxidation of ammonium where ammonium is delivered by aerobic organic nitrogen degradation. Model experiments simulating the effects of ocean deoxygenation and warming show that the productivity of the Peruvian OMZ will increase due to the enhanced release of dissolved iron from shelf and slope sediments. A positive feedback loop rooted in the oxygen-dependent benthic iron release amplifies, both, the productivity rise and oxygen decline in ambient bottom waters. Hence, a 1% decline in oxygen supply reduces oxygen concentrations in sub-surface waters of the continental margin by 22%. The trend towards enhanced productivity and amplified deoxygenation will continue until further phytoplankton growth is limited by the loss of reactive nitrogen. Under nitrogen-limitation, the redox state of the OMZ is stabilized by negative feedbacks. A further increase in productivity and transition to sulfidic conditions is only possible if the rate of nitrogen fixation increases drastically under anoxic conditions. Such a transition would lead to a wide-spread accumulation of toxic sulfide with detrimental consequences for fishery yields in the Peruvian OMZ that currently provides a significant fraction of the global fish catch.

Vol:. (1234567890) twice as high as predicted by ocean models that consider the effects of warming on oxygen uptake in the surface ocean and biological respiration (Oschlies et al. 2018). The gap between observed and predicted trends indicates that physical, biological and biogeochemical processes governing the oxygen content of the ocean are not yet fully understood. It is conceivable that underrepresented respiration processes and nutrient dynamics have contributed to the observed oxygen decline.
Oxygen minimum zones (OMZs) located in the tropical oceans have expanded significantly over the last decades (Stramma et al. 2008). In addition to ocean warming, natural variability induced by e.g. the Pacific Decadal Oscillation has contributed to the observed oxygen depletion (Deutsch et al. 2011;Stramma et al. 2020). Moreover, biogeochemical feedbacks may have amplified the OMZ expansion since marine sediments release dissolved iron and phosphorus into the water column where low oxygen waters interact with the seafloor (Berelson et al. 2003;Dale et al. 2015;Ingall et al. 1993;Ingall and Jahnke 1994;Klar et al. 2018;McManus et al. 1997;Plass et al. 2020;Rapp et al. 2019;Wallmann 2003Wallmann , 2010. The Peruvian OMZ has been studied in detail by e.g. the Collaborative Research Centre 754 (Krahmann et al. 2021) and during cruise GP16 of the GEOTRACES program Moffett and German 2018). It is located in the Eastern Tropical Pacific and marked by low oxygen subsurface waters extending from the Peruvian coast into the open ocean. Marine productivity in this region is constrained by light availability and limited by iron that is released from continental shelf and slope sediments and transported into surface waters by upwelling and vertical mixing Chever et al. 2015;DiTullio et al. 2005;Hutchins et al. 2002;John et al. 2018;Noffke et al. 2012;Plass et al. 2020;Rapp et al. 2020;Schlosser et al. 2018b;Scholz et al. 2016). Dissolved phosphate and nitrate are supplied by the nutrient-rich Peru-Chile Undercurrent which is the main source for the coastal upwelling off Peru (Montes et al. 2014). Reactive nitrogen is consumed by anaerobic oxidation of ammonium (anammox) in low-oxygen subsurface waters  as well as pelagic (Dalsgaard et al. 2012) and benthic denitrification (Bohlen et al. , 2011. The relative importance of pelagic denitrification vs. anammox as sinks for reactive nitrogen is currently debated (Ward 2013). Phosphate is consumed by biological uptake and replenished by anoxic sediments deposited on the shelf and upper slope that release large amounts of dissolved phosphate into ambient bottom waters Noffke et al. 2012).
In the following sections, we provide a detailed budget of nutrient fluxes in the Peruvian OMZ and analyze how the redox state and productivity of this system may respond to global deoxygenation and warming. We selected the Peruvian margin for this study because it is the most productive region in the ocean and provides a significant fraction of the global fish catch (Chavez et al. 2008). We focus on the cycling of nitrogen and iron since these elements serve not only as nutrients facilitating marine productivity but are also involved in a range of redox processes that respond very sensitively to ocean change. Moreover, positive and negative feedbacks that regulate the turnover of these redox-sensitive nutrients may either amplify or mitigate regional and global ocean change. We develop and present a new biogeochemical box model that-for the first time-resolves all known feedbacks in pelagic and benthic nutrient cycling. The large data set available for the Peruvian OMZ is employed to calibrate the new model. After careful calibration, the new model is used to forecast how redox-state, productivity and fishery yields in the Peruvian OMZ may change over the coming decades considering the non-linear response of the system to ongoing ocean warming and deoxygenation. Finally, we discuss how our modeling outcomes may help to explain the fast deoxygenation and OMZ expansion observed in the global ocean.

Study area and model set-up
The model covers the continental margin off Peru between 6° S and 18° S and the offshore area extending to 100° W (Fig. 1). The study area that captures the most intense part of the Peruvian OMZ is divided into 12 subregions. The surface layer (0-50 m), upper (50-300 m) and lower intermediate water  and deep water (1000 m to seafloor) are represented by individual boxes. The horizontal and vertical water fluxes between the resulting 48 boxes and the horizontal fluxes across the outer boundaries of the model are defined using output of a regional ocean model system (ROMS) for the Peruvian OMZ (José et al. 2021;Montes et al. 2014). They represent annual mean volume fluxes averaged over the period 1990-2010 (Supplementary Information, Section S1.1). An extended version of the REDBIO model is used to simulate the pelagic and benthic turnover of oxygen, sulfide, nutrients (nitrate, nitrite, ammonium, phosphate, ferric and ferrous iron) and organic matter (Jurikova et al. 2020;Wallmann et al. 2019). REDBIO was developed to simulate the redox-dependent cycling of nutrients during oceanic anoxic events that occurred in the geological past . Here, REDBIO equations for iron and nitrogen turnover are updated and additional processes are considered to fully resolve nutrient cycling in the Peruvian OMZ. The modeling approach is presented in Sect. 4 and the biogeochemical machinery employed in the model is described in the Supplementary Information (Section S1.2).

Model optimization and sensitivity
Export production and respiration are largely controlled by the turnover of iron since biological productivity is limited by this trace element across the model region (Browning et al. 2018;Bruland et al. 2005). Dissolved iron concentrations are in turn mainly governed by the release of iron from anoxic shelf and slope sediments and by removal through oxidation and scavenging processes. While the benthic iron release can be projected using the rate of organic matter deposition at the seabed and the dissolved oxygen concentrations in bottom waters (Dale et al. 2015), scavenging rates are largely unconstrained. A range of simple and complex models have been used to simulate iron scavenging (Tagliabue et al. 2017). Here, we use a simple scavenging model where scavenging is proportional to the downward flux of biogenic particles and the concentration of dissolved ferric iron that is not complexed by organic ligands. The model has two adjustable parameters: a scavenging constant (k SC ) and a uniform concentration of organic ligands (C L ) forming stable iron complexes that are not removed by scavenging (Gledhill and Buck 2012). The parameter values are estimated by fitting nutrient and oxygen distributions calculated by the model to the corresponding concentrations measured in the study area (Supplementary Information, Section S2). The model is run into steady state applying different k SC -C L combinations covering the parameter space. The scavenging parameter has a strong effect on dissolved nutrients, oxygen and export production since phytoplankton growth in the OMZ is controlled by the availability of dissolved iron (Fig. 2). Dissolved iron concentrations decline when scavenging is enhanced by an increase in the scavenging constant and a decrease in the organic ligand concentration. Export production decreases when iron is removed by scavenging while dissolved oxygen in intermediate waters increases due to a corresponding decline in respiration. The model sensitivity is governed by a positive feedback loop rooted in the benthic iron release and a negative feedback provided by iron scavenging. The best fit to the observations is attained applying k SC = 0.7 µM −1 and C L = 0.8 nM (Fig. 2).
Observations employed for optimization are dissolved oxygen, nitrate and phosphate concentrations taken from the World Ocean Atlas 2018 (mean values averaged over the period 1955-2014) provided by NOAA (World Ocean Atlas 2018 (noaa.gov)) and concentrations of dissolved oxygen, nitrate,  Fig. 1). c Mean concentration of dissolved oxygen in upper intermediate water at the continental margin. d Export of particulate organic carbon (POC) across 50 m water depth integrated over the entire model area. Black dots represent combinations of k SC and C L values applied in the model simulations. The red circle indicates the k SC -C L combination yielding the best fit to the data (k SC = 0.7 µM −1 , C L = 0.8 nM)

3
Vol.: (0123456789) phosphate, nitrite and iron measured during the U.S. GEOTRACES cruise GP16 (Moffett and German 2018;Schlitzer et al. 2018) and SFB 754 cruises RV Meteor 135-138 (Krahmann et al. 2021). A good fit is obtained for all dissolved species applying the best fit k SC and C L values ( Supplementary Information, Fig.  S6). Moreover, all pelagic turnover rates and benthic fluxes obtained in the standard model run are consistent with rates and fluxes that have been measured in the study area (Tables S3 and S4). Tracer distributions obtained in the standard model run are shown in Supplementary Figures S8-S13.
In the following sections, we focus on N and Fe rather than P cycling since the phosphate-rich Peru-Chile Undercurrent and anoxic margin sediments supply large amounts of dissolved phosphate such that P-limitation does not occur in our study area (Fig. S7).

Nitrogen cycling
Steady state nitrogen fluxes obtained in the standard case (k SC = 0.7 µM −1 and C L = 0.8 nM) show that dissolved reactive nitrogen (N R = NO 3 + NO 2 + NH 4 ) is rapidly cycled within the Peruvian OMZ (Fig. 3). All nitrogen fluxes presented in Fig. 3 are dynamically calculated in the model using kinetic rate laws presented in the Supplementary Information (Section S1.2). Only atmospheric N R deposition is externally prescribed employing previous model estimates (Su et al. 2016). Pelagic rates and benthic fluxes calculated in the model are consistent with a large number of measurements that have been carried out in the area (Supplementary Information, Tables S3 and S4). In the following text, numbers in brackets refer to processes listed in Fig. 3. Internal N cycling is driven by net primary production (NPP) of particulate organic nitrogen (PON) by diazotrophs (2, in Fig. 3) and other phytoplankton (1) where the contribution of diazotrophs to NPP is almost neglectable as previously indicated by rate measurements that have been conducted in the study area (Chang et al. 2019;Selden et al. 2019). PON is rapidly degraded in the water column (4). Only 2% of the PON reaches the seabed where most of the PON is subsequently degraded (14) while a small fraction is permanently buried in in Gmol N year −1 . Processes: 1. Uptake of reactive nitrogen (N R = NO 3 + NO 2 + NH 4 ) during primary production of particulate organic nitrogen (PON). 2. Autotrophic nitrogen fixation. 3. Atmospheric deposition of reactive nitrogen. 4. Ammonium release during pelagic PON degradation. 5. Aerobic oxida-tion of ammonium to nitrite. 6. Aerobic oxidation of nitrite to nitrate. 7. First step of pelagic denitrification. 8. Second step of pelagic denitrification. 9. Heterotrophic nitrogen fixation. 10. Anaerobic oxidation of ammonium (anammox). 11. Oxidation of dissolved ferrous iron and sulfide by nitrate and nitrite. 12. Benthic dissimilatory nitrate reduction to ammonium (DNRA). 13. Benthic denitrification. 14. Benthic ammonium release induced by POM degradation in surface sediments. 15. PON burial in sediments. 16. Net flux of reactive N into the model across outer model boundaries sediments (15). The PON turnover is balanced by a total NPP of 12,565 Gmol year −1 (1 + 2), total degradation rate of 12,497 Gmol year −1 (4 + 14) and a PON burial flux of 68 Gmol year −1 (15).
Ammonium is produced by pelagic and benthic PON degradation (4 and 14), heterotrophic nitrogen fixation (9) and benthic dissimilatory nitrate reduction to ammonium (DNRA,12) where pelagic PON degradation is by far the largest ammonium source. Ammonium is consumed by aerobic (5) and anaerobic oxidation (10, anammox) with aerobic oxidation being clearly the dominant process. Nitrite is produced by aerobic ammonium oxidation (5) and by the first step of pelagic denitrification (7). Nitrite is consumed by aerobic oxidation to nitrate (6), reduction to N 2 via anammox (10) and the second step of denitrification (8). Nitrate is formed by aerobic nitrite oxidation (6) and consumed by pelagic and benthic denitrification (7 and 13) and benthic DNRA (12). Nitrate and nitrite turnover are clearly dominated by aerobic pathways.
Input of N R to the model domain is provided by autotropic and heterotrophic nitrogen fixation (2 and 9), atmospheric deposition (3) and the net horizontal influx of N R that is provided by nitrate-rich pelagic water masses entering the model domain across the outer boundaries (16). The lateral N R influx (16) is by far the largest N R source. Heterotrophic nitrogen fixation is included in the model (Supplementary Information, Eq. S11) since the available data indicate that most of the nitrogen fixation in the Peruvian OMZ occurs in sub-surface waters where light is not available (Bonnet et al. 2013;Fernandez et al. 2011). Heterotrophs use organic matter instead of light as energy source (Bonnet et al. 2013;Loescher et al. 2014). In the model, the biomass produced by heterotrophs is not conserved but completely degraded such that heterotrophic nitrogen fixation serves as a source for dissolved ammonium (9). N R is removed by pelagic and benthic denitrification (8 and 13), PON burial (15), pelagic anammox (10) and by the oxidation of dissolved ferrous iron and sulfide where nitrite and nitrate are reduced to N 2 (11). Pelagic N 2 production (8, 10, 11) accounts for 67% of the total N R removal while benthic denitrification (13) and burial of PON (15) contribute the remaining 33%. The external N R cycle is balanced with a total N R removal rate (8 + 10 + 11 + 13 + 15) and N R input (2 + 3 + 9 + 16) of 666 Gmol year −1 .
The rate of annamox (10) clearly exceeds the rate of nitrite removal via denitrification (second step of pelagic denitrification, 8) in our model. Experimental studies show that the second step of pelagic denitrification is less oxygen tolerant than the first step and annamox (Dalsgaard et al. 2014;Kalvelage et al. 2011). This oxygen sensitivity is considered in the model and effectively suppresses the second denitrification step (Supplementary Information, Eqs. S8 and S9, Table S2). Our model is consistent with a large set of rate measurements conducted in the study area that show wide-spread annamox and no detectable nitrite removal via denitrification (Hamersley et al. 2007;Kalvelage et al. 2011Kalvelage et al. , 2013. It should, however, be noted that nitrite reduction to N 2 via denitrification was observed at a few stations in the Peruvian OMZ (Dalsgaard et al. 2012).
The dominance of annamox over nitrite removal via denitrification has been challenged using stoichiometric constraints. According to this perspective, denitrification should dominate over anammox when most organic matter degradation proceeds via denitrification (Babbin et al. 2014;Dalsgaard et al. 2012;Koeve and Kahler 2010;Ward 2013). However, this reasoning does not apply to our study area. The model clearly shows that most of the organic matter produced in the Peruvian OMZ is degraded aerobically. Only the first step of denitrification (NO 3 → NO 2 ) contributes significantly to organic matter degradation while the second step (NO 2 → N 2 ) and dissimilatory sulfate reduction are almost completely suppressed (Table 1). Aerobic PON degradation supplies sufficient ammonium to fuel annamox while nitrite removal via denitrification proceeds at a low rate since the oxygen level in the core of the OMZ does not reach sub-micromolar levels in our model. The ROMS-derived water fluxes applied in our model (Supplementary Information, Figs. S1-S4, Table S1) generate horizontal and vertical oxygen fluxes into the low-oxygen region that are sufficiently high to support aerobic respiration throughout the OMZ.
Numbers indicate depth-integrated rates in mmol PON m −2 year −1 averaged over the respective model areas. Numbers in brackets indicate the contribution of different pathways of PON degradation to total PON degradation at 50-300 m water depth.
A growing set of observations indicates that aerobic respiration is an important pathway for organic matter degradation in many OMZs (Garcia-Robledo et al. 2017;Zakem et al. 2020). Metagenomics and metatranscriptomics suggest widespread aerobic metabolic potential and activity in OMZs (Kalvelage et al. 2015;Stewart et al. 2012) while rate measurements revealed that aerobic degradation of organic matter occurs throughout the Peruvian OMZ, even at nanomolar oxygen concentrations, and accounts for most of the organic matter degradation and ammonium production in our study area (Kalvelage et al. 2015). This perspective, that is based on direct rate measurements and metagenomics data, is supported by our mechanistic model (Table 1).

Iron cycling
Steady state iron fluxes obtained in the standard case are shown in Fig. 4.
Particulate iron inputs by terrigenous sediments (process 1 in Fig. 4) and atmospheric deposition (7) are prescribed employing measurements carried out in the Peruvian OMZ (Kadko et al. 2020;Scholz et al. 2014). All other iron fluxes are dynamically calculated employing the network of biogeochemical processes implemented in the model (Supplementary Information, Section S2.1). Most of the dissolved iron present in the Peruvian OMZ is released from shelf and slope sediments that are covered by oxygen-depleted bottom waters (2) as previously shown by a range of measurements conducted in the Peruvian OMZ Chever et al. 2015;DiTullio et al. 2005;Hutchins et al. 2002;John et al. 2018;Noffke et al. 2012;Plass et al. 2020;Rapp et al. 2020;Schlosser et al. 2018b;Scholz et al. 2016). The deposition of terrigenous iron at the seabed (1) is thus the ultimate source of new dissolved iron in the Peruvian system. The benthic iron release calculated in the model is proportional to the degradation rate of organic matter in surface sediments and increases when ambient bottom waters are deoxygenated (supplementary information, Eq. S20). The resulting benthic iron release (2) is consistent with iron fluxes measured in the study area (Supplementary Information, Table S4). Dust deposition (7) provides additional particulate iron for the surface ocean (Kadko et al. 2020). In the model, 10% of the dust-bound iron is released into surface waters as dissolved ferric iron (8) to support primary production. This high percentage is chosen to allow for a maximum impact of dust deposition on dissolved iron and ocean productivity. However, even at this high value, iron release from dust (8) is too small to compete with the much larger benthic iron release (2). Terrigenous sediments are, hence, by far the most important source for dissolved Fe in the Peruvian OMZ. Ferrous iron released from sediments is rapidly oxidized by oxygen (3), nitrate (4) and nitrite (5). Microorganisms mediate ferrous iron oxidation by nitrate (Scholz et al. 2016) while oxidation by oxygen and nitrite proceeds abiotically (Heller et al. 2017;Millero et al. 1987). Due to the fast oxidation kinetics applied in the model (supplementary information, Table S2), ferrous iron accumulates in the water column only when oxygen, nitrate and nitrite are depleted.
Phytoplankton take up dissolved iron that reaches the surface ocean by upwelling, vertical mixing and dust deposition (6). Most of the dissolved ferric iron is, however, removed by scavenging onto particles (10) in sub-surface waters and sinks down to the seabed before it can be used by phytoplankton. Iron scavenging is assumed to be an irreversible process in the model while iron taken up by phytoplankton is released back into the water column during POM degradation (9). Dissolved ferric iron produced by ferrous iron oxidation is complexed by organic ligands.
The ligand used in the model (uniform concentration of 0.8 nM) is assumed to have a high binding strength such that dissolved ferric iron is complexed by the ligand until the ligand is saturated with iron. Complexed iron is not scavenged but can be taken up by phytoplankton (Gledhill and Buck 2012). The ligand, thus, inhibits scavenging and keeps iron in solution that can be used to support primary production. Modeled iron concentrations in the surface ocean (0-50 m) decrease from 0.76 nM at the continental margin to 0.07 nM at 90-100° W (supplementary information, Fig. S7). Since the ligand concentration exceeds these values, the entire dissolved iron pool in surface waters is organically complexed.
Most of the iron released from sediments and dust particles is ultimately buried in sediments (88%). However, a significant fraction (12%) is exported from the model area to the surrounding regions (12). This model result is supported by iron isotope data indicating that iron released from sediments is partly kept in solution and advected to the open ocean (Chever et al. 2015;John et al. 2018;Scholz et al. 2014). The highest iron concentrations (3.1 nM) are observed at the margin in upper intermediate water (50-300 m, Fig. 2b). However, iron concentrations at these water depths rapidly decline further offshore due to scavenging (supplementary information, Fig.  S13) and reach a mean value of 0.7 nM at 90-100° W such that iron exported across the outer model boundary at 100° W is bound in organic complexes. The horizontally exported iron may contribute to the formation of iron-rich manganese nodules in the adjacent Peru Basin (Bollhöfer et al. 1996) and support marine productivity in regions that are located further offshore at > 100° W Tagliabue et al. 2019).

Model limitations
Our box model has a very coarse spatial resolution and uses constant water fluxes representing the climatological mean averaged over a period of 20 years. It is, thus, not able to resolve interannual and seasonal variations in circulation, productivity and redox state. The distribution of dissolved oxygen predicted by the model is broadly consistent with data reported in WOA 2018 and oxygen measurements conducted during cruise GP16 of the GEOTRACES program and several SFB 754 cruises (supplementary information, Fig. S6). It is, thus, likely that our model adequately represents the mean oxygenation state of the Peruvian OMZ. However, time-series data recorded in the Peruvian OMZ show a strong temporal variability in dissolved oxygen (Cardich et al. 2019;Gutierrez et al. 2008;Ulloa et al. 2012). Hence, dissolved oxygen concentrations at the continental margin have varied between 0.5 and 25 µM at 200 m water depth over the last decades (Cardich et al. 2019). Oxygen concentrations calculated by the model for the margin boxes at 50-300 m water depth (8-17 µM) fall into this range. They are lower than corresponding values reported in WOA 2018 (30-36 µM) but higher than concentrations measured during GP16 and SFB754 cruises in the margin area (0-13 µM). It is thus possible that our model does not appropriately resolve oxygen sensitive processes (e.g. second step of denitrification) during low oxygen periods that prevailed during these cruises. High-resolution ocean models with dynamic circulation would be needed to resolve the timedependent changes in oxygen conditions. The new network of biogeochemical processes developed for our box model could be integrated in such models to allow for a better presentation of the spatial and temporal variability in the Peruvian OMZ.

Response to ocean deoxygenation
The decrease in dissolved oxygen concentrations observed in the eastern equatorial Pacific over the last decades significantly reduces oxygen concentrations in water masses entering the model across its outer boundaries (Schmidtko et al. 2017). To investigate this change and how future ocean deoxygenation may affect the biogeochemical system in the Peruvian OMZ, oxygen concentrations in water masses flowing laterally into the model across the western boundary are diminished by up to 25% and the model is run over a period of 150 years to reach a new steady state. Dynamic tracer concentrations are applied at the southern and northern boundary of the model located at 6° S and 18° S, respectively, since deoxygenation affects the entire eastern Pacific OMZ extending from approximately 40° N to 30° S along the Pacific coast of North and South America (Supplementary Information, Section S1.3). It is presently unknown to what degree autotrophic and heterotrophic nitrogen fixation (NF) may increase when reactive nitrogen (N R ) is depleted by a rise in anammox and denitrification. Hence, the model is run applying different NF sensitivities to N R loss (supplementary information, Eq. S1 and Eq. S11). Autotrophic NF is calculated in the model considering SST and nutrient limitation as well as the N R /PO 4 ratio in surface waters (r NP ) and the parameter e NF (Eq. S1). NF is proportional to where r O NP is the N R /PO 4 ratio observed in surface boxes of the modern OMZ and r M NP is the ratio predicted by the model. With e NF = 0, modeled NF values do not respond to a decline in r M NP , e NF = 1 induces a moderate NF increase when the modeled r NP ratio ( r M NP ) is smaller than r O NP and e NF = 2 induces a drastic NF increase in N R -depleted surface waters such that NF locally compensates for the N-loss under anoxic conditions.
The model experiments show that the productivity of the Peruvian OMZ increases significantly when the lateral oxygen input across the western boundary at 100°W is diminished by 1-2% ( Fig. 5a  and b). This response is caused by an increase in benthic iron release from shelf and slope sediments under low oxygen conditions that is further amplified by the intensification of the biological carbon pump (Fig. 6c, Eq. S20). Oxygen concentrations in intermediate waters (50-1000 m water depth) decrease drastically due to the rise in productivity and respiration. Hence, a 1% decline in oxygen supply across 100°W reduces the mean oxygen concentration at 50-1000 m water depth by 8% (Fig. 5g). An even larger decrease by 22% is observed at the margin (Fig. 6g). These drastic changes show that the oxygen level within the OMZ is highly sensitive with respect to ocean deoxygenation. Considering that the global ocean has lost about 2% of its oxygen content over the last decades (Schmidtko et al. 2017), the strong expansion and intensification of the Peruvian OMZ that has been observed over this period (Stramma et al. 2008) may be attributed to a decline in the lateral oxygen influx that led to an amplified oxygen loss within the OMZ. The positive feedback that induces these drastic OMZ changes is rooted in the increase in benthic iron release from shelf and slope sediments covered by oxygen-depleted bottom waters. Enhanced benthic iron fluxes promote productivity under iron limitation and amplify oxygen consumption via respiration which leads to a further increase in benthic iron release and deoxygenation in a self-amplifying positive feedback loop.
The system response to deoxygenation is less sensitive when the oxygen input declines by more than 2% since nitrogen limitation is initiated and leads to a decrease in productivity ( Fig. 5a and b). While the productivity in the entire model area is iron-limited under present day conditions (Fig. S7), nitrogen becomes the limiting nutrient in the offshore areas at > 250 km distance to the coast. This change first appears at the southern part of the model domain (areas 9-11 in Fig. 1) since the nitrate input across the southern boundary at 18°S is lower than the input across the other model boundaries (Fig. S10). The entire model area is nitrogen limited when the oxygen input is lowered by more than 20% (for e NF = 0). Nitrogen limitation is induced by an increase in anammox in sub-surface waters under low oxygen conditions (Fig. 5f). The first step of denitrification strongly increases and supplies nitrite for anammox (Fig. 5d) while the second step of nitrification proceeds at a much lower rate (Fig. 5e) such that most of the N R depletion is induced by anammox over the entire deoxygenation range explored in the model. Enhanced nitrate consumption by denitrification (first step, NO 3 → NO 2 ) induces a decline in nitrate concentrations (Fig. 5i) and increase in dissolved nitrite (Fig. 5h). Dissolved ammonium is almost completely consumed by aerobic and anaerobic oxidation such that the concentration is kept below 1 µM in the standard case over the entire model area (Fig.  S12) and in all other model runs (data not shown). N R values in surface waters decline drastically (Fig. 5k) since almost the entire N R pool is consumed by phytoplankton under nitrogen-limitation while the N R input by upwelling is diminished by anammox. Total dissolved iron concentrations in surface waters increase since nitrogen limitation inhibits iron drawdown by phytoplankton (Fig. 5l).
Under nitrogen-limitation, the OMZ shows a much lower sensitivity with respect to oxygen supply.
Hence, a decline in the oxygen input across the western boundary by 1% (e.g. from 9 to 10%) induces an oxygen drop of only 0.8 µM under nitrogen limitation while a decline by 4.2 µM is observed under iron-limitation. The diminished sensitivity is induced by negative feedbacks in the nitrogen cycle where deoxygenation results in a decline in productivity and respiration due to the intensification of anammox and other N R consuming processes (Fig. 3). In contrast to the open ocean areas, productivity at the continental margin is maintained at a high level up to a decline in the oxygen supply by 10% since N R is provided by intense coastal upwelling (Fig. 6a). A further decline in oxygen supply leads to a decrease in productivity due to nitrogen limitation. These trends do not change significantly when nitrogen fixation is allowed to increase moderately upon N R depletion (e NF = 1). Dissolved phosphate is maintained at a high concentration level in all model runs (data not shown) due to the large influx of phosphate across the outer boundaries (Fig. S9) and the decline in phosphorus burial under low oxygen conditions (Eq. S19). Hence, the productivity is either limited by nitrogen or iron whereas phosphorus limitation does not occur (Fig.  S7).
Productivity is maintained at a high level over the entire deoxygenation range explored in the model when a strong increase in nitrogen fixation compensates for N R losses (e NF = 2, Figs. 5c and 6b). Dissolved sulfide accumulates in intermediate waters at the continental margin and reaches micromolar concentrations when the oxygen supply is diminished by more than 15% (Fig. 6j). Dissolved sulfide is produced by pelagic and benthic sulfate reduction and consumed in the water column by oxidation using oxygen, nitrite and nitrate as electron acceptors (Table S2). At the continental margin, 40% of the total pelagic PON degradation at 50-1000 m water depth is provided by dissimilatory sulfate reduction when the oxygen supply is diminished by 25% (e NF = 2) whereas this PON degradation pathway contributes less than 1% under present day conditions ( Table 1).
The model runs indicate that the ongoing expansion of the Peruvian OMZ is initiated by a small decline in the lateral oxygen influx that is strongly amplified by benthic iron release within the OMZ. Future deoxygenation induces further oxygen losses in the Peruvian OMZ. These losses may, however, Fig. 5 Effects of ocean deoxygenation. Model results shown in panels a-l are mean values averaged over the entire model area. a Net primary production. b Export production at 50 m water depth. c Nitrogen fixation rate (0-1000 m, autotrophic + heterotrophic). d First step of denitrification, NO 3 → NO 2 . e Second step of denitrification, NO 2 → N 2 . f Anammox. g Oxygen (50-1000 m). h Nitrate (50-1000 m). i Nitrite (50-1000 m). j Iron (50-1000 m). k Reactive nitrogen (0-50 m). l Iron (0-50 m). e NF defines the response of nitrogen fixation to N R depletion ◂ induce nitrogen limitation which enhances the resilience of the system with respect to marine environmental change. Most of the organic matter is degraded via aerobic respiration and denitrification and the oxygen concentrations are maintained at a low micromolar level over a wide range of deoxygenation scenarios (Figs. 5g and 6g) as indicated by previous modeling studies (Canfield et al. 2019). Stabilization of the redox state is caused by negative feedbacks provided by the redox-dependent nitrogen cycle. A switch to sulfidic conditions is only possible if the rate of nitrogen fixation is dramatically enhanced.
Response to ocean warming Sea surface temperatures (SSTs) are enhanced to investigate how the biogeochemical system in the OMZ corresponds to ocean warming. Initial SST values for each of the subregions in Fig. 1 are defined using SST data from the World Ocean Atlas 2018. These values are increased in 0.25 °C steps and the model is run into a new steady state for each incremental SST increase. Since ocean deoxygenation and warming proceed simultaneously, oxygen concentrations in water masses entering the model domain across the western boundary are diminished as in the previous model runs. The oxygen content in the inflowing water is decreased by 4% for each degree of warming applying the historic relationship between ocean warming and deoxygenation observed over the last decades (Schmidtko et al. 2017). It should, however, be noted that ocean warming may induce even stronger oxygen losses on centennial to millennial time scales (Oschlies 2021). As outlined in the supplementary information, net primary production, autotrophic nitrogen fixation and PON degradation in the surface ocean increase with ocean warming (Eq. S1, S2 and S4). Moreover, the solubility of oxygen in surface waters depends on SST (García and Gordon 1992) such that the modeled oxygen exchange with the atmosphere responds to ocean warming . Changes in stratification and vertical mixing that may be induced by ocean warming are, however, not considered in the following model runs since our model is not able to resolve changes in physical circulation and mixing (Sect. 2.5). It is possible that productivity and ventilation will be reduced in the Peruvian OMZ if warming induces a significant decrease in vertical mixing. The biogeochemical consequences of these changes are not considered in the following simulations but could be explored using other more advanced ocean models. According to our model, coupled warming and reduction in the lateral oxygen supply has an even larger effect on dissolved oxygen than ocean deoxygenation alone. Hence, a warming by 0.25 °C accompanied by a decrease in the oxygen influx across the western boundary by 1% reduces the oxygen concentrations at 50-1000 m water depth by 9% averaged over the entire model area (Fig. 7g) and 25% at the continental margin (Fig. 7h). The mean oxygen concentration in water masses entering the model domain at 50-1000 m water depth across the western boundary is reduced by only 0.6 µM when a 1% decline is imposed whereas the oxygen concentrations in the model domain at the corresponding water depths drop by 4.7 µM over the entire model area and 5 µM at the margin ( Fig. 7g and 7h). This amplification is caused by the SST-induced decrease in oxygen solubility, rise in net primary production ( Fig. 7a and 7d) and the positive feedback provided by benthic iron release.
The switch from iron to nitrogen limitation is already initiated when SSTs increase by 0.5 °C. It induces a corresponding drop in primary and export production in the offshore regions ( Fig. 7a and b) whereas strong coastal upwelling provides sufficient N R at the margin up to an SST increase of about 3 °C ( Fig. 7d and e).
A moderate increase in nitrogen fixation (e NF = 1) mitigates the decrease in productivity in the offshore region but amplifies the productivity drop at the margin when warming exceeds 3 °C (12% deoxygenation). This surprising model behavior is caused by enhanced N R consumption in offshore water masses that leads to a decline in the horizontal N R flux to the margin. A strong increase in nitrogen fixation (e NF = 2) provides sufficient N R to compensate for Fig. 6 Effects of ocean deoxygenation at the Peruvian margin. Model results shown in panels a-l are mean values averaged over the margin area (0-250 km distance to the coast) a. Export production at 50 m water depth. b. Nitrogen fixation rate (0-1000 m). c. Benthic iron flux (0-1000 m). d. First step of denitrification, NO 3 → NO 2 . e. Second step of denitrification, NO 2 → N 2 . f. Anammox. g. Oxygen (50-1000 m). h. Nitrate + nitrite (50-1000 m). i. Iron (50-1000 m). j. Sulfide (50-1000 m). k. Reactive nitrogen (0-50 m). l. Iron (0-50 m). e NF defines the response of nitrogen fixation to N R depletion ◂ Fig. 7 Effects of ocean warming and deoxygenation. a Net primary production averaged over the entire area. b Export production at 50 m water depth averaged over the entire area. c Nitrogen fixation (0-1000 m) averaged over the entire area. d Net primary production averaged over the margin area (0-250 km distance to the coast). e Export production at 50 m water depth averaged over the margin area. f Nitrogen fixation rate (0-1000 m) averaged over the margin area. g Oxygen (50-1000 m) averaged over the entire model area. h Oxygen (50-1000 m) averaged over the margin area. i Sulfide (50-1000 m) averaged over the margin area the enhanced N R loss via anammox. It leads to an increase in productivity but induces a rise in dissolved sulfide at the margin where sulfide concentrations exceed 1 µM when a warming of 3 °C is imposed (Fig. 7i).
Ocean deoxygenation and warming may promote sulfidic conditions in sub-surface waters of the continental margin The model experiments indicate that sulfidic conditions are only attained if N R losses via anammox and denitrification are compensated by a massive increase in autotrophic and heterotrophic nitrogen fixation. Field data show that sulfidic conditions occur episodically at the inner shelf region of the Peruvian margin when the oxygen supply is reduced presumably by a decline in physical ventilation (Dugdale et al. 1977;Ohde 2018;Plass et al. 2020;Schlosser et al. 2018b;Scholz et al. 2016;Schunck et al. 2013). These conditions are, however, regionally and temporally restricted, and do not represent the chemistry of the upwelling zone in general (Canfield 2006). Dissolved iron concentrations and nitrogen fixation rates drastically increase during these sulfidic events (Fernandez et al. 2011;Loescher et al. 2014;Schlosser et al. 2018b;Scholz et al. 2016). Dissolved sulfide concentrations in stagnant subsurface water masses on the inner shelf can reach 20 µM while oxygen, nitrite and nitrate are almost completely depleted (Loescher et al. 2014). Depth-integrated nitrogen fixation exceeded 0.3 mol m −2 year −1 during one of these events (Loescher et al. 2014). The method used for this rate measurement underestimates the true value by a factor of up to 6 (Grosskopf et al. 2012;Mohr et al. 2010). The rate may thus have been as high as 2 mol m −2 year −1 which approaches the values predicted by the model for 2 °C warming and e NF = 2 (Fig. 7f). Further support for a pronounced increase in nitrogen fixation in anoxic systems is provided by data recorded in marine sediments deposited during ocean anoxic events (OAEs) that occurred repeatedly in the geological past. They are typically enriched in biomarkers indicative for diazotrophs (Kuypers et al. 2004) and show light PON-δ 15 N values approaching the atmospheric value (Baroni et al. 2015;Luo et al. 2011;Scholz et al. 2019). Moreover, models demonstrate that the high marine productivity that prevails during OAEs can only be maintained when sufficient reactive nitrogen is supplied by elevated rates of nitrogen fixation (Baroni et al. 2014;Canfield 2006;Jurikova et al. 2020;Wallmann et al. 2019). Considering these findings, it seems to be possible that future ocean warming and deoxygenation may ultimately induce a permanent switch to sulfidic conditions across the entire Peruvian margin. However, whether the future Peruvian OMZ will ever attain this state depicted in Fig. 8 remains unclear since pathways and mechanisms of autotrophic and heterotrophic nitrogen fixation are still poorly understood (Zehr and Capone 2020).
If the entire margin would indeed attain a sulfidic state (Fig. 8), such a change would have massive consequences for fisheries off Peru since sulfide is a toxic chemical that significantly reduces the habitat available for fish and may ultimately lead to fish poisoning (Currie et al. 2018;Hamukuaya et al. 1998). The Peruvian margin produces more fish per unit area than any other region in the world ocean and accounts for up to 10% of the global fish catch (Chavez et al. 2008). A sulfide-induced decline in fishery yields could thus have negative impacts on global food security.
Anecdotal evidence suggests that strong sulfidic events that occurred in the past led to hydrogen sulfide release into the atmosphere along the Peruvian coast (Ohde 2018). A further spread and intensification of such events may, thus, affect not only the economic well-being but also the health of coastal communities.

Effects of biogeochemical feedbacks on future ocean change
The modeling results clearly show that the effects of ongoing ocean warming and deoxygenation are strongly amplified by a positive feedback loop (Fig. 9). This loop is active at continental margins where marine productivity is limited by iron and most of the iron is provided by benthic iron release. These areas will likely experience a strong decline in oxygen and rise in productivity over the coming decades driven by enhanced benthic iron release. This trend will continue until the coeval loss in reactive nitrogen induces a transition from iron-to nitrogen-limited conditions. The predicted transition to nitrogen limitation may induce a decline in productivity that mitigates further changes in redox state (Canfield 2006).
However, an increase in nitrogen fixation under low oxygen conditions would weaken this negative feedback such that the positive feedback provided by iron cycling may play out and induce further deoxygenation and productivity gains until the margins attain an euxinic (sulfidic) state (Fig. 9). In this scenario, enhanced benthic fluxes provide dissolved iron to allow for a strong increase in nitrogen fixation. Reactive nitrogen produced by the breakdown of diazotroph biomass then mitigates nitrogen limitation for non-diazotrophs and allows for an increase in overall productivity. These drastic changes may by further promoted by phosphorus cycling that provides positive feedback if P is a limiting nutrient for e.g. nitrogen fixation since benthic phosphate fluxes are enhanced under anoxic conditions (Ingall et al. 1993;Wallmann 2003).
Feedbacks identified in the Peruvian OMZ may also play out at other continental margins. The dissolved iron distribution at the continental margin of Mauretania is largely controlled by benthic iron release even though this margin receives large amounts of dust-bound iron generated in the adjacent Sahara (Rapp et al. 2019). Experimental studies with sediments retrieved at the Mauritanian margin confirm that large amounts of dissolved ferrous iron are released from sediments when they are expose to oxygen-depleted bottom waters (Schroller-Lomnitz et al. 2019). It is, thus, possible that the positive iron feedback is not only driving the expansion and intensification of the Peruvian OMZ but may contribute to the dissolved oxygen decline observed in other eastern boundary upwelling systems and the connected oxygen minimum zones (Stramma et al. 2008) when primary production is limited by iron and/or phosphorus.
We also suggest that the benthic iron feedback may partly account for the gap between observed (Schmidtko et al. 2017) and predicted trends in global ocean deoxygenation (Oschlies et al. 2018). Iron release from continental margin sediments is probably the largest source of dissolved iron in the global ocean (Dale et al. 2015;Somes et al. 2021). Most of the iron is removed by scavenging before it reaches the pelagic realm (Somes et al. 2021). Our model results, however, show that a significant fraction is protected by organic complexation and exported into the open ocean. Organic ligand concentrations may increase with productivity since ligands are produced and released during the mineralization of organic matter (Birchill et al. 2017;Boyd and Ellwood 2010;Schlosser and Croot 2009). An increase in productivity and organic matter degradation would thus diminish iron scavenging and enhance dissolved iron concentrations (Fig. 2) such that iron could be transported further offshore. Organic ligands also affect the redox cycling of iron and are generally found to retard the oxidation rate of dissolved ferrous iron  such that ferrous iron released from shelf sediments may be transported further offshore than predicted by the standard kinetic rate law for ferrous ion oxidation applied in our model (Millero et al. 1987). Iron solubility also increases with ocean acidification that proceeds in parallel with deoxygenation and warming (Zhu et al. 2021).
Another process that may facilitate the offshore transport of iron is the lateral advection of particles ) that can be enriched in iron due to the luxury uptake of iron by phytoplankton on the iron-rich shelf (Schlosser et al. 2018a). Recent studies at the Congo River margin confirm that dissolved iron released from shelf sediments can be transported offshore by rapid lateral advection over distances of up to 1000 km (Vieira et al. 2020). Moreover, meso-scale Fig. 9 Effects of biogeochemical feedbacks on future ocean change. Under iron (Fe) and phosphorus (P) limitation, the external forcing provided by climate change is amplified by a positive feedback loop (+). This feedback will induce a decline in oxygen, rise in productivity (NPP) and transition to nitrogen (N) limitation. Negative feedback provided by the N cycle (−) may induce a decline in NPP mitigating further ocean change (dotted line). This negative feedback may be overcome if nitrogen fixation (NF) rises upon oxygen and reactive nitrogen depletion (broken line). With enhanced NF, positive Fe and P feedbacks drive a transition to anoxic conditions where pelagic organic matter degradation proceeds via anaerobic pathways (nitrogenous: denitrification, euxinic: sulfate reduction) rather than aerobic and microaerobic respiration eddies facilitating the offshore transport of oxygendepleted water masses are frequently formed at continental margins (Callbeck et al. 2018;Löscher et al. 2016;Thomsen et al. 2016). These eddies are widely observed to increase off-shelf iron transport (Aguilar-Islas et al. 2013;Browning et al. 2021;Lippiatt et al. 2011). They also enhance offshore nitrogen fixation rates (Löscher et al. 2016) probably due to elevated iron concentrations in the eddy. All of these processes could weaken the negative feedback provided by scavenging and facilitate iron offshore transport under enhanced productivity and reduced oxygen levels. It may, thus, be possible that a decline in oxygen concentrations of bottom waters at continental margins induces wide-spread iron fertilization and deoxygenation in iron-limited offshore regions of the global ocean if the off-shelf transport is facilitated by ocean eddies and other processes. Since the residence time of dissolved iron the global ocean may be as low as 7.5 years (Somes et al. 2021), changing iron fluxes may induce rapid and global scale changes in productivity and oxygen concentrations. However, a better mechanistic understanding of pelagic iron cycling is needed to evaluate this hypothesis (Tagliabue et al. 2017).

Conclusions, hypothesis and outlook
Our model results support previous studies showing that most of the N R consumption in the Peruvian OMZ occurs via anammox ). Our mechanistic model also confirms the hypothesis that aerobic PON degradation is the dominant PON degradation pathway throughout the OMZ and provides most of the dissolved ammonium that is needed to fuel anammox (Kalvelage et al. 2015). The modeling further supports previous studies showing that the productivity of the Peruvian OMZ region is presently limited by iron where most of the dissolved iron is provided by ferrous iron release from shelf and slope sediments covered by oxygen-depleted bottom waters (Browning et al. 2018;Bruland et al. 2005;Noffke et al. 2012;Plass et al. 2020).
We propose that the expansion and intensification of the Peruvian OMZ observed over the past decades (Stramma et al. 2008) has been induced by the deoxygenation of open ocean water masses that enter the Peruvian OMZ since our model results show that even a small decrease in oxygen concentrations at the western boundary of the OMZ induces a massive decline in dissolved oxygen within the OMZ due to the positive feedback provided by benthic iron release. It is possible that this mechanism also contributes to the deoxygenation observed in other upwelling regions in areas where iron is the limiting nutrient and most of the iron is provided by benthic release. The enhanced benthic iron release may partly account for the gap between observed (Schmidtko et al. 2017) and predicted trends in global ocean deoxygenation (Oschlies et al. 2018) since a significant fraction of the iron released at continental margins may be advected further offshore by ocean eddies and other transport processes and fertilize adjacent areas.
Our modeling strongly suggests that ongoing deoxygenation and warming will induce a switch from iron-to nitrogen-limited conditions in major parts of the Peruvian OMZ, possibly over the coming decades. A transition to nitrogen limitation may also occur at other iron-limited continental margins.
The modeling confirms that oxygen concentrations are stabilized under nitrogen limitation due to the negative feedback provided by anammox and the resulting decline in productivity (Naafs et al. 2019). The transition to nitrogen-limited productivity will enhance the resilience of the system with respect to future ocean change. Negative feedbacks related to the nitrogen cycle are difficult to overcome and our modeling shows that a drastic increase in nitrogen fixation is required to allow for a transition into a sulfidic state. If this transition occurs, it would not only affect biogeochemical cycling but compromise fishery yields in the Peruvian OMZ that currently produces more fish per unit area than any other region in the world ocean.
Our modeling results confirm that a better understanding of nitrogen fixation (Zehr and Capone 2020) and pelagic iron cycling (Tagliabue et al. 2017) is urgently required to improve forecasts on ocean productivity, health and redox state.

Methods
We use the REDBIO model to simulate nutrient, oxygen and sulfur turnover in the Peruvian OMZ. This model has been designed to resolve redox-dependent feedbacks in pelagic and benthic nutrient cycling . The model was extended and revised to consider additional processes in the nitrogen cycle (anammox, DNRA, heterotrophic nitrogen fixation, first and second steps of denitrification) and improve the simulation of pelagic iron turnover. It simulates the distribution of the following dissolved tracers: oxygen, phosphate, nitrate, nitrite, ammonium, ferrous iron, ferric iron and sulfide. It considers the turnover of particulate organic matter in the water column and in surface sediments, a suite of organic matter degradation pathways and a wide range of secondary redox processes. The model employs annual mean water fluxes derived from ROMS simulations of the Peruvian OMZ (José et al. 2021). The model uses a geochemical approach and has, therefore, a number of limitations: It does not simulate the abundance of biota and the transport of particles. It considers the limitation of primary production by phosphorus, reactive nitrogen and iron but does not simulate the effects of light availability on primary production in surface waters. It has a coarse spatial resolution with only four depth levels in the water column (0-50 m, 30-300 m, 300-1000 m, > 1000 m) and 12 areas representing the entire OMZ (Fig. 1). It uses a timeinvariant circulation field and is, therefore, not able to resolve seasonal and interannual variability. A full description of the model is given in Section S1 of the Supplementary Information.
Funding Open Access funding enabled and organized by Projekt DEAL. Funding was provided by the Deutsche Forschungsgemeinschaft (DFG) via the Collaborative Research Centre 754 and the Emmy Noether project ICONOX.
Data availability All data that are used in this publication are taken from publicly available sources.

Code availability
The full code of the model presented in the publication will be provided by the first author on demand.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.