# Evaluation of an online grid-coarsening algorithm in a global eddy-admitting ocean biogeochemical model

Sarah Berthet, Roland Séférian, Clément Bricaud, Matthieu Chevallier, Aurore Voldoire, Christian Ethé

► To cite this version:

Sarah Berthet, Roland Séférian, Clément Bricaud, Matthieu Chevallier, Aurore Voldoire, et al. Evaluation of an online grid-coarsening algorithm in a global eddy-admitting ocean biogeochemical model. Journal of Advances in Modeling Earth Systems, 2019, 11 (6), pp.1759-1783. 10.1029/2019MS001644. hal-03528751

HAL Id: hal-03528751

https://hal.science/hal-03528751

Submitted on 17 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

10.1029/2019MS001644

**Special Section:**

The CNRM Climate and Earth
System Models for CMIP6

**Key Points:**

* The eddy-admitting (0.25°) model replicates observed physical-biogeochemical coupling, whereas the 1° horizontal resolution model does not

* The computation cost of the eddy-admitting model is divided by 2.7 when using a coarse-grained grid for marine biogeochemistry

* The coarse-grained solution inherits the key features of the 0.25° solution, including the physical-biogeochemical coupling

**Supporting Information:**

* Supporting Information S1

**Correspondence to:**

S. Berthet,
sarah.berthet@meteo.fr

**Citation:**

Berthet, S., Séférian, R., Bricaud, C., Chevallier, M., Voldoire, A., & Ethé, C. (2019). Evaluation of an online grid-coarsening algorithm in a global eddy-admitting ocean biogeochemical model. *Journal of Advances in Modeling Earth Systems, 11*, 1759–1783. https://doi.org/10.1029/2019MS001644

Received 1 FEB 2019

Accepted 6 MAY 2019

Accepted article online 13 MAY 2019

Published online 20 JUN 2019

<figure><img src="image_2.png" /></figure>

# Evaluation of an Online Grid-Coarsening Algorithm in a Global Eddy-Admitting Ocean Biogeochemical Model

Sarah Berthet¹ ID, Roland Séférian¹ ID, Clément Bricaud², Matthieu Chevallier¹, Aurore Voldoire¹ ID, and Christian Ethé³

¹Centre National de Recherches Météorologiques (CNRM), Université de Toulouse, Météo-France, CNRS, Toulouse, France, ²Mercator Ocean, Toulouse, France, ³Institut Pierre-Simon Laplace, Paris, France

**Abstract** In order to explore the effects of mesoscale eddies on marine biogeochemistry over climate timescales, global ocean biogeochemical general circulation models (OBGCMs) need at least to be run at a horizontal resolution of a 0.25°, the minimal resolution admitting eddies. However, their use is currently limited because of a prohibitive computational cost and storage requirements. To overcome this problem, an *online coarsening* algorithm is evaluated in the oceanic component (NEMO-GELATO-PISCES) of CNRM-ESM2-1. This algorithm allows to compute biogeochemical processes at a coarse resolution (0.75°) while inheriting most of the dynamical characteristics of the eddy-admitting OBGCM (0.25°). Through the coarse-graining process, the effective resolution of the ocean dynamics *seen* by the biogeochemical model is higher than that which would be obtained from an OBGCM run at 0.75°. In this context, we assess how much the increase from low (1°) to *coarse-grained* horizontal resolution impacts the ocean dynamics and the marine biogeochemistry over long-term climate simulations. The online coarsening reduces the computational cost by 60% with respect to that of the eddy-admitting OBGCM. In addition, it improves the representation of chlorophyll, nutrients, oxygen, and sea-air carbon fluxes over more than half of the open ocean area compared to the 1° OBGCM. Most importantly, the coarse-grained OBGCM captures the physical-biogeochemical coupling between sea-air carbon fluxes and sea surface height and between oxygen minimum zone boundaries and eddies, as produced by the eddy-admitting OBGCM. Such a cost-efficient coarsening algorithm offers a good trade-off to conduct process-based studies over centennial timescales at higher resolution.

## 1. Introduction

Over the recent decades, the oceans have stored ~93% of the anthropogenic heat (Fasullo & Trenberth, 2012) and 22% of the anthropogenic CO₂ emissions (Le Quéré et al., 2018), slowing down the climate change at multidecadal timescales. Recent observations have challenged this large-scale picture of the ocean by adding large fluctuations in ocean fluxes of carbon (Landschützer et al., 2015) and heat (e.g., Liu et al., 2016) at shorter timescales (season to decades). There are now robust evidences that eddies or mesoscale ocean structures also influence a number of biogeochemical processes across temporal and spatial scales (e.g., Bettencourt et al., 2015; Dufour et al., 2015; Gruber et al., 2011; Harrison et al., 2018; Lacour et al., 2017; Lévy & Martin, 2013; Mazloff et al., 2018; Munday et al., 2014; Oschlies, 2002; Resplandy et al., 2009; Resplandy et al., 2013; Sweeney et al., 2003; Villar et al., 2015).

For instance, observational process studies contribute to explain how eddy-induced local vertical motion may impact marine biogeochemistry and carbon cycle. In the Tasman Sea, the presence of a large number of eddies reinforces vertical mixing, which decreases the effective stratification and limits the winter phytoplankton and zooplankton populations (Tilburg et al., 2002). On the other side, Lacour et al. (2017) present observational evidences of unexpected winter phytoplankton blooms in the North Atlantic subpolar gyre. These blooms are triggered by intermittent restratification of the mixed layer in response to the eddy-driven transport of lighter water over denser layers. A similar phenomenon exists in the Southern Ocean, in the Agulhas region, where vertical mixing processes have been shown to alter the plankton transport (Villar et al., 2015). This strong vertical mixing has been suggested to shape the biogeochemical signatures in the Agulhas rings as the rings and associated plankton transit westward, emphasizing the role of mesoscale eddies on marine biogeochemistry.
Despite its importance for both ocean and biogeochemical dynamics, the influence of mesoscale activity on ocean biogeochemical cycles has been identified as a missing process in the current generation of Earth system models as assessed in IPCC AR5 (Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Ciais et al., 2013). For example, Terhaar et al. (2018) recently show that the Arctic Ocean storage of anthropogenic carbon differs when moving from coarse to eddy-admitting resolution in a global ocean circulation-biogeochemistry model. Indeed, theory of flow dynamics suggests that a few kilometers of horizontal resolution is needed to enable global ocean biogeochemical general circulation models (OBGCMs) to resolve eddies and their interactions with submesoscale dynamics. However, most of the current state-of-the-art OBGCMs still resolve both ocean and biogeochemical dynamics at a low horizontal resolution of hundreds of kilometers. Such a low resolution leads to a smoothed ocean dynamics where even huge eddies are poorly represented by large-scale and weak vertical motions. The apparent convergence in horizontal resolution toward 1 to 2° emerging from the literature across a variety of OBGCMs is rather a trade-off between model complexity and high-performance computing capacities than a choice supported by scientific considerations.

Recent theoretical advances in mesoscale dynamics and their modeled interactions with ocean large-scale circulation highlight the importance to consider the effective resolution hidden behind the computational one (Skamarock, 2004). By essence, ocean models produce outputs at a resolution lower than the grid resolution due to a smoothing by diffusive and dissipative processes. As introduced in Lévy et al. (2012), the effective resolution characterizes the smallest size of the structures captured by the model outside the dissipative range and determines at which resolution the large-scale influence of the turbulent dynamics saturates. Because the effective resolution of the ocean dynamics remains always lower than its computational grid resolution, it is more cost-efficient to directly compute the marine biogeochemistry on a grid at lower resolution, close to the effective resolution of the ocean dynamics. In this perspective, an online coarsening algorithm for the NEMO ocean model (Madec, 2008) has been developed based on Aumont et al. (1998) in order to compute biogeochemical processes from a coarse-grained ocean dynamics. With this algorithm, eddy-admitting (0.25°) dynamical fields are coarse-grained down to a lower resolution (0.75°), thus allowing to preserve an important part of their initial fine-scale structures. The resulting coarse-grained dynamics features finer scales than fields that would be obtained from a model run directly on the 0.75° grid. As assessed in section 3, the coarse-grained dynamics preserves most of the characteristics of the original finer-grid dynamics. This algorithm allows to bypass the computational cost bottleneck encountered in centennial-long climate model simulations.

The goal of the present study is twofold. Because the coarse-grained ocean dynamics is derived from an eddy-admitting dynamics, we first aim at evaluating how close the coarse-grained biogeochemical solution is to that produced by a 0.25° resolution. Then, to what extent is it closer to the eddy admitting than the standard state-of-the-art (1°) solution? For that purpose, three simulations have been produced with the same NEMO-GELATO-PISCES OBGCM under three different horizontal resolutions ranging from 1° to 0.25°. Whatever the simulation, our OBGCM has been forced by the same suite of atmospheric forcing and boundary conditions. Although a 0.25° remains clearly not sufficient to get an explicit representation of ocean mesoscale features, both horizontal resolutions of 1° and 0.25° have been selected as they are representative of those currently used in the simulations of the sixth phase of the Coupled Model Intercomparison Project (CMIP6; Eyring et al., 2016). Throughout this paper, the emphasis is given to two major biogeochemical features: the Southern Ocean carbon sink and the tropical oxygen minimum zones (OMZs).

The present manuscript is structured as follows. Section 2 provides details on the ocean and biogeochemical compartments of our NEMO-GELATO-PISCES OBGCM and on numerical setups of the experiments, observations, and analyses used throughout this study. Then two sections evaluate the general performance of our OBGCM against modern observations across the three simulations. Section 3 focuses on the dynamical mean state and energy content of that numerical solutions, and section 4 on their biogeochemical mean states. Section 5 shows how the increase in horizontal resolution changes the representation of the physical-biogeochemical coupling existing over two major key biogeochemical features (in the Southern Ocean and in the Pacific OMZ) and investigates how far biogeochemistry computed from the online coarsening algorithm allows to reproduce this coupling. Finally, section 6 summarizes the added value of the coarse-grained solution and the key findings of this work.
## 2. Experimental Design and Methods

### 2.1. Models

#### 2.1.1. Ocean-Sea Ice Component: NEMOv3.6-GELATOv6

In the following we present results from three configurations named LOW, HIGH, and CRS (see section 2.3.1), which are respectively at 1°, 0.25°, and at a coarse-grained horizontal resolution. LOW corresponds to the configuration of the ocean component of the Earth System Model CNRM-ESM2-1 (including sea ice and marine biogeochemistry), which has been used for CMIP6 and the Ocean Model Intercomparison Project (OMIP). The ocean component uses NEMOv3.6 (Nucleus for European Modelling of the Ocean; Madec, 2008), runs on ORCA tripolar grids (Madec & Imbard, 1996), and has 75 vertical levels. Similar to all three configurations is the lateral diffusion scheme for tracers, which is an isoneutral Laplacian with a constant mixing coefficient (see Table 1). Please note that the effective resolution remains sensitive to the choice of the lateral mixing coefficient; a paper in preparation (Bricaud et al.) will address the sensitivity of the online coarsening algorithm to the physical parameters. Whatever the configuration, tracer advection is treated with a total variance dissipation scheme (Lévy et al., 2001; Zalesak, 1979) and with an additional term coming from a mixed layer eddy-induced velocity parameterization (Fox-Kemper et al., 2011) allowing to enhance deep ocean ventilation. Configuration LOW differs from HIGH and CRS in that it adds an eddy-induced velocity term (Gent & McWilliams, 1990) to the latter ocean velocity. The lateral diffusion of momentum is either a horizontal Laplacian (LOW) or a horizontal bilaplacian (HIGH and CRS) with an eddy viscosity coefficient depending on the horizontal resolution of the experiment (see Table 1). The vertical eddy viscosity and diffusivity coefficients are computed from an improved version of the turbulent kinetic energy (TKE) closure scheme (Blanke & Delecluse, 1993) with a Langmuir cell (Axell, 2002). This improved TKE scheme allows a fraction of surface wind energy to penetrate below the base of the mixed layer, improving the coupling between surface wind and mixed-layer depth. Double-diffusion mixing is enabled, as it contributes to diapycnal mixing in extensive regions of the ocean. The parameterization of tidal mixing follows the general formulation for the vertical eddy diffusivity proposed by St. Laurent et al. (2002). We use the multicategory sea ice model GELATOv6 (Global Experimental Leads and ice for ATmosphere and Ocean; Salas y Mélia, 2002). This model is an elastic-viscous-plastic rheology-5 thickness categories sea ice model. It resolves the sea ice rheology using an elastic-viscous-plastic scheme as proposed by Bouillon et al. (2009), which operates on the same horizontal grid as NEMOv3.6. GELATOv6 describes the sea ice pack using five thickness categories: less than 0.30, 0.3-0.7, 0.7-1.2, 1.2-2, and over 2 m thick. The snow and ice parts of every ice category are respectively split vertically into one and nine layers. The radiative transfer scheme across the snow-sea-ice continuum is resolved using the formulation of Grenfell and Maykut (1977).

#### 2.1.2. Marine Biogeochemistry: PISCESv2-gas

In this work marine biogeochemistry is simulated using the Pelagic Interaction Scheme for Carbon and Ecosystem Studies version 2-gas (PISCESv2-gas; Aumont et al., 2015), which includes two additional trace gases: the dimethylsulfide (DMS) and the nitrous oxide ($N_2O$). As this version derives from the PISCESv2 version, the reader is referred to Aumont et al. (2015) for further details on the general description and evaluation of PISCESv2 and to Masotti et al. (2016) and Martinez-Rey et al. (2015) for the DMS and $N_2O$ modules, respectively. In the following, we only describe its main characteristics relevant for our study.

PISCESv2-gas resolves the evolution of 26 compartments, which enables the resolution of oxygen, carbon, sulfur, and nutrients on the NEMO grid. PISCESv2-gas takes five nutrients into account (from macronutrients nitrate, ammonium, phosphate, and silicate to micronutrient iron) whose external availability limits the growth of two phytoplankton classes (nanophytoplankton and diatoms). Diatoms differ from nanophytoplankton because they need silicon and more iron (Sunda & Huntsman, 1997) and because they have higher half-saturation constants due to their larger mean size. Phytoplankton growth is also limited by light. Redfield ratios between carbon, nitrate, and phosphate are assumed identical for all plankton classes (122:16:1, respectively; Takahashi et al., 1985), while concentrations of oxygen, iron, silicon, and chlorophyll are simulated prognostically, and their ratios to carbon (C) vary dynamically. Chlorophyll (Chl) concentrations for the two phytoplankton types are prognostically computed using the photoadaptative model of Geider et al. (1998), with Chl/C ratio varying as a function of light and nutrient limitation.

The coupling between PISCESv2-gas and NEMO is handled as follows. As mentioned in section 2.1.1, all bio-geochemical tracers of PISCESv2-gas are advected by the velocity fields using a total variance dissipation scheme (Lévy et al., 2001; Zalesak, 1979). The velocity fields used in this scheme consider the eddy-
<table><caption>**Table 1**
Selected Physical Coefficients and Parameters for the Various Configurations as Used in This Study</caption><thead><tr><th rowspan="2">Configuration</th><th rowspan="2">Acronym</th><th colspan="3">Horizontal resolution</th><th colspan="2">Lateral viscosity</th><th colspan="2">Eddy param</th></tr><tr><th>Ocean dynamics</th><th>Marine biogeochemistry</th><th>Coupling type</th><th>Tracers (m<sup>2</sup>/s)</th><th>Momentum: type, value, and units</th><th>GM</th><th>FK</th></tr></thead><tbody><tr><td>Low resolution</td><td>LOW</td><td>1°</td><td>1°</td><td>classic</td><td>1,000</td><td>Horizontal Laplacian: 20,000 m<sup>2</sup>/s</td><td>yes</td><td>yes</td></tr><tr><td>High resolution</td><td>HIGH</td><td>0.25°</td><td>0.25°</td><td>classic</td><td>300</td><td>Horizontal bilaplacian: -1.5e11 m<sup>4</sup>/s</td><td>no</td><td>yes</td></tr><tr><td>High resolution</td><td>CRS</td><td>0.25°</td><td><strong>0.75°</strong></td><td>coarse-grained</td><td>300</td><td>Horizontal bilaplacian: -1.5e11 m<sup>4</sup>/s</td><td>no</td><td>yes</td></tr></tbody></table>

Note. The HIGH simulation at 0.25° of horizontal resolution has a nominal resolution of 27 km at the equator and reaches 12 km in the Arctic Ocean. Eddy param refers to the addition of (left) GM: an eddy induced velocity parameterized following Gent and McWilliams (1990) and/or (right) FK: a mixed layer eddy induced velocity (Fox-Kemper et al., 2011). Coupling type refers to the use (coarse-grained)—or not (classic)—of the online coarsening facilities described in section 2.2.

induced velocities, which are either parameterized or resolved depending on the horizontal grid resolution of the experiment (Table 1). These tracers are impacted by diffusion, which follows an isoneutral Laplacian operator. The vertical transport of biogeochemical tracers except the sink of particles is handled by the vertical motions generated by the mixing scheme of NEMO and the divergence of the water flow. Indeed, particles sink with their own speed, which is modulated by the depth of the mixed-layer depth as described in Aumont et al. (2015). Finally, PISCESv2-gas uses the solar radiation from NEMO to represent the photosynthesis of phytoplankton and other photo-dependent processes. This solar radiation is averaged over a day because PISCESv2-gas does not represent the diurnal cycle of the marine biota.

Note that in a classic coupled configuration with NEMO and PISCES sharing the same horizontal resolution, the use of PISCES multiplies by more than 3 the computational cost compared to NEMO stand-alone simulations and by 5 the data storage needed (as NEMO stand alone involves only five tracers).

## 2.2. Increasing the Effective Resolution of Biogeochemistry: The Online Coarsening Algorithm

In NEMO, the algorithm allowing to compute the marine biogeochemistry on a grid having a horizontal resolution close to the effective resolution of the ocean dynamics is called the online coarsening. It has been adapted from the stand-alone degradation-integration technique described in Aumont et al. (1998).

In the following, we briefly describe the main characteristics and principles of this algorithm. The online coarsening consists in resolving the biogeochemical processes of PISCES on a grid whose horizontal resolution is coarser than the one used to resolve the ocean dynamics of NEMO (Figure 1).

Here the coarse-grained grid of PISCES has a 0.75° of horizontal resolution and hence gathers $3\vec{i}*3\vec{j}$ horizontal grid points of the NEMO mother grid at 0.25° (with $\vec{i}, \vec{j}$ representing the horizontal vector components). At each model time step, the coarse-grained dynamical fields (at 0.75°) are built by gathering a 3 by 3 quad of cells from the ocean dynamical fields resolved at 0.25° (Figure 2). The coarse-grained ocean dynamical fields are then used to compute the biogeochemical tracers' advection, mixing, and diffusion onto a grid at 0.75° of horizontal resolution.

Thanks to the online coarsening algorithm, the coarse-grained grid at 0.75° is expected to inherit a nonnegligible part of the energy spectra represented on the mother grid (Lévy et al., 2012). Consequently, the effective resolution (Skamarock, 2004) of the coarse-grained grid remains higher than the one it would normally have had if dynamics had been computed at the coarse-grained resolution (i.e., 0.75°). Indeed, by including the resolved modes of the ocean dynamics at high resolution, the coarse-grained solution is expected to inherit a better representation of the turbulent dynamical field than the one it would have had if computed directly on the coarser grid (see section 3.3).

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Illustration of NEMO-PISCES grid's interaction through the coarse-graining process.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Coarse-grained (black) and 0.25° parent (dark cyan) Arakawa-C mesh grid structures.</figcaption></figure>

Because our NEMO configuration has a nonlinear free surface moving up/down with convergence/divergence in the first ocean level, the main concern of the online coarsening algorithm is to deal with mass conservation and maintain nondivergence. In practice, five main points have to be considered. *First*, dynamical variables are constructed differently according to their properties (extensive or intensive). Temperature and salinity of the coarse-grained grid result from a weighted volume mean over the nine cells of the parent grid (see Appendix A2.1), while surface and interface fluxes of the coarse-grained grid are computed as a weighted area sum of either the horizontal parent grid cells (see Appendix A2.3) or along one dimension considered as the parent cell border (see Appendix A2.2), respectively. Due to these averaging processes, the coarsening algorithm will smooth local structures (e.g., fronts and eddies) of the dynamical fields compared to those of the 0.25° solution. Consequently, the coarse-grained biogeochemistry is computed from dynamical fields that are close to—but not similar to—those used to compute the biogeochemical quantities in the 0.25° solution. *Second*, we compute two types of vertical scale factors (i.e., vertical cell sizes) per coarse-grained mesh (Figure 3). One of them is used for the divergence operations (see Appendix A3.4) and allows to conserve in the coarse-grained mesh the total volume of water contained within the clustered parent cells. The other preserves the maximum depth of the clustered parent cells as required when applying vertical gradient operators (see Appendix A3.5). *Third*, as no exact theory exists to formalize the coarse-grained vertical diffusion ($K_z$), it is approached by a weighted mean over the parent dynamical grid (0.25° grid). *Last*, to conserve one of the fundamental properties of the ocean (the non-divergence of the velocities), the vertical velocity is not directly coarse-grained from the parent grid mesh as all other dynamical variables, but it is computed from the coarse-grained horizontal divergence, which is in turn derived from the horizontal coarse-grained velocities. Obviously, the two latter points will impact the

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Coarse-grained (black) and 0.25° parent (dark cyan) vertical scale factors construction required either by Appendix A3.4 (divergence) or Appendix A3.5 (gradient) operator.</figcaption></figure>
vertical dynamics and mixing processes and may explain in part why the biogeochemical fields of the coarse-grained solution are not expected to be similar to the 0.25° solution, although they derived from the resolved part of the same ocean dynamics. *Moreover*, a far less obvious difference between our coarse-grained and 0.25° solutions comes from the transport of biogeochemical fields: in the coarsened experiment biogeochemical fields are computed *and* advected over the coarse-grained grid. When numerical schemes are run on the coarse-grained grid, they give solutions for biogeochemistry with truncations that differ from those computed over the 0.25° grid. The truncation and averaging effects also affect the coarse-grained isopycnal slopes, and the directions of diffusion will slightly differ between the coarse-grained and 0.25° solutions.

In summary, it has to be noted that by adding the effects of (i) averaging processes on the coarse-grained dynamical fields (the first and second points above), (ii) coarse-grained vertical dynamics (third and fourth points), and (iii) truncation on biogeochemistry computation and transport (fifth point), the coarse-grained biogeochemical solution is expected to differ from the one computed over the 0.25° grid. For this reason, in the following the coarse-grained solution is evaluated as an interesting cost-efficient solution, which remains clearly distinct from those resulting from the 1° and 0.25° configurations in terms of biogeochemistry. The key aspect of the present evaluation is precisely to determine how close to the 0.25° is the coarse-grained solution and what it provides compared to a 1° solution.

Further details relative to the impact of using initial biogeochemical conditions at different resolutions are discussed in the supporting information. Although biogeochemical sources and sinks have initially been computed over grids at different horizontal resolutions, we show that the differences obtained between the biogeochemical solutions as produced by our coarse-grained and 0.25° solutions occur independently of the biogeochemical inputs and are inherent to the coarsening algorithm (i.e., resulting from the above points).

## 2.3. Experiments and Methods

### 2.3.1. Simulations

All the simulations presented here have been produced with the same setup. Each model configuration is run for 186 years using three cycles of the COREII atmospheric forcing (Griffies et al., 2012), which provides atmospheric drivers of the ocean from 1948 to 2009 (62 years). The atmospheric CO₂ as seen by PISCESv2-gas is kept constant in time and set to its preindustrial level: 284.32 ppm. The first 124 years (i.e., the first two CORE-II cycles) are considered as the spin-up period for the biogeochemical fields. This setup corresponds to the experimental configuration as described in Orr et al. (2017) to simulate the natural variation of the ocean carbon cycle.

In agreement with the OMIP protocol (Griffies et al., 2016; Orr et al., 2017), the ocean component of NEMOv3.6 was initialized at rest with observed climatologies from World Ocean Atlas 2013 observations for salinity (Zweng et al., 2013) and temperature (Locarnini et al., 2013). The sea ice model GELATOv6 was also initialized at rest, but its initial conditions were derived from an analytical solution based on ocean salinity and temperature as detailed in Salas-y-Mélia (2002).

PISCESv2-gas uses different sources of initial conditions. Dissolved inorganic carbon (DIC) and alkalinity (Alk) are initialized using the modern climatology of the Global Ocean Data Analysis Project version 2 (GLODAPv2; Key et al., 2015; Lauvset et al., 2016). Phosphate (PO₄), silicate (Si), and nitrate (NO₃) biogeochemical tracers are initialized from the 2013 World Ocean Atlas (WOA2013; Garcia et al., 2014b). Oxygen (O₂) initial conditions are also obtained from WOA2013 (Garcia et al., 2014a). Dissolved iron is not initialized from observations but from a previous model simulation (Aumont et al., 2015). The other biogeochemical tracers including DMS and N₂O are initialized from global average estimates. For these other biogeochemical tracers, the initialization is less critical (e.g., phytoplankton biomass is restricted to the top 200 m) and equilibrates rapidly.

In addition to these initial conditions, PISCESv2-gas uses several boundary conditions that represent the supply of nutrients from five different sources: atmospheric deposition, rivers, sediment mobilization, sea ice, and hydrothermal vents. This version of PISCESv2-gas does not include iron supply from sea ice melting and hydrothermal vents as described in Aumont et al. (2015). It is important to note that these external boundary conditions are provided to PISCESv2-gas at its grid resolution, meaning that boundary conditions at 0.25° are applied to the eddy-admitting OBGCM, whereas they are at 0.75° for the coarse-grained OBGCM.
Thus, the treatment of iron sediment mobilization differs slightly because of the various model resolutions. Indeed, in order to take subgrid-scale bathymetric variations into account in the iron source as resolved by the model at a given grid resolution, the iron source has been estimated from the high-resolution ETOPO5 data set and gathered up to the model grid resolution (Aumont et al., 2015).

Two types of coupled ocean biogeochemical configurations are discussed in this study: classic or coarse-grained (Table 1). In the classic case, both components (ocean and biogeochemistry) share a unique grid, that is to say either 1° or 0.25° of horizontal resolution at the equator according to the considered experiment. In the coarse-grained case, the algorithm described in section 2.2 is used to compute the spatial distribution of biogeochemical tracers. Consequently, the ocean dynamics of the coarse-grained simulation has 0.25° of horizontal resolution, while the associated biogeochemical fields have been computed on a 0.75° grid. All coupled experiments are listed in Table 1: a set of two classic experiments (LOW and HIGH) is used in the following to evaluate the model biases with increasing horizontal resolution. They are compared to the coarse-grained coupled experiment (namely, CRS) in order to quantify the added value of the coarsening algorithm.

### 2.3.2. Analyses

#### a Model-Data Intercomparison

A first step is dedicated to assess model performance against observations, highlighting the impact of increasing the horizontal resolution on the dynamical and biogeochemical mean states. CRS is compared to LOW and HIGH simulations to quantify to what extent the coarsening algorithm approaches the mean state of these solutions. This assessment relies on a set of standard metrics as recommended by Stow et al. (2009). They consist in the spatial distribution of annual means, standard deviations, and associated biases, projected on a regular 1° grid using a bilinear interpolation. Observations used for these analyses are

1. the monthly SSALTO/DUACS L4 product of sea surface height (SSH) computed by the French space agency (Centre National d'Études Spatiales [CNES]) and Collecte Localisation Satellites (CLS) company. The data set covers the 1993–2011 period and has been remapped down from 0.25° to 1° of horizontal resolution to enable comparison with model outputs.

2. the global monthly climatology of the World Ocean Atlas 2013 for seawater temperature (Locarnini et al., 2013) and salinity (Zweng et al., 2013) at 1° of horizontal resolution with 102 vertical levels between the surface and 5,500-m depth.

3. the monthly climatology of mixed layer depth of Sallée et al. (2010) at a 1° of horizontal resolution.

In order to assess the modeled ocean dynamics, biases are evaluated as a difference between each model solution remapped at 1° of horizontal resolution and these reference datasets at the same horizontal resolution. The global root-mean-square error (rmse) is computed as the square root of the horizontal sum of the field squared bias; the latter squared bias has been weighted by the local area mesh value and normalized by the global area. We exclude the coastal grid points with a depth <2,000 m from the computation of rmse.

In addition, vertical profiles and sections have been used to assess the vertical distribution of hydrographic and biogeochemical tracers. Modeled marine carbon cycle and biogeochemistry is evaluated by the use of

1. the annual climatology of dissolved inorganic carbon and total alkalinity from GLODAPv2 (Key et al., 2015; Lauvset et al., 2016).

2. the observed monthly climatology of nitrate, phosphate, silicate, and $O_2$ from the WOA2013 (Garcia et al., 2014a, 2014b) having comparable resolutions than the World Ocean Atlas dynamics (see above).

3. the climatology of natural (i.e., the component of the flux that already existed in preindustrial times) sea-air carbon fluxes at the surface is estimated from an inverse method using basis functions generated by a suite of 10 different ocean general circulation models and on the basis of ocean interior observations of DIC and other tracers (Mikaloff-Fletcher et al., 2007). The inverse product has 1° of horizontal resolution.

4. European Space Agency (ESA) Ocean Colour Climate Change Initiative (OC-CCI) product for surface chlorophyll (Valente et al., 2016).

Taylor diagrams (Taylor, 2001) are used to summarize the global annual mean distributions of nutrients compared to observations at four depth levels (surface, 150, 500, and 1,000 m), based on normalized standard deviations and spatial correlations between modeled and observed $PO_4$, $NO_3$, Si, and $O_2$.
If not specified otherwise, the period of interest covers the last 10 years from 2000 to 2009 and the analyses are based on monthly outputs.

b Assessment of the Intermodel Distance

In order to quantify how far LOW, CRS, and HIGH stand between each other in terms of simulated fields, we performed a suite of statistical tests. A first set of statistical tests has been employed to estimate how different are the simulated ocean dynamics between the eddy-admitting and coarse-grained solutions. The 10-year mean vertical profiles of temperature and salinity are compared between each other using a confidence inter- val at 95% computed from a *t*-statistics with (2*120-2) degrees of freedom estimated from HIGH. *t*- and *F*-statistics have also been used to evaluate how the annual maximum of mixed-layer depth (MLD) and the annual standard deviation of SSH may significantly differ (at 95%) between the coarse-grained and eddy-admitting solutions, respectively. A second set of tests allows to determine how CRS improves the biogeochemical solution compared to the 1° solution. At each grid point, the intermodel difference between the mean fields of LOW and HIGH has been compared to the difference obtained between CRS and HIGH. The result of this difference is given by superimposing a stippling to the biogeochemical mean maps where the difference CRS — HIGH is strictly lower than the difference LOW — HIGH. This highlights the potential model improvements when moving from LOW to CRS.

c Process-Based Analyses

In order to determine whether CRS conserves the physical-biogeochemical coupling properties emerging in the eddy-admitting HIGH simulation, we developed a set of process-based diagnostics. These latter focus on two major biogeochemical features: the ocean carbon fluxes and the Pacific OMZ. First, the signature of mesoscale activities on ocean dynamics across the various model configurations is deduced from the mean power spectra determined from daily outputs over a single year (2000). Then, two process-based analyses have been performed to evaluate if the coupling between eddies and biogeochemistry exists in our experi- ments as it is suggested by previous studies. Indeed, in order to track the potential impact of eddies on ocean carbon fluxes, we compute the square correlations of deseasonalized monthly anomalies of sea-air carbon fluxes and SSH. Finally, the role of mesoscale structures on oxygen distribution is evaluated from the mean daily O₂ horizontal gradient over year 2000.

3. Ocean Dynamics

Three levels of ocean dynamics evaluation are presented in this section depending on the considered dyna-
mical field and the degree of similarity what we are expecting of it between CRS and HIGH.

**3.1. Temperature and Salinity Profiles**

In order to track the impact of deep convection processes on ocean hydrographical and water masses proper-
ties, the vertical profiles of temperature and salinity are analyzed within the Labrador and the Weddell
Seas (Figure 4).

As detailed in section 2.2 and although extracted from the 0.25° solution, the coarse-grained ocean tempera-
ture and salinity are expected to differ slightly from HIGH due to averaging processes effects of weighted
volume means. In frontal regions, smoothing of local structures will affect regional monthly means. The
more convective the region is, the more pronounced the discrepancies in temperature or salinity between
CRS and HIGH will be. For example, vertical mean profiles of temperature are similar at whatever depth
in HIGH and CRS when considering the mean over the tropical band (5°S-5°N) or the tropical instability
waves (150°W-100°W, 5°S-5°N), the Agulhas region (5°W-40°E, 42°S-27°S), or the Austral band between
70°S and 45°S, but mean vertical profiles over the Irminger (40°W-30°W, 57°N-67°N) or the Gulf Stream
(70°W-40°W, 35°N-50°N) regions may differ in the order of 0.1-0.8 °C (respectively) depending on the size
of the considered region (see supporting information Figures S2 to S7).

In the Labrador Sea, one of the major deep convection sites of the North Atlantic sector, the increase in
horizontal resolution from 1° to 0.25°, tends to reduce differences between model and observations
(Figures 4a and 4b). Indeed, the minimum temperature around 100-m depth (Figure 4a) and the modeled
vertical structure of the thermocline are closer to the climatological observed temperature profile in
HIGH than in LOW. The fact that the CRS temperature mean vertical profile remains within or close
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Mean modeled and observed annual vertical profiles of temperature (°C) and salinity (psu) within (a, b) the Labrador Sea [55°W-45°W, 50°N-60°N] and (c, d) the Weddell Sea [55°W-45°W, 70°S-60°S] over the 2000-2009 period. Gray dashed lines correspond to the confidence interval at 95% of HIGH mean profiles.</figcaption></figure>

to the lower boundary of the confidence interval at 95% of the HIGH mean profile (dashed grey lines) allows to consider that the coarse-grained solution behaves similarly as HIGH. A similar behavior between the various model solutions and observations is also found in the Weddell Sea (Figure 4c), except that all solutions exhibit a warm bias in the thermocline. Deep convection occurs only occasionally in the Weddell Sea, and this may explain the better agreement between the modeled vertical mean profiles from CRS and HIGH.

Noticeable model-data differences in salinity profiles or halocline exist in the Labrador and Weddell Seas. The three configurations broadly capture the shape of the halocline. However, in the Labrador Sea, the halocline is too steep in LOW with respect to the observations (Figure 4b). Contrasting with this configuration, HIGH and CRS capture the shape of the halocline but exhibit a bias of about 0.3 psu higher than the observations below 200-m depth (and higher than 0.5 psu at the surface). In the Weddell Sea (Figure 4d), the three
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Observed (a) and modeled annual maxima of MLD climatological mean (m) over 2000–2009 for (c) LOW, (e) CRS, and (g) HIGH. Associated biases are shown in the right column (d, f, and h, respectively) with their indicative rmse values on the left. Biases of zonal means are presented in (b). Grid points where means of CRS and HIGH differ significantly (95%) have been masked in black on subfigure (f).</figcaption></figure>

solutions display broadly the same mean halocline and exhibit too haline water masses below 50 m with
respect to the observations.

**3.2. Mixing and Large-Scale Circulation**

As MLD depends on vertical dynamics and mixing, coarse-grained MLD is expected to differ from HIGH MLD by adding up the averaging effects of coarse-grained variable building and the use of a coarse-grained vertical dynamics (i.e., based on a vertical velocity computed from the coarse-grained horizontal divergence).

This section focuses on the large-scale structures of maximum of MLD. The maximum of MLD is commonly
used to diagnose how a physical model is able to drive properly ocean biogeochemistry (e.g., Doney, 2004;
Schneider et al., 2008), as it controls nutrient supply within the photic zone (e.g., Behrenfeld, 2010;
Sarmiento & Gruber, 2006).

Figure 5 shows that our three configurations broadly capture the geographical structure of MLD annual maxima with the occurrence of deep convection in the North Atlantic and in the Southern Ocean reaching more than 450 m in average over the year. However, the zonal mean hides large compensating errors present in the 1° resolution (Figures 5b and 5d). When comparing the MLD annual maxima in the Southern Ocean,
<caption>**Table 2**</caption>

Estimation From Observations (e.g., Lherminier et al., 2010; Talley et al., 2003) and Model Solutions of the Subsurface Southward (Between 300 and 2,000 m) and Deep Northward Flow of Water Masses (Below 2,000 m) at 26°N

<table><thead><tr><th>Water masses</th><th>Observations</th><th>LOW</th><th>HIGH</th><th>CRS</th></tr></thead><tbody><tr><td>AMOC</td><td>18 ± 5</td><td>19.1 ± 2.8</td><td>19.1 ± 3.1</td><td>19.1 ± 3.1</td></tr><tr><td>AABW</td><td>0.9 ± 2.1</td><td>1 ± 0.8</td><td>0.9 ± 3.2</td><td>0.9 ± 3.2</td></tr></tbody></table>

Note. These latter are evidence of the Atlantic Meridional Overturning Circulation (AMOC) and the Antarctic bottom waters export (AABW), respectively. All the transports are given in Sverdrup (Sv = $10^6 m^3/s$).

the large dipole of bias present in LOW is weakened in HIGH and in the coarse-grained solution (Figure 5f for CRS and Figure 5h for HIGH). LOW displays too smooth structures (Figure 5c) in the Southern Ocean south of 45°S and in the Kuroshio region compared to observations. In these two regions HIGH (Figure 5g) and CRS (Figure 5e) have more finer scales, which tend to locally reduce biases (Figures 5f and 5h).

A *t*-statistics has been used to evaluate how the means of MLD maxima from CRS and HIGH may be significantly different. In Figure 5f grid points where this difference is significant at 95% have been masked (in black, as the landmask): thus, the *t*-statistics allows to quantify that CRS is similar to HIGH in terms of mean of MLD maxima in most part of the global ocean.

Table 2 shows the southward and northward flows of dense water masses that have an impact on global climate, namely, the Atlantic Meridional Overturning Circulation (AMOC) and the Antarctic bottom waters export (AABW). The three simulations exhibit a magnitude of the AMOC in agreement with observations (e.g., Talley et al., 2003). Such agreement is not surprising with regard to the broad agreement in MLD

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Observed (a) and modeled sea surface height (SSH) annual standard deviation (m) over 2000–2009 for (c) LOW, (e) CRS, and (g) HIGH. Associated biases are shown in the right column (d, f, and h, respectively) with their indicative rmse values on the left. Biases of zonal means of the annual standard deviation are presented in (b). Grid points where standard deviations of CRS and HIGH differ significantly (95%) have been masked in black on subfigure (f).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Mean daily power spectra of vertical velocity (m³/s²) for year 2000 at the surface (a, b) and at 500-m depth (c, d) for the Agulhas [18°E–28°E, 48°S–38°S] (left) and the Kuroshio [150°E–165°E, 30°N–45°N] (right) regions.</figcaption></figure>

maximum in the North Atlantic, which is known to supply the flow of southward dense water masses. In comparison, larger differences are found for the northward transport of AABW, which mirrors model-data differences in buoyancy fluxes (i.e., driven by hydrodynamics properties) in the Southern Ocean.

### 3.3. Eddies

SSH is commonly used to distinguish the ocean fine-scale structures such as eddies and filaments from the large-scale ocean features (e.g., Penduff et al., 2010). Using the annual standard deviation of the SSH, we show that increasing the horizontal resolution improves the representation of the turbulent dynamics in eddy-active regions (Figure 6). Fine scales of SSH variability are captured at middle and high latitudes in HIGH in agreement with the observations (Figure 6g), whereas LOW does not (Figure 6c; see the quite smooth 0.02-m contour of variability). The coarse-grained solution (Figure 6e) inherits the grain texture of the SSH variability of HIGH; this is particularly noticeable around 45°S between 10°E and 90°E, where LOW fails to represent the variability associated to the Agulhas Current. More broadly, the western boundary currents, which are the principal conduits for heat communication between the equatorial regions and the polar regions, contain much more fine structures and have an extension closer to observations in HIGH and CRS (e.g., the Kuroshio Current).

A F-statistics allows to evaluate how the annual standard deviation of the SSH resulting from CRS or HIGH may significantly differ between each other. Grid points where this difference is significant at 95% have been masked (Figure 6f: in black, as the landmask) and are mainly concentrated at midlatitudes (~45°S/N) in the Southern Ocean and over the northern Pacific (from subtropics to the Kuroshio region). Here also the F-statistics highlights how close are CRS and HIGH in most part of the global ocean, with few local exceptions in the varying eddying regions.

Finally, in order to better quantify the impact of the resolution on eddy-active regions, we compute the mean daily power spectra of the vertical velocity in two distinct regions (Figure 7, Agulhas: a, c and Kuroshio: b, d). The spectra have been computed from daily outputs over year 2000. At the surface, the power peak is quite similar whatever the horizontal resolution considered: it is centered on 800 km over the Agulhas region (Figure 7a) and exceeds 1,000 km in the Kuroshio region (Figure 7b). At 500-m depth, the peak of the
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Observed (a) and modeled means of surface chlorophyll (mg Chl m⁻³) over 2000–2009 for (c) LOW, (e) CRS, and (g) HIGH. Associated biases are shown in the right column (d, f, and h, respectively) with their indicative rmse values on the left. Biases of zonal means are presented in (b). Grid points for which the difference CRS — HIGH is strictly lower than the difference LOW — HIGH are stippled on subfigure (f).</figcaption></figure>

vertical velocity power spectra is clearly shifted to finer scales, reaching ~200 km in the Agulhas region (Figure 7c) and 300–400 km in the Kuroshio region (Figure 7d). Such a shift to finer scales illustrates the increase of the effective resolution in HIGH and in CRS compared to LOW. As expected, the vertical velocity of the coarse-grained solution inherits a large part of the energy present in the 0.25° experiment, as shown by the similarity of their spectra extending to the limit of CRS resolution. HIGH and CRS represent the upper part of mesoscale activity, which are much more energetic (e.g., o(10⁻⁵) for the Agulhas in Figure 7c) than the large-scale activity found in LOW (compared to o(10⁻⁷), respectively).

**4. Marine Biogeochemistry**

Due to close but not similar ocean dynamical fields, and different grid truncations on biogeochemistry computation and transport, the coarse-grained biogeochemical solution (CRS) is expected to differ from the one computed over the 0.25° grid in HIGH (see section 2.2). In the following we will explore the potential of CRS to serve as a cost-efficient emulator for HIGH: the key aspect of the present evaluation is then to determine how close to HIGH is CRS compared to LOW and what does it bring more than LOW.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Idem Figure 8 for sea-air carbon fluxes (g C m⁻²/year).</figcaption></figure>

**4.1. Surface Chlorophyll**

This section focusses on the evaluation of simulated surface chlorophyll against modern observations. This field is routinely used to assess skill of ocean biogeochemical models due to the availability of global-scale remote sensing data product since 1997 (e.g., Carr et al., 2006; Siegel et al., 2002). Besides, surface chlorophyll provides evidence of nutrient and light limitation controlled by the ocean mixed-layer dynamics (e.g., Behrenfeld, 2010; Behrenfeld et al., 2005, 2006).

No major difference in large-scale geographical structure of surface chlorophyll concentrations is found between the various model configurations (Figure 8). All model solutions tend to underestimate the surface chlorophyll concentrations by more than 0.5 mg Chl m⁻³ in the north between 40°N and 80°N in zonal mean (Figure 8b). In comparison, all configurations overestimate surface chlorophyll concentrations by more than 0.15 mg Chl m⁻³ with respect to the observations between 40°S and 60°S in the Southern Ocean.

However, the three solutions exhibit noticeable differences at regional scale, reflecting the impact of the mesoscale eddies on the marine biogeochemistry. Indeed, an increase in horizontal resolution from 1° to 0.25° slightly reduces the too abundant modeled chlorophyll in eddy-active regions (Kuroshio, Agulhas, and Gulf Stream) even though it does not increase the low amounts of chlorophyll seen in the Atlantic upwelling regions. While values of rmse remain very close from one solution to another (0.13 in LOW, 0.12 in CRS, versus 0.11 in HIGH), increasing the horizontal resolution from 1° to 0.25° reduces the bias
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Schematic illustration of the impact of increasing horizontal resolution on ocean dynamics and marine biogeochemistry. Here the uppermost panel depicts the main processes as simulated in LOW (1°), whereas middle and bottom panels highlight the main processes as simulated in HIGH (0.25°) or CRS (coarse-grained solution), respectively.</figcaption></figure>

associated with western boundary currents in HIGH (Figure 8h), but not in CRS (Figure 8f). Monthly snapshots over the Kuroshio region suggest that the nutrients front is more extended in the latitudinal direction in the CRS solution than in HIGH (not shown). This latitudinally over extended fronts of NO₃ are consistent with a more latitudinally extended bias of chlorophyll in the western boundary currents (Figure 8f), although no imprint of these overextension is detected in the SSH annual variability (Figure 6f).

At each grid point of the regular 1° projected means (Figures 8c, 8e, and 8g), the difference between CRS and HIGH is compared to the one between LOW and HIGH. Grid points for which the difference CRS — HIGH is strictly lower than the difference LOW — HIGH are stippled (Figure 8f), and a weighted area percentage is computed. As a result, CRS improves 50.21% of the global chlorophyll area compared to LOW. It has to be noted that the areas where CRS is the closest to HIGH are mostly key areas in terms of physical processes and biophysical interac tions (Figure 8f). This is the case for most of the upwelling coastal areas (e.g., California, Peru-Chili, Benguela, Senegal, and Australia upwelling systems) and for the Southern Ocean where deep mixed layers occur affecting light availability and chlorophyll productivity. Such areas depend on interrelated biogeochemical (e.g., iron and light limitation and size of the phytoplankton species) and physical processes (e.g., MLD) and so are sensitive to the physical improvements inherited from HIGH.

## 4.2. Sea-Air Carbon Fluxes

The impact of refining the horizontal resolution is now evaluated on sea-air carbon fluxes because of its relevance for climate and carbon cycle interactions (e.g., Friedlingstein et al., 2014). Figure 9 shows that all the model configurations broadly capture the regions of carbon uptake (North Atlantic and southern midlatitudes) and outgassing (equatorial Pacific and western Arabian Sea). However, noticeable differences in simulated sea-air carbon fluxes exist across the various model configura tions. While modeled carbon fluxes exceed the observations by almost 20 g C m⁻²/year in zonal mean south of 45°S in LOW (Figures 9b–9d), this bias is largely reduced in HIGH (Figures 9b, 9g, and 9h) and even more in CRS (Figures 9b, 9e, and 9f). Such an outgassing reduction at high resolution is thought to take place because the eddy-induced dynamics counteracts the Ekman transport and reduces the horizontal divergence. Consequently, the vertical upwelling of deep carbon-rich water masses is reduced in HIGH and CRS compared to LOW (Figure 10).

The different bias magnitudes found in CRS and HIGH (particularly in the Agulhas region) illustrate that CRS is a stand-alone solution, as

explained in section 2.2. However, 54.97% of the global mean sea-air carbon fluxes obtained in HIGH are still closer to CRS than to LOW (stippling in Figure 9f).

The chemical capacity for the ocean to take up atmospheric CO₂ is reflected by the carbonate ion concentra tion. The concentration difference between Alk and DIC is recognized as a good proxy for the latter buffer capacity of seawater (Sarmiento & Gruber, 2006). Indeed, CRS positive bias of sea-air carbon fluxes in the Agulhas reflects a lower uptake capacity of the buffer (which is even reversed to CO₂ outgassing) and is explained by a local lack of Alk-DIC concentration (i.e., an excess of DIC) at the surface (Figures 11e and 11f) compared to HIGH and LOW. This local lack is driven by a northward shift of the concentration front at 350 mmol/m³. As detailed in section 5.1, the structures reflecting the physical-biogeochemical coupling are shifted northward in the Agulhas region in CRS attesting that both biogeochemical and dynamical
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Idem Figure 8 for Alk-DIC concentration (mmol/m³), except that no stippling has been added on subfigure (f). Black contours highlight fronts of concentration equal to 300 and 350 mmol/m³.</figcaption></figure>

processes are affected by this movement. This shift is a direct consequence of the averaging processes used in the coarsening algorithm that smooth local structures (e.g., fronts and eddies) of the dynamical fields compared to those of the 0.25° solution (see section 2.2). In spite of this local northward shift of the Alk-DIC fronts, CRS still improves more than the half global ocean compared to LOW.

### 4.3. Nutrients Distributions

The nutrients distribution is a useful metric to assess the global performance of simulated marine biogeochemistry (e.g., Duteil et al., 2012; Ilyina et al., 2013; Séférian et al., 2013).

Taylor diagrams summarize how nutrient distributions behave as a function of depth, by quantifying model-data global averaged spatial correlations and normalized standard deviations between modeled and observed nutrients (PO₄, NO₃, Si, and O₂) at various depth levels (Figure 12). Similar calculations made over the tropics or the Southern Ocean (not shown) confirm the main characteristics of the model as discussed hereafter. As expected, HIGH outbes the other configurations with spatial correlation coefficient and centered standard deviation close to 0.97 and 0.9 (i.e., the barycenter of the shaded grey triangle in Figure 12c). LOW and CRS display results comparable to HIGH in terms of spatial correlation coefficients (0.96) but exhibit greater variations in centered standard deviation (Figures 12a and 12b, respectively; in this case the
<figure><img src="image_16.png" /></figure>

**Figure 12.** Taylor diagrams showing the correspondence of global annual mean distributions of PO₄ (red), NO₃ (blue), Si (green), and O₂ (cyan) between observations and (a) LOW, (b) CRS, and (c) HIGH. Shapes and colors of the symbols are indicative of the considered depth level and tracer, respectively. The shaded grey triangle represents the mean skill score characteristics of the considered experiment. Its coordinates on the angular axis are the maximal correlations, the correlation at the maximal standard deviation, and the correlation at the minimal standard deviation. On the radial axis, its coordinates are the standard deviation at the maximal correlation, the maximal standard deviation, and the minimal standard deviation. The blue (black hatched) polygon superimposed on each Taylor diagram is indicative of the skill of LOW (HIGH, respectively).

shaded grey polygon represents either LOW skill or CRS skill, while the black hatched one displays HIGH skill). LOW tends to underestimate the spatial magnitudes of nutrients concentrations at the surface, whereas it overestimates them at 500 m. Interestingly, CRS (Figure 12b) reproduces more properly the geographical structure of O₂ distribution obtained in HIGH; the skill score is especially close to HIGH in the first levels (i.e., surface and 150 m). However, CRS is less good at replicating the NO₃ and PO₄ distributions at the surface: although both CRS and HIGH share the same ocean dynamics, their vertical dynamics slightly differs (see section 2.2) affecting differently these two tracer concentrations particularly sensitive to vertical motions. The fact that the hatched and grey polygons of HIGH and of CRS match indicates that these two model solutions simulate close nutrient distributions over the 2000–2009 period and hence a quite similar mean state for marine biogeochemistry in global average.

### 4.4. Vertical Distribution of Oxygen

This section focuses on the oxygen vertical distribution, which is relevant for numbers of biological and biogeochemical processes in the ocean, as discussed in Bopp et al. (2013) and Cocco et al. (2013).

The oxygen distribution is shown at 150-m depth, where the signature of the OMZ is particularly visible. Figure 13 highlights that the annual meridional mean structure of the Pacific OMZ is clearly improved when the horizontal resolution is refined. Pacific OMZs in CRS (Figure 13f) and HIGH (Figure 13h) are more constrained along the East Pacific coastal area in the 10°-45° latitudinal bands compared to LOW (Figure 13d), which shows large negative biases exceeding -80 mmol/m³. This improvement is supported by lower values of rmse for HIGH (19.88 mmol/m³) and CRS (19.18 mmol/m³) than for LOW (22.33 mmol/m³). While some tropical biases persist in the 150-m O₂ mean whatever the horizontal resolution considered, the CRS solution is shown to improve 61.58% of the global ocean area compared to LOW (stippling in Figure 13f).

We now assess the vertical profile of oxygen in the tropical area [5°S-5°N]. HIGH and CRS show vertically overextended OMZ structures over the East Pacific and Atlantic basins compared to the climatological observations (Figures 14b and 14c versus 14d). However, the OMZ intensity (i.e., O₂ minimum values in the core of the OMZ) is improved at eddy-admitting resolution and also in the coarse-grained solution: that core strengthening with the horizontal resolution affects the East Pacific, the tropical Atlantic, and the Indian Oceans (Figure 14). In LOW, the O₂ minimum value of the Pacific OMZ remains too high and never goes below 40 mmol/m³ (Figure 14a), whereas it lies between 14 and 27 mmol/m³ in CRS (Figure 14b) in agreement with HIGH (Figure 14c) and the observations (Figure 14d). Moreover, the vertical slope of the O₂ isolines between 1,000 and 2,000 m is close to zero in LOW (Figure 14a) over the East Pacific and gets closer to observations when the horizontal resolution is refined: CRS shows more realistic slopes than LOW.

The mean horizontal extension of the OMZ is also affected by the grid resolution. Over the East Pacific, the horizontal extension of the OMZ is improved with the increase of horizontal resolution (Figure 14). Indeed, in HIGH and CRS the 66-mmol/m³ O₂-contour extends westward of
<figure><img src="image_17.png" /><figcaption>**Figure 13.** Idem Figure 8 for O₂ (mmol/m³) at 150 m.</figcaption></figure>

160°W in agreement with the observations, while it does not reach 150°W in LOW. Similarly, the 118-mmol/m³ O₂-isoline delimiting the OMZ over the Atlantic basin does not extend westward enough in LOW, while it reaches 30°W in all other experiments and in the observations.

Over the Indian Ocean, increasing the horizontal resolution to 0.25° allows to establish a local vertical structure with two minima (around 700 m and in subsurface) closer to the observed structure. While the second maximum disappears in CRS, the coarse-grained solution is closest to the observed O₂ spatial variability in this basin than is LOW. Finally, none of the model solutions success to fully catch the O₂ vertical structure over the tropical West Pacific. All of the model configurations underestimate O₂ concentrations between 500 and 1,500-m depth. LOW displays continuous quasi-zonal O₂ layers. West of 180°E, CRS is in better agreement with the observations than the other model configurations, because it displays the highest values of O₂ in the core of the vertical gradient.

**5. Physical-Biogeochemical Coupling**

The recent literature suggests that major biogeochemical features such as the Southern Ocean carbon fluxes and OMZ could be impacted by mesoscale eddies because of a fine-scale physical-biogeochemical coupling. In this section we assess how far this physical-biogeochemical coupling is represented in our various model configurations.
<figure><img src="image_18.png" /><figcaption>**Figure 14.** O₂ zonal section averaged between [5°S–5°N] for the 2000–2009 period.</figcaption></figure>

**5.1. Carbon Uptake in the Southern Ocean**

In order to quantify the covariability of mesoscale activity with carbon fluxes across our three numerical solutions, the square correlations ($r^2$) between monthly anomalies (with respect to the seasonal cycle) of sea-air carbon fluxes and SSH have been performed at each grid point of the southern ocean (Figure 15) over the whole third cycle (1948–2009). At higher horizontal resolution (HIGH; Figure 15c) significant values of $r^2$ largely exceed 0.5 in the 70°W–180°E band of the Southern Ocean, while LOW only has a small signal higher than 0.35 confined east of South America. Although weaker, a spatial pattern of $r^2$ similar to that of HIGH is detected in CRS (Figure 15b) with, however, a slight northward shift in the Agulhas region south

<figure><img src="image_19.png" /><figcaption>**Figure 15.** Monthly square correlation between modeled anomalies (with respect to the seasonal cycle) of sea-air carbon fluxes (g C m⁻²/year) and SSH (m) over 1948–2009 for (a) LOW, (b) CRS, and (c) HIGH. Shaded areas are significant at 99% following a Student's t test.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 16.** Mean annual daily O₂ horizontal gradient (10⁻⁵ * mmol/m²) averaged between 150 and 500 m. O₂ contour at 80 mmol/m³ (black) is superimposed and depicts the mean OMZ area over the layer. The grey contour on (b) and (c) replicates the 80-mmol/m³ O₂-isoline associated to LOW in (a).</figcaption></figure>

of Africa (east of 13°E) compared to HIGH. In agreement with this north-
ward shift, monthly snapshots over eddy-active regions (i.e., the Agulhas
and Kuroshio regions) suggest that the tortuous NO₃ front obtained in
both HIGH and CRS is broadened in the latitudinal direction in the
CRS solution (not shown). This latitudinally extended fronts of NO₃ are
linked to the sea-air carbon fluxes through the buffered capacity of the
solubility pump. Note that this NO₃ gradient remains extremely smooth
in LOW due to the absence of mesoscale activity at 1°.

Both spatial patterns of r² are consistent with the spatial patterns of the SSH intrinsic variability as highlighted in Penduff et al. (2011). It suggests a potential effect of SSH intrinsic variability on the variability of sea-air carbon fluxes. More generally, the r² pattern suggests that CRS and HIGH generate a part of the intrinsic variability of the ocean dynamics, which in turn affects biogeochemistry. Sérazin et al. (2017) distinguish the intrinsic variability from the forced atmospheric variability. As the signal generated by the intrinsic variability is included in observations, it needs to be represented in model solutions in order to catch some resolution-dependent processes.

**5.2. Oxygen Minimum Zones**

A comparable physical-biogeochemical coupling is found in the tropics for the oxygen distribution. A recent study from Bettencourt et al. (2015) links the OMZ shape to mesoscale activity at the OMZ boundaries. The authors show that the mean positions and paths of mesoscale structure delimit and maintain the OMZ boundaries. Here we diagnose this physical-biogeochemical coupling from the daily O₂ horizontal gradient averaged over a year (2000) and over a layer between 150 and 500 m. When applying this diagnostics to the outputs of the various model configurations, we highlight much more finer scales in HIGH and CRS coinciding with the averaged OMZ contour (marked by the O₂ isolines at 80 mmol/m³ in Figure 16). As expected, the smooth O₂ horizontal gradient found in LOW reveals that this physical-biogeochemical coupling is missing at 1° of horizontal resolution (Figure 16a). As a consequence, there is a better agreement between the Bettencourt et al. (2015) findings and the physical-biogeochemical coupling as simulated in high-resolution configurations (HIGH and CRS) than in the LOW configuration. Mechanisms of OMZ shaping are expected to be much more realistics in the 0.25° solution and also in CRS than in LOW.

In addition, while oxygen regulates global marine cycles of major nutri-
ents and carbon, its long-term changes remain an important source of
uncertainties for the future of marine ecosystems (Breitburg et al.,
2018); its content in the open ocean and coastal waters has been declining
for at least the past half century in response to ocean global warming and
increased global temperatures and nutrients discharged to coastal waters.
In that context, evaluating precisely the OMZs area and volume appears of
major concern (Resplandy, 2018). Interestingly, comparing the averaged
OMZ contour at 80 mmol/m³ obtained in HIGH and CRS with the one
obtained in LOW (grey contour) highlights the nonnegligible error done
for the estimation of the OMZ area at 1° (which is of order 11 to 22% of
the HIGH OMZ volume delimited by O₂ ≤ 80 mmol/m³). This result appears consistent with Dietze and
Loeptien (2013) that emphasize the reduction of the phosphate bias at 0.25° and its potential effect on spur-
ious suboxic areas.
<caption>**Table 3**</caption>

Mean Computational Performances of an Oceanic-Only Configuration at 0.25° of Horizontal Resolution (NEMOV3.6-GELATOV6 Without PISCESv2-Gas), a Classic 0.25° Coupled Simulation (HIGH) and the Coarse-Grained Solution (CRS)

<table><thead><tr><th>Configuration Type</th><th>Coupling Type</th><th>Resolution</th><th>Performance (957 processors)</th></tr></thead><tbody><tr><td>NEMOV3.6-GELATOV6 without PISCESv2-gas</td><td>-</td><td>0.25°</td><td>1 month ~35 min; 1 year ~7 hr</td></tr><tr><td>HIGH</td><td>Classic</td><td>0.25°</td><td>1 month ~2 hr 25 min; 1 year ~29 hr</td></tr><tr><td>CRS</td><td>Coarse-grained</td><td>Dynamics: 0.25° Biogeochemistry: 0.75°</td><td>1 month ~53 min; 1 year ~10 hr 40 min</td></tr></tbody></table>

Note. See section 2.3.1 for details about the two latter simulations. Given durations correspond to the mean values over the whole third CORE cycle (1948–2009).

## 6. Summary and Conclusions

This study explores the potential of a coarse-graining approach to serve as a cost-effective eddy-admitting emulator for global climate-long simulations. Indeed, the use of the online coarsening algorithm allows to increase the horizontal resolution of a coupled ocean biogeochemical simulation above that of its computational grid. In particular, we evaluated how close the coarse-grained biogeochemical solution is to that produced by an eddy-admitting (0.25°) resolution; and finally, to what extent is it closer to the eddy-admitting than is the coarse state-of-the-art (1°) solution?

Increasing the horizontal resolution by the mean of the coarse-graining approach leads to a mean state that is significantly more realistic than a traditional 1° solution, for both ocean dynamics and marine biogeochemistry. This study has shown that the ocean dynamics of our coarse-grained solution (CRS) is close to that of the eddy-admitting solution (HIGH). This has been evidenced by the appearance of finer scales in dynamical fields that are missing at 1° of horizontal resolution (LOW). Indeed, sections 3.2 and 4.1 show that CRS inherits a better representation of eddy-active regions, which tends to improve locally the mean states of the MLD annual maximum. As a result, surface chlorophyll was improved in more than half of the world ocean by our CRS solution compared to LOW, ensuring a better agreement with observations. Our results also suggest that the representation of eddies tends to reduce the magnitude of the sea-air carbon outgassing in the Southern Ocean in HIGH but also in CRS (section 4.2). Here again CRS has been shown to improve more than half of sea-air carbon fluxes of the global ocean compared to LOW.

The analysis of the OMZ vertical mean state (section 4.4) demonstrates that it is not straightforward to determine which horizontal resolution is needed for a configuration to outbest the others. However, grouping several criteria has shown that the CRS solution improves important climate features of the OMZ mean state. Indeed, if LOW may be considered as having a better vertical extension for Pacific and Atlantic OMZ, the coarse-grained solution improves the meridional OMZ structure at 150 m over the Pacific and the Atlantic. In addition, it is characterized by better Pacific OMZ horizontal extensions and vertical structures, as well as more realistic core intensities and Indian OMZ structure.

Further, we highlight (section 5) that one of the major impacts of increasing the horizontal resolution on simulated marine biogeochemistry lies in the introduction of physical-biogeochemical coupling elements. This coupling is in better agreement with the observations and published findings of process-based studies (e.g., Lacour et al., 2017; Villar et al., 2015). Although the biogeochemical solution of CRS is not strictly similar to that of HIGH, the coarse-grained solution has been shown to inherit either comparable or better dynamical and biogeochemical mean states than LOW (e.g., in terms of sea surface height variability, OMZ mean state, and distribution of nutrients), as well as most of the improvements linked to mesoscale activity at eddy-admitting resolution including the physical-biogeochemical coupling (sections 5.1 and 5.2). The latter aspect is decisive to better characterize and understand how physical-biogeochemical interactions affect global climate at multidecadal timescales in a process-based analysis perspective. Moreover, it suggests that the suite of processes as simulated in the coarse-grained configuration are more reliable than in the low-resolution configuration.

The coarse-grained configuration reduces by 60% the computational burden of the equivalent eddy-admitting (0.25°) coupled configuration (Table 3). Thanks to the online coarsening, adding biogeochemistry only represents an extra charge of 42% compared to an oceanic-only configuration with the same number of
processors. The online coarsening algorithm is thus assumed to offer a good trade-off between model resolution and computational cost.

In conclusion, when considering these results in the context of future climate change simulations, or inter-comparison projects, we recommend comparisons with high-resolution or coarse-grained models to discuss findings inferred from coarse-resolution ocean biogeochemical models.

# Appendix A

## A1. Land-Sea Mask Construction

If the $\overline{3i*3j}$ parent grid mesh includes at least one cell of ocean, the coarse-grained land-sea mask mesh is set to ocean. In the same way, for interface fluxes, if at least one cell of the parent mesh boundary contains sea-water, the coarse-grained land-sea mask is set to ocean at the boundary between the two coarse-grained meshes.

## A2. Dynamical Variables Construction at a Given Vertical Depth z

1. Temperature and salinity are coarse-grained down to a lower resolution following a weighted volume mean over the nine cells of the parent grid:

$$X_{\text{CRS},z} = \frac{\sum_{n}^{9} (dx_{n} * dy_{n} * dz_{n} * X_{n,z})}{\sum_{n}^{9} (dx_{n} * dy_{n} * dz_{n})}$$

where n stands for the cell number (as illustrated by Figure 2, each biogeochemical grid cell is built by gathering a 3 by 3 quad of cells from the dynamical parent grid), $X_n$ the value of the dynamical variable within this nth cell, $X_{CRS,z}$ its resulting value over the coarse-grained cell, and $dx_n$, $dy_n$, $dz_n$ the scale factors associated to each direction of the nth cell.

2. Horizontal transport of fluxes (u and v velocities) at the interface between two parent cells are coarse-grained down to a lower resolution following a weighted area sum along one dimension considered as the parent cell border:

$$U_{\text{CRS},z,i+1} = \sum_{j-1}^{j+1} (dx_{i+1} * U_{i+1,z})$$

where $U_{\text{CRS},z,i+1}$ is the resulting value of the velocity field transported through the three cells (between j - 1 and j + 1) interface (the interface here is the i + 1st boundary) at vertical depth level z.

3. Surface fluxes are coarse-grained down to a lower resolution following a two-dimensional weighted area sum over the horizontal:

$$X_{\text{CRS},z} = \sum_{n}^{9} (dx_{n} * dy_{n} * X_{n,z})$$

## A3. Vertical Scale Factors Construction

The vertical scale factors ($df_n$, where f stands for the direction and n for the cell number) of the coarse grid are approached differently depending on the applied operator:

4. Divergence operation must conserve the total volume of water contained within the nine constituent grid cells of the parent mesh and uses a vertical scale factor conserving the total volume of the parent mesh (Figure 3):

$$dz_{\text{CRS}} = \frac{\sum_{n}^{9} (dx_{n} * dy_{n} * dz_{n})}{dx_{\text{CRS}} * dy_{\text{CRS}}}$$

where the numerator represents the sum of the volumes of the nine cells and the denominator the area of the coarse-grained mesh.

5. Gradient operation has to preserve the maximum depth of the parent mesh (Figure 3):

$$dz_{\text{CRS}} = \text{MAX}_{n}^{9}(dz_{n})$$
Acknowledgments

The authors are grateful to the anonymous reviewers and the Associate Editor for their constructive comments. This work was supported by the H2020 project CRESCENDO Coordinated Research in Earth Systems and Climate: Experiments, kKnowledge, Dissemination and Outreach, which received funding from the European Union's Horizon 2020 research and innovation program under grant agreement 641816. We particularly acknowledge the support of the team in charge of the CNRM-CM climate model. Supercomputing time was provided by the Météo-France/DSI supercomputing center. A number of analyzing tools developed at CNRM, or in collaboration with CNRM scientists, is available as Open Source code (see https://opensource.cnrm-game.meteo.fr/). Model simulations used in this paper are available upon request on https://climatedata.cnrm-game.meteo.fr/public/cnrm/. Part of the development of the coarse-graining algorithm in NEMO was funded by the European Commission's 7th Framework Programme, under grant agreement 282672, EMBRACE project. The Ssalto/Duacs altimeter products were produced and distributed by the Copernicus Marine and Environment Monitoring Service (http://www.marine.copernicus.eu).

References

Aumont, O., Ethé, C., Tagliabue, A., Bopp, L., & Gehlen, M. (2015). PISCES-v2: An ocean biogeochemical model for carbon and ecosystem studies. *Geoscientific Model Development*, *8*, 2465–2513.

Aumont, O., Orr, J. C., Jamous, D., Monfray, P., Marti, O., & Madec, G. (1998). A degradation approach to accelerate simulations to steady-state in a 3-D tracer transport model of the global ocean. *Climate Dynamics*, *14*, 101–116.

Axell, L. (2002). Wind-driven internal waves and Langmuir circulations in a numerical ocean model of the southern Baltic Sea. *Journal of Geophysical Research*, *107*, 3204. https://doi.org/10.1029/2001JC000922

Behrenfeld, M. J. (2010). Abandoning Sverdrup's Critical Depth Hypothesis on phytoplankton blooms. *Ecology*, *91*(4), 977–989.

Behrenfeld, M. J., Boss, E., Siegel, D. A., & Shea, D. M. (2005). Carbon-based ocean productivity and phytoplankton physiology from space. *Global Biogeochemical Cycles*, *19*, GB1006. https://doi.org/10.1029/2004GB002299 (also cited as Editor's Choice in Science 307: 646)

Behrenfeld, M. J., O'Malley, R., Siegel, D., McClain, C., Sarmiento, J., Feldman, G., et al. (2006). Climate-driven trends in contemporary ocean productivity. *Nature*, *444*, 752–755.

Bettencourt, J. H., López, C., Hernández-García, E., Montes, I., Sudre, J., Dewitte, B., et al. (2015). Boundaries of the Peruvian oxygen minimum zone shaped by coherent mesoscale dynamics. *Nature Geoscience*, *8*, 937–940. https://doi.org/10.1038/ngeo2570

Blanke, B., & Delecluse, P. (1993). Variability of the tropical Atlantic Ocean simulated by a general circulation model with two different mixed layer physics. *Journal of Physical Oceanography*, *23*(7), 1363–1388.

Bopp, L., Resplandy, L., Orr, J. C., Doney, S. C., Dunne, J. P., Gehlen, M., et al. (2013). Multiple stressors of ocean ecosystems in the 21st century: Projections with CMIP5 models. *Biogeosciences*, *10*, 6225–6245. https://doi.org/10.5194/bg-10-6225-2013

Bouillon, S., Maqueda, M. A. M., Legat, V., & Fichefet, T. (2009). An elastic–viscous–plastic sea ice model formulated on Arakawa B and C grids. *Ocean Modelling*. 27(3), 174–184.

Breitburg, D., Levin, L. A., Oschlies, A., Grégoire, M., Chavez, F. P., Conley, D. J., et al. (2018). Declining oxygen in the global ocean and coastal waters. *Science*, *359*(6371), eaam7240. https://doi.org/10.1126/science.aam7240

Carr, M.-E., Friedrichs, M. A. M., Schmeltz, M., Noguchi Aita, M., Antoine, D., Arrigo, K. R., et al. (2006). A comparison of global estimates of marine primary production from ocean color. *Deep-Sea Research Part II*, *53*(5-7), 741–770. https://doi.org/10.1016/j.dsr2.2006.01.028

Ciais, P., Sabine, C., Bala, G., Bopp, L., Brovkin, V., Canadell, J., et al. (2013). Carbon and other biogeochemical cycles. In T. F. Stocker, D. Qin, G.-K. Plattner, M. Tignor, S. K. Allen, J. Boschung, A. Nauels, Y. Xia, V. Bex, & P. M. Midgley (Eds.), *Climate change 2013: The physical science basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change*. Cambridge, United Kingdom and New York, NY, USA: Cambridge University Press.

Cocco, V., Joos, F., Steinacher, M., Frölicher, T. L., Bopp, L., Dunne, J., et al. (2013). Oxygen and indicators of stress for marine life in multi-model global warming projections. *Biogeosciences*, *10*, 1849–1868. https://doi.org/10.5194/bg-10-1849-2013

Dietze, H., & Loeptien, U. (2013). Revisiting “nutrient trapping” in global coupled biogeochemical ocean circulation models. *Global Biogeochemical Cycles*, *27*, 265–284. https://doi.org/10.1002/gbc.20029

Doney, S. C. (2004). Evaluating global ocean carbon models: the importance of realistic physics. *Global Biogeochemical Cycles*, *18*, GB3017. https://doi.org/10.1029/2003GB002150

Dufour, C. O., Griffies, S. M., de Souza, G. F., Frenger, I., Morrison, A. K., Palter, J. B., et al. (2015). Role of mesoscale eddies in cross-frontal transport of heat and biogeochemical tracers in the Southern Ocean. *Journal of Physical Oceanography*, *45*, 3057–3081.

Duteil, O., Koeve, W., Oschlies, A., Aumont, O., Bianchi, D., Bopp, L., et al. (2012). Preformed and regenerated phosphate in ocean general circulation models: can right total concentrations be wrong? *Biogeosciences*, *9*, 1797–1807. https://doi.org/10.5194/bg-9-1797-2012

Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B., Stouffer, R. J., & Taylor, K. E. (2016). Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization. *Geoscientific Model Development*, *9*, 1937–1958. https://doi.org/10.5194/gmd-9-1937-2016

Fasullo, J. T., & Trenberth, K. E. (2012). A less cloudy future: The role of subtropical subsidence in climate sensitivity. *Science*, *338*(6108), 792–794. https://doi.org/10.1126/science.1227465

Fox-Kemper, B., Danabasoglu, G., Ferrari, R., Griffies, S. M., Hallberg, R. W., Holland, M. M., et al. (2011). Parameterization of mixed layer eddies. III: Implementation and impact in global ocean climate simulations. *Ocean Modelling*, *39*(1-2), 61–78. https://doi.org/10.1016/j.ocemod.2010.09.002

Friedlingstein, P., Meinshausen, M., Arora, V. K., Jones, C. D., Anav, A., Liddicoat, S. K., & Knutti, R. (2014). Uncertainties in CMIP5 Climate Projections due to carbon cycle feedbacks. *Journal of Climate*, *27*, 511–526. https://doi.org/10.1175/JCLI-D-12-00579.1

Garcia, H. E., Locarnini, R. A., Boyer, T. P., Antonov J. I., Baranova, O. K., Zweng, M. M., et al. (2014a). World Ocean Atlas 2013, Volume 3: Dissolved oxygen, apparent oxygen utilization, and oxygen saturation. S. Levitus, Ed., A. Mishonov Technical Ed.; NOAA Atlas NESDIS 75, 27 pp.

Garcia, H. E., Locarnini, R. A., Boyer, T. P., Antonov, J. I., Baranova, O. K., Zweng, M. M., et al. (2014b). World Ocean Atlas 2013, Volume 4: Dissolved inorganic nutrients (phosphate, nitrate, silicate). S. Levitus, Ed., A. Mishonov Technical Ed.; NOAA Atlas NESDIS 76, 25 pp.

Geider, R., MacIntyre, H., & Kana, T. (1998). A dynamic regulatory model of phytoplanktonic acclimation to light, nutrients, and temperature. *Limnology and Oceanography*, *43*(4), 679–694.

Gent, P. R., & Mcwilliams, J. C. (1990). Isopycnal mixing in ocean circulation models. *Journal of Physical Oceanography*, *20*(1), 150–155.

Grenfell, T. C., & Maykut, G. A. (1977). The optical properties of ice and snow in the Arctic Basin. *Journal of Glaciology*, *18*(80), 445–463.

Griffies, S. M., Danabasoglu, G., Durack, P. J., Adcroft, A. J., Balaji, V., Böning, C. W., et al. (2016). OMIP contribution to CMIP6: Experimental and diagnostic protocol for the physical component of the Ocean Model Intercomparison Project. *Geoscientific Model Development*, *9*, 3231–3296. https://doi.org/10.5194/gmd-9-3231-2016

Griffies, S. M., Winton, M., Samuels, B., Danabasoglu, G., Yeager, S., Marsland, S., et al. (2012). Datasets and protocol for the CLIVAR WGOMD Coordinated Ocean-sea ice Reference Experiments (COREs). WCRP Report, 21, pp. 21

Gruber, N., Lackhar, Z., Frenzel, H., Marchesiello, P., Münnich, M., McWilliams, J. C., et al. (2011). Eddy-induced reduction of biological production in eastern boundary upwelling systems. *Nature Geoscience*, *4*, 787–792. https://doi.org/10.1038/NGEO1273

Harrison, C. S., Long, M. C., Lovenduski, N. S., & Moore, J. K. (2018). Mesoscale effects on carbon export: A global perspective. *Global Biogeochemical Cycles*, *32*, 680–703. https://doi.org/10.1002/2017GB005751

Ilyina, T., Six, K. D., Segschneider, J., Maier-Reimer, E., Li, H., & Núñez-Riboni, I. (2013). Global ocean biogeochemistry model HAMOCC: Model architecture and performance as component of the MPI-Earth system model in different CMIP5 experimental realizations. *Journal of Advances in Modeling Earth Systems*, *5*, 287–315. https://doi.org/10.1029/2012MS000178
Key, R.M., A. Olsen, S. van Heuven, S. K. Lauvset, A. Velo, X. Lin, et al. (2015). Global Ocean Data Analysis Project, Version 2 (GLODAPv2), ORNL/CDIAC-162, ND-P093. Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, US Department of Energy, Oak Ridge, Tennessee. https://doi.org/10.3334/CDIAC/OTG.NDP093_GLODAPv2

Lacour, L., Ardyna, M., Stec, K. F., Claustre, H., Prieur, L., Poteau, A., et al. (2017). Unexpected winter phytoplankton blooms in the North Atlantic subpolar gyre. *Nature Geoscience*, *10*, 836–839. https://doi.org/10.1038/ngeo3035

Landschützer, P., Gruber, N. & Bakker, D. C. E. (2015). A 30 year observation-based global monthly gridded sea surface pCO2 product from 1982 through 2011, http://cdiac.ornl.gov/ftp/oceans/spco2_1982_2011 ETH SOM-FFN. Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, US Department of Energy, Oak Ridge, Tennessee. doi: https://doi.org/10.3334/CDIAC/OTG.SPCO2_1982_2011 ETH SOM-FFN

Lauvset, S. K., Key, R. M., Olsen, A., van Heuven, S., Velo, A., Lin, X., et al. (2016). A new global interior ocean mapped climatology: The 1°x1° GLODAP version 2. *Earth System Science Data*, *8*(2), 325–340. https://doi.org/10.5194/essd-8-325-2016

Le Quéré, C., & 70 others (2018). Global carbon budget 2018. *Earth System Science Data*, *10*, 405–448.

Lévy, M., Estublier, A., & Madec, G. (2001). Choice of an advection scheme for biogeochemical models. *Geophysical Research Letters*, *28*, 3725–3728. https://doi.org/10.1029/2001GL012947

Lévy, M., & Martin, A. P. (2013). The influence of mesoscale and submesoscale heterogeneity on ocean biogeochemical reactions. *Global Biogeochemical Cycles*, *27*, 1139–1150. https://doi.org/10.1002/2012GB004518

Lévy, M., Resplandy, L., Klein, P., Capet, X., Iovino, D., & Ethé, C. (2012). Grid degradation of submesoscale resolving ocean models: Benefits for offline passive tracer transport. *Ocean Modelling*, *48*, 1–9. https://doi.org/10.1016/j.oceemod.2012.02.004

Lherminier, P., Mercier, H., Huck, T., Gourcuff, C., Perez, F. F., Morin, P., et al. (2010). The Atlantic Meridional Overturning Circulation and the subpolar gyre observed at the A25-OVIDE section in June 2002 and 2004. *Deep Sea Research Part I*, *57*(11), 1374–1391. https://doi.org/10.1016/j.dsr.2010.07.009

Liu, W., Xie, S.-P., & Liu, J. (2016). Tracking ocean heat uptake during the surface warming hiatus. *Nature Communications*, *7*, 10,926.

Locarnini, R. A., A. V. Mishonov, J. I. Antonov, T. P. Boyer, H. E. Garcia, O. K. Baranova, et al. (2013). World Ocean Atlas 2013, Volume 1: Temperature. S. Levitus, Ed., A. Mishonov Technical Ed.; NOAA Atlas NESDIS 73, 40 pp.

Madec, G. (2008). Nemo ocean engine. *Note du Pôle de modélisation de l'Institut Pierre-Simon Laplace*, *27*, 1–217.

Madec, G., & Imbard, M. (1996). A global ocean mesh to overcome the North Pole singularity. *Climate Dynamics*, *12*, 381–388.

Martinez-Rey, J., Bopp, L., Gehlen, M., Tagliabue, A., & Gruber, N. (2015). Projections of oceanic N₂O emissions in the 21st century using the IPSL Earth system model. *Biogeosciences*, *12*, 4133–4148.

Masotti, I., Belviso, S., Bopp, L., Tagliabue, A., & Bucciarelli, E. (2016). Effects of light and phosphorus on summer DMS dynamics in subtropical waters using a global ocean biogeochemical model. *Environmental Chemistry*, CSIRO Publishing, *13*(2), 379–389.

Mazloff, M. R., Cornuelle, B. D., Gille, S. T., & Verdy, A. (2018). Correlation lengths for estimating the large-scale carbon and heat content of the Southern Ocean. *Journal of Geophysical Research: Oceans*, *123*, 883–901. https://doi.org/10.1002/2017JC013408

Mélia, S. D. (2002). A global coupled sea ice-ocean model. *Ocean Modelling*, *4*, 137–172.

Mikaloff-Fletcher, S., Gruber, N., Jacobson, A. R., Gloor, M., Doney, S. C., Dutkiewicz, S., et al. (2007). Inverse estimates of the oceanic sources and sinks of natural CO₂ and the implied oceanic carbon transport. *Global Biogeochemical Cycles*, *21*, GB1010. https://doi.org/10.1029/2006GB002751

Munday, D. R., Johnson, H. L., & Marshall, D. P. (2014). Impacts and effects of mesoscale ocean eddies on ocean carbon storage and atmospheric pCO₂. *Global Biogeochemical Cycles*, *28*, 877–896. https://doi.org/10.1002/2014GB004836

Orr, J. C., Najjar, R. G., Aumont, O., Bopp, L., Bullister, J. L., Danabasoglu, G., et al. (2017). Biogeochemical protocols and diagnostics for the CMIP6 Ocean Model Intercomparison Project (OMIP). *Geoscientific Model Development*, *10*, 2169–2199. https://doi.org/10.5194/gmd-10-2169-2017

Oschlies, A. (2002). Can eddies make ocean deserts bloom? *Global Biogeochemical Cycles*, *16*(4), 1106. https://doi.org/10.1029/2001GB001830

Penduff, T., Juza, M., Barnier, B., Zika, J., Dewar, W. K., Treguier, A., et al. (2011). Sea level expression of intrinsic and forced ocean variabilities at interannual time scales. *Journal of Climate*, *24*, 5652–5670. https://doi.org/10.1175/JCLI-D-11-00077.1

Penduff, T., Juza, M., Brodeau, L., Smith, G. C., Barnier, B., Molines, J.-M., et al. (2010). Impact of global ocean model resolution on sea-level variability with emphasis on interannual time scales. *Ocean Science*, *6*, 269–284.

Resplandy, L. (2018). Will ocean zones with low oxygen levels expand or shrink? *Nature*, *557*, 314–315.

Resplandy, L., Bopp, L., Orr, J. C., & Dunne, J. P. (2013). Role of mode and intermediate waters in future ocean acidification: Analysis of CMIP5 models. *Geophysical Research Letters*, *40*, 3091–3095. https://doi.org/10.1002/grl.50414

Resplandy, L., Lévy, M., d'Ovidio, F., & Merlivat, L. (2009). Impact of submesoscale variability in estimating the air-sea CO₂ exchange: Results from a model study of the POMME experiment. *Global Biogeochemical Cycles*, *23*, GB1017. https://doi.org/10.1029/2008GB003239

Sallée, J.-B., Speer, K., & Rintoul, S. R. (2010). Zonnally asymmetric response of the Southern Ocean mixed-layer depth to the Southern Annular Mode. *Nature Geoscience*, *3*, 273–279.

Sarmiento, J. L., & Gruber, N. (2006). *Ocean biogeochemical dynamics*. Princeton: Princeton University Press.

Schneider, B., Bopp, L., Gehlen, M., Segschneider, J., Frölicher, T. L., Cadule, P., et al. (2008). Climate-induced interannual variability of marine primary and export production in three global coupled climate carbon cycle models. *Biogeosciences*, *5*(2), 597–614. https://doi.org/10.5194/bg-5-597-2008

Séférian, R., Bopp, L., Gehlen, M., Orr, J. C., Ethé, C., Cadule, P., et al. (2013). Skill assessment of three Earth system models with common marine biogeochemistry. *Climate Dynamics*. https://doi.org/10.1007/s00382-012-1362-8

Sérazin, G., Jaymond, A., Leroux, S., Penduff, T., Bessières, L., Llovel, W., et al. (2017). A global probabilistic study of the ocean heat content low-frequency variability: Atmospheric forcing versus oceanic chaos. *Geophysical Research Letters*, *44*, 5580–5589. https://doi.org/10.1002/2017GL073026

Siegel, D. A., Doney, S. C., & Yoder, J. A. (2002). The North Atlantic spring phytoplankton bloom and the Sverdrup's critical depth hypothesis. *Science*, *296*, 730–733. https://doi.org/10.1126/science.1069174

Skamarock, W. (2004). Evaluating mesoscale NWP models using kinetic energy spectra. *Monthly Weather Review*, *132*, 3019–3032.

St. Laurent, L., Simmons, H., & Jayne, S. (2002). Estimating tidally driven mixing in the deep ocean. *Geophysical Research Letters*, *29*(23), 2106. https://doi.org/10.1029/2002GL015633

Stow, C. A., Jolliff, J., McGillicuddy, D. J. J., Doney, S. C., Allen, J. I., Friedrichs, M. A. M., et al. (2009). Skill assessment for coupled biological/physical models of marine systems. *Journal of Marine Systems*, *76*, 4–15. https://doi.org/10.1016/j.jmarsys.2008.03.011

Sunda, W., & Huntsman, S. (1997). Interrelated influence of iron, light and cell size on marine phytoplankton growth. *Nature*, *390*(6658), 389–392.
Sweeney, E. N., McGillicuddy Jr., D. J., & Buesseler, K. O. (2003). Biogeochemical impacts due to mesoscale eddy activity in the Sargasso Sea as measured at the Bermuda Atlantic Time-series Study (BATS). *Deep-Sea Research II*, *50*, 3017-3039. https://doi.org/10.1016/j.dsr2.2003.07.008

Takahashi, T., Broecker, W., & Langer, S. (1985). Redfield ratio based on chemical-data from isopycnal surfaces. *Journal of Geophysical Research*, *90*, 6907-6924.

Talley, L. D., Reid, J. L., & Robbins, P. E. (2003). Data-based meridional overturning stream functions for the global ocean. *Journal of Climate*, *16*(19), 3213–3226. https://doi.org/10.1175/1520-0442(2003)016

Taylor, K. E. (2001). Summarizing multiple aspects of model performance in a single diagram. *Journal of Geophysical Research*, 106(D7), 7183–7192. https://doi.org/10.1029/2000JD900719

Terhaar, J., Orr, J. C., Gehlen, M., Ethé, C., & Bopp, L. (2018). Model constraints on the anthropogenic carbon budget of the Arctic Ocean. *Biogeosciences*. https://doi.org/10.5194/bg-2018-283 in review

Tilburg, C. E., Subrahmanyam, B., & O'Brien, J. J. (2002). Ocean color variability in the Tasman Sea. *Geophysical Research Letters*, 29(10), 1487. https://doi.org/10.1029/2001GL014071

Valente, A., Sathyendranath, S., Brotas, V., Groom, S., Grant, M., Taberner, M., et al. (2016). A compilation of global bio-optical in situ data for ocean-colour satellite applications. *Earth System Science Data*, 8, 235–252. https://doi.org/10.5194/essd-8-235-2016

Villar, E., Farrant, G. K., Follows, M., Garczarek, L., Speich, S., Audic, S., et al. (2015). Environmental characteristics of Agulhas rings affect interocean plankton transport. *Science*, 348(6237). https://doi.org/10.1126/science.1261447

Zalesak, S. T. (1979). Fully multidimensional flux corrected transport algorithms for fluids. *Journal of Computational Physics*, 31.

Zweng, M. M, Reagan, J. R., Antonov, J. I., Locarnini, R. A., Mishonov, A. V., Boyer, T. P., et al. (2013). World Ocean Atlas 2013, Volume 2: Salinity. S. Levitus, Ed., A. Mishonov Technical Ed.; NOAA Atlas NESDIS 74, 39 pp.