Domain Editor-in-Chief: Jody W. Deming; University of Washington, US

## 1 Introduction

The global ocean nitrogen inventory depends on the balance between dinitrogen fixation (DNF) and denitrification, and provides an important control on ocean uptake of atmospheric CO2 (Galbraith et al., 2008; Matear et al., 2010). Rapid climate change associated with the anthropogenic enhanced greenhouse effect could change the balance between these two biological processes so as to either enhance or attenuate the ocean sink for CO2. Climate change in the coming century will directly affect ocean ecosystems and temperature-dependent rates of biogeochemical processes, but what the net effect on ocean CO2 uptake will be is not known (Gruber, 2011).

At present, the open ocean (bottom depth > 180 m) accounts for the largest share of global photosynthetic biomass production at 83% of the ocean fraction; the remaining 17% includes coastal zones and upwelling regions (Karl et al., 1996). Field et al. (1998) estimated about equal marine and terrestrial contributions, implying an open ocean contribution of >40% of the global (ocean + terrestrial) total. In low-latitude subtropical oceans, as much as half of primary production is supported by DNF as its principal nitrogen source (Karl et al., 1997). At higher latitudes, recent observations in the Canadian Arctic (Blais et al., 2012) suggest a small (~4%) contribution of DNF to new primary production. Enhanced stratification under anthropogenic warming can be expected to result in reduced primary production via reduced vertical nutrient supply (e.g., Steinacher et al., 2010; Gruber, 2011; Roxy et al., 2016; Yasunaka et al., 2016). Yet enhanced stratification can also potentially increase DNF if adequate sources of phosphorus or iron are available, with opposing effects on primary production (Karl et al., 1995, 1997; Karl, 1999). As diazotrophic phytoplankton thrive in warm, stratified waters (Sohm et al., 2011), investigating the potentially enhanced role of diazotrophs in tropical and subtropical ocean ecosystems in a warming climate is important. Proxy δ15N measurements of the DNF contribution to the N inventory in the subtropical Pacific suggest that the recent apparent enhancement discussed by Karl (1999) is a long-term trend throughout the 20th century, possibly associated with anthropogenic warming (Sherwood et al., 2014).

We have examined trends in DNF in a variety of Earth System Models associated with the CMIP5 project (Taylor et al., 2012) to determine which if any show future trends and what the controlling processes may be. The simulations considered include present day conditions and future warming scenarios such as RCP 8.5 (a ‘no-mitigation’ scenario in which atmospheric CO2 concentration approaches 1000 ppm by 2100). The available models present a diversity of approaches to modelling DNF and therefore encompass a range of possible responses to future climate changes that are consistent with present knowledge of ocean DNF.

## 2 Data and methods

### 2.1 CMIP5 models and experiments

The CMIP5 data repository contains five models that include ocean DNF as an output field: CanESM2 (Zahariev et al., 2008; Arora et al., 2011), CESM1-BGC (Moore et al., 2002, 2004, 2013), GFDL-ESM2M (Dunne et al., 2013), IPSL-CM5A-LR (Dufresne et al., 2013; Aumont et al., 2015) and MPI-ESM-LR (Ilyina et al., 2013). In addition we included output from the intermediate-complexity UVic Earth System Climate Model (UVicESCM; Weaver et al., 2001), which includes a prognostic diazotroph group within its ocean biogeochemistry module (Schmittner et al., 2005, 2008), and has been run with radiative forcing from most of the CMIP5 core experiments (M. Eby, pers. comm.).

The experiments considered were the Historical and RCP 8.5 experiments. The Historical experiment (1850–2005) uses radiative forcing based on observations, which include concentrations or emissions of natural and anthropogenic aerosols, solar forcing, atmospheric composition and land use (Taylor et al., 2012), while the RCPs are hypothetical future scenarios over 2006–2100 (Moss et al., 2010; Taylor et al., 2012). We examined the differences between a “present day” climate taken as a monthly climatology of 1986–2005 of the Historical run, and the most comprehensive data base of DNF observations available (see Section 2.3). In examining the time series stations HOT (22.75°N, 158°W) and BATS (31.67°N, 64.17°W), we also considered the full 250 years of Historical + RCP8.5 (1850–2100). In the case of UVicESCM, only two periods of 20 years (1980–1999 and 2080–2099) were available to us. We refer to the individual models as CanESM, CESM, GFDL, IPSL, MPI and UVic to signify CanESM2, CESM1-BGC, GFDL-ESM2M, IPSL-CM5A-LR, MPI-ESM-LR, and UVicESCM.

Since understanding of ocean DNF is still very limited (Diez et al., 2008; Sohm et al., 2011; Luo et al., 2012, 2014; Voss et al., 2013), the models are quite diverse in their approaches to modelling DNF. Most models represent DNF as a function of light, temperature, and nutrient concentrations, or some combination of these. Table 1 summarizes the major characteristics of the models, in particular diazotroph requirements or parameterizations of dependence on nutrients and environmental conditions. Some models have a prognostic diazotroph group (e.g., CESM), while others do not. Most use parameterizations based on observations of and experiments with Trichodesmium spp., although many different pelagic diazotrophs are now known, and others may be regionally more important than Trichodesmium in terms of total contribution to DNF (Zehr and Ward, 2002; Karl and Church, 2014; Zehr and Bombar, 2015). A brief description of the various parameterizations follows.

Table 1

Controls on marine dinitrogen fixation in the different models. DOI: https://doi.org/10.1525/elementa.277.t1

Model DNF dependencea References

T PAR N P Fe O2 Biomass

CanESM +b + c ndd nd nd nd Zahariev et al., 2008
CESM + + + + nd + Holl and Montoya, 2005; Moore et al., 2002, 2007, 2013
GFDL + + + + + Holl and Montoya, 2005; Dunne et al., 2013
IPSL + + + + nd nd Séférian et al., 2013; Aumont et al., 2015
MPI nd nd e +e nd nd nd Ilyina et al., 2013
UVic + + nd + nd nd + Schmittner et al., 2008

a Except for MPI, all models depend on temperature and PAR (see Section 2.2).

b Indicates an increase of DNF with increase of the environmental factor (e.g., P limitation); under Biomass, + implies prognostic diazotrophs, which are also subject to the other limitations or inhibitions.

c Indicates a decline in DNF with increase of the environmental factor (e.g., N inhibition).

d No dependence on the environmental factor.

e In MPI, N and P dependencies are linearly combined; i.e., DNF depends on –N* = 16P – N (see Section 2.2).

### 2.2 Brief description of DNF models

The simplest parameterization, as exemplified by MPI (Equation 1), scales DNF to availability of nitrogen relative to phosphorus. This scaling is similar to the concept of the quasi-conservative tracer ${N}^{\ast }\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\left[{\text{NO}}_{3}^{-}\right]-{R}_{N:P}\left[{\text{PO}}_{4}^{3-}\right]$, where RN:P = 16 is the N:P Redfield ratio, and N* is a proxy of DNF when N* > 0 (Gruber and Sarmiento, 1997). In MPI, DNF is assumed to occur, at rate λMPI, in the surface layer wherever there is an excess of P with respect to N.

(1)

where λ0 is a reference rate in d–1 (Table 2).

Table 2

Model parameters. DOI: https://doi.org/10.1525/elementa.277.t2

Model Parameter Symbol Value Unit

MPI reference DNF rate λ0 5 × 10–3 d–1
Redfield nitrogen to phosphorus ratio RN:P 16 mol N (mol P)–1
CanESM reference DNF rate λ0 7.2 × 10–5 µmol trichome–1d–1
nitrogen half-saturation constant KN 0.1 µmol N L–1
maximum reference diazotroph concentration C1 500 trichomes L–1
surface reference diazotroph concentration C0 50 trichomes L–1
inverse depth of concentration maximum aNf 0.1 m–1
minimum temperature Tmin 20 °C
maximum temperature Tmax 30 °C
maximum reference surface irradiance Imax 350 W m–2
light attenuation coefficient for water kw 0.04 m–1
light attenuation coefficient for chlorophyll kchl 0.03 m–1 (mg chl m–3)–1
IPSL reference DNF rate λ0 0.013 µmol L–1
temperature-dependent growth coefficient µ(T) 0.6 × 1.066T d–1
minimum growth rate µ0 2.15 d–1
photosynthetic parameter Ifix 50 W m–2
half-saturation constant for P limitation KP 0.8 µmol P L–1
half-saturation constant for Fe limitation KFe 0.1 nmol Fe L–1
minimum half-saturation constants for N limitation KminNO3 0.13 µmol N L–1
KminNH4+ 0.013 µmol N L–1
UVic temperature growth coefficient µ0 0.055 d–1
reference temperature Tref 15.65 °C
minimum temperature Tmin 15 °C
half-saturation constant for P limitation KP ~0.044 µmol P L–1
initial slope of P-I curve α 0.1 (W m–2)–1 d–1
CESM reference DNF growth coefficient λ0 0.4 d–1
Arrhenius coefficient AT 4000 K
reference temperature Tref 303.15 K
minimum and maximum N:C cell quotas Qmin, Qmax 0.034, 0.17 mol N mol C–1
GFDL reference growth coefficient µ0 0.65 d–1
half-saturation constant for N inhibition KN ~7 µmol N L–1
dissolved oxygen reference concentration O2ref ~300 µmol O2 L–1
temperature sensitivity coefficient α 0.063 °C–1
maximum P:N cell quota [QP:N]max 0.167 mol P mol N–1
maximum Fe:N cell quota [QFe:N]max 4.4 × 10–3 mol Fe mol N–1
Fe:N half saturation constant KFe:N 2.4 × 10–4 mol Fe mol N–1

An intermediate approach has a more complex parameterization of DNF, but still parameterizes DNF as introduction of ‘new’ inorganic N without having a prognostic diazotroph group. CanESM prescribes a vertical profile of diazotroph abundance, linear light and temperature dependence, and inhibition by reactive nitrogen (Equation 2).

(2)

where C0, C1 and aNf are coefficients determining the vertical profile of abundance of Trichodesmium spp., T is the sea-surface temperature in Kelvin, KN and Imax are coefficients determining the dependence of DNF on inorganic N concentration and irradiance, Io is surface PAR, and kx are light attenuation coefficients for water and (surface) chlorophyll (chl0). CanESM has no limitation by other nutrients such as P or Fe. Parameters for all of the models are defined in Table 2.

IPSL has the most complex diagnostic parameterization (Equations 3a–g). In particular, inhibition by reactive nitrogen depends on both ammonium and nitrate. Temperature dependence is exponential for T ≥ 20°C, and factors for iron and phosphorus limitation are included (Equation 3d). Note that, in IPSL, the light dependence of DNF is an exponential function of PAR, represented as I(z) here.

(3a)
(3b)
(3c)
${L}_{\text{\hspace{0.17em}Fe}}=\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\frac{\text{[Fe]}}{{K}_{\text{Fe}}\text{\hspace{0.17em}}+\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\left[\text{Fe]}},\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{L}_{\text{\hspace{0.17em}}{\text{PO}}_{4}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\frac{{\text{PO}}_{4}^{3-}}{{K}_{P}+\text{\hspace{0.17em}\hspace{0.17em}}\left[{\text{PO}}_{4}^{3-}\right]}$
(3d)

Kx, where x is either ammonium, nitrate, phosphate or iron, is also a function of phytoplankton biomass for ammonium and nitrate (Aumont et al., 2015, their equations 7a–c). Aumont et al. (2015) explain that “half-saturation constants increase with the size of the phytoplankton cell as a consequence of a smaller surface-to-volume ratio”:

(3e)
(3f)
${K}_{X}=\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{K}_{X}^{\mathit{min}}\text{\hspace{0.17em}\hspace{0.17em}}\frac{{P}_{1}\text{\hspace{0.17em}\hspace{0.17em}}+\text{\hspace{0.17em}\hspace{0.17em}}{S}_{\mathit{sat}}{P}_{2}}{{P}_{1}\text{\hspace{0.17em}\hspace{0.17em}}+\text{\hspace{0.17em}\hspace{0.17em}}{P}_{2}}$
(3g)

where P is the biomass, Pmin a minimum phytoplankton biomass, Kxmin the minimum half-saturation constant (for either ${\text{NO}}_{3}^{-}$ or ${\text{NH}}_{4}^{+}$; Table 2) and Ssat the size ratio of phytoplankton (set to 3).

UVic (Equation 4) exemplifies the general form of models with a prognostic diazotroph group. Diazotroph biomass (DN in mmol N m–3) is combined with a temperature-dependent growth rate and a nutrient limitation factor, in this case phosphate limitation (Schmittner et al., 2008). In the UVic model, the rate of DNF is determined by the most limiting factor, either nutrient or light at depth z (Schmittner et al., 2008).

${\lambda }_{\mathit{UVic}}={\lambda }_{\mathit{max}}\text{\hspace{0.17em}\hspace{0.17em}}{D}_{N}\text{\hspace{0.17em}\hspace{0.17em}}\mathit{min}\text{\hspace{0.17em}\hspace{0.17em}}\left(\frac{\left[{\text{PO}}_{4}^{3-}\right]}{{K}_{P}+\text{\hspace{0.17em}\hspace{0.17em}}\left[{\text{PO}}_{4}^{3-}\right]}\text{\hspace{0.17em}\hspace{0.17em}},\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\frac{\mathit{\alpha I}\left(z\right)}{\sqrt{{\lambda }_{\mathit{max}}^{2}+\text{\hspace{0.17em}\hspace{0.17em}}{\left(\mathit{\alpha I}\left(z\right)\right)}^{2}}}\right)$
(4a)
(4b)

Diazotrophs can use nitrate where available, but are not growth-limited by nitrogen. Equation A3 in Schmittner et al. (2008) defined ${\text{NO}}_{3}^{-}$ uptake as follows:

(4c)

CESM (Moore et al., 2013) is mostly based on the earlier BEC model by Moore et al. (2002). The CESM DNF implementation differs in several ways from the original BEC DNF implementation. For instance, DNF depends on the carbon biomass, growth rate, and a biomass-specific DNF rate. The modelled diazotrophs can utilize both fixed dinitrogen and inorganic nitrogen as a nutrient source following Holl and Montoya (2005); nitrate and ammonium half-saturation constants have been set so that diazotrophs do not strongly compete with the other phytoplankton (0.1 μM, about 5 to 10 times larger than for the other groups; Moore et al., 2013). Below (Equation 5) we reproduce the BEC DNF (Moore et al., 2002), which depends on cell N quota (Q), diazotroph biomass and temperature. Following Geider et al. (1998), as for the other phytoplankton groups, the modelled inorganic nitrogen uptake declines as Q approaches Qmax. Temperature-dependence is parameterized according to the Arrhenius equation, with a minimum temperature for diazotroph growth of 16°C.

(5a)
(5b)

DN (DC) is diazotroph biomass in moles N (C) m–3; Q, Qmin, and Qmax are the N:C cell quota and minimum and maximum quotas, respectively; and T is the temperature expressed in Kelvin. Diazotroph biomass is modelled prognostically according to modified Geider et al. (1998) iron- and phosphorus-limited growth (5 nM and 90 pM, phosphate and iron half-saturation constants, respectively, in Moore et al., 2007).

GFDL (Equation 6) is the most complex and complete prognostic model, with micro- and macronutrient limitations based on iron and phosphorus cell quotas, an inorganic nitrogen inhibition factor, and an oxygen inhibition factor. Temperature dependence is exponential and not limited to the tropical range as in the four previous models (no temperature dependence in the case of MPI).

(6a)
(6b)
$\begin{array}{l}{f}_{P}=\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\frac{{Q}_{P:N}}{{\left[{Q}_{P:N}\right]}_{\mathit{max}}},\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{f}_{\mathit{Fe}}=\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\frac{{Q}_{\mathit{Fe}:N}^{2}}{{K}_{\mathit{Fe}:N}^{2}+\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{Q}_{\mathit{Fe}:N}^{2}}\text{,\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\\ \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{Q}_{X:N}=\mathit{min}\left(\frac{{D}_{x}}{{D}_{N}},\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{\left[\frac{{D}_{x}}{{D}_{N}}\right]}_{\mathit{max}}\right)\end{array}$
(6c)

The term e1–f(PAR) represents the light limitation following Geider et al. (1997), where f(PAR) depends on the initial slope of the PI curve, the carbon to chlorophyll ratio and temperature. DN (DP, DFe) is the diazotroph biomass in moles N (P, Fe). Q and Qmax are the phosphorus or iron cell quota (relative to nitrogen) and maximum quota, respectively. O2 and O2ref are the dissolved oxygen concentration and inhibition concentration (Dunne et al., 2013); KN is a half-saturation constant for reactive nitrogen inhibition. The modelled diazotroph biomass follows Geider et al. (1997) and can assimilate dissolved inorganic nitrogen following Holl and Montoya (2005).

### 2.3 Observations

Luo et al. (2012) provide the largest up-to-date dataset of observed DNF rates. We have used a subset of the Luo et al. (2012) dataset: 15N2 assimilation measurements over the whole euphotic zone. These are believed to be the most accurate estimates of water column DNF rate, due to consistent use of a modern and accurate method (Montoya et al., 1996).

A preliminary examination of IPSL and UVic rates suggested that they have extremely high DNF rates (total global 657.2 and 315.9 Tg N y–1, respectively), and their maximum areal rates are outside the range of the other CMIP5 models. In IPSL the long-term nitrogen budget was balanced by scaling DNF to denitrification (Séférian et al., 2013), and denitrification is high because volume of the oxygen minimum zone (OMZ) in the Pacific is overestimated (Séférian et al., 2013, Cabré et al., 2015). For illustrative purposes, IPSL and UVic DNF rates have been rescaled in Figure 1 by 0.17 and 0.40, respectively. UVic was scaled to the global total DNF rate estimate of Gruber and Sarmiento (1997); IPSL was scaled to its value in an ocean-only simulation with the same biogeochemistry model (Aumont et al., 2015). For IPSL only, this rescaling is retained in all of the figures and tables, and IPSL is excluded from any discussion of absolute magnitude.

Figure 1

Global depictions of mean depth-integrated DNF rate from the six models examined. Values as estimated by a) CanESM, b) CESM, c) IPSL, d) UVic, e) GFDL, and f) MPI, for 1986–2005 of the Historical experiment, except UVic which is based on the period 1980–1999. All models have been regridded to a common 2° × 2° grid. A common colour scale is used to indicate DNF rate; black areas correspond to values >400 μmol m–2 d–1. DOI: https://doi.org/10.1525/elementa.277.f1

Most of the available observations in Luo et al. (2012) are from subtropical regions. We limited our analyses to latitudes between 45°S and 45°N. Since the observations are mostly in the open ocean and the models have coarse resolution, we also limited our analysis to the open ocean, excluding coastal waters and some of the marginal seas. In the Atlantic, we included the Gulf of Mexico and the Caribbean Sea, as these are large areas of warm, stratified ocean known from historical observations to have active diazotroph populations. In the Indian Ocean, we included the Arabian Sea but excluded the Bay of Bengal (north of 5°N), because coarse resolution models create a nonphysical “nutrient trap” there, and because observations suggest that DNF is much greater in the Arabian Sea (Naqvi, 2008). All Pacific marginal seas were excluded, as they are poorly resolved by coarse-resolution global models; DNF is known to occur in the South China Sea, but at much lower rates than in open ocean waters of the low-latitude Pacific (Chen et al., 2003, 2014).

Modelled DNF rates were regridded to a common regular 2° × 2° grid. Observed DNF was projected onto the same grid and available observations pooled and averaged within each grid cell. Total DNF was calculated for the models and the observations using the regridded and averaged maps. Individual-basin totals are estimated for 0–45°N or S. The Indian Ocean total includes both the northern and southern hemispheres (north of 45°S), excluding the Bay of Bengal north of 5°N. Estimates of the total observed DNF rates are based on Luo et al. (2012) within the same areas. Mean observed rates are calculated (weighted by grid cell area) within each region and multiplied by the area of the region.

We also compared the modelled DNF and primary production (PP) at the benchmark time series stations HOT (22.75°N, 158°W) and BATS (31.67°N, 64.17°W). DNF data at these locations were taken from Luo et al. (2012) while PP rates are available over the period 1989–2014 (Bates et al., 2001; Lomas et al., 2002; Brix et al., 2006; Karl and Church, 2014) represented by 247 and 320 approximately monthly profiles at HOT and BATS, respectively. For DNF and PP time series comparisons, a single model grid cell containing the stations was chosen, and data were averaged to annual means. At BATS only a very limited number of observed DNF rates was available (see Section 3). At BATS we also considered a geochemical estimate of DNF based on long-term nutrient concentration time series (Singh et al., 2013) and modelled DNF estimates by Hood et al. (2001). For modelled DNF and PP at these locations, a single 2° × 2° grid cell containing the stations was chosen, and model outputs were averaged to annual means.

Linear trends were calculated for modelled DNF and PP at HOT and BATS for the period 1986–2100, except in the case of UVic where only the 1980–1999 and 2080–2099 time periods were available. In this case, estimating the trend might be more appropriately done as simply the difference between means for these two periods, but calculating regression coefficients for all available data gives essentially the same answer. PP trends were calculated for HOT and BATS for the time period of data availability. Not enough data were available to estimate trends in observed DNF. Regression coefficients were tested for significance using a standard t-test, assuming a decorrelation time scale of 3 years for the models and 3 months for the observations. This approach is conservative in all but a few cases, but it makes no difference to the results. All significant trends were significant regardless of what assumptions were made regarding autocorrelation, with most p values at 0.01 or 0.001.

## 3 Results

Table 3 shows the range of total DNF rates and their distribution among ocean basins. In the Pacific, most of the estimated basin totals show good agreement with observations. UVic systematically overestimates DNF, but the distribution among basins is similar to other models. In the Atlantic, the total DNF rates tend to be larger than observations, especially in the South Atlantic. South Atlantic rates are low in the Luo et al. (2012) data base, but similar rates have been reported elsewhere (Fernandez et al., 2010). CanESM, CESM, and GFDL represent the lower end of the range of modelled rates of DNF. Global total DNF for CanESM, CESM, and GFDL fall within 100–200 TgN y–1, with MPI just slightly higher; observation-based estimates range from 103 to 177 TgN y–1 (Großkopf et al., 2012; Luo et al., 2012). As mentioned earlier, IPSL CMIP5 data were rescaled to the rate of denitrification, resulting in a global total about five times larger than the observed range; ocean-only simulations with the same biogeochemistry model produce rates within this range (Aumont et al., 2015). UVic represents the higher end showing areal values as high as 1135 μmol m–2 d–1 (see Table S1 for statistics of observed and modelled distributions) and a global total of 316 TgN y–1. Latitudes less than 45° account for the vast majority of the global totals (>90%, except for GFDL with 80% and MPI with 73%; Table 3). GFDL has substantial DNF in the mid-latitudes, particularly in the Southern Hemisphere (Figure 1), along the southern boundaries of the ocean basins where there is potential for episodic equatorward transport of surface nutrients from the HNLC region. The Pacific Ocean represents the largest share across all models with 40–57% of the global total (Table 3), which is mainly due to the larger surface area of the Pacific Ocean as areal DNF rates are usually within the same range across ocean basins. Totals for latitudes less than 45° tend to be higher than the observations. The Atlantic and Indian Ocean fractions (Table 3) are of similar magnitude, particularly in CanESM, CESM, and UVic. The modelled total for the Indian Ocean ranges between 9 and 27% (mean of 18.0%) of the global total; Atlantic Ocean fractions are slightly higher (13–26%, mean of 19.9%). In all models, the totals for the northern and southern hemispheres in the Pacific and Atlantic Oceans are similar.

Table 3

Total DNF rates (Tg N y–1) by region. DOI: https://doi.org/10.1525/elementa.277.t3

Model Global total Regional totalsa

N Pacific S Pacific N Atlantic S Atlantic Indian Ocean Latitude < 45°

CanESM 144.4 40.2 30.4 19.1 11.5 31 132.3
CESM 175.1 38.5 32.1 26.4 19.2 47.4 163.6
IPSLb 111.8 34.1 29.0 10.2 11.8 15.4 100.6
UVic 315.9 87.3 92.7 34.5 18.1 65.2 297.8
GFDL 158.8 29.4 38.2 15.7 20.8 23.6 127.7
MPI 203.7 48.1 56.6 11.8 13.9 18.3 148.7
Observedc nad 34.3 49.5 10.9 2.3 na 97

a Sum of individual basin/hemisphere totals to latitude < 45°.

b IPSL rescaled by 0.17 (see Section 2.3).

c Values estimated from data compiled by Luo et al. (2012).

d Indicates not available.

Figure 1 shows the geographic distribution of DNF for all six models. Note that UVic has been rescaled for illustrative purposes here (see Section 2.3). Models can be grouped according to spatial distribution. CanESM, CESM, IPSL and UVic all have high DNF in oligotrophic regions on the western boundaries, i.e., maximal outside of the coastal and equatorial zones with high concentrations of upwelled inorganic nitrogen. This distribution is not surprising as the models are parameterized to have inhibition of DNF by inorganic nitrogen (see Section 2.2). In these models, DNF is also dependent on warm temperatures, with minimum temperatures from 15°C (UVic) to 20°C (CanESM/IPSL).

GFDL and MPI exhibit very different distributions from the other four models, with the largest DNF in the coastal and equatorial upwelling zones (Figure 1). In these models, DNF occurs in waters with excess P relative to N. In MPI, it is a simple linear relationship of DNF to –N* (Equation 1), and DNF is not inhibited by low temperature. In GFDL, nitrogen inhibition is combined with oxygen inhibition (Equation 6) which occurs at surface concentrations of oxygen (~300 μM; Table 2). GFDL nitrogen inhibition has a high half-saturation constant (7 μM) compared to the other models, and is restricted mostly to subsurface waters or upwelling regions. In these two models, there is also high DNF at mid-latitudes where temperature limits DNF activity in the other models: GFDL in the Southern Ocean and MPI in the North Pacific.

Spatial pattern correlations among models also suggest two groups (Table S2). MPI and GFDL are weakly correlated with the other models (r < 0.18), and are relatively weakly correlated with each other (r = 0.32). The other four models tend to have stronger spatial correlation with each other (r > 0.5), with the exception of CESM and UVic (r = 0.31).

In the Indian Ocean, most of the models agree south of 15–20°S. Most models show a steep gradient from high to very low DNF rates, except for GFDL and MPI which have enhanced rates at higher latitudes (Figure 1). Rates are low near Western Australia, except in CESM (along most of the coast) and GFDL. In the interior of the Indian Ocean, CanESM and CESM have uniformly high DNF rates, while MPI is uniformly low and the other models are spatially variable. In CanESM, CESM, IPSL and UVic, high DNF areas spread across the tropical and subtropical oceans with a large area of low rates in the eastern tropical Pacific Ocean. IPSL tends to have high DNF on the eastern side of the Indian Ocean and low DNF in the southeast Pacific Ocean.

Figure 2 shows comparisons of modelled DNF rates to observations compiled by Luo et al. (2012, 2014). Figure 2a shows the spatial distribution of measurements; almost all of the data are located in the Pacific and Atlantic Oceans and during the 2000–2010 time period (Luo et al., 2012). Within the same region, the data can range over several orders of magnitude. For example, in the North Atlantic the range is 1 to 248 μmol m–2 d–1.

Figure 2

Observed depth-integrated DNF and scatter plots of modelled vs. observed values. a) Observed depth-integrated DNF rate (monthly climatology) from the MAREMIP database, and (b–g) comparison between mean modelled (y-axis) and observed (MAREMIP, x-axis) rates: b) CanESM, c) CESM, d) IPSL (rescaled by 0.17), e) UVic, f) GFDL, and g) MPI. The black dashed line indicates 1:1. An outlier observed value (1189 μmol m–2 d–1) has been omitted. Note different time period for UVic. DOI: https://doi.org/10.1525/elementa.277.f2

The models overestimate DNF rates, especially in the south Atlantic Ocean (Figure 2b–g). The small number of rates exceeding 200 μmol m–2 d–1 are not consistently reproduced by any model. Counting the points above and below the 1:1 line, all of the models overestimate DNF, with IPSL (rescaled by 0.17) having the highest fraction of locations overestimated (95 of 133) and GFDL the lowest (79 of 131). The relative overestimation is usually larger in the lower part of the observed range, e.g., 0–100 μmol m–2 d–1; 77% of the observations fall in this range. CanESM, CESM, GFDL and MPI (excluding an outlier in the latter) have an overall range lower than the observed rates: observed values go up to ~600 μmol m–2 d–1 (one outlier at 1189 μmol m–2 d–1 has been excluded for visualization purposes) while modelled values are less than 300 μmol m–2 d–1. GFDL and MPI mostly cluster close to the 1:1 line at the lower end of the range. The spatial correlation between models and observations is not statistically significant in most cases. Correlations are significant only for CESM (r = 0.21) and IPSL (r = 0.23) for untransformed data, and for IPSL (r = 0.32) and UVic (r = 0.21) for log-transformed data (Table S2).

Figure 3 shows frequency histograms of DNF rates. Observed rates vary from 0 to 600 μmol m–2 d–1 (excluding one outlier as above). The distribution is strongly positively skewed with a large proportion of low rates and a tail of high rates. Histograms of modelled DNF rates at the locations of the observations show that IPSL, UVic, and GFDL agree on the shape of the distribution. MPI mostly agrees with the observations, although it has a larger frequency of high DNF rates and lower frequency of low rates. In the North and South Pacific and South Atlantic, MPI exhibits a mode around 75 μmol m–2 d–1. CanESM and CESM differ from the observations in having a more Gaussian distribution with a strong mode at intermediate DNF rates.

Figure 3

Histograms of DNF rates (monthly climatology) at locations of observations. a) Mean observed depth-integrated DNF rate from the MAREMIP database, and (b–g) modelled DNF rates: b) CanESM, c) CESM, d) IPSL (rescaled by 0.17), e) UVic, f) GFDL, and g) MPI. Modelled values are monthly, for a climatology calculated over the 20-year period indicated. An outlier observed value (1189 μmol m–2 d–1) has been omitted. Note different scales and time period for UVic. DOI: https://doi.org/10.1525/elementa.277.f3

Figure 4 shows frequency histograms similar to Figure 3 but considering all locations, not just those where observations exist. In some cases, the shapes of the histograms change substantially. Overall the shapes of the histograms are closer to that of the observations, particularly for CanESM and CESM. IPSL and UVic change the least, and overall have a shape closest to the observations (rescaling does not affect this result). All of the models approximate the observed pattern of a positively skewed distribution (with varying degrees of skew) but several (CanESM, GFDL, MPI) have modes at intermediate rates that are not present in the observations. Statistics of the distributions are provided in Table S1. GFDL is the closest to the observed distribution in terms of maximum, standard deviation and skewness. All exhibit positive skew ranging from 0.2 (CanESM) to 4.8 (MPI); even the lowest values are significantly different from 0 (P < 0.001).

Figure 4

Histograms of DNF rates for all locations. a) Mean observed depth-integrated DNF rate from the MAREMIP database (same as Figure 3a, reproduced here for ease of comparison), and (b–g) modelled DNF rates for all locations: b) CanESM, c) CESM, d) IPSL (rescaled by 0.17), e) UVic, f) GFDL, and g) MPI. An outlier observed value (1189 μmol m–2 d–1) has been omitted. Note different scales and time period for UVic. Histogram construction and colour codes are identical to Figure 3. DOI: https://doi.org/10.1525/elementa.277.f4

Figure 5 compares observations and models at the benchmark subtropical stations HOT and BATS. For DNF, observations at HOT cover a longer time period than at BATS and thus capture more of the variance in actual rates. At HOT, observed primary production (PP) shows an increasing trend (Figure 5a; Table 4). For the models, trends are calculated over a much longer period (1986–2100), with only CanESM and UVic showing positive trends in PP. While the absolute rate of PP is very low in CanESM, the relative magnitude of the trend is similar to that in the observations, suggesting approximately a doubling of PP over a century. At BATS (Figure 5b), CanESM shows a similar trend, but the observations do not show a statistically significant trend. Based on 1980–1999 and 2080–2099, UVic PP increases slightly over time at both stations, but declines at HOT after 2080. All of the other models show a decrease of PP with climate warming. In CanESM, DNF supports about ~45% of PP at HOT and 3–30% at BATS, assuming a Redfield scaling.

Figure 5

Time series of modelled and observed DNF and primary production at stations HOT and BATS. Model outputs represent the 2° × 2° grid point nearest the stations. a) Primary production at HOT, b) primary production at BATS, c) DNF rate at HOT, d) DNF rate at BATS. Observed DNF rates are depth-integrated 15N2 assimilation measurements (Luo et al., 2012); primary production is from 14C uptake as per HOT and BATS protocols. The upper bound for DNF is 222 μmol N m–2 d–1 for strong mixing and 518 μmol N m–2 d–1 for weaker mixing (Hood et al., 2001). The time period of the observations is 1989–2014. The horizontal positions from Hood et al. (2001) for strong mixing, weaker mixing, and from Singh et al. (2013) are located in the centre of the associated time periods, 1989–1994, 1995–1996, 1988–2009, respectively (grey-scaled bars). IPSL DNF rate has been rescaled by 0.17. DOI: https://doi.org/10.1525/elementa.277.f5

Table 4

Linear trends in DNF and primary production (PP) over 100 years, expressed as percentage of the mean value in 1986a, at the time-series stations HOT and BATS. DOI: https://doi.org/10.1525/elementa.277.t4

Observations or models HOT (22.75°N, 158°W) BATS (31.67°N, 64.17°W)

PP DNF PP DNF

Observations 108.3 nab 22.2 (NS) na
CanESM 79.0 78.7 201.2 93.2
CESM –8.6 21.2 –11.4 –63.9
IPSL –25.4 –103.0 –12.5 (NSc) 224.8
UVic 4.8 –73.6 4.0 –89.9
GFDL –15.5 35.1 –1.9 (NS) –12.7
MPI –14.7 –16.4 –32.1 –45.5

a Except for UVic (1980) and observed PP (1988).

b Not available due to insufficient data.

c Not significant; all other trends are statistically significant (α = 0.05 or less).

Modelled DNF rates (Figure 5c) all fall within the observed range at the HOT station (16–628 μmol m–2 d–1), although CESM is clearly at the high end. At BATS (Figure 5c), only two models (UVic, MPI) fall within the lower observed range of 15–39 μmol m–2 d–1 (which is based on only 6 data points) in the present climate. The other three models fall within the higher range of 111–129 μmol m–2 d–1 estimated from long-term time series of nutrient concentrations at BATS (Singh et al., 2013). When comparing modelled DNF to earlier estimates based on models constrained by observations of DIC and N export (Hood et al., 2001), most of the CMIP5 models fall within the lower half of this range. CESM has higher rates in the present climate but declines into the observed range after about 2040. IPSL (rescaled by 0.17) and UVic rates are the lowest at this location. Rates in the tropical-subtropical Atlantic decline with latitude in CESM and CanESM, but this gradient is much more abrupt in IPSL and UVic (Figure 1a–d).

DNF rates at HOT show a range of trends (Figure 5c), with half of the models projecting a long-term increase (CanESM, CESM and GFDL) while the others show a decline. At BATS (Figure 5d), there is similar divergence of trends with two models (CanESM and IPSL) showing a long-term increase and the other four declining. Only half of the models show trends of the same sign at both stations, and all four possible combinations (+ +, + –, – +, and – –) are represented among the six models (Tables 4 and 5).

Table 5

Signsa of the modelled trends of DNF and environmental factors at the time-series stations HOT and BATS. DOI: https://doi.org/10.1525/elementa.277.t5

Model DNF SSTb NO3 NH4 Phosphorus Iron Biomassc

HOT BATS HOT BATS HOT BATS HOT BATS HOT BATS HOT BATS HOT BATS

CanESM + + + + +d + nd nd nd nd nd nd nd nd
CESM + + + d +d d + +
IPSL + + + + + + + + + nd nd
UVic + + nd nd nd nd nd nd
GFDL + + + d d +d + + + +
MPI nd nd nd nd nd nd nd nd

a + and – signs indicate non-zero trends; nd, no dependence on the environmental factor.

b Sea surface temperature.

c When DNF depends on biomass (CESM, GFDL, and UVic), the other dependencies (N, P, Fe) constrain a biomass-specific DNF rate.

d Indicates trend is not statistically significant.

The physical processes underlying the positive and negative trends vary. Trends in environmental factors are summarized in Table 5. Figure S1 shows the corresponding time series of the environmental controls for CESM, IPSL, UVic, and GFDL. Table S3 shows the temporal correlation between DNF rates and environmental controls across all the models. In CanESM, DNF increases at both stations despite an increase in the concentration of dissolved inorganic N (DIN) at BATS (Table 5). The concentration increase is not enough to inhibit DNF, and DNF appears to be driving DIN rather than the reverse. In MPI, DNF declines at both stations, as the supply of P declines relative to that of N (Table 5). Note that in CESM and IPSL, inhibition of DNF by ammonium is stronger than by nitrate (Table 2). In CESM, GFDL, and IPSL, processes leading to nutrient limitation differ at BATS and HOT, but P limitation generally prevails over iron limitation. Dissolved iron concentration is high; e.g., the minimum iron concentration at HOT in IPSL is ~150 pM (Figure S1e), whereas the half-saturation constant is 100 pM (Table 2). P decreases over time in all but one case: GFDL at BATS (Table 5). In IPSL, P concentration at BATS is higher than at HOT, but the rate of decline is greater at HOT (Figure S1e–f); P is completely depleted by 2016. In UVic, enhanced stratification limits the resupply of inorganic nutrients; in particular, surface PO4 steadily declines (Table 5, Figure S1g–h). As a result, both diazotroph biomass and per-cell DNF decline at both stations (Table 5; Figure S1g–h).

In Table S4, we compare the trends of PP with the same environmental factors. As noted above, because iron concentration (CESM, IPSL, GFDL) is mostly above the half-saturation constant, the PP decline is likely to depend on N or P. Except for CanESM (and UVic at BATS), enhanced thermal surface stratification reduces PP despite higher metabolic rates in warming waters. In UVic the long-term trend in PP at HOT is small (Figure 5). In CanESM, PP correlates significantly with surface nitrate at BATS; i.e., increasing PP is driven by increasing DNF. In CESM and GFDL, declining PP (Figure 5) follows declining N and P concentrations (Table S4) at both HOT and BATS; i.e., the trends that emerge from the simple parameterizations in CanESM are not present in the more complex models with prognostic diazotrophs and multiple nutrient limitations.

## 4 Discussion

### 4.1 How well do the models compare with the observations? Does one model emerge as best representing the observed spatial pattern? How do the models compare with each other?

The spatial pattern correlation between models and observations at the locations of the observations (Figure 2) suggests that to meaningfully distinguish among these models using these observations is not possible. The highest spatial correlation (log-transformed) was 0.32, and four of the six models showed no significant correlation. Using the observations to distinguish between the two model groups was also difficult (Figure 2; Table S2), as observations are mostly lacking in eastern boundary waters (Figure 2a). This finding indicates a need for future measurement campaigns in undersampled regions (Figure 2a). Figure 2b–g further show that models systematically overestimate (underestimate) in the low (high) range of the observations. In particular, the observed DNF rate in the South Atlantic is low (0–10 μmol m–2 d–1), while the modelled range is much broader (up to 300 μmol m–2 d–1). In the South Atlantic, MPI has the least spread and the least positive bias with values mostly within 30–100 μmol m–2 d–1; these rates are the closest to the observations in this particular region (Figure 2a).

Waters upwelling along eastern boundaries contain an excess of dissolved inorganic P with respect to N requirements (high negative N*). High P:N ratios tend to favour diazotrophs over other phytoplankton groups, which may result in strong spatial coupling between surface DNF and subsurface denitrification (Tyrrell et al., 1999; Deutsch et al., 2007; Landolfi et al., 2013). Denitrification removes DIN and potentially enhances DNF over OMZs in five of the six models (CESM, GFDL, IPSL, MPI, and UVic), all of which have P dependence of DNF. The spatial coupling is obvious in GFDL and MPI in equatorial and coastal upwelled waters (Figure 1). Model biases in ocean circulation and thermocline ventilation may result in expanded OMZs (Dunne et al., 2013; Cabré et al., 2015), so that the models may overestimate denitrification and upwelling of low N/P waters. In GFDL this bias is compounded by O2 inhibition of DNF, particularly in the Pacific.

At HOT and BATS, only CanESM shows a positive trend in PP (mainly attributable to increasing DNF). These trends are consistent with previous observations (Karl, 1999; Sherwood et al., 2014). The other models encompass all possible combinations of PP and DNF trends. Increasing P limitation appears to be the most common environmental control among the models (Table 5; Figure S1).

### 4.2 Is the distribution of DNF Gaussian? Are the available observations sufficient to know? Do the models reproduce the observed distribution?

Figure 3a suggests a strongly positively skewed distribution of observed DNF rates (statistics are summarized in Table S1). This result is based only on 15N2 assimilation, but the distribution is similar to what was found using both 15N2 assimilation and C2H2-reduction assay measurements (Luo et al., 2012). Whether the positive skew of modelled DNF is due to the same processes as in the observations is not clear, but a positive skew of even 0.2 does not occur by chance (P < 0.0005 for N > 2000). Using all model grid points (Figure 4) gives a distribution of modelled DNF rates that more closely resembles the observed distribution than one considering only locations where observations exist (Figure 3). CanESM, CESM, GFDL, and MPI show a much more positively skewed distribution, while IPSL and UVic distributions change little and fit the observed distribution fairly well.

Luo et al. (2012) argued that the currently available observations did not cover the eastern equatorial regions for historical reasons (e.g., effort was focused on nutrient-depleted oligotrophic regions), even though these regions are important to understand the connection of DNF to processes like denitrification, upwelling, and O2 deficiency. The diversity of parameterizations in the CMIP5 models makes this data set potentially useful for understanding these processes as new observations emerge, although model error in processes such as ventilation and denitrification may limit their usefulness in some cases.

### 4.3 What are the important dependencies on environmental conditions? Are the models likely to respond to future changes in environmental conditions in a realistic way?

Besides temperature and PAR dependence, DNF parameterizations may assume iron and/or phosphorus limitation, and inhibition by DIN and O2. DIN inhibition is a dependence suggested by studies on Trichodesmium spp. (Capone et al., 1997; Holl and Montoya, 2005; Zehr, 2011); many of the models considered here do not have complete inhibition at high DIN concentrations.

Recent laboratory experiments have qualified nitrogen inhibition by showing that DNF can occur when diazotrophs are acclimated to chronic high DIN concentrations, although at a lower rate than at low DIN concentration (Knapp, 2012; Knapp et al., 2012). DIN inhibition may be offset by phosphorus supply in iron-replete conditions as observed in laboratory experiments (Knapp et al., 2012) and field studies (Fernandez et al., 2011, 2015). Knapp et al. (2012) suggest that enhanced abundance of diazotrophs in P- and Fe-rich waters offsets reduced DNF due to nitrogen inhibition, resulting in DNF rates similar to low DIN conditions. In CESM, however, iron deficiency limits DNF in the southeastern tropical Pacific so that spatial coupling between DNF and denitrification is weak (Moore et al., 2007).

Additional supply of dissolved organic phosphorus in oligotrophic regions is a mechanism that could potentially explain some discrepancies between observations and models (Anderson et al., 2015). However, there is no particular reason for the supply of dissolved organic phosphorus to increase in a warming climate, and even in the present climate there are net sources only in specific regions fed by transport from more productive regions (Anderson et al., 2015).

It is now accepted that a large diversity of diazotrophs exists, including heterotrophic species (Sohm et al., 2011; Zehr et al., 2011). In particular, some recent studies strongly suggest that DNF occurs in nutrient-rich but low N/P waters (e.g., the Peru and Chile upwelling regions) with rates similar to those observed in oligotrophic surface waters, and may be mainly heterotrophic (Fernandez et al., 2011, 2015). High DNF in eastern boundary current environments is consistent with the parameterizations used by GFDL and MPI, although photoautotrophy is implied because of light (GFDL) or depth (MPI) dependence.

Diazotroph adaptation to variable N/P environments can be represented by a dependence of diazotroph abundance on P concentration and DNF rate on N concentration or N/P ratio (Knapp et al., 2012). A prognostic approach (CESM, GFDL, UVic) may more realistically represent these effects than a diagnostic approach (CanESM, IPSL, MPI). In addition, a prognostic approach takes into account competition between diazotrophs and other phytoplankton groups for common resources such as P or Fe.

Another potentially important control on DNF that has barely been explored in models is CO2. Experiments on Trichodesmium sp. by Levitan et al. (2007) showed that DNF, growth, biomass and cell C content were enhanced when Trichodesmium sp. cultures were acclimated to high pCO2 (900 ppm). Fowler et al. (2015) estimated trends of DNF under higher pCO2 (750–1000 ppm) based on similar experiments with Trichodesmium sp. and Crocosphaera watsonii. Hutchins et al. (2013) estimated DNF dependence on CO2 in isolates from the Atlantic and Pacific Oceans; these results suggest an increased DNF contribution to global new N supply with high CO2, an increase that is potentially irreversible (Hutchins et al., 2013, 2015). Gradoville et al. (2014), by contrast, found that DNF activity in Trichodesmium was insensitive to pCO2 enrichment in natural seawater conditions. These authors suggested that long-term experiments or field studies would help determine whether or not high-pCO2 species will eventually dominate the diazotroph community. It is important to resolve this question and develop adequate parameterizations for models, because CO2 will probably change at least as much as any of the environmental factors considered here.

### 4.4 What are the implications of CMIP5 results for DNF and primary production in an enhanced greenhouse climate?

None of the CMIP5 models has shown that a particular set of environmental controls on DNF satisfactorily explains the spatial pattern in the observations. Furthermore, the models divide into two distinct groups, one of which associates DNF with relatively nutrient-rich regions (GFDL, MPI), while the other favours stratified, oligotrophic regions (Figure 1; Table S2). The lack of observations in the eastern boundary current waters and the generally weak spatial correlation between models and observations prevents any firm statement that one or the other group of models is more realistic. However, the limited metrics available suggest that GFDL and MPI are as skillful as any of the others, despite the bias in the distribution of observations towards habitats that the other group assumes are most important.

Warmer, more stratified oceans and a poleward expansion of such regions are robust consequences of increased greenhouse gas forcing (e.g., Xu et al., 2012). Four of the six models have DNF occurring primarily in this type of environments (Figure 1). Among these models, only CanESM shows a strong increasing trend in DNF over the 20th and 21st centuries at HOT and BATS (Figure 5c–d and Table S3). If we assume that the trend is real (Sherwood et al., 2014), we must then ask whether the differences between CanESM and the other models indicate that CanESM has captured something essentially correct in the mechanism, or whether the model having reproduced this trend is coincidental. In CanESM, the secular trend is a result of DNF being primarily a function of temperature (Bissett et al., 1999), without limitation by other environmental factors except for inhibition at high concentrations of DIN. In the other models other factors, especially phosphorus, are limiting as the area of habitat that appears generally favourable expands. The Bissett et al. (1999) model was based on the idea that other phytoplankton groups would outcompete diazotrophs in nutrient-replete waters and that diazotrophs would outcompete other phytoplankton in oligotrophic regions (Capone et al., 1997). However, it was not intended to represent the sort of long-term trend seen in CanESM under anthropogenic global warming, which would require either great metabolic flexibility with regard to N/P ratios or additional sources of P or both.

CanESM is the only model with no nutrient (P, Fe) limitation. The trends at HOT suggest that either diazotrophs are accessing new sources of P or depleting existing stocks. The ability of diazotrophs to efficiently use dissolved organic phosphorus represents a competitive advantage over non-diazotrophic phytoplankton (Dyhrman et al., 2006; Landolfi et al., 2015). However, if the dissolved organic pool is a source of ‘new’ phosphorus in a novel climate that favours diazotrophy, that pool could be depleted over time. The changes to date have been small relative to those expected in the future, and whether the kind of trends seen in CanESM can be sustained is not known.

The negative correlation between PP and sea-surface temperature (Table S4) highlights the effect of enhanced thermal stratification. Enhanced stratification combined with declining DNF in most models suggests that PP (and in some cases DNF) is increasingly limited by the supply of nutrients. In several models the trend is toward greater limitation of DNF by P, especially at BATS (Table S3).

At HOT, PP shows an upward trend over fully 25 years of observations (Figure 5; Table 4). Although this trend could result from interdecadal climate variability, its time scale approaches the range of scales for which this explanation is unlikely (e.g., Henson et al., 2010; Christian, 2014; Cummins and Masson, 2014; Rodgers et al., 2015). The PP trend has been associated with DNF, presumed to result from climate variability (Karl et al., 1995; Karl, 1999); more recently, Sherwood et al. (2014) identified enhanced DNF in this region as a long term trend. There are many weaknesses in the CanESM solution and no a priori reason to believe that it is qualitatively ‘better’ than the other models, but the increasing trend in DNF in CanESM represents a potentially large negative feedback on atmospheric CO2 growth and deserves further investigation. It is also important to consider that none of the models address the effect of CO2 on diazotrophy. In a world with much higher CO2, direct effects of CO2 on DNF may be as large as climate change effects on temperature and nutrient supply.

## 5 Summary and conclusions

Comparisons among models, and between models and observations, show that models are qualitatively consistent with the observations in terms of relative spatial distribution among the Pacific, Atlantic, and Indian Oceans (Table 3). However, the magnitude of DNF and the spatial distribution vary widely among models, and the spatial pattern correlation between modelled and observed DNF rates is weak, and in most cases not statistically significant. The frequency distribution of DNF rates (Figures 3 and 4) shows a positively skewed distribution which the models reproduce to varying degrees. The model distributions are generally narrower and less skewed than appears to be the case in the observations; i.e., the modelled space/time distribution of DNF is less variable and episodic than in the real world.

The available observational data base is an impressive achievement, as prior to about 2000 there were almost no such data, yet spatial and temporal coverage is still very incomplete. As a result, these data are insufficient to distinguish among the examined models in terms of their relative skill. In the present generation of climate models, the first order question is the global aggregate rate of DNF, and the spatial distribution is secondary. If the process is as episodic as it appears in this compilation of observations, then modelling the underlying mechanisms correctly is likely to be extremely difficult. Continuing collection of high quality observations is necessary, and should expand further into areas not traditionally considered primary areas of diazotrophy; models can continue to be tested against these observations in a variety of local, regional and global frameworks.

## Data Accessibility Statement

UVic output products are available through Michael Eby at the School of Earth and Ocean Sciences, University of Victoria. A subset of these data used for the paper is available at ftp://ftp.cccma.ec.gc.ca/pub/jchristian/DNF/.

The CMIP5 data are publically available through one of the CMIP web nodes (for instance, the Lawrence Livermore National Laboratory node at http://pcmdi9.llnl.gov), part of the Earth System Grid Federation peer-to-peer enterprise system. The Luo et al. (2012) dataset is publically available through the PANGAEA portal at https://doi.pangaea.de/10.1594/PANGAEA.774851.