How does the phytoplankton-light feedback affect the marine N 2 O inventory? 1

Abstract

to explain the sign of the final heat perturbation: either an indirect dynamical response (Murtugudde et al., 2002;Löptien et al., 2009) or a direct thermal effect (Mignot et al., 2013;Hernandez et al., 2017).Hernandez et al. (2017) further distinguished a local from a remote thermal effect by highlighting the important role played by the advection of offshore CHLinduced cold anomalies in the Benguela upwelling waters.The interplay of these mechanisms is regionally variable (Park et al., 2014).Despite the diversity of modelled responses, a consensus emerges on the first order effect of PLF on the ocean physics, which is to perturb the ocean thermal structure (Nakamoto et al., 2001;Murtuggude et al., 2002;Oschlies, 2004;Manizza et al., 2005Manizza et al., , 2008;;Anderson et al., 2007;Lengaigne et al., 2007;Gnanadesikan and Anderson, 2009;Löptien et al., 2009;Patara et al., 2012;Mignot et al., 2013;Hernandez et al., 2017).By trapping more heat at the ocean surface in eutrophic regions, such as coastal or equatorial upwellings areas, the presence of phytoplankton initially increases the surface warming.Confining heat at the surface leads to less heat penetrating in subsurface.In some cases, the advection and upwelling of subsurface cold anomalies can lead to remote cooling effects (Hernandez et al., 2017;Echevin et al., 2022).Dynamical readjustment in response to perturbations in thermal structure has also been shown to have a cooling effect, by increasing upwelling of cold water to the ocean surface (Manizza et al. 2005;Marzeion et al., 2005;Nakamoto et al., 2001;Löptien et al., 2009;Lengaigne et al., 2007;Park et al., 2014).Because these effects depend on upper ocean stratification, an important role is attributed to modelled seasonal deepening of the mixed layer as it determines the intensity of the underlying temperature anomaly and its vertical movement to the surface.In other terms, whatever the temporality of the causal chain, changes in the PLF representation are expected to both perturb the ocean heat uptake, and trigger perturbations of both the water column stratification and associated ocean dynamics.b) This study: implications for N 2 O budget uncertainties Nitrous oxide (N2O) is a major ozone-depleting substance (Ravishankara et al., 2009;Freing et al., 2012) and a potent greenhouse gas, whose global warming potential is 265-298 times that of CO2 for a 100-year timescale (Myhre et al., 2013).The spatial coherence between marine productive areas and observed hot-spots of N2O production leads to question the impact of an incomplete representation of the PLF on the simulated N2O inventory.Recent observational studies highlight that N2O production is high in low-oxygen tropical regions and cold upwelling waters (Arévalo-Martinez et al. 2018;2020;Yang et al., 2020;Wilson et al., 2020).N2O becomes increasingly saturated in surface waters of equatorial upwelling regions due to the upward advection of N2O-rich waters (Arévalo-Martínez et al., 2017).Regions known to account for the most productive areas of the ocean spatially coincide with highest N2O production: 64% of the annual N2O flux occurs in the tropics, and 20% in coastal upwelling systems that occupy less than 3% of the ocean area (Yang et al., 2020).Despite recent advances, a large range of uncertainties still surrounds oceanic N2O emissions as large areas of both the open and coastal ocean remain undersampled by observations (Wilson et al., 2020).In particular, the paucity of observational data over key source regions contributes to increase uncertainties.The recent global budget of Tian et al. (2020) estimates natural sources from soils and oceans to contribute with up to 57% to the total N2O emissions between 2007 and 2016, with the ocean flux reaching 3.4 (2.5−4.3)Tg N yr -1 .A large uncertainty range is associated to the ocean flux estimate, as it is based on outputs from only a small number of global ocean-biogeochemical models.Very few climate models, even in the current CMIP6 generation, include emissions (and beforehand a complete representation of N cycling) of N2O fluxes: only 4 out of the 26 Earth system models considered in Séférian et al. (2020) simulate marine N2O emissions.The last generation of Earth system models projects an enhanced ocean warming in response to climate change, which is in turn expected to increase upper-ocean stratification (Sallée et al., 2021) and to contribute to greater reductions in upper-ocean nitrate and subsurface oxygen ventilation (Kwiatkowski et al., 2020).Ocean warming and deoxygenation constitute two triggers of high-probability high-impact climate tipping points (Heinze et al., 2021) and are identified as two of the main environmental factors influencing marine N2O distributions (IPCC, 2019;Hutchins and Capone, 2022).Through its expected impacts on the upper ocean stratification, the PLF representation could further change the oceanic N2O source by modulating the mixing between N2O-rich water and intermediate depths, perturbing the way N2O-rich water reaches the air-sea interface (Freing et al., 2012).
Here we investigate how an incomplete representation of the PLF leads to uncertainties in N2O projection in an up-to-date global ocean-biogeochemical model making up the current generation of Earth system models.Section 2 describes the numerical model and the set of simulations, as well as the existing options to consider CHL modulations of the incoming shortwave radiation.Section 3 presents the effect of an interactive PLF on the ocean heat content, associated ocean stratification and dynamics, and its feedback on marine N2O inventory.Finally, Section 4 summarizes the main results, addresses their broader implications, and discusses the future work motivated by this study.

Methodology a) Configuration of the global ocean-biogeochemical model
Recent projections of future N2O emissions contributing to intercomparison projects like CMIP6 are still based on Earth system models with a low spatial resolution (Séférian et al., 2020).For sake of coherence with CMIP biogeochemical modelling efforts, in the following we use a global ocean-biogeochemical configuration of the NEMO-PISCESv2 model (Madec, 2008;Aumont et al., 2015) at 1° of horizontal resolution.This model corresponds to the oceanic component of CNRM-ESM2-1 (Séférian et al., 2019) and is one of the few CMIP6-class models that contributed to the Global N2O budget (Tian et al., 2020).Our modelled ocean has 75 vertical levels and the first level is at 0.5 meter depth.Vertical levels are unevenly spaced with 35 levels being in the first 300 meters of depth.Atmospheric forcings of momentum, incoming radiation, temperature, humidity, and freshwater are provided to the ocean surface by bulk formulae following Large and Yeager (2009).Details on physical configuration are given in Berthet et al. (2019).Using an ocean-only configuration allows to isolate the local response induced by the PLF by not confounding it with potential inter-basin feedbacks acting through the atmosphere.JRA55-do atmospheric reanalysis (Tsujino et al., 2018;Tsujino et al., 2020) provided the atmospheric forcings of the ocean.The global domain was first spun-up under preindustrial conditions during several hundred years ensuring that all fields approached a quasi-steady state.The historical evolution of atmospheric CO2 and N2O concentrations was prescribed since 1850.To avoid the warming jump between the end of the spin-up and the onset of the reanalyses in 1958, the first 5 years of JRA55-do forcings were cycled, followed by the complete period of JRA55-do atmospheric forcing from 1958 to 2018.b) Experimental design: three representations of the PLF The control simulation (hereafter REF) together with the spin-up both account for a fully interactive PLF: the penetration of shortwave radiation into the ocean surface is constrained by the CHL concentration ([CHL]) produced by the PISCESv2 biogeochemical component (Figure S1, REF).PISCESv2 (Pelagic Interactions Scheme for Carbon and Ecosystem Studies v2) is a 3D biogeochemical model which simulates the lower trophic levels of marine ecosystems (nanophytoplankton, diatoms, microzooplankton and mesozooplankton), the biogeochemical cycles of carbon and of the main nutrients (phosphate, nitrogen, iron, and silicate) along the 75 levels of our numerical ocean.A comprehensive presentation of the model is found in Aumont et al. (2015).PISCESv2 simulates prognostic 3D distributions of nanophytoplankton and diatom concentrations.The evolution of phytoplankton biomasses is the net outcome of growth, mortality, aggregation and grazing by zooplankton.Growth rate of phytoplankton mainly depends on the length of the day, depths of the mixed layer and of the euphotic zone, the mean residence time of the cells within the unlit part of the mixed layer and includes a generic temperature dependency (Eppley, 1972).Nanophytoplankton growth depends on the external nutrient concentrations in nitrogen and phosphate (Monod-like parameterizations of N and P limitations), and on Fe limitation which is modeled according to a classical quota approach.The production terms for diatoms are defined as for nanophytoplankton, except that the limitation terms also include silicate.Light absorption by phytoplankton depends on the waveband and on the species (Bricaud et al., 1995).A simplified formulation of light absorption by the ocean is used in our experiments to calculate both the phytoplankton light limitation in PISCESv2 and the oceanic heating rate (Lengaigne et al., 2007).In this formulation, visible light is split into three wavebands: blue (400-500 nm), green (500-600 nm) and red (600-700 nm); for each waveband, the CHLdependent attenuation coefficients,  " ,  $ and  & , are derived from the formulation proposed in Morel and Maritorena (2001): (1) where WLB means the wavelength band associated to red (R), green (G) or blue (B), and bounded by the wavelengths  ; and  < as detailed above.() is the attenuation coefficient for optically pure sea water.() and () are fitted coefficients which allows to determine the attenuation coefficients due to chlorophyll pigments in sea water (Morel and Maritorena, 2001).Figure 1: Modelled tropical [35°S-35°N] heat content of upper 300 m (OHC300; in ZJ) for each simulation described in Table 1: REF (black; empty stars), climZCST (green; full downward triangles) and climZVAR (blue; empty rightward triangles).In (a) final part of the spin-up has been added in gray to illustrate the branching protocol in year 1999, and OHC300 anomalies have been computed with respect to year 1999.Subplot (b) zooms over the Argo period to compare modelled tropical OHC300 anomalies with 3 in situ-based products (see section 2c).
At year 1999 two sensitivity experiments were branched off (Figure 1).Both simulations climZCST and climZVAR account for an incomplete and external PLF, as they consider an observed climatology of surface [CHL] from ESACCI (Valente et al., 2016) in order to compute the light penetration into sea water (Equation 1; Figure S1).These two simulations differ from each other by the "realism" of the vertical profile derived in each grid point from the surface value of the ESACCI CHL climatology to the level of light extinction (Table 1).climZCST uses constant profiles of CHL spreading uniformly in the vertical direction (  Model results are compared with available observational-based gridded temperature and salinity datasets.Ocean heat content (OHC) of the upper 0-300 meter layer was inferred from three different products: i) the global objective analysis of subsurface temperature EN4 (Good et al., 2013), ii) the SIO product of the Scripps Institution of Oceanography (Roemmich anf Gilson, 2009), and iii) the ISAS20 optimal interpolation product released by Ifremer (Kolodziejczyk et al., 2019;Kolodziejczyk et al., 2021).While the SIO and ISAS20 products consider only Argo temperature and salinity profiles, the EN4 dataset considers all types of in situ profiles providing temperature and salinity (when available).These three in situ-based datasets are considered since 2005, the year the Argo coverage became sufficient to characterize the global ocean.Details on OHC computation are given in Llovel and Terray (2016) and Llovel et al. (2022).The authors also refer to cross-validations of OHC of deeper layers (0-700 m and 0-2000 m) against OHC anomalies from World Ocean Atlas 2009 (Levitus et al., 2012).A monthly climatology  of oceanic temperature from World Ocean Atlas 2013 version 2 (Locarnini et al., 2013) was used to evaluate modelled temperatures.Modelled O2 was compared to the annual climatology of O2 from World Ocean Atlas 2013 (Garcia et al., 2014) and modelled CHL was compared to the 3D monthly climatological global product estimated from merged satellite and hydrological data of Uitz et al. (2006).Modelled N2O partial pressure difference across the air-sea interface (Dpn2o) was compared to the recent dataset of Dpn2o observations compiled by Yang et al. (2020).
3. Results a) Impact of PLF on the upper ocean heat content and dynamics Meridional sections reveal that heat perturbations in response to changing CHL fields interacting with light are limited to the top 0-300 m layer of the ocean and predominantly affect the tropical area (Figure 3 and Figure S2, c-d).
The largest temperature anomalies are observed near the thermocline depth and reflect upper ocean warming and deepening of the thermocline in climZCST (Figure 3c), and cooling and shallowing of the thermocline in climZVAR (Figure 3d).In climZCST the ocean warming reflects large-scale patterns of a tropical CHL deficit compared to REF (Figure 2, b).Temperature differences are lower in the near-surface layer (0-50 m) than in the 50-300 m layer.This is expected as a result from weak stratification but also from simulations run with a forced atmosphere in which the temperature of the ocean surface layer is constrained by the atmospheric prescribed state.
When using an incomplete representation of the PLF, two contrasting trends of the upper ocean heat content (OHC) emerge compared to our control run REF (Figure 1a).Over the Argo period (2005-present) EN4 estimates of tropical OHC300 are in very good agreement with our warmest simulation climZCST (Figure 1b), while the two other dataproducts SIO and ISAS20 are in better agreement with our control run REF and with climZVAR.The good accordance between modeled OHC300 and observations is not a systematic feature of model-data comparisons (Cheng et al., 2016;Liao et al., 2022).Moreover, non-negligible differences exist among OHC dataproducts which are generally particularly strong in the upper 0-300 m layer (Lyman et al., 2010;Liang et al., 2021).The spread between these products at the end of the 2005-2018 period (12.1 10 21 J) is comparable to that of our numerical set (13.6 10 21 J).tropical domain in climZVAR compared to REF (Figure 2, c).As a result, the energy associated with the incoming radiation is caught in surface waters without being distributed over the water column.
In both climZCST and climZVAR the subsurface temperature anomaly deepens progressively over the first six years of simulation as a result of vertical mixing (Figure S3).This evolution indicates that part of the OHC300 differences between simulations comes from the adjustment of climZCST and climZVAR to the spin-up mean state yielded by an interactive PLF.It can be expected that experiments having spin-ups run with different representations of the PLF, would give even stronger sensitivities than those highlighted in this study.The sensitivities of OHC300 to the PLF formulation evaluated here should be considered at the lower end of estimate of OHC discrepancies that may emerge from changing the PLF representation.Prescribing a constant vertical profile of CHL (climZCST) to compute the penetration of the radiation into the ocean increases the OHC300 by more than 20 10 21 J during the last two decades (1999-2018) compared to REF (Figure 1).This rise of OHC300 decreases the verticallyintegrated tropical potential density of the upper 300 m at the end of the simulated period by 5 kg/m 2 compared to REF (Figure S4).The opposite trend (a reduced OHC300 compared to REF) is simulated with the same state-of-the-art CMIP6 ocean-biogeochemical model when considering a variable vertical profile of CHL (climZVAR).However Figure 1 highlights that the simulation using a consistent CHL for interacting with both incoming shortwave radiation and biogeochemical cyclings (REF) does not amplify one of these two trends, as climZCST and climZVAR surround REF.Average ranges of uncertainties associated with the PLF representation over the extended tropical domain (35°S-35°N) exceed 40 10 21 J in terms of OHC300 (Figure 1), 4 meters for the thermocline depth and more than 9 kg/m 2 for the vertically-integrated potential density perturbation (Figure S4).Similar to OHC300, ranges of uncertainty for the OHC estimates of deeper layers (0-700 m and 0-2000 m) also slightly exceed 40 10 21 J.Such uncertainty ranges are quite important as they are obtained by only changing the PLF representation in a single ocean-biogeochemical model.By comparison and in the context of OMIP protocols, Tsujino et al. (2020) give spreads between CMIP model estimates of the order of 50 10 21 J for the OHC of the upper 700m after 20 years (please refer to their Figure 24, a-b).Regarding the OHC integrated over the 0-2000m layer, they report an inter-model spread between 50 and 100 10 21 J, depending on the OMIP protocol considered (see their Figure 24, d-e).The OHC300 uncertainty of 40 10 21 J triggered by the representation of the PLF in our set of simulations has a comparable order of magnitude than the current spread of multi-model estimations of OHC.The present study suggests that part of the OHC multi-model uncertainty in current climate models may be due to different representations of the phytoplankton-light interaction.
The heat and associated density perturbations also cause dynamical modifications of upper ocean currents (Figure 4).Absolute differences in upper ocean velocities (average between 0 and 300m depth) are between |0.05| and |0.6| cm/s with strongest differences along the equator revealing perturbations of the equatorial undercurrent (Figure 4, b and c).Circulation around the subtropical gyres is also impacted, in particular for the South Pacific subtropical gyre.These modifications of zonal and meridional dynamics spread over the entire tropical latitudes, from 30°S to 30°N, strongly supporting the idea that heat perturbations induced by different interactions between CHL and incoming shortwave cause non-negligible modifications of the equatorial and tropical ocean dynamics.c).As stressed by Sweeney et al. (2005), small changes in CHL concentration (Figure S5) may have important effects on the mixed layer depth in these subtropical gyres due to low local wind speeds and low mixing conditions.This is thought to explain the large sensitivity we observe in terms of pycnocline depth (Figure 5) and ocean heat content in these regions.In line with their results, our set of simulations highlights that small CHL changes in low productivity regions trigger a vertical redistribution of density anomalies affecting the stratification.The relationship between N2O concentration and OHC300 in the Tropical Ocean is derived from a linear regression for each of the three 20-years simulations (Figure 6).The resulting slopes allow to identify three distinct tropical N2O production pathways along time as a function of the oceanic heat uptake: from 0.3 µmolN m -2 per ZJ for the most simplified PLF scenario climZCST, to 1 µmolN m -2 per ZJ for climZVAR.The slope of the simulation with the higher level of realism in terms of interactivity (REF) appears a solution between the two previous extremes, as it increases its N2O production by 0.8 µmolN m -2 per ZJ.Each of these N2O production pathways will translate into a different temporal evolution of the N2O budget and hence future climate.This result stresses the importance of having an interactive PLF in order to neither overestimate nor underestimate the N2O production projections due to a simplified representation of the PLF. Figure 6: Annual N2O inventory vertically integrated over the first 300 meters depth (µmolN/m 2 ) as a function of the annual OHC300 (ZJ) and averaged over an extended tropical domain (35°S-35°N).All points reflect anomalies compared to year 1999.c) Impacts on oceanic N2O emissions By perturbing the OHC, the ocean dynamics and the N2O production, the way PLF is modelled has non-negligible consequences on Dpn2o and thus on N2O emissions at the air-sea interface (Figure 7).Because the atmospheric partial pressure of N2O is identical among simulations, differences in Dpn2o are driven by changes in surface oceanic N2O concentration normalized by those in N2O solubility.Since solubility is mainly driven by temperature and because surface temperature anomalies are very weak (Figure S3, c and d), we do not expect solubility perturbations close to the surface.It results that spatial patterns of Dpn2o anomalies (Figure 7) reflect differences in surface oceanic N2O concentration.
Compared to a scenario considering a fully interactive PLF (REF), an incomplete representation of the PLF underestimates Dpn2o in all oxygen minimum zones of the northern hemisphere, which are strong emission zones (Figure 7, c and d).Large Dpn2o anomalies of -2.5 natm encompasses northern parts of the oxygen minimum zones of the Indian, Pacific and Atlantic oceans and anomalies reach up to -5 natm locally.Consequently, climZCST and climZVAR underestimate N2O fluxes by more than 12% in these oxygen minimum regions compared to REF.This result highlights that the representation of the PLF can be an important source of uncertainty in modelling N2O fluxes.As a matter of fact, the oceanic contribution to the recent global N2O budget by Tian et al. (2020) is based on only five global ocean-biogeochemical models (as still only few models simulate marine N2O emissions).These models have different configurations of the PLF which adds considerable uncertainty to simulated marine N2O emissions.
In subtropical gyres, the strong and direct effect of temperature (Figure S2, c and d) on in-depth N2O concentration (Figure 5, e and f) is in line with Yang et al. (2020) who demonstrate that the seasonality of Dpn2o in that regions is driven by a solubility regime.Both climZCST and climZVAR overestimate Dpn2o in subtropical gyres of the South Pacific and South Atlantic (Figure 7, c and d).This leads to an overestimation of the regional N2O fluxes by 24% compared to a simulation having a complete and interactive PLF representation (REF).Analyses are based on differences between a control run with an interactive PLF (REF) and two experiments with an incomplete PLF (climZCST and climZVAR) using a prescribed CHL climatology to interact with the incoming solar radiation.Changing the approach to compute how CHL filters the light penetrating into the ocean highlights the consequences of using an interactive PLF.
Our results demonstrate that the approach commonly used to account for the impact of the phytoplankton on light penetration significantly interfers with upper ocean heat uptake (Figure 1), the associated dynamics (Figure 4) and stratification in the tropics (Figure 5, a-c).Our set of forced ocean-biogeochemical simulations reveals that marine production of nitrous oxide (N2O) is sensitive to the representation of the PLF (Figure 5, d-f).The heat perturbations add to the uncertainty of modelled oceanic N2O production and result in three N2O production trajectories along time (Figure 6) that in turn trigger regional differences of Dpn2o and sea-air N2O fluxes (Figure 7).In forced ocean simulations, atmospheric forcings constrain surface temperature, salinity and thus solubility.However, the N2O concentration integrated over the upper 300 meters depth of the water column (Figure 5, e-f) showed differences with the control run that follow those of the in-depth temperature (Figure S2, c-d): in climZCST (climZVAR), a warmer (colder) tropical ocean leads to a decreased (an increased) N2O concentration.Because higher marine greenhouse gas emissions will increase the temperature of the coupled atmosphere-ocean system, adding an interactive atmospheric component is expected to amplify the PLF-induced mean changes in marine N2O concentration highlighted in this ocean-only numerical set (Park et al., 2014;Asselot et al., 2022).
Our results also question the reliability of current modelled estimates of the area and volume of oxygen minimum zones, as well as their trends in a future climate.The expansion rate of O2depleted waters still remains unclear and its controlling mechanisms are not yet fully understood and represented in today's models.Observation based assessments suggested that the ocean has already lost around 2% of the global marine oxygen since 1960 (Schmidtko et al., 2017).The expansion of oxygen minimum zones is expected to result in an increase of the volume of suboxic water and to have an impact on the production and decomposition of N2O (Freing et al., 2012).Our set of simulations highlights that an incomplete representation of the PLF underestimates the expansion of oxygen-depleted waters over the 20 years of simulation in comparison to REF.In climZCST and climZVAR the global volume (0-1000 m) of hypoxic water with [O2] under 50 mmol m -3 is up to 2.3 10 14 m 3 lower in 2018 compared to that of the control run REF.Thus an incomplete representation of the PLF might lead to an underestimation by 1.2 % of the modelled tropical volume of low-oxygenated waters after 20 years.
Recent regional studies demonstrated that the interactive PLF strongly affects upwelling systems of the South Pacific and Atlantic oceans (Hernandez et al., 2017;Echevin et al., 2021).Coastal upwellings are known to be sites of high N2O production with an annual N2O flux amounting to approximately 20% of the global fluxes while these systems occupy less than 3% of the ocean area (Yang et al., 2020).However, in the present study main modelled perturbations are rather localized over oxygen minimum zones or subtropical gyres (Figure 5; Figure 7).While the latter regional studies were performed at horizontal resolutions compatible with the complex dynamics of coastal upwellings (from 10 km to about 28 km), the resolution of climate models (~1-degree of horizontal resolution) does not allow to resolve these dynamics.A step further would be to evaluate how the sensitivity of N2O emission to the representation of the PLF depends on the horizontal resolution by running simulations at higher resolution with the same climate model.This would help to better determine how coastal upwelling systems may impact the modelled N2O inventory through different PLF representations, as well as the associated modelled range of uncertainty.
Figure1: Modelled tropical [35°S-35°N] heat content of upper 300 m (OHC300; in ZJ) for each simulation described in Table1: REF (black; empty stars), climZCST (green; full downward triangles) and climZVAR (blue; empty rightward triangles).In (a) final part of the spin-up has been added in gray to illustrate the branching protocol in year 1999, and OHC300 anomalies have been computed with respect to year 1999.Subplot (b) zooms over the Argo period to compare modelled tropical OHC300 anomalies with 3 in situ-based products (see section 2c).At year 1999 two sensitivity experiments were branched off (Figure1).Both simulations climZCST and climZVAR account for an incomplete and external PLF, as they consider an observed climatology of surface [CHL] from ESACCI(Valente et al., 2016) in order to compute the light penetration into sea water (Equation1; FigureS1).These two simulations differ from each other by the "realism" of the vertical profile derived in each grid point from the surface value of the ESACCI CHL climatology to the level of light extinction (Table1).climZCST uses constant profiles of CHL spreading uniformly in the vertical direction (Figure2, b and d-f).climZVAR uses variable vertical profiles computed followingMorel and Berthon (1989) (Figure2, c and d-f).This set of simulations is representative of the several configurations used in the case of CMIP intercomparison project.Table 1: Experimental set-up.Simulation Which CHL fields to interact with incoming shortwave radiation?PLF nature

Figure 2 :
Figure2: CHL concentration (mgCHL m -3 ) interacting with the incoming shortwave radiation for each numerical experiment (Table1).Maps a-c) show annual means of the vertical sum over 0-6000 m, a) as modelled over the 2009-2018 period for REF, and its differences with the external CHL prescribed for b) climZCST and c) climZVAR experiments.Labels P1 to P3 on subplot a) locate vertical profiles shown on subplots d-f).

Figure 4 :
Figure 4: Annual mean speed (color; cm/s) and streamlines of oceanic currents between 0-300 m over the 2009-2018 period for a) REF, and its differences with b) climZCST and c) climZVAR.In b-c) streamlines are colored when absolute speed are larger than 0.05 cm/s.
Figure 5: a-c) Depth of annual pycnocline (m) for 2009-2018 computed as the annual mean depth of the maximum of the Brunt-Väisälä frequency N 2 (T, S) over the water column (Maes and O Kane, 2014).d-f) Mean [N2O] (µmolN/m 3 ) over the first 300 meters depth.For REF (upper panel) and its mean-state differences with climZCST (middle panel) and climZVAR (bottom panel).Anomalies of N2O concentration integrated over the first 300 meters of the water column (Figure5, e and f) are in good agreement with patterns of pycnocline anomalies over the tropics (Figure5, b and c).These comparable spatial structures attest that N2O anomalies are driven by perturbations of stratification in large parts of the tropical domain.In the South Pacific subtropical gyre, the concomitance of i) an increased temperature (FigureS2, c and d), ii) a reinforced transport (Figure4, b and c) and iii) a weakened stratification illustrated by a local deepening of the pycnocline (Figure5, b and c), contributes to decrease the N2O concentration in both climZCST and climZVAR (Figure5, e and f).In contrast, in the South Indian Ocean and North tropical Atlantic the increase of N2O concentration seems to be mainly driven by the mean shoaling of the local pycnocline, as both regions exhibit contrasted perturbations in terms of transport and temperature.Finally, in the North Pacific oxygen minimum zone, the strong N2O deficits in both climZCST and climZVAR compared to REF do not respond to stratification and transport anomalies but are rather driven by a local rise of O2 concentration (FigureS6).Considering an incomplete PLF contributes to overestimate the oxygen concentration in this oxygen minimum zone and leads to a lack of local N2O production.The relationship between N2O concentration and OHC300 in the Tropical Ocean is derived from a linear regression for each of the three 20-years simulations (Figure6).The resulting slopes

Figure 7 :
Figure 7: Mean sea-to-air Dpn2o (natm) computed from a) observations, b) REF over the 2009-2018 period, and its differences with c) climZCST and d) climZVAR compared to REF.
Compared to an ocean model using a fully interactive PLF (REF), an incomplete PLF results in an overestimation of N2O fluxes by up to 24% in the South Pacific and South Atlantic subtropical gyres, and a reduction by up to 12% in oxygen minimum zones of the northern hemisphere.Our results based on a model at CMIP6 state-of-the-art emphasize an overlooked important source of uncertainty in climate projections of marine N2O production and in current estimations of the marine nitrous oxide budget.In subtropical gyres of the southern Hemisphere which are regions of low productivity, small CHL changes have a strong and direct effect on temperature (FigureS2, c and d), transport (Figure4, b and c) and local stratification (Figure5, b and c).These concomittant effects result in a local decrease of N2O concentrations in both experiments with a simplified PLF representation (climZCST and climZVAR).