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 5 a mid-latitude double-gyre circulation, and perform global warming simulations under increasing horizontal resolution (from 1◦ to 1/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 resolution : 1◦ :::::::: resolution :::::::::: simulations, the simulated decline in primary production in the subpolar 10 gyre is halved at the finest eddy-resolving resolution (−12 % at 1/27◦ vs −26 % at 1◦ : ) : at the end of the 70 years long global warming simulations). This difference stems from the high sensitivity of the sub-surface nutrient transport to model resolution. Our ::::::: 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 on nutrient transport below the seasonal mixed-layer to better constrain the future evolution of marine biomass and fish catch potential. 15

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., 2020a) :::::::::::::::: (Bahl et al., 2020b). Here, we assess the difference that eddy 40 resolution makes to shaping the response of NPP to future global warming. To circumvent the computational constrains, 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 Kuroshoïo :::::::: 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 45 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 (1/9 • and 1/27 • ) model configurations ( Fig.   1 and 2), we show that the projected NPP decline under global warming is halved with the eddy-resolving configurations.

Methodology
To overcome the computational constraints associated with the challenge of running eddy-resolving simulations of climate 65 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.
The basin is initialised at rest with vertical profiles of temperature and salinity uniformly applied to the whole domain and homogeneous concentrations of the biogeochemical tracers. A 2000 years spin-up with a : an : eddy-parameterized coarse 105 resolution configuration (1 • , k gm =k redi ::: iso =1000 m 2 .s −1 ) is then performed under preindustrial forcings (wind, air temperature and fresh water 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 simulations are carried out from this initial coarse resolution spin-up state in 2 steps ( Fig. 1 and A1): 1. 100 years of spin-up under the same preindustrial control forcings, in order to reach an equilibrium for each configuration.

110
These seven additional spin-up (one for each configuration) are initialized with the annual mean coarse resolution spinup state. Notable is : It :: is :::::::: important :: to :::: note that the final dynamical and biogeochemical equilibrium are different for each configuration ( Fig. 1, 2, A1, A2, A4). The following section briefly describes these equilibrium and their differences.
In total, 14 simulations of 70 years were run: 7 configurations × 2 scenarios. In the following we will compare the differences between the eddy-resolving simulations (1/9 • and 1/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 120 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 choice. For each simulation, the impact of climate change is evaluated as the difference between the climate change 125 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 130 has achieved an equilibrium states that comes close to the eddy-resolving configurations (Fig. A1).

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 set up in ); Lévy and Martin (2013) :::::::::::::::::::::::: ; Lé are briefly discussed in the following. It can be noted that the differences between the coarse resolution equilibrium states on 135 one hand, and the fine resolution equilibrium states on the other had :::: hand, are much larger than the differences within the five (resp. two) coarse resolution (resp. fine resolution) spin-up ( 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 ( Fig. 2 and A2). On 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 150 at fine resolution than at 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 1/9 • and 1/27 • (Fig. 2) with a slightly larger sub-surface nitrate concentration at 1/9 • compared to 1/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.

Nitrate budgets
In our simulation, NPP is constrained by nitrate supplies from 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 ( Fig. 2 and 3). This subpolar box is shown by the black dashed horizontal lines in Fig. 3); it is bounded by the change of sign in the wind stress curl forcing 160 in the South and excludes the northern boundary. The local nitrate budget can be expressed as: is lateral ::::::: 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 parametrization Gent and McWilliams (1990a)) at coarse resolution.
Integrating this equation over the subpolar box and depth D leads to: dz Lateral mixingIsopycnal mixing ::::::::::: + Vertical mixing The bracket stands for the horizontal integral over the subpolar box. The first term on the left hand side is the integral of the isopycnal mixing, which is of second order (Fig. A5).

175
In the preindustrial control simulations, above ∼ 50 metres depth, nitrate supplies to the subpolar box are dominated by vertical mixing whereas deeper than ∼ 50 metres 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 1/9 • and 1/27 • resolution than at 1 • resolution (Fig. 4c, d).
The separation between mean and eddy transport in this budget was achieved following the Reynolds averaging method. u 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/spatial scale larger than the typical eddy scales. With this oper-185 ator, we decomposed the total advection term into a mean part including advection by mean velocities ( u cot N ds :::::::: u · N ds) and an eddy part ( u · N ds), that 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 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 CT L ) and nitrate from the climate change simulations (N CC ) or 195 vice versa (u CC and N CT L ). More precisely, we used the following decomposition of the change in nitrate advective fluxes, with: 3 Results

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 205 temperature increase at the end of 21st century in the SSP5-8.5 scenario (Tokarska et al., 2020a) ::::::::::::::::::: (Tokarska et al., 2020b)), 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  representing −26 ± 1 %, −14 % and −12 %, with horizontal resolution set at 1 • , 1/9 • and 1/27 • , respectively (Fig. 3). Thus 210 they differ markedly between each other, with the decline in NPP halved at the finest resolution compared with the decline projected at coarse resolution. Locally, the intensity of the decline in the coarse resolution model's versions exceeds −1.5 mmolN.m −2 .d −1 . The intensity and distribution of these local changes slightly varies among the five coarse resolution configurations. But with finer grid resolution, they are strongly reduced (Fig. 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 largescale 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) ::::::::::::::::::: Lévy and Martin (2013), with a spatio-temporal annual-mean 250 averaging over boxes of 1 • . In the coarse resolution simulations, eddy transport is parameterized by a bolus advection term (Gent and McWilliams, 1990a) :::::::::::::::::::::::::: (Gent and McWilliams, 1990b) and isopycnal mixing.
This analysis reveals the dominant role of the mean advection in fuelling the sub-surface nitrate reservoir at equilibrium in the two eddy-resolving models (preindustrial control simulations, blue dashed lines in Fig. 4c, d). Eddy advection supplies a minor addition of nutrient in the subpolar box (mostly at 1/27 • ), essentially through the eddy meridional advection (Fig. A5).

255
In the eddy-parameterized models, the eddy parametrization strongly opposes the resolved transport, removing nutrient from the surface of the subpolar box mainly through the vertical bolus advection term ( Fig. 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 ∼ 1/3 of the total transport decrease at 200 metres (Fig. 4). As for the preindustrial control simulation, it is dominated by changes in 260 the vertical bolus advection term (Fig. A6).

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 265 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 (grey 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 270 for most of the change. In fact at 400 metres 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 subsurface 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 model resolution increases.

275
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 sub-surface (between 50 and 400 m) (Fig. 6a, c). With warming, these upward and subsurface northward transport both decline (Fig. 6b, d). The mean transport are similar in the equilibrium states between the coarse and fine resolution simulations, but they respond differently to warming. At coarse resolution, the decline of the 280 upward transport is particularly strong, −0.5 Sv at ∼ 100 metres depth (Fig. 6b). This is in phase with the strong decline in nitrate vertical advection at similar depth discussed above (blue lines in Fig. A6). At finer resolution, the decline in the vertical flow is close to zero at ∼ 100 metres, but increases with depth reaching an amplitude similar to the coarse resolution simulations at ∼ 400 metres, ∼ −0.5 Sv. In the finer resolution simulations, the sub-surface northward flow decline is also weaker, ∼ −0.15 Sv against ∼ −0.3 Sv at coarse resolution below 100 metres depth (Fig. 6d).

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 310 at eddy resolution: −12 ::::::: −13 ± 1 % versus −26 ± 1 % at 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 r edi and k g m :: iso :::: and ::: . This result is due to the much weaker decline in nutrient transport at fine resolution, −29 ::::::: −27 ± 2 % versus −83 ± 5 % at coarse resolution, combined with a weaker decrease in vertical mixing, −18 ::::::: −18 ± 1 : % vs. −48 ± 4 % (Fig. 7).
Another possible reason for the high sensitivity to resolution of projected ocean circulation changes could be derived from the "eddy saturation" and "eddy compensation" phenomenon which are actively studied in the southern ocean  They stand respectively for the low sensitivity of stratification and overturning circulation to an increase in wind forcing in The GM parameterization has improved the representation of circulation in ocean models (Gent, 2011) ::::::::::::::::::::::::: (Gent and McWilliams, 1990b), but its ability to fully represent the impact of eddies is still questioned (Eden and Greatbatch, 2008a;Marshall et al., 2012a;Zanna et al., 201 The parameterization represents a net sink of energy as the potential energy extracted from the iso-neutral slopes is lost.

385
Such advances would be significant for the representation of the global carbon cycle and consequently for climate projections.