Impact of intensifying nitrogen limitation on ocean net primary production is fingerprinted by nitrogen isotopes

Modelling approach

We used the PISCES-v2 biogeochemical model, attached to the Nucleus for European Modelling of the Ocean version 4.0 (NEMO-v4) general ocean circulation model²⁹. PISCES-v2 includes five nutrients pools (nitrate, ammonium, phosphate, silicic acid and dissolved iron), dissolved oxygen, the full carbon system and accounts for two phytoplankton (nanophytoplankton and diatoms) and two zooplankton types (microzooplankton and mesozooplankton). Bioavailable nitrogen in our simulations is considered to be the combination of nitrate and ammonium. Its nitrogen cycle includes nitrogen fixation, nitrification, burial, denitrification in both the water column and sediments, and coupled nitrification–denitrification. Nitrogen isotopes were integrated within PISCES-v2 for the purposes of this study, using nine new tracers (Supplementary Note 1). Horizontal model resolution varied between ~0.5° at the equator and poles, and 2° in the subtropics, whereas vertical resolution varied between 10 and 500 m thickness over 31 levels.

We conducted simulations under both preindustrial control and climate change scenarios. The preindustrial control scenario from 1801 to 2100 maintained preindustrial greenhouse gas concentrations and only included internal modes of variability. The climate change simulation from 1851 to 2100 included natural variability, prescribed changes in land use, as well as historical changes in concentrations of greenhouse gases and aerosols until 2005, after which future concentrations associated with RCP8.5 were imposed³⁰. The biogeochemical model (PISCES-v2) was run offline from the physical model (NEMO-v4) using monthly transports and other physical conditions generated by the low resolution version of the IPSL-CM5A ESM⁵⁷.

Experiments were initialized from biogeochemical fields created from an extensive spin-up of 5000 years under repeat physical forcing, followed by a 300-year simulation under the preindustrial control scenario. The preindustrial control simulation used in analysis was therefore the final 300 years of a 5600-year spin-up involving two repeat simulations of the preindustrial control scenario. We utilized a global compilation of δ¹⁵N_NO3²⁰ supplemented with recent data to assess the isotopic routines in the model and conducted a thorough model-data skill assessment at replicating observed patterns in space (Supplementary Note 2 and Supplementary Figs. 1–3).

Anthropogenic nitrogen deposition

The effect of increasing aeolian deposition of nitrogen was assessed in our simulations. Preindustrial nitrogen deposition was prescribed as the preindustrial estimate at 1850, whereas the historical to future deposition was created by linear interpolation between preindustrial (1850) and modern/future fields (2000, 2030, 2050 and 2100). These fields were provided by Hauglustaine et al.⁸. However, the rapid rise between 1950 and 2000 was maintained, such that 60% of the increase between the preindustrial and modern fields occurred after 1950 (Supplementary Fig. 4).

The historical rise in anthropogenic nitrogen deposition was assessed by including it in additional simulations under both preindustrial control and climate change scenarios. Four initial experiments were therefore conducted: preindustrial control; preindustrial control plus anthropogenic nitrogen deposition; climate change; and climate change plus anthropogenic nitrogen deposition.

Global model experiments

We undertook four initial simulations to quantify the impacts of anthropogenic climate change and nitrogen deposition: a preindustrial control simulation from 1801 to 2100; a full anthropogenic scenario from 1851 to 2100; a climate change-only scenario without the increase in anthropogenic nitrogen deposition from 1851 to 2100; and a nitrogen deposition scenario without anthropogenic climate change from 1851 to 2100. Anthropogenic effects to nitrogen cycling were quantified by comparing mean conditions over the final 20 years of the twenty-first century (2081–2100) with mean conditions over the final 20 years of the preindustrial control simulation, whereas effects on nitrogen isotopes were quantified by comparing mean conditions over the final 20 years of the twenty-first century (2081–2100) with mean conditions over the historical period (1986–2005) from the same simulation.

To understand the direct and indirect effects of climate change, we undertook two additional idealized simulations. First, we imposed temperature changes on biogeochemical rates, while maintaining ocean circulation associated with the preindustrial control scenario, to assess the direct effects of warming on biogeochemical processes. Second, we imposed the preindustrial control temperature field on biogeochemical processes, while altering the circulation in line with the climate change scenario, to assess the indirect effects of climate change (i.e., how changing circulation alters substrate supply to biogeochemical reactions). Each experiment was run from 1851 to 2100 and without the anthropogenic increase in atmospheric nitrogen deposition, parallel with the full climate change simulation.

Agreement between the climate change simulation without anthropogenic nitrogen deposition was quantified using a pixel-by-pixel correlation analysis using Spearman’s rank correlation based on the non-parametric nature of the two-dimensional fields used for comparison. Fields were euphotic zone nitrate, twilight zone δ¹⁵N_NO3, euphotic zone δ¹⁵N_POM, and vertically integrated NPP, zooplankton grazing, nitrogen fixation, water column denitrification and sedimentary denitrification.

Depth zones

We assessed changes in biogeochemical variables related to nitrogen cycling in two depth zones defined by light. The euphotic zone was defined by depths between the surface and 0.1% of incident irradiance as recommended by Buesseler et al.⁴². The twilight zone was also defined using light, as advocated by Kaartvedt et al.⁵⁸. Depths between 0.1% and 0.0001% of incident irradiance defined the twilight zone. These definitions typically returned euphotic zone thicknesses of 137 ± 23 m (mean ± SD), and twilight zone thicknesses of 233 ± 37 m. The boundary between these depth zones were deepest in oligotrophic tropical and subtropical waters, and were shallowest in equatorial and temperate waters (Supplementary Fig. 7).

Time of emergence

ToE calculations determined when anthropogenic, anomalous trends emerged from the noise of background variability. ToE was calculated at each grid cell within both the euphotic and twilight zones (depth-averaged) and using annually averaged fields of ocean tracers. We therefore ignored temporal trends and variability at seasonal and sub-seasonal scales. Raw time series were first detrended and normalized using the linear slope and mean of the preindustrial control experiment, such that the preindustrial control time series varied about zero, while anomalous trends in experiments with climate change and/or nitrogen deposition deviated from zero. These detrended and normalized time series were smoothed using a boxcar (flat) moving average with a window of 11 years to filter decadal variability (Supplementary Fig. 12). Differences with the preindustrial control experiment were then computed.

To determine whether the differences with the preindustrial control experiment were anomalous, we calculated a measure of noise from the raw, inter-annual time series of the preindustrial control experiment (1801–2100). A signal emerged from the noise if it exceeded 2 SDs, a threshold that represents with 95% confidence that a value was anomalous and is therefore a conservative envelope to distinguish normality from anomaly¹⁶.

Furthermore, we required that anomalous values must consistently exceed the noise of the preindustrial control experiment until the end of the simulation (2100) to be registered as having emerged. Temporary emergences were therefore rejected, making our ToE estimates more conservative. A graphical representation of this process is shown in Supplementary Fig. 12.

Isolating biogeochemical ¹⁵NO₃ fluxes

We analysed the biogeochemical fluxes of ¹⁵NO₃ and NO₃ into and out of each model grid cell within the twilight zone, to determine whether the trends in δ¹⁵N_NO3 were related to biogeochemical or physical changes. Fluxes of ¹⁵NO₃ and NO₃ included a net source from nitrification (NO₃^nitr) and net sinks due to new production (NO₃^new) and denitrification (NO₃^den). Although nitrification did not directly alter the ¹⁵N : ¹⁴N ratio in our simulations, the release of ¹⁵NO₃ and NO₃ by nitrification conveyed an isotopic signature determined by prior fractionation processes that produce ammonium (NH₄). These processes include remineralization of particulate and dissolved organic matter, excretion by zooplankton and nitrogen fixation. The isotopic signatures of these processes were thus included implicitly in NO₃^nitr. For each grid cell, we calculated the biogeochemical tendency to alter δ¹⁵N_NO3 based on the ratio of inputs minus outputs:

This calculation excluded any upstream biological changes and circulation changes that might have altered δ¹⁵N_NO3.

0D water parcel model

We simulated the nitrogen isotope dynamics in a recently upwelled water parcel during transit to the subtropics by building a 0D model. The model simulates state variables of dissolved inorganic nitrogen (DIN), particulate organic nitrogen (PON) and exported particulate nitrogen (ExpN), as well as their heavy isotopes (DI¹⁵N, PO¹⁵N and Exp¹⁵N) in units of mmol N m⁻³ over 100 days given initial conditions and constants listed in Supplementary Table 1.

First, the model calculates maximum potential growth rate of phytoplankton (μ_max) in units of day⁻¹ (Eq. 8) using temperature and then finds nitrogen uptake (N_uptake, Eq. 10) using PON and limitation terms for nitrogen (N_lim, Eq. 9), light (L_lim, Supplementary Table 1) and iron (Fe_lim, Supplementary Table 1).

At a constant temperature of 18 °C, μ_max is equal to ~1.9 day⁻¹. Limitation terms for light and iron are set as constant and are used to prevent unrealistically high nitrogen uptake when nitrogen is high, such as occurs immediately following upwelling in the high-nutrient low-chlorophyll regions of the tropics. Fractionation by phytoplankton is calculated assuming an open system²¹, in this case where nitrogen can be lost through export of organic matter. To calculate the fractionation associated with uptake (¹⁵N_uptake, Eq. 11), we multiply the total nitrogen uptake (N_uptake, Eq. 10) by the heavy to light isotope ratio (({r}_{{{{{{rm{DIN}}}}}}}^{15}), Eq. 12) and the fractionation factor (ε_phy, Supplementary Table 1), which is converted from units of per mil (‰) to a fraction relative to one. This fractionation factor (ε_phy) is constant at 5‰ but is decreased towards 0‰ by the nitrogen limitation term (N_lim, Eq. 9), such that when nitrogen is limiting to growth, the fractionation during uptake decreases (last term on the right-hand side approaches 1).

At each timestep, a fraction of the PON pool becomes detritus (Eq. 15) and this detritus is instantaneously recycled back to DIN or exported to ExpN and removed from the water parcel. The amount of detritus produced per timestep is calculated as the sum of linear respiration (Eq. 13) and quadratic mortality (Eq. 14) terms, where P_resp (units of day⁻¹), K_resp (units of mmol N m⁻³) and P_mort (units of (mmol N m⁻³)⁻¹ day⁻¹) are constants (Supplementary Table 1).

Once we know the fraction of PON that becomes detritus at any given timestep, we must solve for the fraction of that detritus that becomes DIN through recycling (Eq. 17), and that which becomes ExpN through export (Eq. 18). The fraction of detritus that is recycled back into DIN is temperature dependent (Eq. 16), with higher temperatures increasing rates of recycling above a minimum fraction set by f_recmin (Supplementary Table 1). The relationship with temperature is exponential, similar to phytoplankton maximum growth (μ_max), but the degree of increase associated with warming is scaled down by a constant factor equal to T_rec (Supplementary Table 1). The fraction that is exported to ExpN is the remainder (Eq. 18).

The major fluxes of N_uptake, N_recycled and N_exported are now solved for. All that remains is to calculate the isotopic signatures of the recycling (Eq. 19) and export (Eq. 20) fluxes. These, similar to ¹⁵N_uptake (Eq. 11), are solved by multiplying against a standard ratio of heavy to light isotope (({r}_{{{{{{rm{PON}}}}}}}^{15}), Eq. 21).

Finally, we calculate the δ¹⁵N values of the major pools in the model (DIN, PON and ExpN) as output (Eqs. 22–24). We assume in this model that the major pools of DIN, PON and ExpN represent the total amount of the light isotope (¹⁴N), whereas the DI¹⁵N, PO¹⁵N and Exp¹⁵N pools represent the relative enrichment in ¹⁵N compared to a standard ratio. For simplicity, we make the standard ratio equal to 1. Therefore, taking the ratio of the DI¹⁵N to DIN pools and subtracting one returns the isotopic signature. Multiplying this by 1000 converts this signature to per mil units (‰).

Source: Ecology - nature.com