Oceanic primary production decline halved in eddy-resolving simulations of global warming

The decline in ocean primary production is one of the most alarming consequences of anthropogenic climate change. This decline could indeed lead to a decrease in marine biomass and fish catch, as highlighted by recent policy-relevant reports. Because of computational constraints, current Earth system models used to project ocean primary production under global warming scenarios have to parameterize flows occurring below the resolution of their computational grid (typically 1^∘). To overcome these computational constraints, we use an ocean biogeochemical model in an idealized configuration representing a mid-latitude double-gyre circulation and perform global warming simulations under an increasing horizontal resolution (from 1 to $\mathrm{1}/\mathrm{27}$^∘) and under a large range of parameter values for the eddy parameterization employed in the coarse-resolution configuration. In line with projections from Earth system models, all our simulations project a marked decline in net primary production in response to the global warming forcing. Whereas this decline is only weakly sensitive to the eddy parameters in the eddy-parameterized coarse 1^∘ resolution simulations, the simulated decline in primary production in the subpolar gyre is halved at the finest eddy-resolving resolution (−12 % at 1/27^∘ vs. −26 % at 1^∘) at the end of the 70-year-long global warming simulations. This difference stems from the high sensitivity of the sub-surface nutrient transport to model resolution. Although being only one piece of a much broader and more complicated response of the ocean to climate change, our results call for improved representation of the role of eddies in nutrient transport below the seasonal mixed layer to better constrain the future evolution of marine biomass and fish catch potential.

1 Introduction

A decrease in global marine animal biomass and in fisheries' catch potential is projected over the 21st century under all emission scenarios, and this decrease is mostly driven by the projected decline in phytoplankton primary production in response to anthropogenic climate change (Bindoff et al., 2019). Phytoplankton are microorganisms essential to the Earth system through their influence on the global carbon cycle and biological sequestration of atmospheric CO₂ (Volk and Hoffert, 1985; Falkowski et al., 1998; Field et al., 1998). They are also at the base of most oceanic food webs and as such ultimately constrain fish biomass and fish catch potential in the ocean (Pauly and Christensen, 1995; Chassot et al., 2010; Friedland et al., 2012; Stock et al., 2017). For these reasons, net primary production (NPP) projections are at the heart of policy-relevant reports (IPCC, 2019; IPBES, 2019), whose conclusions on the critical role of the microbial world in addressing climate change and accomplishing the United Nations Sustainable Development Goals have been recently reaffirmed in “Scientists' warning to humanity”, a consensus statement (Cavicchioli et al., 2019).

Earth system models (ESMs) have become essential tools for projecting how greenhouse gas emissions will affect global biogeochemical cycling in the future (Bonan and Doney, 2018). Multi-model ensemble mean projections show that under a wide range of emission scenarios, global ocean NPP will decline in the 21st century and beyond, principally due to the reduced supply of inorganic nutrients from the sub-surface where they are abundant to the sunlit ocean where phytoplankton photosynthesis occurs (Steinacher et al., 2010; Bopp et al., 2013; Moore et al., 2018; Bindoff et al., 2019; Kwiatkowski et al., 2020). Due to strong computational constraints, an important limitation of ESMs is that transient eddies and jet-like flows of horizontal scales of 100 km and less fall below the size of their computational ocean grid and are crudely parameterized (Gent and McWilliams, 1990). Eddy-parameterized coarse-resolution (1^∘ or coarser) ESMs capture the main contrasts in nutrient distribution at the global scale (Séférian et al., 2020). But they fail to capture the full complexity of these distributions at the regional scale related to the action of eddies and fine-scale flows (McGillicuddy, 2016; Mahadevan, 2016; Lévy et al., 2018), despite ongoing effort in the development of parameterizations (Fox-Kemper et al., 2008).

Current projections of NPP decline are based on eddy-parameterized models and thus are potentially biased by their inadequate representation of sub-grid processes (Bahl et al., 2020). Here, we assess the difference that eddy resolution makes to shaping the response of NPP to future global warming. To circumvent the computational constraints, we use a biophysical model in a size-reduced setting representing a double-gyre circulation (Fig. 1). Such a configuration is a crude approximation of the mid-latitude North Atlantic or North Pacific circulations comprising subpolar and subtropical gyres separated by the Gulf Stream or the Kuroshio. We perform transient climate change simulations under an idealized global warming scenario in the form of a linear increase in atmospheric temperature, corresponding to a typical high-carbon-emission scenario in the 21st century (see the Methodology section and Fig. 1). By comparing the response to the same scenario for five eddy-parameterized coarse-resolution (1^∘) and two eddy-resolving fine-resolution ($\mathrm{1}/\mathrm{9}$ and $\mathrm{1}/\mathrm{27}$^∘) model configurations (Figs. 1 and 2), we show that the projected NPP decline under global warming is halved with the eddy-resolving configurations.

The traditional paradigm attributes declining NPP to increased vertical stratification that diminishes the vertical supply of nutrients by diapycnal mixing (Doney, 2006; Steinacher et al., 2010; Fu et al., 2016; Bindoff et al., 2019; Kwiatkowski et al., 2020). More recent understanding has shed light on the significant contribution of lateral advective supplies (Palter et al., 2005; Lévy, 2005; Lozier et al., 2011; Letscher et al., 2016; McKinley et al., 2018). This lateral nutrient transport is known as the nutrient stream (Williams et al., 2006; Palter and Lozier, 2008; Williams et al., 2011), a strong advective nutrient flux carried in sub-surface waters and associated with western boundary currents (WBCs) and Meridional Overturning Circulation (MOC). Recent results based on ESM projections have suggested that the weakening of the MOC under global warming would reduce the nutrient stream and act in concert with the reduced vertical mixing to induce a decline in NPP (Whitt, 2019; Whitt and Jansen, 2020; Tagklis et al., 2020). An explicit resolution of eddies in ocean models is known to change the position of WBCs (Chassignet and Marshall, 2008; Lévy et al., 2010; Chassignet and Xu, 2017), change the strength of the MOC (Hirschi et al., 2020; Roberts et al., 2020), and increase both stratification (Chanut et al., 2008; Lévy et al., 2010; Karleskind et al., 2011) and the advective supply of nutrients (Lévy et al., 2012a; Uchiyama et al., 2017; Uchida et al., 2020). These effects from resolving eddies have just recently begun to be explored in global warming scenarios (van Westen and Dijkstra, 2021), with model resolutions not finer than $\mathrm{1}/\mathrm{10}$^∘ at most and to our knowledge not in terms of their impacts on projected NPP.

2 Methodology

To overcome the computational constraints associated with the challenge of running eddy-resolving simulations of climate change, we adopt an idealized framework using a size-reduced setting configuration with an idealized scenario of global warming. To investigate the processes driving the decline in NPP, we base our analysis on nitrate budgets. The following sections describe (1) the experimental design, (2) the mean equilibrium states and (3) the diagnostics used in this study.

2.1 Models, configurations, simulations and experimental design

Ocean physics are simulated using the primitive-equation ocean model NEMO (Madec et al., 2017). Density is defined from salinity and temperature using a bilinear state equation. Vertical mixing is computed with a turbulent kinetic energy closure scheme, with a background value of 10⁻⁵ m² s⁻¹ and local enhancement to 100 m² s⁻¹ in the case of convection (unstable density profile). For biogeochemistry, we used the LOBSTER model (Lévy, 2005; Lévy et al., 2012b), implemented in NEMO, which uses nitrogen as a currency and computes the evolution of six tracers: phytoplankton, zooplankton, dissolved organic material, nitrate, ammonium and sinking detritus. Note that in the LOBSTER model, none of the biological rates are directly dependent on temperature, allowing us to focus our analysis on effects due to changes in physical transport.

Figure 1 (a) Barotropic circulation over the model domain (average of the five 1^∘ resolution preindustrial control simulations; see also Fig. A2 for individual configurations) and analytical external forcings as a function of latitude: (b) wind stress (N m⁻²), and (c) pseudo-atmospheric temperature (^∘C). The forcings vary seasonally in a sinusoidal manner between winter (blue) and summer (orange) extrema. In the climate change simulations, the pseudo-atmospheric temperature increases linearly in time between the dashed and solid lines (c). After 70 years of simulation, the linear atmospheric temperature increase reaches +2.8 ^∘C. (d) Time series of annual ocean temperature (^∘C) averaged between 0 and 700 m depth over the model domain. Years −100 to 0 correspond to the spin-up period under seasonal preindustrial forcing. The preindustrial control and climate change simulations (CTL and CC, respectively) discussed in this paper start at year 0. The 1^∘ resolution time series shows the average of the five 1^∘ configurations. Shading indicates ±1 standard deviation. See Fig. A1 for individual 1^∘ model configurations. The small peaks in the time series of the $\mathrm{1}/\mathrm{27}$^∘ resolution spin-up around years −80 and −60 are due to a problem during the chaining of the spin-up jobs (which were set by periods of 2 years). These errors did not influence the equilibrium states, and we decided not to rerun the spin-up simulations to save computational resources.

The model is run in a configuration adapted from previous studies (Krémeur et al., 2009; Lévy et al., 2012b; Resplandy et al., 2019). The domain geometry is a closed square basin on the β plane centered at ∼30^∘ N, 3180 km wide and long and 4 km deep. The domain is bounded by vertical walls and by a flat bottom with a free slip boundary condition. The sea surface is kept free. The circulation is forced by analytical zonal distributions of wind stress, net heat flux and freshwater flux. The wind and buoyancy forcings are such that a double-gyre circulation is set up (Figs. 1 and A2). The wind stress is zonal and varies latitudinally in a sinusoidal manner between the extrema at the edges and the middle of the domain. The net heat flux takes the form of a restoration toward a zonal apparent air temperature profile. A portion of the net heat flux comes from solar radiation and is allowed to penetrate within the water column. A freshwater flux is also prescribed and varies zonally. It is determined such that, at each time step, the basin-integrated flux is zero.

We use three horizontal resolutions: 106 km (1^∘), 12 km ($\mathrm{1}/\mathrm{9}$^∘) and 4 km ($\mathrm{1}/\mathrm{27}$^∘). For each resolution, time steps, numerical schemes and isopycnal/horizontal diffusion are adapted (see Table 1). Note that diffusion is along the isopycnals only for the 1^∘ coarse resolution. In the following, we will nevertheless use the term isopycnal diffusion regardless of the resolution for simplicity and because the diffusion is null or minimal at finer resolutions. On the vertical, 30 z levels are used whose thickness varies from 10 m at the surface to 500 m at depth. The upper 120 m contains 10 vertical levels; the upper 500 m contains 20 vertical levels. For the 1^∘ resolution configurations, we used the Gent and McWilliams (1990, GM hereafter) eddy parameterization. This parameterization relies on two coefficients, an isopycnal diffusion coefficient (k_iso) and a GM coefficient (k_gm) that translates into a bolus velocity that acts to flatten isopycnals. For testing the sensitivity to the GM parameterization, we used five combinations of the isopycnal diffusion and GM coefficients: (1) 500 m² s⁻¹, (2) 1000 m² s⁻¹ and (3) 2000 m² s⁻¹ for both parameters and (4) 500 m² s⁻¹ and (5) 2000 m² s⁻¹ for the isopycnal diffusion parameter but keeping the GM coefficient at 1000 m² s⁻¹. We thus end up with seven different configurations: five eddy-parameterized at a coarse resolution (1^∘) and two eddy-resolving at fine resolutions ($\mathrm{1}/\mathrm{9}$ and $\mathrm{1}/\mathrm{27}$^∘). In the following, results from the eddy-parameterized coarse-resolution configurations are synthesized by showing the average ±1 standard deviation across the five different configurations. Contrary to the $\mathrm{1}/\mathrm{27}$^∘ configuration, the qualifier “eddy-permitting” is probably more appropriate for the $\mathrm{1}/\mathrm{9}$^∘ configuration. Nevertheless to simplify and as the emphasis is put on the differences between the 1^∘ resolution and the finer ones, we use the term eddy-resolving for both.

Table 1 Resolution-dependent model features and parameters. For details on numerical schemes see Madec et al. (2017) and Lévy et al. (2001). Note that we added a minimum level of bi-Laplacian tracer diffusivity at $\mathrm{1}/\mathrm{27}$^∘ to insure numerical stability, and this was not needed at the $\mathrm{1}/\mathrm{9}$^∘ resolution. This bi-Laplacian diffusion acts on the horizontal unlike the Laplacian diffusion acting along isopycnals in the 1^∘ resolution simulations.

The basin is initialized at rest with vertical profiles of temperature and salinity uniformly applied to the whole domain and homogeneous concentrations of the biogeochemical tracers. A 2000-year spin-up with an eddy-parameterized coarse-resolution configuration (1^∘; $k_{\mathrm{gm}}=k_{\mathrm{iso}}=\mathrm{1000}$ m² s⁻¹) is then performed under preindustrial forcings (wind, air temperature and freshwater flux) varying seasonally in a sinusoidal manner between winter and summer extrema and repeated every year. Then, for each configuration, a climate change and a preindustrial control simulation are carried out from this initial coarse-resolution spin-up state in two steps (Figs. 1 and A1):

1. Step 1 comprises 100 years of spin-up under the same preindustrial control forcings, in order to reach an equilibrium for each configuration. These seven additional spin-ups (one for each configuration) are initialized with the annual mean coarse-resolution spin-up state. It is important to note that the final dynamical and biogeochemical equilibrium is different for each configuration (Figs. 1, 2, A1, A2, A4). The following section briefly describes these equilibrium and their differences.

2. Step 2 comprises a 70-year simulation for each configuration, following two scenarios. The preindustrial control scenario is the continuation of the spin-up; forcings are kept seasonal. In the climate change scenario, for air temperature, in addition to seasonal variations, we impose a linear trend of +0.04 ^∘C yr⁻¹ ($\sim +\mathrm{2.8}$ ^∘C after 70 years). This trend roughly corresponds to the trend in surface air temperature of the simulations following the SSP5-8.5 scenario (Tokarska et al., 2020).

   In total, 14 simulations of 70 years were run: seven configurations × two scenarios. In the following we will compare the differences between the eddy-resolving simulations ($\mathrm{1}/\mathrm{9}$ and $\mathrm{1}/\mathrm{27}$^∘ resolutions) and the mean of the eddy-parameterized simulations (1^∘ resolution). A notable point is that after year 0, the residual drifts (CTL) are much smaller than the climate change trends (CC), which allows us to draw conclusions on the effect of the climate change forcing despite the fact that the spin-up is not yet fully achieved.

It is noteworthy that this experimental design implies that each climate change simulation starts from a different initial state. This strategy has the advantage that each initial state represents an equilibrated control state for the given set of resolution and parameter choices. For each simulation, the impact of climate change is evaluated as the difference between the climate change and preindustrial control simulations, and this difference depends on model choices but may also depend on the differences in the control simulations. Another option would have been to start the climate change simulations from the same initial conditions, but this would have had the consequences of a strong drift at the start of the eddy-resolving climate change simulations (as can be seen in the first 40 years of the spin-up, Fig. 1). Our choice is also consistent with the fact that shaping different equilibrium states is part of the role of the eddies. We can also note that none of the eddy-parameterized configurations that we explored has achieved an equilibrium state that comes close to the eddy-resolving configurations (Fig. A1).

Figure 2 (a, b, c) Surface relative vorticity (% of planetary vorticity), (d, e, f) vertically integrated net primary production (NPP, $\mathrm{mmol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) and (g, h, i) vertical distribution of nitrate (zonal mean, mmol N m⁻³) in the preindustrial control simulations at (a, d, g) 1^∘, (b, e, h) $\mathrm{1}/\mathrm{9}$^∘ and (c, f, i) $\mathrm{1}/\mathrm{27}$^∘ resolution. Vorticity is computed with 2 d averaged velocity fields. NPP and nitrate distribution are averaged over the last 5 years of the simulations (years 66 to 70). The 1^∘ resolution distribution is the mean of the five 1^∘ configurations (see Figs. A4 and A2 for individual distributions).

2.2 Equilibrium states

As noted above, the dynamical and biogeochemical equilibrium states differ for each model configuration. These differences have been more extensively discussed in a similar setup in Lévy et al. (2010, 2012b) and Lévy and Martin (2013) and are briefly discussed in the following. It can be noted that the differences between the coarse-resolution equilibrium states on one hand and the fine-resolution equilibrium states on the other hand are much larger than the differences within the five (two) coarse-resolution (fine-resolution) spin-ups (Fig. 1).

The main features of the model's solution comprise a western boundary current separating an oligotrophic subtropical gyre in the south of the domain from a productive subpolar gyre in the north (Figs. 2 and A2). As resolution increases, mesoscale eddies and filamentary structures emerge in the relative surface vorticity field (Fig. 2a, b, c). In all simulations, the barotropic circulation is characterized by an anticyclonic circulation in the south and a cyclonic circulation in the north (Fig. A2). But as resolution increases, the non-linear effect become more important and small recirculation gyres appear close to the domain's northern and southern boundaries.

The divergence of the horizontal Ekman flux induces downwelling and a depressed nitracline in the subtropical gyre, as well as upwelling and a raised nitracline in the subpolar gyre (Figs. 2 and A4). The model's nitracline follows the thermocline very closely (Fig. A10). Nitrate supplies to the euphotic layer are seasonal and driven by the seasonal deepening of the mixed layer, which is strongest in the north of the domain (Fig. A9). This seasonal supply leads to an intense spring bloom in the subpolar gyre and to a weaker winter bloom in the subtropical gyre, where nitrate supplies are lower (Figs. 2 and A2). Over an annual average, this leads to a strong contrast between the subpolar gyre, where annual mean NPP is large, and the subtropical gyre characterized by low levels of annual mean NPP. We can also note that nitrate gradients in the thermocline are weaker at fine resolutions than at the coarse resolution (Fig. 2), in accordance with weaker stratification (Fig. A10). This is associated with lower sub-surface nitrate concentrations. We can note the opposed effect between $\mathrm{1}/\mathrm{9}$ and $\mathrm{1}/\mathrm{27}$^∘ (Fig. 2) with a slightly larger sub-surface nitrate concentration at $\mathrm{1}/\mathrm{9}$^∘ compared to $\mathrm{1}/\mathrm{27}$^∘. These differences in equilibrated nitrate distribution lead to differences in equilibrated NPP between our different models runs (Fig. 2), but importantly the strong contrast between the two gyres is present is all simulations.

2.3 Nitrate budgets

In our simulation, NPP is constrained by nitrate supplies from the sub-surface into the sunlit surface layers where photosynthesis can occur. Nitrate budgets are thus used to investigate the processes driving changes in NPP. We focus our analysis on a predefined region within the model's subpolar gyre where the strongest NPP decline is observed (Figs. 2 and 3). This subpolar box is shown by the dashed horizontal black lines in Fig. 3; it is bounded by the change of sign in the wind stress curl forcing in the south and excludes the northern boundary. The local nitrate budget can be expressed as $-\mathbf{\nabla }\cdot (\mathit{u}\cdot N)+L\left(N\right)+{\partial }_z(k\cdot {\partial }_{z}N)+S\left(N\right)={\partial }_{t}N$, where $\mathbf{\nabla }\cdot (\mathit{u}\cdot N)$ is the divergence of the advective fluxes, ${\partial }_z(k\cdot {\partial }_{z}N)$ is the vertical diffusion term, L(N) is isopycnal diffusion, and S(N) represents the biological nitrate sources (regeneration) and sinks (NPP). u is the total velocity and includes the bolus velocity of the GM parameterization at the coarse resolution. Integrating this equation over the subpolar box and depth D leads to

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& -\stackrel{\text{Advective fluxes}}{\overbrace{\oint \mathit{u}\cdot N\mathrm{d}s}}+\stackrel{\text{Isopycnal mixing}}{\overbrace{\underset{\mathrm{0}}{\overset{D}{\int }}\langle L\left(N\right)\rangle \mathrm{d}z}}+\stackrel{\text{Vertical mixing}}{\overbrace{\langle k\cdot {\partial }_{z}N\rangle }}\\ & +\stackrel{\text{Biological sources/sinks}}{\overbrace{\underset{\mathrm{0}}{\overset{D}{\int }}\langle S\left(N\right)\rangle \mathrm{d}z}}={\partial }_t\left(\underset{\mathrm{0}}{\overset{D}{\int }}\langle N\rangle \mathrm{d}z\right).\end{array}\end{array}$$

The angle brackets stand for the horizontal integral over the subpolar box. The first term on the left hand side is the integral of the advective fluxes entering/exiting the subpolar box between the surface and depth D. It is obtained from the volume integral of the advective fluxes' divergence using Stoke's theorem. The subpolar box being closed at the east and west, it is computed as the sum of three nitrate advective fluxes: (1) the meridional flux entering through the southern boundary of the box between the surface and depth D, (2) the vertical flux entering at depth D, and (3) the meridional flux going out through the northern boundary of the box between the surface and depth D. Of interest here are the advective fluxes supplying nitrate and thus the fluxes entering the subpolar box, i.e., the vertical flux and the meridional flux through the southern boundary (Fig. A5). The total transport term shown in Fig. 4 represents these entering fluxes and contains the sum of the vertical flux and of the meridional flux at the southern boundary. It also contains isopycnal mixing, which is of second order (Fig. A5).

In the preindustrial control simulations, above ∼50 m depth, nitrate supplies to the subpolar box are dominated by vertical mixing, whereas deeper than ∼50 m they are dominated by the transport term (Fig. 4b). As resolution increases, this separation continues to hold although the respective intensities are modified: advective fluxes are stronger at the $\mathrm{1}/\mathrm{9}$ and $\mathrm{1}/\mathrm{27}$^∘ resolutions than at the 1^∘ resolution (Fig. 4c, d).

The separation between mean and eddy transport in this budget was achieved following the Reynolds averaging method. u and N are broken into a mean component ($\stackrel{\mathrm{‾}}{u}$ and $\stackrel{\mathrm{‾}}{N}$) and a fluctuating eddy component (u^′ and N^′ with $\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }}=\mathrm{0}$ and $\stackrel{\mathrm{‾}}{N^{\prime }}=\mathrm{0}$), leading to

$$\begin{array}{}\text{(2)}& {\displaystyle }\oint \stackrel{\mathrm{‾}}{\mathit{u}\cdot N}\mathrm{d}s=\underset{\text{Mean}}{\underbrace{\oint \stackrel{\mathrm{‾}}{\mathit{u}}\cdot \stackrel{\mathrm{‾}}{N}\mathrm{d}s}}+\underset{\text{Eddy}}{\underbrace{\oint \stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }\cdot N^{\prime }}\mathrm{d}s}}.\end{array}$$

The overbar is an averaging operator (in this study a spatio-temporal annual-mean averaging over boxes of 1^∘) used to separate the mean from the fluctuating component over a temporal and/or spatial scale larger than the typical eddy scales. With this operator, we decomposed the total advection term into a mean part including advection by mean velocities ($\oint \stackrel{\mathrm{‾}}{\mathit{u}}\cdot \stackrel{\mathrm{‾}}{N}\mathrm{d}s$) and an eddy part ($\oint \stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }\cdot N^{\prime }}\mathrm{d}s$), which we computed as the residual between total and mean advection.

This eddy–mean separation is by nature imperfect and leads to mean and eddy components that remain heterogeneous (Fig. A11) but is nevertheless useful to obtain an estimate of the contribution of each component. At a 1^∘ resolution, the explicit (Reynolds) eddy term is close to zero and we estimated the parameterized eddy component as the sum between the two terms included in the GM parameterization: advection by the bolus velocity and isopycnal mixing (Figs. 4, A5 and A6).

Finally, the change in nitrate transport in response to global warming is due to changes in nitrate distribution as well as changes in circulation. In order to evaluate these two contributions, we computed offline diagnostics of nitrate advection, using velocities from the preindustrial control simulations (u_CTL) and nitrate from the climate change simulations (N_CC) or vice versa (u_CC and N_CTL). More precisely, we used the following decomposition of the change in nitrate advective fluxes, $\mathrm{\Delta }(\mathit{u}\cdot N)=\oint {\mathit{u}}_{\mathrm{CC}}\cdot N_{\mathrm{CC}}\mathrm{d}s-\oint {\mathit{u}}_{\mathrm{CTL}}\cdot N_{\mathrm{CTL}}\mathrm{d}s$:

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }(\mathit{u}\cdot N)=\stackrel{\mathrm{\Delta }\text{Nitrate}}{\overbrace{{\mathit{u}}_{\mathrm{CTL}}\cdot \mathrm{\Delta }N}}+\underset{\mathrm{\Delta }\text{Circulation}}{\underbrace{\mathrm{\Delta }\mathit{u}\cdot N_{\mathrm{CTL}}}}+\stackrel{\text{Non-linear}\mathrm{\Delta }}{\overbrace{\mathrm{\Delta }\mathit{u}\cdot \mathrm{\Delta }N}},\end{array}$$

with

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathit{u}}_{\mathrm{CTL}}\cdot \mathrm{\Delta }N=\oint {\mathit{u}}_{\mathrm{CTL}}\cdot N_{\mathrm{CC}}\mathrm{d}s-\oint {\mathit{u}}_{\mathrm{CTL}}\cdot N_{\mathrm{CTL}}\mathrm{d}s,\text{(5)}& {\displaystyle }\mathrm{\Delta }\mathit{u}\cdot N_{\mathrm{CTL}}=\oint {\mathit{u}}_{\mathrm{CC}}\cdot N_{\mathrm{CTL}}\mathrm{d}s-\oint {\mathit{u}}_{\mathrm{CTL}}\cdot N_{\mathrm{CTL}}\mathrm{d}s,\text{(6)}& {\displaystyle }\mathrm{\Delta }\mathit{u}\cdot \mathrm{\Delta }N=\mathrm{\Delta }(\mathit{u}\cdot N)-\mathrm{\Delta }\mathit{u}\cdot N_{\mathrm{CTL}}-{\mathit{u}}_{\mathrm{CTL}}\cdot \mathrm{\Delta }N.\end{array}$$

3 Results

3.1 Sensitivity of the decline in net primary production to model resolution

After 70 years of forcing by a linear increase in atmospheric temperature (reaching +2.8 ^∘C, similar in magnitude to the global temperature increase at the end of 21st century in the SSP5-8.5 scenario; Tokarska et al., 2020), NPP has declined in all model simulations (Fig. 3). The simulated declines are particularly strong in the model's subpolar gyre, a region with elevated initial levels of NPP (Fig. 2) and where we will now focus the analysis. The integrated NPP over the model's subpolar box decreases by $-\mathrm{0.83}\pm \mathrm{0.03}$, −0.51 and −0.44 $\mathrm{mmol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ in our three model resolution versions, representing $-\mathrm{26}\pm \mathrm{1}$ %, −14 % and −12 %, with the horizontal resolution set at 1, $\mathrm{1}/\mathrm{9}$ and $\mathrm{1}/\mathrm{27}$^∘, respectively (Fig. 3). Thus the NPP declines differ markedly between each other depending on the resolution, with the decline in NPP halved at the finest resolution compared with the decline projected at the coarse resolution. Locally, the intensity of the decline in the coarse-resolution model's versions exceeds −1.5 $\mathrm{mmol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. The intensity and distribution of these local changes slightly vary among the five coarse-resolution configurations. But with a finer grid resolution, they are strongly reduced (Fig. 3).

Figure 3 Projected decrease in net primary production (NPP, $\mathrm{mmol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) under a high-emission scenario in (a, b, c, d, e) five 1^∘, (f) one $\mathrm{1}/\mathrm{9}$^∘ and (g) one $\mathrm{1}/\mathrm{27}$^∘ model simulations. The NPP change is vertically integrated over the entire water column and taken as the difference between the last 5 years of the climate change and preindustrial control simulations (years 66 to 70). The dashed horizontal black lines delineate the model's subpolar box over which diagnostics for (h) and Figs. 4, 5, 6 and 7 are performed. This is where simulated NPP is strongest (Fig. 2) and where the simulated NPP decline is also strongest. (h) Time evolution of the NPP decrease in the model's subpolar box in response to climate change, for the three model resolutions. The decrease is computed relative to the preindustrial control simulations. The 1^∘ resolution time series shows the average of the five 1^∘ configurations. Shading indicates ±1 inter-model standard deviation. See Fig. A1 for individual model simulations.

Figure 4 Nitrate concentrations (mmol N m⁻³) and nitrate supplies ($\mathrm{mmol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) averaged over the model's subpolar box (dashed horizontal black lines in Fig. 3). (a) Nitrate profiles in the preindustrial control simulations, for the three resolutions. (b, c, d) Dynamical nitrate supplies, cumulated between the surface and depth D in the preindustrial control simulations, against D. (e, f, g, h) Changes in these profiles between the climate change and preindustrial control simulations. In the coarse-resolution simulations (CTL1/CC1), the profiles show the average of the five 1^∘ configurations. Shading indicates ±1 standard deviation. The total transport term (solid blue line) comprises the sum of vertical advection across depth D and meridional advection through the southern border of the model's subpolar box; in the case of the 1^∘ simulations, total transport also includes eddy parameterizations (including advection by bolus velocity and isopycnal mixing). Total transport is the sum of mean and eddy transport. In the fine-resolution simulations (CTL9/CC9, CTL27/CC27), the mean and eddy transport terms are derived from Reynolds decomposition (see the Methodology section). In the coarse-resolution simulations (CTL1/CC1), the mean transport represents the resolved transport by the coarse-resolution velocity field, and the eddy transport includes the two eddy parameterization terms. See also Figs. A5 and A6 for individual fluxes. All profiles are averaged over years 66 to 70.

3.2 Decline in nutrient supply pathways

The decline in NPP in the model's subpolar box stems from the decline in the supply of nutrients supporting phytoplankton growth, which is associated with the decline in nitrate concentrations in the thermocline (Figs. 4a, e and A4). About half of the NPP decrease results from a decline in NPP supported by newly supplied nutrients, which are mainly in the form of nitrate, with the remaining half arising from decreased NPP supported by locally regenerated nutrients (mostly in the form of ammonium, Fig. A3). The mean f ratio on the subpolar box is about 0.43 all along the simulations (Fig. A3, fourth column) albeit with a very small decline (from ∼0.43 to ∼0.40) in the climate change simulations. This slight decline is consistent with the lower primary production regime having a lower f ratio (Sarmiento and Gruber, 2006, p. 166). We should also note that changes in the f ratio in response to warming might be underestimated here because the direct impact of increasing temperature on biological rates, which has been shown to affect the response of biology to climate change (Olonscheck et al., 2013; Lewandowska et al., 2014; Laufkötter et al., 2015), is not accounted for in this study.

A closer examination of nitrate trends within the model's subpolar box provides insights into the physical mechanisms controlling the declining nitrate concentrations in the thermocline in our coarse-resolution simulations (Fig. 4). The traditional paradigm (Steinacher et al., 2010; Fu et al., 2016; Bindoff et al., 2019; Kwiatkowski et al., 2020) for the decrease in NPP is confirmed. Indeed in the euphotic layer (top 150 m), there is a strong reduction in nutrient supply by vertical mixing (Fig. 4f, orange line). This decrease is associated with increasing stratification (Fig. A10). However, similarly to what has been shown for the North Atlantic in coarse-resolution ESM projections (Whitt, 2019; Whitt and Jansen, 2020; Tagklis et al., 2020), the decrease in nutrient supply is also controlled by the sub-surface nutrient reservoir and the fluxes fueling it. Nitrate concentrations in the top 400 m strongly decrease under our global warming scenario (Fig. 4e) because of a reduction in nutrient advection. Under global warming, the net advective supply of nitrate in the sub-surface (between 100 and 400 m) is strongly reduced (Fig. 4f, blue line; see also each component of the advective flux in Fig. A6).

As the resolution increases from 1 to $\mathrm{1}/\mathrm{9}$ and $\mathrm{1}/\mathrm{27}$^∘, the mechanisms described previously still hold. However, despite similar mean states for the three resolutions at equilibrium (preindustrial control simulations; Fig. 4b, c, d), the climate-change-induced decline in vertical mixing to the euphotic layer and in sub-surface advection are both dampened when resolution increases (Fig. 4f, g, h). At the coarse resolution (1^∘), vertical mixing decreases by 48±4 %, whereas it only decreases by 18±1 % in the finer-resolution models (Fig. 7). Likewise, nutrient advection decreases by 83±5 % at the coarse resolution compared to only 27±2 % at the finer resolutions. A large part of the attenuation occurs when the resolution is changed from 1 to $\mathrm{1}/\mathrm{9}$^∘, but the attenuation is much more moderate between $\mathrm{1}/\mathrm{9}$ and $\mathrm{1}/\mathrm{27}$^∘. This is consistent with the expectation that resolution of mesoscale dynamics should be nearly achieved at $\mathrm{1}/\mathrm{9}$^∘ in the latitudinal range considered here but not fully (see also Figs. 1 and A1). Moreover, we expect further differences to emerge with grids finer than $\mathrm{1}/\mathrm{27}$^∘ which would allow resolving sub-mesoscale dynamics (Lévy et al., 2010; Gula et al., 2016; McWilliams, 2016; Uchida et al., 2019).

3.3 On the role of local eddy processes

A legitimate question is whether the differences between our coarse-resolution and fine-resolution models are due to local small-scale eddy processes not adequately captured by the eddy parameterization or to different responses of the mean large-scale state to global warming. To address this question, we separated the mean and eddy transport of nitrate, in the two eddy-resolving simulations following Lévy and Martin (2013), with a spatio-temporal annual-mean averaging over boxes of 1^∘. In the coarse-resolution simulations, eddy transport is parameterized by a bolus advection term (Gent and McWilliams, 1990) and isopycnal mixing.

This analysis reveals the dominant role of the mean advection in fueling the sub-surface nitrate reservoir at equilibrium in the two eddy-resolving models (preindustrial control simulations; dashed blue lines in Fig. 4c, d). Eddy advection supplies a minor addition of nutrients in the subpolar box (mostly at $\mathrm{1}/\mathrm{27}$^∘), essentially through the eddy meridional advection (Fig. A5). In the eddy-parameterized models, the eddy parameterization strongly opposes the resolved transport, removing nutrients from the surface of the subpolar box mainly through the vertical bolus advection term (Figs. 4b and A5).

In the eddy-resolving climate change simulations, most of the decline in nutrient input is explained by a reduction in mean nitrate transport (Fig. 4g, h). In the eddy-parameterized simulations, the eddy component plays a larger role and explains approximately one-third of the total transport decrease at 200 m (Fig. 4). As for the preindustrial control simulation, it is dominated by changes in the vertical bolus advection term (Fig. A6).

Figure 5 Role of circulation change versus nitrate distribution change in explaining the decline in the nitrate advection flux ($\mathrm{mmol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) in response to global warming in the model's subpolar box, cumulated between the surface and depth D, for the (a) 1^∘, (b) $\mathrm{1}/\mathrm{9}$^∘ and (c) $\mathrm{1}/\mathrm{27}$^∘ resolution configurations. The black lines show the simulated change, which is decomposed into three components: the red lines show the component related to the change in nitrate distribution; the blue lines show the component related to the change in circulation; the gray lines show the residual due to non-linear changes. See Eqs. (3) to (6) for details. In (a), the profiles show the average of the five 1^∘ configurations. Shading indicates ±1 standard deviation.

3.4 On the role of circulation changes

Another question is how changes in nitrate transport (black lines in Fig. 5) are related to changes in circulation. To examine this question, we first computed in offline mode the advective transport of nitrate keeping the nitrate distribution fixed to its preindustrial state. By doing so, we accounted only for changes in circulation and neglected changes in nitrate distribution (blue lines in Fig. 5). In a similar manner, we kept the circulation constant and accounted only for the changes in nitrate distribution (red lines in Fig. 5). Due to non-linearities, the sum of circulation-induced (blue lines) and nitrate-induced (red lines) changes is not equal to the total change (black lines) and a residual remains (gray lines in Fig. 5).

In the euphotic layer (upper 0–100 m), the change in nitrate transport is mostly attributable to a change in nitrate concentrations, at all resolutions. But deeper in the water column, the situation is reversed and the change in circulation is responsible for most of the change. In fact at 400 m depth, the change in nitrate transport is largely explained by circulation changes in all of our simulations. This result highlights the crucial role of circulation changes, which directly affect the resupply of the sub-surface nitrate reservoir (between 100–400 m), with the consequence of affecting the mean nitrate gradients, particularly in the upper layers (0–100 m, Fig. 4e) and thus the nitrate fluxes closer to the surface. We can note a stronger contribution of non-linearities in this simple analysis as the model resolution increases.

Figure 6 Vertical profiles of (a) the northward and (c) upward water transport (in Sv) entering the model's subpolar box (dashed horizontal black lines in Fig. 3) in the preindustrial control simulations for each resolution. (b, d) Changes in these profiles between the climate change and preindustrial control simulations. The thick black lines show the average of the five 1^∘ configurations. Shading indicates ±1 standard deviation.

To further link our results to circulation changes, we then examined how the water transport through the subpolar box was modified in our climate change simulations (Fig. 6). In our model, the main transport pathway to the subpolar box originates from upwelling and from northward advection at the sub-surface (between 50 and 400 m) (Fig. 6a, c). With warming, these upward and sub-surface northward transports both decline (Fig. 6b, d). The mean transports are similar in the equilibrium states between the coarse- and fine-resolution simulations, but they respond differently to warming. At the coarse resolution, the decline in the upward transport is particularly strong, −0.5 Sv at ∼100 m depth (Fig. 6b). This is in phase with the strong decline in nitrate vertical advection at a similar depth discussed above (blue lines in Fig. A6). At finer resolutions, the decline in the vertical flow is close to zero at ∼100 m, but it increases with depth, reaching an amplitude similar to the coarse-resolution simulations at ∼400 m, $\sim -\mathrm{0.5}$ Sv. In the finer-resolution simulations, the sub-surface northward flow decline is also weaker, $\sim -\mathrm{0.15}$ Sv against $\sim -\mathrm{0.3}$ Sv at the coarse resolution below 100 m depth (Fig. 6d).

Previous coarse-resolution modeling studies (Whitt, 2019; Whitt and Jansen, 2020; Tagklis et al., 2020) have suggested that the climate-change-driven decline in North Atlantic NPP is strongly linked to changes in the Meridional Overturning Circulation (MOC). This link presumes that the slowdown in the upper branch of the MOC is associated with a decrease in the strength of the portion of the Gulf Stream that conveys nutrients northward (Whitt, 2019). In our simplified model configuration, the MOC (here diagnosed by the meridional stream function in the z coordinate, Fig. A8) transports waters from the subtropical gyre to the north at a depth between ∼100 and ∼400 m: 1.75±0.16 Sv at the 1^∘ resolution, 3.14 Sv at $\mathrm{1}/\mathrm{9}$^∘ and 2.94 Sv at $\mathrm{1}/\mathrm{27}$^∘ (flow through the vertical black lines in Fig. A8). Once these waters reach the northern boundary of the domain, they are entrained downward by convection to depths of 1000 m and more (Fig. A9) before they return southward toward the southern boundary where they are upwelled back to the surface. In our coarse-resolution simulations, this overturning cell is strongly slowed down under global warming, inducing in particular a strong reduction, $-\mathrm{1.79}\pm \mathrm{0.22}$ Sv ($-\mathrm{102}\pm \mathrm{13}$ %), in the water transport between the two gyres (Fig. A8). The slowing of the overturning circulation is associated with a decrease in the maximum depth of the mixed layer (MLD) in the northern part of the domain where convection occurs (Fig. A9). The maximum MLD shoals from over 2000 to ∼300 m after 70 years of climate change forcing (e.g., Fig. A9c). At finer resolutions, the MOC is also slowed down, but this slowdown is less pronounced (Fig. A8). The flow between 100 and 400 m between the two gyres is reduced by −0.58 and −0.85 Sv (−18 % and −29 %) at $\mathrm{1}/\mathrm{9}$ and $\mathrm{1}/\mathrm{27}$^∘ resolution, respectively. It is also exemplified by the maximum depth of the mixed layer in the northern part of the domain which shoals less than it does at the coarse resolution, from ∼1000 to ∼700 m (Fig. A9r, u).

The different responses of the MOC between our coarse- and fine-resolution simulations can be partly related to the differences in stratification in the north of the domain. The fine-resolution equilibrium states are initially less stratified than the coarse-resolution ones (Fig. A10, first column). As a consequence, surface heating during the transient climate warming for the fine-resolution simulations penetrates deeper into the interior and leads to a weaker increase in stratification (Fig. A10, last column) and thus a weaker decrease in convection and in meridional circulation.

Figure 7 Schematic representation of the projected response of net primary production (NPP) and of the nutrient fluxes (advection and vertical mixing) supporting it to warming, with our (a) eddy-parameterized coarse-resolution (1^∘) and (b) eddy-resolving fine-resolution ($\mathrm{1}/\mathrm{9}$, $\mathrm{1}/\mathrm{27}$^∘) simulations, in the model's subpolar gyre. Decline in nutrient advection is assessed from changes between the surface and 200 m depth (Tot. trp. in Fig. 4). Decline in vertical mixing is evaluated at 25 m depth, which is where the decrease is maximal. NPP change is vertically integrated over the entire water column. Changes at the coarse resolution are the average of the five 1^∘ configurations ±1 standard deviation. Changes at the fine resolutions are the average of the $\mathrm{1}/\mathrm{9}$ and $\mathrm{1}/\mathrm{27}$^∘ configurations ± half of the difference between the two. The background diagrams in gray represent the initial control situation, while the colored diagrams represent the situation after 70 years of warming.

4 Conclusions

Based on a wind- and buoyancy-driven double-gyre model configuration with an idealistic scenario of global warming, we have shown that the projected decline in NPP is strongly sensitive to horizontal grid resolution, with a decline which is halved at an eddy resolution: $-\mathrm{13}\pm \mathrm{1}$ % versus $-\mathrm{26}\pm \mathrm{1}$ % at the coarse resolution (Fig. 7). Moreover, we have also found that this sensitivity is much stronger than the sensitivity to a large set of eddy parameterization coefficients (k_iso and k_gm). This result is due to the much weaker decline in nutrient transport at fine resolutions, $-\mathrm{27}\pm \mathrm{2}$ % versus $-\mathrm{83}\pm \mathrm{5}$ % at the coarse resolution, combined with a weaker decrease in vertical mixing, $-\mathrm{18}\pm \mathrm{1}$ % vs. $-\mathrm{48}\pm \mathrm{4}$ % (Fig. 7).

A key result of our study is that the resolution-related uncertainties in the NPP projections stem from the high sensitivity of mean ocean transport projections to resolution. This sensitivity is currently being debated. Regarding the Atlantic MOC (AMOC), some studies show a weaker decline as resolution increases (Roberts et al., 2004; Delworth et al., 2012; as with our model MOC), whereas other studies show an amplified response (Weijer et al., 2012; Spence et al., 2013; Jackson et al., 2020; Roberts et al., 2020). It has been suggested that the intensity of the AMOC decline may be related to its intensity in the control climate: the stronger the initial AMOC, the stronger the decline (Gregory et al., 2005; Winton et al., 2014; Jackson et al., 2020; Roberts et al., 2020). The opposite is found here: at fine resolutions, the model's MOC is stronger in the control climate but the response is dampened in the global warming scenario. The AMOC intensity in the control climate itself has been shown to be sensitive to resolution but with no consensus yet on the sign of this sensitivity (Winton et al., 2014; Hirschi et al., 2020). The reason for the sensitivity of control and projected AMOC to resolution remains unclear and differs between studies. The potential causes mentioned are stronger air–sea interactions at fine resolution (Hirschi et al., 2020), different spatial distributions of perturbations by the eddies (Spence et al., 2013) or the introduction of biases (Delworth et al., 2012).

Another possible reason for the high sensitivity to resolution of projected ocean circulation changes could be derived from the “eddy saturation” and “eddy compensation” phenomena which are actively studied in the Southern Ocean (Farneti et al., 2010; Farneti and Delworth, 2010; Meredith et al., 2012; Munday et al., 2013; Mak et al., 2017). They stand for the low sensitivity of stratification and overturning circulation, respectively, to an increase in wind forcing in eddy-resolving models. The connection between the two phenomena is subject to ongoing research (Meredith et al., 2012). Eddy saturation results from the balance between the steepening of the isopycnal by the wind stress and the flattening by eddies (Mak et al., 2017). One could imagine that similar competing processes are at work in our eddy-resolving simulations of climate change, warming stratifying the ocean versus eddies destratifying it. One clue in that direction is that, as in Munday et al. (2013), surface eddy kinetic energy increases with global warming in the eddy-resolving simulations. Eddy saturation is poorly represented by the GM parameterization (Farneti et al., 2010; Munday et al., 2013, 2014); however simple modification of this parameterization can lead to the emergence of the phenomenon (Mak et al., 2017).

The GM parameterization has improved the representation of circulation in ocean models (Gent and McWilliams, 1990), but its ability to fully represent the impact of eddies is still questioned (Eden and Greatbatch, 2008; Marshall et al., 2012; Zanna et al., 2017; Mak et al., 2017). The parameterization represents a net sink of energy as the potential energy extracted from the iso-neutral slopes is lost. Moreover, even coupled with Laplacian viscosities, it neglects the eddy effects on the momentum equation, the so-called eddy Reynolds stresses (Marshall et al., 2012; Zanna et al., 2017; Mak et al., 2018). These deficiencies in the representation of the eddies have consequences for the mean flows and the transport of tracers (Li et al., 2016; Tamarin et al., 2016; Zanna et al., 2017; Klocker, 2018). In our simulations, the mean equilibrium state in the eddy-parameterized simulations deviates from the eddy-resolving ones, with little sensitivity to the choice of the GM and isopycnal diffusion coefficients. And importantly, the global-warming-driven changes in circulation and the subsequent NPP decline greatly differ from the simulated changes in the eddy-resolving simulations (Fig. 7). This suggests that strong biases in global warming NPP projections could come from inadequate eddy parameterizations. On a positive note, most of the biases stem from differences in the response of the mean flow to warming, which suggests that a better representations of the eddy impact on the mean flow may correct for it. Ongoing activities in that direction open optimistic perspectives (Marshall et al., 2012; Zanna et al., 2017; Mak et al., 2018).

The difficulty of attributing recently observed changes in NPP to climate change and even of detecting long-term trends (Wernand et al., 2013; Beaulieu et al., 2013; Boyce et al., 2010; Henson et al., 2016) make ESMs essential for anticipating anthropogenic changes in NPP. In this regard, it is key to understand the origins and consequences of ESM biases and/or projection uncertainties. One such bias, identified here, is the high sensitivity of NPP projections to model resolution. It is questionable how the bias highlighted here in an idealized regional setting can be extrapolated to global ESMs. The oceanic regime closest to our configurations is found in the North Atlantic. Our model captures important characteristics of the North Atlantic NPP system, in particular the reduction in NPP in response to global warming (Kwiatkowski et al., 2020) as well as the leading role of the nutrient advection in refueling nutrients at the sub-surface in the productive subpolar region (Williams et al., 2011). With our simplified model, we were able to investigate the consequences of a decline in the MOC driven by changes in air–sea heat fluxes and associated with the reduction in deep water formation of a unique water mass. However, changes in the AMOC may also be driven by freshwater input (Bras et al., 2021), driven by changes in wind stress pattern (Polovina et al., 2008; Yang et al., 2020), or related to changes in adjacent regions and involving the formation of different water masses (Delworth and Zeng, 2008; Bronselaer et al., 2016; Lique and Thomas, 2018; Bras et al., 2021). Moreover, the link between the reduction in deep water formation at high latitudes in the North Atlantic and the slowing of the AMOC may be more tenuous than previously thought (Lozier, 2012). Investigation of these different aspects would require more realistic configurations and more complex climate change scenarios. Global ESMs and CMIP6 framework scenarios incorporate such additional features, as well as the temperature sensitivity of biological processes that has been revealed to also affect the global-warming-driven response of biology (Olonscheck et al., 2013; Lewandowska et al., 2014; Laufkötter et al., 2015). All these elements are likely to further influence the sensitivity of NPP projections to model resolution.

Constraining uncertainties in ESM NPP projections is crucial because these projections are essential to assess the future evolution of marine biomass and the potential implications for fish catch potential and food security (Lotze et al., 2019). Although being only one piece of a much broader and more complicated response of the ocean to climate change, our results suggest that the uncertainty in NPP decline intensity under global warming may have been underestimated in recent policy-relevant reports such as IPCC (2019) and IPBES (2019). A natural next step to better assist in political decision-making is now to conduct climate change simulation with ESMs better representing oceanic eddies. Two conceivable strategies to conduct such simulations have already been explored for ocean dynamics: simulations with ESMs at finer ocean resolutions (Haarsma et al., 2016; Gutjahr et al., 2019; Roberts et al., 2004; van Westen and Dijkstra, 2021) and simulations with ESMs using improved eddy parameterizations (Zanna et al., 2017; Mak et al., 2018). Our study calls for pursuing research in that direction. Work is in progress to examine the sensitivity of carbon fluxes to eddy resolution. Such advances would be significant for the representation of the global carbon cycle and consequently for climate projections.

Appendix A

Figure A1 Time series of (a) annual ocean temperature (^∘C) and (b) annual nitrate concentration (mmol N m⁻³) averaged between 0 and 700 m depth and (c) vertically integrated net primary production ($\mathrm{mmol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). All variables are averaged over the whole domain. Years −100 to 0 correspond to the spin-up period under seasonal preindustrial forcing. The preindustrial control and climate change simulations discussed in this paper start at year 0. Line thickness stands for the different resolutions, and the various shades of gray stand for the various eddy parameterizations. The small peaks in the time series of the $\mathrm{1}/\mathrm{27}$^∘ resolution spin-up around years −80 and −60 are due to a problem during the chaining of the spin-up jobs (which were set by periods of 10 years). These errors did not influence the equilibrium states, and we decided not to rerun the spin-up simulations to save computational resources.

Figure A2 The left column shows barotropic circulation (Sv), and the right column shows vertically integrated net primary production (NPP, $\mathrm{mmol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) in preindustrial control simulations. Panels (a)–(j) show the five 1^∘ resolution configurations, and (k)–(n) show the $\mathrm{1}/\mathrm{9}$ and $\mathrm{1}/\mathrm{27}$^∘ resolution configurations. All figures show the average over the last 5 years of the simulations (years 66 to 70).

Figure A3 Percentage of local change in vertically integrated net primary production (NPP) due to a change in new NPP (New NPP, first column) and regenerated NPP (Reg NPP, second column). Respective share of the NPP decrease from new NPP (black) and regenerated NPP (gray), averaged over the subpolar box (third column, $\mathrm{mmol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). Evolution of the f ratio averaged over the subpolar box in the preindustrial control and climate change simulations (fourth column, dashed and solid lines, respectively). Vertical integration up to 100 m depth where most of the NPP occurs. New NPP is the NPP supported by nitrate minus nitrification, and regenerated NPP is NPP supported by ammonium plus nitrification. The nitrification is handled in this way so that locally produced nitrate is counted as a source of regenerated production and not as a source of new production. The f ratio is computed as the ratio between new NPP and total NPP. Panels (a)–(t) show the five 1^∘ resolution configurations, and (u)–(ab) show the $\mathrm{1}/\mathrm{9}$ and $\mathrm{1}/\mathrm{27}$^∘ resolution configurations. The first two columns show the average over the last 5 years of the simulations (years 66 to 70).

Figure A4 Vertical distribution of nitrate across the domain (zonal mean, mmol N m⁻³) in the preindustrial control simulations (first column) and the climate change simulations (second column) and the difference between the two (third column). Panels (a)–(o) show the five 1^∘ resolution configurations, and (p)–(u) show the $\mathrm{1}/\mathrm{9}$ and $\mathrm{1}/\mathrm{27}$^∘ resolution configurations. All figures show the average over the last 5 years of the simulations (years 66 to 70).

Figure A5 Vertical profiles of dynamical nitrate supplies ($\mathrm{mmol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) in the model's subpolar box (dashed black lines in Fig. 3), cumulated between the surface and depth D in the preindustrial control simulations, against D. Total transport (gray lines here, solid blue lines in Fig. 4) is decomposed into three components: in orange, the meridional advective flux at the southern border; in blue, the vertical advective flux at depth D; and, in green, the isopycnal mixing. For each resolution, the total advective fluxes (first column) are broken into mean and eddy transport terms (the second and third columns, respectively). In the fine-resolution simulations ($\mathrm{1}/\mathrm{9}$, $\mathrm{1}/\mathrm{27}$^∘), they are derived from Reynolds decomposition (see the Methodology section). In the coarse-resolution simulations (CTL1/CC1), the mean transport represents the resolved transport by the coarse-resolution velocity field, and the eddy transport includes the two eddy parameterizations. Panels (a)–(o) show the five 1^∘ resolution configurations, and (p)–(u) show the $\mathrm{1}/\mathrm{9}$ and $\mathrm{1}/\mathrm{27}$^∘ resolution configurations. All figures show the average over the last 5 years of the simulations (years 66 to 70). Positive values indicate an inflow into the subpolar box.

Figure A6 Changes in the dynamical nitrate supplies ($\mathrm{mmol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) between the climate change and preindustrial control simulations. See Fig. A5 for details. The change is relative to the preindustrial control simulations. Positive values indicate an increase in nitrate input either because an inflow increases or because an outflow decreases.

Figure A7 Vertical profiles of changes in advective nitrate supplies (black lines, $\mathrm{mmol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) broken into the nitrate distribution change term (red lines), circulation change term (blue lines) and residual non-linear term (gray lines). Changes are cumulated between the surface and depth D in the model's subpolar box (dashed black lines in Fig. 3) and plotted against D. Change is the difference between climate change and preindustrial control simulations For details, see the Methodology section, Eqs. (3) to (6). The first column shows total advective fluxes which are the sum of the vertical advective fluxes at depth D (second column) and the meridional advective fluxes at the southern border (third column). Panels (a)–(o) show the five 1^∘ resolution configurations, and (p)–(u) show the $\mathrm{1}/\mathrm{9}$ and $\mathrm{1}/\mathrm{27}$^∘ resolution configurations. All figures show the average over the last 5 years of the simulations (years 66 to 70).

Figure A8 Meridional stream functions (Sv) in the preindustrial control simulation (first column) and climate change simulations (second column) and the difference between the two (third column). Meridional stream functions are computed in z coordinates. The vertical black lines range from 100 to 400 m depth and show the area through which we diagnose the water flow between the two gyres. Panels (a)–(o) show the five 1^∘ resolution configurations, and (p)–(u) show the $\mathrm{1}/\mathrm{9}$ and $\mathrm{1}/\mathrm{27}$^∘ resolution configurations. All figures show the average over the last 5 years of the simulations (years 66 to 70).

Figure A9 Maximum depth of the mixed layer (meters) reached along the seasonal cycle (monthly averages) in the preindustrial control simulations (first column) and differences in the climate change simulations (second column). The third column shows distributions of maximum mixed-layer depth north of the model's subpolar box (dashed black lines in Fig. 3) in the preindustrial control simulations (blue bars) and the climate change simulations (red bars). Panels (a)–(o) show the five 1^∘ resolution configurations, and (p)–(u) show the $\mathrm{1}/\mathrm{9}$ and $\mathrm{1}/\mathrm{27}$^∘ resolution configurations. All figures show the average over the last 5 years of the simulations (years 66 to 70).

Figure A10 Vertical distribution of the Brunt–Vaisala frequency across the domain (zonal mean, ${\mathrm{10}}^{-\mathrm{5}}{\mathrm{s}}^{-\mathrm{1}}$) in the preindustrial control simulations (first column) and the climate change simulations (second column) and the difference between the two (third column). Panels (a)–(o) show the five 1^∘ resolution configurations, and (p)–(u) show the $\mathrm{1}/\mathrm{9}$ and $\mathrm{1}/\mathrm{27}$^∘ resolution configurations. All figures show the average over the last 5 years of the simulations (years 66 to 70).

Figure A11 Decomposition of the local divergence of the nitrate advective fluxes ($\mathrm{mmol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) into a mean component and a fluctuating eddy component, for the $\mathrm{1}/\mathrm{27}$^∘ simulation. (a, d, g) Total local divergence, (b, e, h) mean component and (c, f, i) residual eddy component. (a–c) The preindustrial control simulations, (d–f) the climate change simulations and (g–i) the difference between the two. The total component is computed using 2 d averages; the mean component is computed using a spatio-temporal (1^∘ and 1 year) average. The eddy component is the residual between the two. The decomposition is integrated vertically between the surface and 400 m depth and averaged over the last 5 years of each simulation (years 66 to 70). Details of the decomposition are in the Methodology section. Positive values indicate an input of nitrate.