Nitrogen isotopic constraints on nutrient transport to the upper ocean

Ocean circulation supplies the surface ocean with the nutrients that fuel global ocean productivity. However, the mechanisms and rates of water and nutrient transport from the deep ocean to the upper ocean are poorly known. Here, we use the nitrogen isotopic composition of nitrate to place observational constraints on nutrient transport from the Southern Ocean surface into the global pycnocline (roughly the upper 1.2 km), as opposed to directly from the deep ocean. We estimate that 62 ± 5% of the pycnocline nitrate and phosphate originate from the Southern Ocean. Mixing, as opposed to advection, accounts for most of the gross nutrient input to the pycnocline. However, in net, mixing carries nutrients away from the pycnocline. Despite the quantitative dominance of mixing in the gross nutrient transport, the nutrient richness of the pycnocline relies on the large-scale advective flow, through which nutrient-rich water is converted to nutrient-poor surface water that eventually flows to the North Atlantic. Much of the nutrient transport from the deep ocean into the ocean’s upper water column occurs through the Southern Ocean, with mixing and advection playing complementary roles, according to a box model analysis of the isotopic composition of ocean nitrate.

circulation is now understood to be the residual of the wind-driven Ekman transport and the largely opposing response of eddies to the baroclinic instability (that is, eddy advection) that results from it 14 . The contribution of the Southern Ocean surface to the global pycnocline water is widely thought to be substantial 11,14 , but its quantitative estimation remains subject to considerable uncertainty. In large part, this uncertainty arises from a lack of strong observational constraints, especially when only traditional physical and biogeochemical parameters are considered 3,13 .
Much of the uncertainty in the ocean's circulation relates to mixing. While mixing (m) is thought to have first-order importance in the transport of heat, salt and biogeochemical properties, its contribution remains poorly constrained by observations [15][16][17] . Isopycnal mixing (m S in Fig. 1a) by mesoscale eddies between the pycnocline and its Southern Ocean ventilating region (that is, the Polar Frontal Zone (PFZ) and Subantarctic Zone (SAZ)) induces tracer fluxes by stirring properties along interior isopycnal surfaces (that is, eddy diffusion) 16,17 . In the ocean interior, diapycnal mixing (m D in Fig. 1a) leads to exchanges of water and buoyancy between the pycnocline and the underlying deep ocean [4][5][6] .
With regard to large-scale advective and mixing fluxes, a distinction must be made between mass (that is, bulk seawater) and seawater properties (for example, the concentrations of nutrients). There is no net transfer of mass associated with mixing; however, when in the presence of property gradients, mixing induces both gross and net transports of properties which will, in turn, alter the density structure of the pycnocline and thus affect large-scale advective fluxes. Such transports occur because the water parcels experiencing mixing-associated motions in one direction have a systematically different property concentration from the parcels experiencing motions in the opposite direction. Integrated over large areas and long time periods, the net transport by mixing is down-gradient 17 . In our treatment of nutrients, the capital letters Ω and M refer to the nutrient transport rates by the large-scale advective water transport, ω, and by mixing, m, respectively. Gross nutrient transport by mixing is given by M in and M out for the transport into and out of the pycnocline, respectively. For example, M S-in refers to the gross nutrient input to the pycnocline from the Southern Ocean. In contrast, large-scale advection only transports nutrients in the direction of the water flow. For example, Ω S refers to the advective nutrient input to the pycnocline from the Southern Ocean, whereas Ω N refers to the advective nutrient loss from the pycnocline to the high-latitude North Atlantic (Fig. 1).
Here, we use the nitrogen isotopic ( 15 N/ 14 N) ratio of nitrate (NO 3 − ) as a new tracer with which to reconstruct the sources of water and nutrients to the pycnocline. The sensitivity of this tracer for this application stems from the elevation in NO 3 − 15 N/ 14 N in Southern Ocean surface waters due to partial consumption of NO 3 − by phytoplankton, which makes the Southern Ocean's contribution of NO 3 − to the low-latitude pycnocline isotopically distinct from the NO 3 − entering from the underlying deep ocean. We calculate the mean pycnocline and extra-pycnocline values for the concentrations of oxygen (O 2 ), phosphate (PO 4 3− ) and NO 3 − and for the δ 15 (Fig. 1a). The model is constrained by the calculated mean pycnocline properties to estimate the relative importance of the Southern Ocean surface versus the deep ocean in the gross nutrient supply to the pycnocline (Fig. 1b). For this purpose, we define the term 'pycnocline recipe' as the proportion of nutrient in the pycnocline deriving from the Southern Ocean surface: (Ω S + M S-in )/(Ω S + M S-in + Ω D + M D-in ). We also consider an advection-only, water-specific analogue, which we term the 'overturning ratio': ω S /(ω S + ω D ). When combined with other information, the model also provides novel information on the respective roles of the large-scale advective transport (ω) and mixing (m) in nutrient transport (Fig. 1c).

Nitrate 15 N/ 14 N as a tracer of nutrient fluxes
Southern Ocean processes cause NO 3 − δ 15 N in the ocean's pycnocline to be higher than that in the deep ocean [18][19][20] . Partial NO 3 − assimilation in the Southern Ocean raises the δ 15 N of NO 3 − in surface waters because of the preferential extraction of 14 N for incorporation into new biomass 21 . In the PFZ and SAZ, surface waters enter the ocean interior as AAIW and SAMW, spreading the elevated δ 15 N signal northward into the global pycnocline while lowering the NO 3 − δ 15 N of the deep ocean through the regeneration of the low-δ 15 N sinking organic matter produced in the Southern Ocean surface 19,20 (Fig. 2d). In the low-and mid-latitude ocean, the NO 3 − supply to the surface mixed layer is nearly completely consumed by phytoplankton assimilation, such that sinking nitrogen δ 15 N is similar to that of the NO 3 − supply to the mixed layer 22,23 . Since nitrification in the ocean interior typically competes with no other process, the δ 15 N of the NO 3 − added to the ocean interior by remineralization is largely determined by the δ 15 N of the sinking N from the surface   Straight arrows represent water exchange terms, while curved arrows show export production and remineralization. b, Nitrogen fluxes, including the gross inputs of nitrate to the pycnocline (that is, Ω + M in , in Tg N per year) and export production and remineralization fluxes (in Tg N per year) (mean ± s.d.). c, Deconvolution of advective and mixing components, Ω and M (mean ± s.d.), derived from the model best fits and prior estimates of the advective flows 38,39 (ω, in Sv and shown in blue). OAZ, Open Antarctic Zone; Low-lat. ML, low-latitude mixed layer.
ocean 20,24 . The internal cycling of N in the low-and mid-latitude ocean therefore preserves the high-δ 15 N imprint of the Southern Ocean on pycnocline NO 3 − (refs. 19,20,24 ). Our model corroborates a positive relationship between the pycnocline recipe and the mean difference in NO 3 − δ 15 N between the pycnocline and the deep ocean (Fig. 3a). That is, the δ 15 N difference increases as the relative contribution of SAMW + AAIW to the pycnocline nutrient budget increases. The upper bound produced by the model is indistinguishable from the conservative mixing line between the two end-members of the PFZ-SAZ source and the deep ocean source (Table 1 and black line in Fig. 3a).
The export of low-δ 15 N organic N from the PFZ-SAZ and its associated remineralization in the pycnocline counteract the increase in pycnocline NO 3 − δ 15 N from the subduction of 15 N-enriched NO 3 − . Increasing the fraction of the exported organic N from the PFZ-SAZ that is remineralized within the pycnocline decreases the δ 15 N difference between the pycnocline and the deep ocean and the pycnocline relative to the conservative mixing line (Fig. 3a). Given the relatively small contribution of the PFZ-SAZ to total remineralization in the global pycnocline (for example, <10% in the model best fits), this effect is limited.
Sources and sinks of fixed N also affect the mean pycnocline NO 3 − δ 15 N (ref. 25 ). Being the main source of fixed N, N 2 fixation (and the subsequent nitrification of the newly fixed N) is largely responsible for introducing new NO 3 − to the ocean, with a δ 15 N of approximately -2‰ to 0‰ (refs. 26,27 ). Denitrification and related anaerobic pathways (for example, anaerobic ammonium oxidation) are the dominant mechanisms of oceanic fixed N loss. Denitrification in the water column preferentially removes low-δ 15 N NO 3 − and is thus instrumental in raising the δ 15 N of ocean NO 3 − above that of N 2 fixation 25,28,29 . The isotopic discrimination associated with N loss in sediment porewaters is typically much weaker than for water-column N loss 30,31 . In an oceanic fixed N budget that is balanced for N and its isotopes, the δ 15 N of the fixed N loss equals the δ 15 N of the fixed N input 25,28 . Ocean NO 3 − δ 15 N will reflect this budget, hereafter referred to as the 'budget-driven' NO 3 − δ 15 N. To consider the effects of the fixed N budget on the δ 15 N difference between the pycnocline and the deep ocean, we develop a prognostic, time-dependent analogue of our inverse model. To capture the relevant pycnocline dynamics, it consists of an ocean subdivided into five boxes, and it allows a balance to be reached between N 2 fixation and denitrification (Supplementary Information and Extended Data Fig. 1).
We first consider the model scenario where N 2 fixation and denitrification are entirely restricted to the base of the pycnocline and above. This is likely an accurate representation of the real ocean: (i) N 2 fixation occurs mostly in the low-latitude mixed layer 32,33 and introduces newly fixed NO 3 − to the pycnocline largely through low-latitude export production and its remineralization within the pycnocline 34 ; (ii) most water column denitrification occurs in the oxygen-deficient zones, which are also located within the pycnocline 35 ; and (iii) the bulk of benthic denitrification occurs on the seafloor at depths that fall within the pycnocline or the overlying surface mixed layer 35 . In this scenario and independent of the model parameters, the N 2 fixation-denitrification balance forces the pycnocline NO 3 − δ 15 N towards the budget-driven NO 3 − δ 15 N. However, the difference in NO 3 − δ 15 N between the pycnocline and the deep ocean is indistinguishable from a model scenario without N 2 fixation and denitrification (Extended Data Fig. 2). Partial assimilation in the Southern Ocean followed by the sinking of organic N into the deep Southern Ocean and the subduction of AAIW + SAMW into the pycnocline lowers the NO 3 − δ 15 N of the deep ocean relative to the pycnocline. This pycnocline-deep ocean NO 3 − δ 15 N difference is proportional to the pycnocline recipe, regardless of whether N 2 fixation and denitrification are included in the model. However, fixed N inputs and losses may not be entirely restricted to the base of the pycnocline and above. In particular, a fraction of benthic denitrification is likely to occur outside the pycnocline, either at high latitudes or below the depth of the pycnocline 35 . In the prognostic model, we explore the sensitivity of the pycnoclinedeep ocean NO 3 − δ 15 N difference to the inclusion of some benthic denitrification occurring outside the pycnocline (Supplementary Information). This scenario increases the mean difference in NO 3 − δ 15 N between the pycnocline and the deep ocean by 0.10 ± 0.05‰ (mean ± s.d.) in comparison with the model scenario without N 2 fixation and denitrification (Extended Data Fig. 2). This is a small effect, within the uncertainties assumed in the data fitting activity   (Table 1). Colour indicates the pycnocline PFZ-SAZ remineralization/gross nutrient transport ratio (in the currency of N) into the pycnocline (that is, R SAZ /{Ω + M in }). Black dots are model best fits. b, Density function of the model best fits for the pycnocline recipe. The black line is with the constraints given by NO 3 − isotopes, the dashed grey line is without the constraints given by NO 3 − isotopes, and the blue line is with the constraints given by NO 3 − isotopes and the additional constraints given by prior estimates of the advective flows 38,39 . with the inverse, one-box model (Methods). In summary, while the ocean's fixed N input-output budget affects the NO 3 − δ 15 N of both the pycnocline and the deep ocean, it has remarkably little effect on the pycnocline-deep ocean NO 3 − δ 15 N difference.

Nutrient sources to the pycnocline
The central goal of this study is to exploit the sensitivities of NO 3 − δ 15 N to characterize the relative importance of the supply of water and nutrients to the low-latitude pycnocline from the Southern Ocean versus the underlying deep ocean. To do so, we solve the one-box model equations by varying the model parameters over a range well beyond literature estimates (Extended Data Table 1), targeting the best agreement between the observations and the model counterparts (Methods). Temperature and salinity are not used as constraints in the standard version of the model inversion because these parameters (and temperature in particular) are characterized by intense lateral variation in the overlying surface mixed layer, such that correlated lateral variation in water fluxes could yield substantial systematic errors. This concern is validated by the observation that the model best fit predicts the most accurate pycnocline temperature when the temperature for the low-latitude mixed layer is taken from the higher-latitude margins of the low-latitude surface mixed layer (that is, the low-latitude regions poleward of 30° N and 30° S; Supplementary Information). These regions are known for their ventilation of the pycnocline 36,37 . Regardless, including temperature and salinity constraints (using either the values of the bulk low-latitude mixed layer or those of the mixed layer poleward of 30° N and 30° S) has an insignificant impact on the model best fit, with the pycnocline recipe being indistinguishable from the estimate produced without temperature and salinity constraints (Supplementary Information).
Despite the large range that is tested for the model parameters, the pycnocline recipe ({Ω S + M S-in }/{Ω S + M S-in + Ω D + M D-in }) is well constrained by the observations (Fig. 3b), with a median for the best fits of 0.61 (0.56-0.68, 10th to 90th percentiles). This implies that 62 ± 5% (mean ± s.d.) of the NO 3 − in the global pycnocline originates from the Southern Ocean. The same values are recovered if the pycnocline recipe is defined in terms of PO 4 3− rather than NO 3 − . Our estimate for the pycnocline recipe is not significantly different if we allow that the mean difference in NO 3 − δ 15 N between the pycnocline and the deep ocean is also impacted by a fraction of benthic denitrification occurring outside the pycnocline (Extended Data Fig. 2), which yields a median for the best fits of 0.59 (0.53-0.66, 10th to 90th percentiles).
Our estimate for the pycnocline recipe is modestly lower than the estimate of 75% obtained using a nutrient-restoring approach coupled to an ocean general circulation model 12 . Our findings call for more equal contributions of nutrients from Southern Ocean surface and deep waters than the previous study but still indicate that the Southern Ocean supplies the greater part of the pycnocline nutrient reservoir. General circulation models parameterized with different physical terms (that is, subgrid-scale mixing terms) yield divergent pycnocline recipes, but these different models can nevertheless be tuned to produce similar physical and biogeochemical tracer fields 2,3 . Using a suite of such models, Palter et al. 3 estimated a Southern Ocean contribution to the low-latitude pycnocline of 33-75% for PO 4 3− . Excluding the constraints imposed by the NO 3 − isotope data, our calculations yield a similarly broad range for the pycnocline recipe (Fig. 3b), illustrating the benefit of including the NO 3 − isotopes as a model constraint. The decomposition of the fluxes into their advective and mixing components is less constrained by the observations. Nevertheless, our estimates are consistent with the prevalent view of the Southern Ocean's wind-driven upwelling as greater than deep upwelling in the advective input of water to the pycnocline 11-14 , with a median overturning ratio of 0.75 (0.56-0.94, 10th to 90th percentiles). In addition, a significant proportion of gross nutrient transport into the pycnocline is due to mixing (Fig. 4b), with a median M in /(Ω + M in ) ratio of 0.52 (0.33-0.67, 10th to 90th percentiles).

importance of mixing for the exchange of nutrients
To further decompose the nutrient transports from both the Southern Ocean and the deep ocean into their advective (Ω) and mixing (M in ) components, we make use of prior literature estimates 38,39 for large-scale advective fluxes (ω). By combining altimetric and gridded hydrographic data, Karsten and Marshall 38 estimated a subduction rate (that is, residual overturning) from the PFZ-SAZ ventilating area (ω S ) of 30 ± 5 Sv. Based on the energy conversion from tidal and geostrophic motions into internal waves combined with a turbulent mixing parameterization, Nikurashin and Ferrari 39 estimated an interior upwelling rate across the 27.5 g kg −1 isopycnal (ω D ) of 3.5 ± 0.8 Sv. These estimates yield an overturning ratio of 0.90 ± 0.20, at the high end of the estimates given by our model best fits (0.75 ± 0.15).
The model best fits for which ω S and ω D are in the ranges of prior literature estimates 38,39 yield a median (10th to 90th percentiles) for the pycnocline recipe of 0.63 (0.57-0.75), indistinguishable from the pycnocline recipe estimated without these additional constraints (Fig. 3b). The nutrient supplies by advection (Ω) and mixing (M in ) are well constrained by the observations (Fig. 1c), with a median M in /(Ω + M in ) ratio of 0.69 (0.55-0.75, 10th to 90th percentiles) (Fig. 4b): 0.59 (0.43-0.67) for the Southern Ocean and 0.86 (0.69-0.92) for the deep ocean. The mixing components (m D and m S ) correspond to diffusivity coefficients (κ in m 2 s −1 ) that match prior literature estimates 8,40 (Supplementary Information). The results imply that mixing contributes ~70% of the gross nutrient supply to the pycnocline and therefore plays a substantial role in sustaining low-latitude biological productivity. As predicted, our model yields a positive linear relationship between the gross nutrient import to the pycnocline (Ω + M in ) and the export production in the low-latitude areas (Fig. 4a). With the additional constraints of the large-scale advective fluxes 38,39 , export production for the low-latitude areas is estimated to be 629 ± 160 Tg N per year (mean ± s.d.), in the lower range of previous estimates (620-1290 Tg N per year) [41][42][43][44] .
While both advection and mixing drive large gross nutrient fluxes into the pycnocline, mixing also carries nutrients out of the pycnocline (M out ), both towards the Southern Ocean ventilation region and downwards into the deep ocean. In particular, net mixing-associated nutrient transport is from the pycnocline to the Southern Ocean, almost completely countering the advective nutrient input from the Southern Ocean (Fig. 1c). Thus, from the perspective of net fluxes, the deep ocean and the Southern Ocean contribute equally to the reservoir of nutrients in the pycnocline. One would not have inferred this situation from a purely advective view of ocean water fluxes.
At the same time, the mixing-associated nutrient effluxes from the pycnocline can be thought of as a consequence of the pycnocline's nutrient richness, which itself relies on the advective nutrient inputs. Advection drives the conversion of nutrient-rich water to predominantly nutrient-poor low-latitude surface waters that eventually flow to the North Atlantic: 82 ± 15 % of the advective outflow from the pycnocline passes through the low-latitude surface ocean (Fig. 1c). The essentially complete nutrient consumption in low-latitude surface waters requires that its nutrient supply is mostly returned to the pycnocline with export production and remineralization, effectively 'trapping' nutrients in the pycnocline. Over the last 30 years, the oceanographic community has come to appreciate the quantitative importance of eddies and their mixing in fuelling biological productivity 45 . While it is true that mixing dwarfs advection in the gross nutrient fluxes of the upper ocean, the advective circulation is ultimately responsible for the nutrient richness of the pycnocline.

online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41561-021-00836-8.

Data. New and previously published NO 3
− δ 15 N measurements are reported. New hydrographic sections (n = 147 hydrocasts) are GOSHIP I08S (February 2016) in the eastern Indian sector, P18S (January 2017) in the eastern Pacific sector, GOSHIP P02 (March to June 2013) in the subtropical north Pacific, and GOSHIP P06 (November 2009 to February 2010) in the subtropical south Pacific for which the methodology has already been described 47 . These new datasets are compiled in combination with a database of NO 3 − δ 15 N in the ocean, which includes most of the measurements published to date ( Fig. 2a; see Supplementary Material for a description of the database and references). All data are available in the PANGAEA database 48 , except the GOSHIP P02 and P06 sections, for which we report average profiles versus depth.
We compare the observations with the output of a one-box isotope model representing the pycnocline (Fig. 1a). We delineate the lower boundary of the pycnocline with the 27.5 g kg −1 isopycnal, located at ∼1,000-1,500 m depth in the low-latitude ocean and outcropping near the Antarctic Polar Front (APF) in the Southern Ocean and near the Subarctic Front in the North Atlantic (Fig. 2a). This boundary corresponds to the base of AAIW, which is well represented by the mid-depth salinity minimum spreading northwards in the different oceanic basins (Fig. 2b). In the North Pacific and Indian basins, this isopycnal does not outcrop to the north and remains relatively deep (~1,000 m depth, Fig. 2a). The northern limit in these basins is either the continental boundaries or the continental shelves (for example, the Bering Shelf). Between the APF and the Subtropical Front (STF), the upper pycnocline boundary is the base of the winter mixed layer (∼400 m depth) 46 . Above this boundary and up to the surface is the region we define as the combined PFZ and SAZ ventilating area, where AAIW and SAMW are subducted into the pycnocline at the Antarctic Convergence. North of the STF, the upper pycnocline boundary is the depth of the 26.0 g kg −1 isopycnal. Above this boundary and up to the surface is defined here as the low-latitude mixed layer, which exchanges gases directly with the atmosphere and in which nutrients (including NO 3 − ) are nearly completely consumed (Fig. 2c). The deep ocean is the region below the 27.5 g kg −1 isopycnal down to the seafloor, and between the southern and northern limits of the pycnocline defined above.
We used the World Ocean database 2013 (WOA13; annual properties with a resolution of 1°) 49,50 to compute the weighted averages for O 2 , PO 4 3− and NO 3 − concentrations for each domain according to the definitions outlined above. We first estimated the integrated properties at the scale of each domain (for example, moles of NO 3 − in the pycnocline) and then divided the integrated properties by the volume of the corresponding domain (Table 1). For the low-latitude mixed layer, we set the NO 3 − and PO 4 3− concentrations to ∼0.0 µmol L −1 and the O 2 concentration at saturation with the atmosphere, using the solubility as described in ref. 51 . For NO 3 − δ 15 N, the NO 3 − isotope database was used (Fig. 2a). The weighted averages were first computed for each station by taking the depth-integrated values divided by the representative thickness (see, for example, Extended Data Fig. 3 for the pycnocline). Given the relatively small amount of available data, there are inevitable biases in the dataset in terms of spatial coverage. For the PFZ-SAZ ventilating area and deep ocean, the median of the station's weighted averages is assumed to be representative of the mean properties of each domain. This approach appears to be robust because similar NO 3 − concentrations are reported for WOA13 (that is, annual properties) and the NO 3 − isotope database (Table 1). For the pycnocline, significant isotopic gradients are observed (Extended Data Fig. 3; see Supplementary Information for a discussion of the large-scale distribution of nitrate isotopes in the pycnocline), which could induce a sampling bias in the assessment of the mean pycnocline properties. These gradients are mostly driven by variable NO 3 − isotope budget-related fluxes and are restricted to the areas influenced by the North Atlantic Ocean and the Oxygen Minimum Zones 19 (OMZs, contributing to ~24% of the pycnocline NO 3 − pool). Relatively constant nitrate δ 15 N is observed outside these areas (contributing to 76% of the pycnocline NO 3 − pool and with isotopic gradients of <0.2‰ for δ 15 N between oceanic basins). We thus sorted the stations into three different categories, reflecting these variable isotope budgets: OMZ-influenced stations, North Atlantic-influenced stations and remaining pycnocline stations (Extended Data  Fig. 3). The median properties were computed for each category and then weighted by the volume contribution of each category to the total pycnocline volume. The contributions of the North Atlantic Ocean and OMZ-influenced areas to the NO 3 − pool in the pycnocline are relatively small and, taken together, increase the mean pycnocline NO 3 − δ 15 N by 0.2‰. In addition, the computed NO 3 − concentration for the pycnocline is indistinguishable from the weighted average estimated from the WOA13 (Table 1), giving us confidence in our approach.
One-box model. Nutrients and dissolved O 2 are supplied into the pycnocline from the PFZ-SAZ ventilating area and the deep ocean (Fig. 1). This occurs by large-scale advective fluxes, that is, the subduction of AAIW and SAMW (ω S , the residual overturning circulation) and interior upwelling (ω D ), and by mixing, that is, isopycnal mixing (eddy diffusion) between the PFZ-SAZ ventilating region and the pycnocline (m S ) and diapycnal mixing between the deep ocean and the pycnocline (m D ). The removal of nutrients and dissolved O 2 occurs via a large-scale advective flux, that is, sinking (ω N = ω S + ω D ) in the North Atlantic, which is partitioned between surface (ω N-LL ) and intermediate (ω N-P ) outflows, with y denoting the ω N-LL /ω N ratio, and by mixing (m S and m D ). In the model, ω N-LL also includes the component of upwelling into the low-latitude surface that flows towards the PFZ-SAZ ventilating region 38 (X* in Fig. 1a). Mixing between the pycnocline and low-latitude mixed layer (m LL ) is also represented.
For nitrogen (N), export production (ϕ) is controlled by the NO 3 − supply into the surface layer and the degree of NO 3 − consumption (f ϕ ). NO 3 − assimilation is assumed to proceed with an isotope effect of 5‰ for both the N and O isotopes 17,20,47 . A fraction (f R ) of the export production is regenerated in the pycnocline (R; hereafter the subscripts 'LL' and 'SAZ' refer to low-latitude mixed layer and combined PFZ-SAZ ventilating region, respectively), including from ϕ SAZ . The δ 15 N of the NO 3 − added to the ocean interior by remineralization is largely determined by the δ 15 N of the N exported from the surface ocean (that is, R δ 15 N = ϕ δ 15 N) 23 . The same rationale is applied to PO 4 3− , using the same values for ϕ, f ϕ and f R as for NO 3 − . In the pycnocline, the consumption of O 2 during remineralization is based on PO 4 3− and the O 2 /PO 4 3− oxidation ratio of 140 associated with aerobic respiration 52 . Complete nutrient consumption in low-latitude surface waters (that is, f ϕ-LL = 1) yields an N/P ratio for export production equivalent to the ratio of the nutrient supply (the pycnocline N:P ratio of 14.1/1). This N/P ratio is lower than the average N/P ratio of marine phytoplankton (16/1) and than the still higher ratio of exported organic matter in the subtropical ocean 53 . The low N/P ratio calculated for low-latitude export production is a consequence of the absence of the fixed N budget in the one-box model. In the real ocean, fixed N inputs and losses are such that the NO 3 − /PO 4 3− ratio of the global mean ocean (and the pycnocline) is less than 16/1 (ref. 34 ). In the PFZ-SAZ area, we define export production (ϕ SAZ ) for N (0-250 Tg N per year), using an N:P ratio of 16/1 to determine export production for P. In the Supplementary Information, we test the sensitivity of the one-box model to variable N/P ratios. The outcome is not significantly different than in the standard run.
We solve the steady-state one-box model equations by varying the parameters over the sensitivity range given in Extended Data Table 1, for O 2 concentration, PO 4 3− concentration, NO 3 − concentration and NO 3 − δ 15 N, such that the inputs and outputs are balanced (that is, dx/dt = 0):  Table 1, except the following: V is the volume of the pycnocline (2.8 × 10 20 L) and (O/P) R is the respiratory quotient (140) (ref. 52 ).
Our assumption of a steady state for the pycnocline is justified by the residence time of its nutrients calculated from the model best fits with the additional constraints given by prior estimates of the advective flows 38,39 , that is, 121 ± 40 years. This residence time is beyond the observed timescale for variations in the ventilation of the pycnocline by the Southern Ocean over recent decades 54 , and the fluxes investigated here are unlikely to change dramatically within the Holocene, including the early stages of anthropogenic global warming. Thus, while the assumption of steady state deserves further investigation, it is justified for this first application of the NO 3 − isotopes to infer nutrient fluxes at the scale of the pycnocline.

Model optimization.
To estimate the best fit with the observations, we solve the model equations by varying the parameters over the sensitivity ranges given in Extended Data  15 N}. For each model parameter, we generated 30 × 10 6 random numbers in the parameter sensitivity range (Extended Data Table 1), giving the same number of model scenarios. For each model scenario, we assess the discrepancy between observations and model counterparts, using the standardized residual (SR), as follows 55 : with X i representing the observation, X m the model counterpart and σ i the standard error for the corresponding observation. Normalizing the standard errors solves the problem of scaling between different model variables. In WOA13, the mean standard error is 0.5 µM for NO 3 − concentration, 0.05 µM for PO 4 3− concentration and 3.4 µM for O 2 concentration 49,50 . In the NO 3 − isotope database, the standard errors on the mean weighted properties are 0.05‰ for NO 3 − δ 15 N. The model best fits are then estimated using the residual sum of squares (RSS), as follows: (6) and correspond to the lowest RSS values, which represent the best agreement between observations and model counterparts. All model scenarios with an RSS value lower than 7 are considered a model best fit. The model best fits reproduce satisfactory O 2 concentration, PO 4 3− concentration, NO 3 − concentration and NO 3 − δ 15 N: with SRs symmetrically distributed with a mean not significantly different from zero and within two standard errors of the observations (Extended Data Fig. 4a). Accordingly, there is no evidence of systematic errors in the model optimization scheme. We also tested an optimization scheme in which the model best fits are the model scenarios falling within two standard errors for each model variable (Extended Data Fig. 4b). No significant differences are reported between the two optimization schemes (see Extended Data Fig. 5 for the pycnocline recipe), which supports our RSS approach.
To assess the influence of the values selected for the standard errors, we run the optimization with different values for the standard errors, that is, two times lower and higher than the values reported in WOA13 and in the NO 3 − isotope database. Extended Data Fig. 6 shows the distribution of the model best fits for O 2 concentration, NO 3 − concentration and NO 3 − δ 15 N. Lower (higher) standard errors decrease (increase) the spread of the model best fits without significantly changing the mean values, including for the pycnocline recipe (Extended Data Fig. 5). Fig. 1 | Prognostic five-box ocean model. Schematic of the meridional overturning circulation showing the reservoirs and processes included in the prognostic five-box ocean model (see Supplementary Material). This model is used to test the effects of the fixed N budget on the NO 3 − δ 15 N difference between the deep ocean and the pycnocline. The symbols in the diagram follow those in the main text and Fig. 1. Straight arrows represent water exchange terms, curved arrows show export production and remineralization, and dashed arrows indicate the addition and removal of nitrate by N 2 fixation (N 2 fix.) and denitrification (Deni.), respectively, occurring either in the pycnocline or in the deep ocean. The five boxes are the deep ocean, the pycnocline, the AZ mixed layer, the PFZ-SAZ ventilating area, and the low-latitude mixed layer. The values used for the model parameters are varied over the same sensitivity ranges as in the one-box model (Extended Data Table 1); except for m AZ and the degree of NO 3 − consumption in both the AZ mixed layer and PFZ-SAZ ventilating area, which are held constant (at 20 Sv, 0.12 and 0.16, respectively), and denitrification, which varies from 100 to 250 Tg N yr −1 . Fig. 4 | Standardized residuals for model variables. Standardized residuals for model variables (equation 5) using the optimization scheme with RSS < 7 (A) and by the requirement that model best fits are within two standard errors for all model variables (B). On each box, the red line indicates the median, and the bottom and top edges of the box indicate the 25 th and 75 th percentiles, respectively. The whiskers extend to the most extreme data points not considered outliers. The y-axis indicates the deviation between the model outputs and the observations: if equal to 0, there is no difference between the model outputs and the observations; if equal to 2, the difference between the model outputs and the observations is two standard errors (SE). The values at the top of the panels correspond to the assigned SE for each parameter.  (Table 1).