DINO: a diabatic model of pole-to-pole ocean dynamics to assess subgrid parameterizations across horizontal scales

Climate models are limited in resolution by computational constraints. The ocean component is currently resolved at spatial scales between approximately 10 to 100 km, which is too coarse to adequately capture the mesoscale. Eddies at these scales play a major role in the global energy cycle, and therefore it is crucial that they are accurately parameterized. In this context, we propose DINO (DIabatic Neverworld Ocean), an ocean-only model configuration of intermediate complexity designed as a test protocol for eddy parameterizations across a range of horizontal scales. It allows for affordable simulations, even at very high resolution, while crucial aspects of the global ocean such as the meridional overturning circulation (MOC), subtropical and subpolar gyres, or the Antarctic circumpolar current (ACC) are maintained. We compare key metrics across eddy-resolving ($\mathrm{1}/\mathrm{16}\mathit{°}$), eddy-permitting ($\mathrm{1}/\mathrm{4}\mathit{°}$), and eddy parameterizing (1°) simulations to showcase the evaluation of eddy parameterizations in two ways: by testing their impact on the mean state and by directly diagnosing the missing eddy fluxes from coarse-grained high-resolution experiments.

1 Introduction

The vast range of spatial and temporal scales of Earth's climate system and the underlying processes involved make numerical climate simulations a computationally costly endeavor: it requires representing the effect of small scales (<100 km) in long simulations (>500 years). Limited available computational resources therefore impose constraints on the horizontal resolution of future climate projections. At the time of writing this article, the ocean component of such models ranges in grid size between approximately 10 and 100 km (Hewitt et al., 2020). These grid scales coincide with the horizontal scale of geostrophic turbulence (Chelton et al., 1998). A prominent feature associated with turbulence at these scales is the formation of mesoscale eddies. They are the dominant reservoir of kinetic energy (KE) in the ocean and play a key role in its energy cycle (Ferrari and Wunsch, 2009).

Winds inject KE at the surface and, along with heterogeneous buoyancy forcing, they sustain a reservoir of potential energy (PE) at large scales. The thereby excited baroclinic modes and nonlinear interactions between them lead to energy transfers across scales (Charney, 1971). Upon reaching scales close to the deformation radius of the first mode, baroclinic instabilities drive the formation of eddies, thereby converting available PE (APE) to eddy KE (EKE). Turbulence theory predicts an inverse energy cascade, which then transfers KE back to larger scales (Rhines, 1977).

When numerical models do not resolve mesoscale eddies explicitly nor parameterize them accurately, APE is not sufficiently converted to EKE, which in turn is not transferred to the larger scale KE budget. This leads to biases in the model's mean state and variability, adding uncertainty to the future response of climate projections (Hewitt et al., 2020). As computational cost limits the explicit representation of eddies in the global oceans, substantial efforts of the ocean modeling community have been dedicated to the development of eddy parameterizations over the last decades.

Gent and McWilliams (1990) (GM) introduced a parameterization for tracer transport driven by baroclinic instabilities. It diffuses the APE of the large-scale flow, leading to the flattening of isopycnal surfaces, consistent with the effects of resolved baroclinic instabilities. However, it does not account for the injection nor inverse cascade of KE. To address this, KE backscatter parameterizations have been developed to represent the inverse KE cascade from the subgrid scale when eddies are only partially resolved (Jansen and Held, 2014; Juricke et al., 2019; Eden and Greatbatch, 2008; Bachman, 2019).

Recent studies have used machine learning techniques to learn the missing eddy fluxes of momentum and buoyancy directly from high-resolution model data (Zanna and Bolton, 2020; Guillaumin and Zanna, 2021). Once a parameterization is implemented into an ocean general circulation model, it is not straightforward to evaluate the effect on the model solution given the limited availability of high-resolution reference simulations or observations. To address this challenge, we propose an idealized configuration of pole-to-pole ocean dynamics, a DIabatic Neverworld Ocean (DINO). As the name suggests, it is broadly based on the Neverworld2 configuration (Marques et al., 2022; in the following NW2), but it accounts for diabatic processes such as convection, diapycnal mixing, and dense water formation.

Neverworld2 is strictly adiabatic, and for good reasons: the adjustment time is much shorter; the calculation of APE is straightforward; and ultimately the circulation is nearly adiabatic in most parts of the ocean in any case (Pedlosky, 1996). However, these benefits are not without trade-offs. The model does not include explicit tracer equations for temperature and salinity but a single tracer with a linear equation of state. The volume of water masses represented by each layer remains constant after initialization. Important metrics of the climate system, such as the meridional heat transport or the meridional overturning circulation (MOC), are not captured accurately with such a model. The introduced DINO configuration aims to complement Neverworld2 in these aspects but retains its overall objective of serving as a test case for the development and assessment of eddy parameterizations.

In Sect. 2, we describe the DINO configuration in detail, followed by its numerical implementation in Sect. 3. In Sect. 4, we outline the experiment design and present the results in Sects. 5 and 6: a comparison of the mean state at eddy-parameterizing (1°) and eddy-permitting ($\mathrm{1}/\mathrm{4}\mathit{°}$) resolutions in Sect. 5, and an analysis of the missing subgrid momentum fluxes at eddy-permitting resolutions in Sect. 6. We discuss the findings in Sect. 7 and conclude in Sect. 8.

2 The DINO configuration

2.1 Model equations

We solve the primitive equations, which describe an incompressible fluid in a rotating system under the Boussinesq and hydrostatic approximation. The primary prognostic equations are given by the momentum balance for the zonal (u) and meridional (v) components of the velocity vector field $\mathit{u}={\mathit{u}}_{\mathrm{h}}+\mathit{k}w$,

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\partial }_t{\mathit{u}}_{\mathrm{h}}& +{\left[\left(\mathbf{\nabla }\times \mathit{u}\right)\times \mathit{u}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathbf{\nabla }{\mathit{u}}^{\mathrm{2}}\right]}_{\mathrm{h}}\\ & +f\mathit{k}\times {\mathit{u}}_{\mathrm{h}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\mathbf{\nabla }}_{\mathrm{h}}p={\mathit{D}}^{\mathit{u}}+{\mathit{F}}^{\mathit{u}}\end{array}\end{array}$$

and the conservation equations for conservative temperature Θ and absolute salinity S_A

$$\begin{array}{}\text{(2)}& {\displaystyle }{\partial }_t\mathrm{\Theta }+\mathrm{\nabla }\cdot \left(\mathrm{\Theta }\mathit{u}\right)& {\displaystyle }=D^{\mathrm{\Theta }}+F^{\mathrm{\Theta }},\text{(3)}& {\displaystyle }{\partial }_t{S}_{\mathrm{A}}+\mathrm{\nabla }\cdot \left(S_{\mathrm{A}}\mathit{u}\right)& {\displaystyle }=D^{S_{\mathrm{A}}}+F^{S_{\mathrm{A}}},\end{array}$$

where p is the pressure, k is the unit vector in the vertical direction, $F^{\mathit{u}/\mathrm{\Theta }/S_{\mathrm{A}}}$ represent forcing terms and $D^{\mathit{u}/\mathrm{\Theta }/S_{\mathrm{A}}}$ potential parameterizations of small-scale physics for the corresponding prognostic variable. The subscript (…)_h denotes the local vector component in the horizontal plane. The remaining parameters are the Coriolis frequency f and a reference density ρ₀. The vertical velocity component (w) and the pressure (p) follow diagnostic equations assuming that the fluid is incompressible and in hydrostatic equilibrium, that is,

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathbf{\nabla }\cdot \mathit{u}& {\displaystyle }=\mathrm{0},\text{(5)}& {\displaystyle }{\partial }_{z}p& {\displaystyle }=-\mathit{\rho }g,\end{array}$$

where g is the gravitational acceleration and the density ρ is a function of Θ, S_A, and the depth z, defined by the equation of state (EOS) for seawater. We follow the formulation proposed by Roquet et al. (2015a) and approximate a simplified EOS as

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\mathit{\rho }(\mathrm{\Theta },S_{\mathrm{A}},z)=& {\mathit{\rho }}_{\text{ref}}\left(z\right)-\left(a_{\mathrm{0}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{C}_b{\mathrm{\Theta }}_a+T_{\mathrm{h}}z\right){\mathrm{\Theta }}_a\\ & +b_{\mathrm{0}}{S}_a,\end{array}\end{array}$$

where ${\mathrm{\Theta }}_a=\mathrm{\Theta }-\mathrm{10}$ °C and $S_a=S_{\mathrm{A}}-\mathrm{35}$ g kg⁻¹. We use the linear thermal expansion coefficient a₀, thermal cabbeling coefficient C_b, thermobaric coefficient T_h, and linear haline expansion coefficient b₀ following Caneill et al. (2022). The reference density profile ρ_ref is arbitrary for a Boussinesq fluid, as only horizontal density gradients enter the governing equations through the pressure gradient term (Roquet et al., 2015b). We use characteristic values of the in situ density where absolute values are needed, for example ${\mathit{\rho }}_{\text{ref}}(z=\mathrm{0})=\mathrm{1026}$ kg m⁻³ and ${\mathit{\rho }}_{\text{ref}}(z=\mathrm{2000})=\mathrm{1035}$ kg m⁻³ to compute the potential density referenced to the surface or 2000 m depth. Numerical values of all constant parameters and their physical meaning are listed in Table 1.

Table 1 Numerical values, physical units, and a short description of the constant parameters of the model.

2.2 Model domain

The DINO configuration is an idealized Atlantic sector model (Fig. 1). The domain is broadly adapted from the NW2 configuration. It spans 50° in longitude and approximately 70° north and south of the Equator in latitude. The boundaries are closed everywhere except for a zonally periodic re-entrant channel between 45 and 65° S.

Figure 1 Bathymetry of the DINO configuration. The outlines of the continents are displayed in the background for reference only, in order to provide a comparison to the actual Atlantic Ocean. Lateral boundaries are closed everywhere except a periodic re-entrant channel between 45 and 65° S.

The bathymetry of DINO is a two-hemisphere extension of the configuration by Caneill et al. (2022), with a smooth exponential slope at boundaries and channel walls. The slope steepens from a maximum depth of 4000 m toward a minimum depth of 2000 m, above which the domain is bounded by vertical walls. A semicircular ridge at the western edge of the channel represents an idealized Scottia ridge. It introduces a horizontal pressure gradient to the channel at depth, which reduces the momentum of the circumpolar flow and retroflects the deep western boundary currents. We do not add a mid-Atlantic ridge to the bathymetry, unlike in NW2. We found that this prohibits a coherent subpolar gyre and introduces an undesirable separation into two basins with respect to dense water formation and meridional overturning. While such a separation is indeed observed for the recirculation at depth, the idealized geometry and buoyancy forcing of DINO cannot capture this appropriately. The complete analytical formulation of the bathymetry can be found in Appendix E. For the momentum equations, we apply free-slip boundary conditions at the coast and a nonlinear friction term at the bottom. The surface boundary conditions are described in the following.

2.3 Surface forcing

The surface boundary conditions for the prognostic equations are analytically prescribed by zonally uniform profiles (Fig. 2). Regarding momentum, we follow Marques et al. (2022) using a purely zonal wind stress profile τ_u constructed by piecewise cubic interpolation between fixed values of 0, 0.2, −0.1, −0.02, −0.1, 0.1, and 0 N m⁻², respectively, for the latitudes −70, −45, −15, 0, 15, 45, and 70° N (panel a of Fig. 2). The wind stress applies to the momentum equations as a Neumann boundary condition of vertical momentum diffusion, but can also be interpreted as a constant source term to the topmost fluid layer of thickness Δz₀.

$$\begin{array}{}\text{(7)}& {\mathit{F}}^{\mathit{u}}={\displaystyle \frac{{\mathit{\tau }}_u}{{\mathit{\rho }}_{\mathrm{0}}\mathrm{\Delta }{z}_{\mathrm{0}}}}\mathit{i}\end{array}$$

For the tracer equations, we apply Haney-type boundary conditions (Haney, 1971). Temperature and salinity are restored to meridional profiles Θ^⋆ and S^⋆, adapted from Munday et al. (2013) (panel c and d of Fig. 2). As for momentum, they enter the tracer equations as a source term to the topmost fluid layer.

$$\begin{array}{}\text{(8)}& F_{\text{ns}}^{\mathrm{\Theta }}={\displaystyle \frac{\mathrm{1}}{c_p{\mathit{\rho }}_{\mathrm{0}}\mathrm{\Delta }{z}_{\mathrm{0}}}}\left(A_{\mathrm{\Theta }}\right({\mathrm{\Theta }}^{*}-\mathrm{\Theta })-Q_{\text{sr}})\end{array}$$

and

$$\begin{array}{}\text{(9)}& F^{S_{\mathrm{A}}}={\displaystyle \frac{A_{\mathrm{S}}}{{\mathit{\rho }}_{\mathrm{0}}\mathrm{\Delta }{z}_{\mathrm{0}}}}(S^{*}-S_{\mathrm{A}}),\end{array}$$

where we use the same temperature and salinity restoring coefficient A_Θ=40 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\mathrm{°}{\mathrm{C}}^{-\mathrm{1}}$ and $A_{\mathrm{S}}=\mathrm{3.858}\times {\mathrm{10}}^{\mathrm{3}}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ as Caneill et al. (2022) and the specific heat capacity c_p=3991.86 $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$. Note that Eq. (8) only represents the non-solar heat flux, where we subtract the solar heat flux from short-wave radiation Q_sr. This radiative flux does not only enter the topmost fluid layer, but can penetrate the water column. Its absorption corresponds to oligotrophic type I water in the classification of optical properties of seawater by Jerlov (1968), with two wavebands of e-folding scale ζ₀=0.35 m and ζ₁=23 m. The resulting source term in Eq. (2) is given by

$$\begin{array}{}\text{(10)}& F_{\text{sr}}^{\mathrm{\Theta }}={\displaystyle \frac{Q_{\text{sr}}}{c_p{\mathit{\rho }}_{\mathrm{0}}}}{\partial }_z\left[\mathrm{0.58}{e}^{-\frac{z}{{\mathit{\zeta }}_{\mathrm{0}}}}+\mathrm{0.42}{e}^{-\frac{z}{{\mathit{\zeta }}_{\mathrm{1}}}}\right].\end{array}$$

The solar heat flux follows a seasonal cycle (indicated by the blue shading in Fig. 2). The seasonal cycle of the temperature restoring profile lags by 1 month, replicating an atmospheric response to insolation. Its amplitude is larger in the Northern Hemisphere. From Θ^⋆ and S^⋆ we can diagnose an effective density restoring ρ^⋆ (panel e in Fig. 2). Three key observations in ρ^⋆ motivate our choice of restoring profiles:

1. They ensure that water forming at the southern boundary is always denser than the water forming at the northern boundary. This represents the model analog of North Atlantic Deep Water (NADW) located above denser Antarctic Bottom Water (AABW).

2. The meridional density gradient is positive and approximately uniform throughout most of the periodic channel. This sets the strength of the Antarctic circumpolar current (ACC) through thermal wind balance.

3. Water forming at the northern boundary shares isopycnals outcropping within the periodic channel. This condition was found necessary to support a pole-to-pole component of the MOC (Wolfe and Cessi, 2011).

Note that the model has no sea-ice component. This limits the model's ability to accurately represent the formation of AABW. Complete analytical formulations of Θ^⋆, S^⋆, and Q_sr are provided in Appendix B.

Figure 2 Meridional profiles of the surface forcing fields: zonal wind stress (a), solar heat flux (b), temperature (c), salinity (d) and effective density (e) restoring. The density restoring is diagnosed from Θ^⋆ and S^⋆ through the equation of state (6). All profiles are uniform along the zonal direction and depict the yearly average. Blue shading indicates the minimum and maximum value due to the seasonal cycle.

2.4 Diagnosing subgrid fluxes through coarse-graining

Since we aim to compare solutions of the described configuration across different horizontal resolutions, it is useful to introduce a coarse-graining formalism that allows us to average high-resolution results onto a lower-resolution grid. We follow the approach of Mana and Zanna (2014) and use this method to diagnose eddy fluxes from the subfilter scale of the coarse-grained model solution, reflecting the unresolved subgrid scale of coarse model simulations. The idea is that a low-resolution model with ideal parameterizations of these unresolved processes should behave as the coarse-grained high-resolution model. For simplicity, we demonstrate the concept for the momentum equation (Eq. 1), but it also holds for the tracer equations.

The coarse-graining operator, denoted by an overbar $\stackrel{\mathrm{‾}}{\left(\mathrm{\dots }\right)}$, is a Reynolds average weighted by the local grid size. We assume that it only acts in the horizontal direction and commutes with spatial derivatives. While this assumption is not strictly valid for nonlinear dynamics discretized on a sphere (Aluie, 2019), we consider it sufficient for the purposes of this study. Applying the operator to Eq. (1) then yields

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\partial }_t{\stackrel{\mathrm{‾}}{\mathit{u}}}_{\mathrm{h}}& +{\stackrel{\mathrm{‾}}{\left[\left(\mathbf{\nabla }\times \mathit{u}\right)\times \mathit{u}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathbf{\nabla }{\mathit{u}}^{\mathrm{2}}\right]}}_{\mathrm{h}}\\ & +f\mathit{k}\times {\stackrel{\mathrm{‾}}{\mathit{u}}}_{\mathrm{h}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\stackrel{\mathrm{‾}}{\mathbf{\nabla }}}_{\mathrm{h}}\stackrel{\mathrm{‾}}{p}={\mathit{D}}^{\stackrel{\mathrm{‾}}{\mathit{u}}}+{\mathit{F}}^{\stackrel{\mathrm{‾}}{\mathit{u}}}.\end{array}\end{array}$$

By adding the nonlinear advection term of the coarse-grained velocity fields to both sides and rearranging the terms, we arrive at the coarse-grained momentum equation:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\partial }_t{\stackrel{\mathrm{‾}}{\mathit{u}}}_{\mathrm{h}}& +{\left[(\stackrel{\mathrm{‾}}{\mathbf{\nabla }}\times \stackrel{\mathrm{‾}}{\mathit{u}})\times \stackrel{\mathrm{‾}}{\mathit{u}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\stackrel{\mathrm{‾}}{\mathbf{\nabla }}{\stackrel{\mathrm{‾}}{\mathit{u}}}^{\mathrm{2}}\right]}_{\mathrm{h}}\\ & +f\mathit{k}\times {\stackrel{\mathrm{‾}}{\mathit{u}}}_{\mathrm{h}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\stackrel{\mathrm{‾}}{\mathbf{\nabla }}}_{\mathrm{h}}\stackrel{\mathrm{‾}}{p}={\mathit{D}}^{\stackrel{\mathrm{‾}}{\mathit{u}}}+{\mathit{F}}^{\stackrel{\mathrm{‾}}{\mathit{u}}}+{\mathcal{S}}_{\text{sgs}}^{\mathit{u}},\end{array}\end{array}$$

with

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\mathcal{S}}_{\text{sgs}}^{\mathit{u}}=& {\left[(\stackrel{\mathrm{‾}}{\mathbf{\nabla }}\times \stackrel{\mathrm{‾}}{\mathit{u}})\times \stackrel{\mathrm{‾}}{\mathit{u}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\stackrel{\mathrm{‾}}{\mathbf{\nabla }}{\stackrel{\mathrm{‾}}{\mathit{u}}}^{\mathrm{2}}\right]}_{\mathrm{h}}\\ & -{\stackrel{\mathrm{‾}}{\left[\left(\mathbf{\nabla }\times \mathit{u}\right)\times \mathit{u}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathbf{\nabla }{\mathit{u}}^{\mathrm{2}}\right]}}_{\mathrm{h}}.\end{array}\end{array}$$

Equation (12) resembles Eq. (1) for the coarse velocity fields plus an additional eddy source term ${\mathcal{S}}_{\text{sgs}}^{\mathit{u}}$. It arises from the nonlinearity of the advection term and represents interactions of the resolved flow with the unresolved spatial scales (Aluie, 2019). From Eq. (12) we can derive the behavior of the KE budget due to the subgrid terms through the chain rule as follows:

$$\begin{array}{}\text{(14)}& {\partial }_t{\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\stackrel{\mathrm{‾}}{\mathit{u}}}^{\mathrm{2}}\right]}_{\text{sgs}}=\stackrel{\mathrm{‾}}{\mathit{u}}\cdot {\left[{\partial }_t\stackrel{\mathrm{‾}}{\mathit{u}}\right]}_{\text{sgs}}=\mathit{u}\cdot {\mathcal{S}}_{\text{sgs}}^{\mathit{u}}.\end{array}$$

Similar subgrid-scale interaction terms can be derived for the tracer equations as well, namely, ${\mathcal{S}}_{\text{sgs}}^{\mathrm{\Theta }}$ and ${\mathcal{S}}_{\text{sgs}}^S$. Finding good approximations to these terms to close the equations is ultimately what developing parameterizations of subgrid-scale processes is all about.

3 Numerical implementation in NEMO

In the previous section, we introduced the ingredients defining the DINO configuration. Now, we focus on some numerical choices, essential to arrive at the results presented here. These should be seen as user choices for the purposes of this study, but can be adapted for other use cases. We use the well-established NEMO framework, version 4.2.1. For a more detailed description of its numerical methods, we refer to its documentation (Madec et al., 2023).

We use 36 vertical levels with a grid spacing increasing with depth, adapted from Lévy et al. (2010). The levels are time- and space-dependent as they are referenced to the sea surface height obtained from the nonlinear free-surface boundary (quasi-Eulerian, or z^⋆ coordinates in Adcroft and Campin, 2004). In the horizontal direction, we chose an isotropic Mercator grid, so that the horizontal grid spacing Δx=Δy decreases poleward. This is convenient because the first baroclinic deformation radius L_D decreases with latitude similarly, and explicitly resolving mesoscale eddies requires that L_D is captured by at least two grid points (Hallberg, 2013).

We show results for DINO within three classes of models, categorized by their horizontal resolution following Hewitt et al. (2020): eddy-parameterizing with 1°, eddy-permitting with $\mathrm{1}/\mathrm{4}\mathit{°}$, and eddy-resolving with $\mathrm{1}/\mathrm{16}\mathit{°}$ zonal grid spacing. The corresponding reference experiments are denoted R1, R4, and R16, respectively. A detailed description of the spatial discretization is provided in Appendix C.

For the Coriolis and advection term in the momentum equations, we chose an energy- and enstrophy-conserving advection scheme in vector invariant form (see Chap. 5.2 of Madec et al., 2023). In the tracer equations, the flux-corrected transport scheme is employed.

The pressure gradient term is computed directly by a cumulative sum of horizontal density gradients from the surface to the bottom, where the surface pressure follows from the nonlinear free-surface formulation of NEMO. The forcing terms are introduced in Sect. 2.3 and formulated purely analytically. They are simply computed on the model grid and applied as source terms to the right-hand sides. The remaining terms in Eqs. (1)–(3), denoted as $D^{\mathit{u}/\mathrm{\Theta }/S_{\mathrm{A}}}$, represent subgrid parameterizations and can vary depending on the specific experiment and the resolution of the model. We present DINO as a testbed for horizontal closures, so parameterizations for vertical mixing do not vary between experiments. We use the turbulent KE closure adapted from Blanke and Delecluse (1993) with a background vertical eddy viscosity of $\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻¹ and diffusivity of $\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹. Where the model becomes statically unstable, the vertical eddy mixing coefficients are drastically increased to 100 m² s⁻¹ to replicate a fast convective adjustment.

For R1, momentum and tracers are dissipated by Laplacian friction along isopycnal surfaces. The diffusivity and viscosity parameters scale linearly with the grid spacing. For higher horizontal resolution (R4 and R16), we replace the Laplacian friction with a bilaplacian operator. The diffusivity parameter then scales cubically with the grid spacing, following Willebrand et al. (2001). For viscosity, we chose a grid- and flow-dependent model proposed by Smagorinsky (1963). Our choices were motivated by the need to maintain numerical stability while avoiding excessive dissipation of the eddy field. At eddy-parameterizing resolution, we employ the GM parameterization mentioned in Sect. 1 to account for the unresolved tracer transport by mesoscale eddies. The GM coefficient is chosen to be dependent on the local growth rate of baroclinic instability, as suggested by Tréguier et al. (1997). In eddy-permitting and eddy-resolving horizontal resolution, we assume this process to be at least partially resolved and consequently omit the GM parameterization. A summary of the different parameterizations along with their numerical parameter values is listed in Table 2.

Smagorinsky (1963)Tréguier et al. (1997)

Table 2 Differences between DINO reference experiments at eddy-parameterizing, eddy-permitting and eddy-resolving horizontal resolutions. The indicated hCPU are only estimates for 1 simulated year on the Jean Zay supercomputer used for this study and cannot be generalized.

4 Spin-up and model experiments

DINO incorporates diabatic processes, such as dense water formation and diapycnal mixing. These are slow mechanisms that adjust the stratification to perturbations over centuries. Consequently, the spin-up of DINO takes a few thousand years to equilibrate. This is computationally too expensive for eddy-resolving horizontal resolution. Therefore, we compute the spin-up with 1° zonal grid spacing. We start from rest with idealized initial temperature and salinity fields. Their complete analytical description is provided in Appendix D. The model is integrated for 3000 years, until the tracer and velocity fields have reached a quasi-equilibrated state (see Appendix E for details on equilibration). We then use the spun-up state as new initial conditions for all other experiments (Fig. 3). For experiments with a higher horizontal resolution, the tracer fields are first interpolated onto the respective grid using the Python package Xesmf (Zhuang et al., 2023). The interpolation method we used was designed for scalar fields and cannot ensure conservation of some key properties of vector fields, such as divergence or vorticity. Since the velocity fields spin up rather quickly, we chose to initialize all experiments from rest, after interpolating only the tracer and sea surface height fields.

Figure 3 Schematic representation of the experiment design. R1, R4 and R16 denote experiments at 1°, $\mathrm{1}/\mathrm{4}\mathit{°}$ and $\mathrm{1}/\mathrm{16}\mathit{°}$ horizontal resolution, respectively. C16 is the R16 output, but coarse-grained to the same resolution as R4, following Sect. 2.4. The thick black line denotes the initial state for all reference experiments, interpolated from the R1 spin-up to higher horizontal resolutions where necessary. Gray boxes indicate the last 50 years of the 400-year simulations and the last 10 years of the 30-year simulations used for data collection. They also mark the experiments used for the eddy-permitting/eddy-parameterizing and eddy-resolving/eddy-permitting comparison in Sects. 5 and 6.

In the following sections, we present two intercomparisons of DINO at different resolutions. In Sect. 5, we compare R1 and R4 to illustrate parameterized and partially resolved eddy-induced tracer transport across the two resolution regimes. We analyze the effect on the mean circulation and stratification fields. Hence, we afford 400 years of integration time for both experiments to allow for the slow adjustment of the mean state to the change in resolution.

At $\mathrm{1}/\mathrm{16}\mathit{°}$ horizontal resolution we cannot afford this long adjustment. Consequently, we compare only 30 years of the R4 and R16 experiments in Sect. 6. This integration time was found to be sufficient for faster adjusting metrics, such as KE and its transfer across scales. Fortunately, it is precisely these fast adjusting processes in the energy cycle which are addressed by parameterizations for the eddy-permitting regime. However, we cannot make quantitative statements about their effect on the mean state with the short experiments presented here.

Furthermore, we coarse-grain the output of the high-resolution data to the same resolution as R4, denoted as C16. As described in Sect. 2.4, we use an area-weighted Reynolds average on snapshots of the model output for coarse-graining. The convolution is computed at the cell center for all fields, so the velocity components are first interpolated to T-points. This approach treats the velocity components as scalar fields and does not strictly preserve the conservation properties of the three-dimensional vector field. However, it is the simplest approach and considered sufficient for the purposes presented in this study. Before coarse-graining, we smoothed the velocity fields with a simple fixed-factor filter of the GCM-filters Python package (Grooms et al., 2021; Loose et al., 2022), to suppress high-amplitude noise at the grid scale (Perezhogin et al., 2023b).

In the following, we climb the ocean equivalent of Charney's famous ladder (Balaji, 2021) in two steps: from the eddy-parameterizing into the eddy-permitting regime in Sect. 5 and further into the eddy-resolving regime in Sect. 6.

5 The eddy-parameterizing regime

DINO simulates an idealized Atlantic Ocean circulation and reproduces essential features such as the ACC, MOC, subtropical and subpolar gyres, western boundary currents, and dense water formation. In this section, we compare the mean state circulation of DINO for the eddy parameterizing and eddy-permitting reference experiments R1 and R4, respectively.

The barotropic transport stream function reflects the volume transport of the vertically averaged horizontal flow field. Figure 4 shows the time average over the last 50 years of data collection. The horizontal flow depicted in Fig. 4 is similar for both experiments: a system of alternating subpolar, subtropical, and tropical gyres and the intense circumpolar current at the zonally periodic re-entrant channel. For R1, the flow is mostly laminar, while it is disturbed by partially resolved geostrophic turbulence for R4. The western flank of the subtropical gyres and its zonal extension along 40° N strengthen with increased horizontal resolution. In contrast, ACC transport is weaker for the R4 case.

Figure 4 Barotropic transport stream function for (a) 1° reference (R1), (b) $\mathrm{1}/\mathrm{4}\mathit{°}$ reference (R4). Time average over the last 50 years of simulation.

Figure 5 complements the barotropic transport in Fig. 4 with the vertical structure of the circulation. It shows the MOC transport stream function in potential density space, referenced to 2000 m. For the lightest water masses above approximately 32 kg m⁻³ we find two shallow tropical cells, antisymmetric about the Equator. Just below and further poleward extends a second pair of overturning cells, visible between 32 and 34 kg m⁻³. These cells are an imprint of the tropical and subtropical gyres on the meridional overturning, transporting light, warm waters poleward and colder, denser water equatorward. Below approximately 34 kg m⁻³ the antisymmetry about the Equator breaks. There is a distinct overturning cell associated with the subpolar gyre apparent in the Northern Hemisphere, but its return flow partially crosses the Equator along 36 kg m⁻³ iso-lines as part of a pole-to-pole overturning circulation. A deep anti-clockwise overturning cell develops for the densest water masses below 36 kg m⁻³ in R4, but not in R1. We can associate the overturning to model analogues of water masses by having a look at zonally averaged potential density contours, referenced to 2000 m (Fig. 6). The deep cell just below the iso-lines of 36 kg m⁻³ transports the model analogue of AABW, formed by convection close to the southern boundary. The flow of the upper cell is along isopycnals, which outcrop at the northern boundary and the southern end of the periodic channel. This allows interhemispheric, adiabatic transport of NADW formed by convection at the northern boundary. The circulation and associated stratification are broadly similar for both horizontal resolutions and fit the theoretical framework discussed in previous studies (Vallis, 2017; Wolfe and Cessi, 2011; Munday et al., 2013). However, the transport of the diabatic cell in Fig. 5 is larger for R1 compared to R4. And the deep cell in the Southern Ocean almost vanishes for the R1 experiment. This occurs alongside a volume reduction of AABW and a volume increase of NADW, with steeper isopycnals in the Southern Ocean channel.

Figure 5 Meridional overturning circulation (MOC). Contours of the zonally averaged transport streamfunction in density space for the (a) R1 and (b) R4 reference experiment. Solid contours with red shading represent clockwise, dashed contours with blue shading anti-clockwise circulation. The potential density anomalies are referenced to 2000 m depth. Time average over the last 50 years of simulation.

Figure 6 Zonally averaged contours of potential density referenced to 2000 m depth (σ₂) for the (a) R1 and (b) R4 experiment. Dashed vertical lines mark the channel position. Time average over the last 50 years of simulation.

The Atlantic MOC plays an essential role in regulating Earth's climate system, primarily by transporting heat from the tropics to the North Atlantic. This is illustrated in Fig. 7, which shows the zonally and vertically integrated meridional heat transport for the R1 and R4 experiments. The mean transport (solid lines) for R1 is oriented northward everywhere, except for a small band close to the southern boundary where the transport is southward. It is also overall more northward than for the R4 experiment, where we find wider bands of southward transport in the Southern Hemisphere. This partially reflects the changes in circulation described above. A weaker, cross-equatorial meridional overturning in the R4 experiment leads to less northward heat transport. In addition, a strengthened circulation of AABW leads to more southward heat transport in the Southern Hemisphere. We separate the eddy component of meridional heat transport, indicated by the dashed lines in Fig. 7. It is computed as the residual of an average taken over the last 50 years of data. For R1, the eddy-induced advection of heat is almost entirely due to the GM parameterization. The eddy-induced heat transport is remarkably similar for the partially resolved eddy field in R4 and the parameterized one in R1. The most significant difference occurs in the tropics, where large eddies can be explicitly resolved at $\mathrm{1}/\mathrm{4}\mathit{°}$ resolution and partially resolved at 1° resolution. Here, GM is deactivated, causing a very small eddy-induced heat transport. It follows that the striking differences in meridional heat transport between R1 and R4 are not a direct consequence of the partially resolved or parameterized eddy transport, but rather an indirect result of changes in the mean circulation.

Figure 7 Meridional heat transport as a function of latitude in R1 and R4, produced by the mean currents (straight lines) and by mainly the parameterized (R1) or partially resolved (R4) eddies (dashed lines). Time average over the last 50 years of simulation.

To better understand this process we have a closer look at the role of eddies in the Southern Ocean. The strong eastward winds, combined with the absence of zonal boundaries, drive the ACC (Fig. 4). In thermal wind balance, this intense current is characterized by steeply sloping isopycnals. At $\mathrm{1}/\mathrm{4}\mathit{°}$ horizontal resolution, we start to resolve large eddies. In their formation through baroclinic instabilities, they extract APE and thereby act to flatten these steep isopycnals, while decelerating the ACC. We cannot explicitly resolve this physical process at 1° horizontal resolution in high latitudes. Instead, we parameterize it with the GM scheme. Figure 6 suggests that GM with the parameter choice presented here is less efficient at extracting APE than the partially resolved eddy field of R4, leading to steeper isopycnals and a stronger ACC. In the R1 experiment, we observe an ACC transport of 206.0 Sv, whereas it is reduced to 149.7 Sv in the R4 experiment (see Fig. 4). It is through this delicate balance that the parameterized eddies of R1 or the partially resolved eddies of R4 are setting the global stratification and hence circulation.

The presented results show that the idealized DINO configuration is capable of simulating the primary mechanisms through which horizontal resolution influences the mean circulation. We illustrate its utility in testing the impact of eddy parameterizations on the mean state, using the GM parameterization as a case study. Note that the GM coefficient could be tuned to better match the R4 experiments regarding the shown metrics. We chose to use standard values of a commonly used scheme in NEMO to demonstrate the test protocol with a widely recognized subgrid parameterization.

At eddy-permitting horizontal resolution, only the largest scales of geostrophic turbulence (approx. 50 to 200 km) are explicitly resolved, but not its entire spectrum. The unresolved small scales have a major impact on the energy cycle. In the following, we proceed to eddy-resolving horizontal resolution to showcase the processes still missed in the eddy-permitting regime.

6 The eddy-permitting regime

Figure 8 shows a time series of KE integrated over the whole domain for the reference experiments R4, R16, and the coarse-grained high-resolution experiment C16 for the 30-year integration period shown in Fig. 3. We use the last 10 years of simulation for data collection (indicated by the blue shading in Fig. 8). The KE exhibits a slow drift resulting from changes in the mean circulation, discussed in the previous section. However, it occurs on much slower time scales than the processes studied in this section and is considered negligible for our purposes. The total KE in the R16 experiment is more than two and a half times higher compared to the R4 experiment, that is, 2.1×10¹⁸ J as compared to 0.8×10¹⁸ J, respectively. This drastic increase with horizontal resolution is only partially captured by the added range of smaller scales resolved in R16. After removing them by coarse-graining to the same resolution as R4 (C16), most of the total KE is retained. This suggests that a large portion of KE introduced at smaller scales is transferred to larger scales, which are resolved by R4. The process behind this phenomenon is an upscale energy transfer known as the inverse energy cascade in turbulence theory. In the eddy-permitting regime, this source of KE is missing and needs to be parameterized.

Figure 8 Total KE for R4 (pink), R16 (black) and C16 (dashed black). The light-blue shading indicates the time span used for data collection.

With the coarse-graining approach applied to the high-resolution experiment, we can directly diagnose the subgrid eddy momentum fluxes, following Eq. (13). Figure 9 shows surface KE snapshots for R4, R16, and C16 (panel a to c), as well as the instantaneous KE tendency due to the subgrid forcing, as derived in Eq. (14) (panel d). The surface flow is more energetic for the R16 experiment throughout the domain. A fine mesh of small-scale turbulence is revealed at higher resolution, which is not resolved in the R4 experiment. In addition, large-scale filaments show a notable increase of KE, mainly in the Gulf Stream region (upper white box), Southern Ocean region (lower white box), and around the Equator. A large part of the increase is retained after coarse-graining. Again, this demonstrates the inverse energy cascade that transfers KE from smaller to larger scales. It is instantaneously captured by the subgrid forcing of KE in panel (d). However, the simultaneous injection and extraction of KE at various scales obscures any evidence for such an upscale transfer. To better quantify the phenomenon, we analyze the KE spectra across two distinct dynamical regimes: the Gulf Stream region and the Southern Ocean (see Fig. 10). We compute the spectra with a two-dimensional Fourier transform of the instantaneous velocity fields, interpolated to an equidistant Cartesian grid. We use the Python library xrft (Zhuang et al., 2023) and apply a linear detrending along with a Hann window. By integrating along circles of constant wave numbers in the spectral x–y plane, we obtain the isotropic power density spectra.

Figure 9 Surface KE snapshots for the (a) R4, (b) R16 and (c) C16 reference experiment. Panel (d) shows the KE tendency per unit mass due to the subfilter eddy source term of the C16 experiment. White boxes indicate the Gulf-Stream (top) and Southern Ocean (bottom) region. All snapshots are computed for the last time step after 30 years of simulation.

Figure 10 The upper panels (a, b) depict the KE spectrum for the R4, R16, and C16 experiments in the (a) Southern Ocean and (b) Gulf Stream regions. The gray dashed line represents a $k^{-\mathrm{5}/\mathrm{3}}$ power law. The lower panels, (c, d), show the corresponding KE transfer from the sub-filter scale, computed as a cross-spectrum between the sub-grid forcing S_sgs and the coarse-grained velocities. All spectra are derived from surface velocity fields within the white boxes in Fig. 9, interpolated onto an equidistant Cartesian grid.

Panel (a) of Fig. 10 shows power density spectra of surface KE in the Southern Ocean region, indicated by the lower white rectangle in Fig. 9. We include a logarithmic slope $k^{-\frac{\mathrm{5}}{\mathrm{3}}}$, which is predicted for the inverse energy cascade by theory (Graham and Ringler, 2013). The R16 spectrum (in black) follows the slope of $k^{-\frac{\mathrm{5}}{\mathrm{3}}}$ from length scales of approximately 40 to 200 km, where it has its maximum. Below length scales of 40 km it steepens rapidly towards the grid scale where KE is dissipated by numerical and explicit diffusion. For coarser horizontal resolution (R4), the KE is lower across all scales, but has a broadly similar shape to that of R16. The spectrum of C16 resembles that of R4 close to the grid-/filter scale. Towards larger scales, it approaches R16, as these scales remain unaffected by the coarse-graining operator. Panel (c) shows the KE transfer due to the sub-grid flux ${\mathcal{S}}_{\text{sgs}}^{\mathit{u}}$ in C16 as given by Eq. (14). On average, the sub-grid fluxes extract KE for length scales smaller than 70 km and inject it at larger scales. This agrees well with the range of scales predicted by the $k^{-\frac{\mathrm{5}}{\mathrm{3}}}$ power law fitted to the R16 spectrum. The KE transfer has a maximum of $\mathrm{3.0}\times {\mathrm{10}}^{-\mathrm{9}}$ m² s⁻³ injected at scales just above 200 km and a minimum of $-\mathrm{1.2}\times {\mathrm{10}}^{-\mathrm{9}}$ m² s⁻³ extracted at scales around 40 km. When we integrate the mean KE tendency per unit mass across the entire region, we find that the subgrid flux leads to a net increase of KE by $\mathrm{9.6}\times {\mathrm{10}}^{-\mathrm{12}}$ m² s⁻³.

For the Gulf Stream region (see panel b) all power density spectra are considerably lower in KE than those for the Southern Ocean region, across all scales. They are shifted slightly towards larger length scales, due to the increased grid spacing for lower latitudes. The same holds true for the KE transfer shown in panel (d). It is approximately one order of magnitude smaller in the Gulf Stream region. Here, KE is injected above and extracted below horizontal scales of 80 km. The maxima and minima of the transfer are found at similar horizontal scales as in the Southern Ocean, that is, just above 200 km and around 50 km with KE transfers of $\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{10}}$ and $-\mathrm{0.4}\times {\mathrm{10}}^{-\mathrm{10}}$ m² s⁻³, respectively. The net KE tendency per unit mass averaged across the Gulf Stream region is $\mathrm{1.8}\times {\mathrm{10}}^{-\mathrm{12}}$ m² s⁻³. As for the Southern Ocean, this represents a net increase of total KE in the Gulf Stream region. Once again, this illustrates the cascade of KE from the subfilter scale to the resolved flow at eddy-permitting resolution, explaining the weakened total KE of R4, where these scales remain unresolved (see Fig. 8).

7 Discussion

DINO is an idealized model configuration, and as such, comes with inherent limitations. In this section, we discuss some of these limitations and point to possible future improvements.

The model is ocean-only and does not account for any feedbacks from the atmosphere, biosphere, or cryosphere. The latter two are entirely absent, while the atmospheric boundary is only represented through simplified, zonally uniform profiles. DINO also has an idealized bathymetry with no continental shelves, no mid-Atlantic ridge, a closed northern boundary, and an overall reduced volume when compared to the Atlantic. All these choices make it difficult to compare DINO output directly to observations. For example, the strength of the ACC can be modified by making small adjustments to the meridional surface density gradient through the temperature and salinity restoring (not shown). The meridional overturning and associated heat transport in DINO could only align with those of the real Atlantic by coincidence, as the volume and thermodynamic properties of the transported water masses are not the same. However, the lack of realism comes with the advantage of isolating the main processes influenced by horizontal resolution from the overwhelming complexity of Earth's climate system. It also enables simulations at high horizontal resolution for comparatively low computational cost.

Despite the low cost compared to global ocean models, we still encounter limitations regarding the horizontal resolution for this study. We afford 30 years of simulation for the R16 experiment, clearly not enough time for the mean state to fully adjust to the interpolation to higher resolution. Hence, we have no ground truth of the mean state of DINO. Nevertheless, we demonstrate in Sect. 6 that 30 years is adequate to study the subgrid momentum fluxes missed at eddy-permitting resolution. Moreover, the net effect of potential subgrid parameterizations on the mean state can still be evaluated, as demonstrated in Sect. 5.

For R16 we reach a horizontal grid spacing of a few kilometers in the polar oceans of DINO. While this is sufficient for the purposes of this study, we admit that the simulated outputs could possibly differ if spatial resolution were further increased. We could not afford to resolve the submesoscale regime in DINO, which plays a key role in the vertical restratification of the upper ocean and the forward KE cascade towards dissipation (Capet et al., 2008). We did not include any submesoscale parameterization, to avoid interference with the mesoscale parameterizations we tested with DINO. That being said, it would be straightforward to extend the experiment design proposed in Sect. 4 with even higher-resolution experiments when the computational resources are available.

In a step beyond the assessment of mesoscale parameterizations, DINO could be used to train machine-learning-based parameterizations as was done previously for simpler, more idealized models (Zanna and Bolton, 2020; Guillaumin and Zanna, 2021; Perezhogin et al., 2023a). The high-resolution data of R16 captures information on resolved geostrophic turbulence in the Southern Ocean, Gulf Stream extension, boundary currents, equatorial currents, and deep currents along topography. This richness in different dynamic regimes is an advantage over the previously mentioned studies.

8 Conclusions

We introduce the DINO configuration, a diabatic Atlantic basin model of intermediate complexity, designed for the testing of parameterizations in the eddy-parameterizing and eddy-permitting regime. For this purpose, we present results of DINO at horizontal resolutions of 1° (R1), $\mathrm{1}/\mathrm{4}\mathit{°}$ (R4) and $\mathrm{1}/\mathrm{16}\mathit{°}$ (R16).

We compare the mean state of R1 and R4 based on key metrics relevant to the climate system. We use the Gent and McWilliams (1990) parameterization as an illustrative example in our test configuration and showcase that, as implemented, it does not sufficiently extract APE. Particularly in the Southern Ocean channel, this results in steeper isopycnals, subsequently accelerating the ACC through thermal wind balance from 149.7 to 206.0 Sv. It also inhibits the formation of AABW and promotes the formation of NADW. Both these effects are reflected by the MOC. Compared to R4, the circulation of R1 shows a weaker deep overturning cell associated with the transport of AABW and a stronger diabatic cell associated with the cross-equatorial transport of NADW. This also has consequences for the mean meridional heat transport, which we find to be overall more northward in the R1 experiment. This does not necessarily indicate a flaw in the Gent and McWilliams (1990) parameterization, as improved tuning could yield better results. Rather, it demonstrates that the DINO configuration captures useful metrics for the evaluation, development, tuning, or training of novel parameterizations at coarse resolution.

In addition, we compare instantaneous fields of the R4 and R16 experiment to quantify the KE transfers not resolved in the eddy-permitting regime. DINO's total KE at eddy-resolving resolution is approximately two and a half times higher than at eddy-permitting resolution. Using a simple coarse-graining approach, we provide evidence that most of the increase of KE in R16 is retained at scales resolved by R4, which we attribute to an inverse cascade of KE. We diagnose the missing subgrid fluxes of KE at eddy-permitting resolution from the coarse-grained high-resolution data. On average, these fluxes extract KE below and inject KE above scales of around 70 km in the Southern Ocean and 80 km in the Gulf Stream region. The upscale transfer is almost one order of magnitude larger in the Southern Ocean. As for the eddy-parameterizing resolution, these results serve as benchmark metrics for developing and testing novel parameterizations in the eddy-permitting regime. In particular, they hold significant value for an emerging class of data-driven, machine-learning-based parameterizations.

Appendix A: Bathymetry

The shape of the bathymetry is defined through a normalized and tapered exponential function:

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}& g(x,x_{\mathrm{1}},x_{\mathrm{2}},s,d)=\\ & \left\{\begin{array}{ll}\mathrm{1}-\frac{e^{-s(x-x_{\mathrm{1}})}}{\mathrm{1}+e^{-s\mathrm{\Delta }\mathit{\lambda }}}\cdot & \\ (\mathrm{1}-S(x,x_{\mathrm{1}},x_{\mathrm{1}}+d\left)\right),& x_{\mathrm{1}}\le x\le x_{\mathrm{1}}+d,\\ \mathrm{1},& x_{\mathrm{1}}+d\le x\le x_{\mathrm{2}}-d,\\ \mathrm{1}-\frac{e^{s(x-x_{\mathrm{2}})}}{\mathrm{1}+e^{-s\mathrm{\Delta }\mathit{\lambda }}}S(x,x_{\mathrm{2}}-d,x_{\mathrm{2}}),& x_{\mathrm{2}}-d\le x\le x_{\mathrm{2}},\end{array}\right.\end{array}\end{array}$$

where d is the tapering distance, s is a parameter controlling the slope, Δλ is the longitudinal width of the basin, x is one of the horizontal coordinates (λ,φ) with its respective boundary values x₁<x₂, and S is the smooth step function for tapering:

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}& S(x,a,b)=\\ & \left\{\begin{array}{ll}\mathrm{0},& x<a,\\ \mathrm{6}{\left(\frac{x-a}{b-a}\right)}^{\mathrm{5}}-\mathrm{15}{\left(\frac{x-a}{b-a}\right)}^{\mathrm{4}}+\mathrm{10}{\left(\frac{x-a}{b-a}\right)}^{\mathrm{3}},& a\le x\le b,\\ \mathrm{1},& x>b.\end{array}\right.\end{array}\end{array}$$

The tapering ensures that the slope reaches the sea floor without discontinuities after a maximum distance of d from the boundary. This is particularly important to avoid discontinuities in the periodic channel, where opposing lateral boundaries are closest. Consequently, we chose half the channel width Δφ_c as a tapering distance and arrived at the bathymetry b as

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}b(\mathit{\lambda },\mathit{\phi })=& \underset{g_{\mathit{\phi }}}{\underbrace{g\left(\mathit{\phi },{\mathit{\phi }}_{\mathrm{1}},{\mathit{\phi }}_{\mathrm{2}},s_{\mathit{\phi }},{\displaystyle \frac{\mathrm{\Delta }{\mathit{\phi }}_{\mathrm{c}}}{\mathrm{2}}}\right)}}\cdot \\ & \underset{g_{\mathit{\lambda }}}{\underbrace{g\left(\mathit{\lambda },{\mathit{\lambda }}_{\mathrm{1}},{\mathit{\lambda }}_{\mathrm{2}},s_{\mathit{\lambda }},{\displaystyle \frac{\mathrm{\Delta }{\mathit{\phi }}_{\mathrm{c}}}{\mathrm{2}}}\right)}}(H_{max}-H_{min})\\ & +H_{min},\end{array}\end{array}$$

where ${\mathit{\phi }}_{\mathrm{1}}\approx -\mathrm{70}$° N, φ₂≈70° N, ${\mathit{\lambda }}_{\mathrm{1}}=-\mathrm{50}$° E, λ₁=0° E define the domain extent and H_max=2000 m, H_min=4000 m are the maximum and minimum depth of the bathymetry, respectively. A minimum depth of 2000 m is arguably too deep for an accurate representation of the continental bathymetry. Deep western boundary currents flow above 2000 m and their interactions with the continental slope impact the separation of the Gulf Stream (Zhang and Vallis, 2007). This mechanism is not captured here. We chose 2000 m to enable a hybrid vertical coordinate, which employs z-levels in the upper 1000 m and terrain-following σ-levels below; z-levels near the surface are necessary to prevent large errors in the pressure gradient term, particularly near the Equator, where it is not balanced by the Coriolis effect. Ultimately, we want DINO to be usable by ocean model developers in any vertical coordinate. Therefore, we did not revisit the bathymetry for this study, where only pure z-levels are employed. The slope parameter is chosen as $s_{\mathit{\lambda }}=\frac{\mathrm{1}}{\mathrm{3}\mathit{°}}$ and $s_{\mathit{\phi }}=\cos \left(\frac{\mathit{\pi }{\mathit{\phi }}_{max}}{\mathrm{180}}\right)s_{\mathit{\lambda }}$ to account for the grid deformation of the Mercator projection. The channel width is chosen as Δφ_c=20°, ranging from ${\mathit{\phi }}_{{\mathrm{c}}_{\mathrm{1}}}=-\mathrm{65}$° N to ${\mathit{\phi }}_{{\mathrm{c}}_{\mathrm{1}}}=-\mathrm{45}$° N. It is added to the bathymetry when periodic boundary conditions are chosen, by modifying g_λ as

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{g}}_{\mathit{\lambda }}=& g\left(\mathit{\lambda },{\mathit{\lambda }}_{\mathrm{1}},{\mathit{\lambda }}_{\mathrm{2}},s_{\mathit{\lambda }},{\displaystyle \frac{\mathrm{\Delta }{\mathit{\phi }}_{\mathrm{c}}}{\mathrm{2}}}\right)\cdot \\ & \left[\mathrm{1}-g\left(\mathit{\phi },{\mathit{\phi }}_{{\mathrm{c}}_{\mathrm{1}}},{\mathit{\phi }}_{{\mathrm{c}}_{\mathrm{1}}},s_{\mathit{\lambda }},{\displaystyle \frac{\mathrm{\Delta }{\mathit{\phi }}_{\mathrm{c}}}{\mathrm{2}}}\right)\right]\\ & +g\left(\mathit{\phi },{\mathit{\phi }}_{{\mathrm{c}}_{\mathrm{1}}},{\mathit{\phi }}_{{\mathrm{c}}_{\mathrm{1}}},s_{\mathit{\lambda }},{\displaystyle \frac{\mathrm{\Delta }{\mathit{\phi }}_{\mathrm{c}}}{\mathrm{2}}}\right).\end{array}\end{array}$$

Finally, the semi-circular sill is added in the form of a Gaussian ring centered around the midpoint of the Drake passage $({\mathit{\lambda }}_{\mathrm{m}},{\mathit{\phi }}_{\mathrm{m}})=(-\mathrm{50}\mathit{°}\mathrm{E},-\mathrm{55}\mathit{°}\mathrm{N})$, yielding the modified bathymetry as

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}& \stackrel{\mathrm{̃}}{b}(\mathit{\lambda },\mathit{\phi })=\\ & \left\{\begin{array}{ll}b(\mathit{\lambda },\mathit{\phi })+(H_{\text{sill}}-b(\mathit{\lambda },\mathit{\phi }\left)\right)\cdot & \\ \exp \left(\frac{-{\left(\sqrt{(\mathit{\lambda }-{\mathit{\lambda }}_{\mathrm{m}}{)}^{\mathrm{2}}+(\mathit{\phi }-{\mathit{\phi }}_{\mathrm{m}}{)}^{\mathrm{2}}}-\frac{{\mathit{\phi }}_{\mathrm{c}}}{\mathrm{2}}\right)}^{\mathrm{2}}}{s^{\mathrm{2}}}\right)\cdot & \\ S(\mathit{\lambda },{\mathit{\lambda }}_{\mathrm{m}},{\mathit{\lambda }}_{\mathrm{m}}+s),& \mathrm{\text{if }}b(\mathit{\lambda },\mathit{\phi })\le H_{\text{sill}}\\ b(\mathit{\lambda },\mathit{\phi }),& \mathrm{\text{otherwise}},\end{array}\right.\end{array}\end{array}$$

where the depth of the sill is chosen as H_sill=2500 m.

Appendix B: Surface forcing

The temperature and salinity restoring profiles are given by

$$\begin{array}{}\text{(B1)}& {\displaystyle }{\mathrm{\Theta }}^{*}(t,\mathit{\phi })={\mathrm{\Theta }}_{\text{n/s}}^{*}\left(t\right)+({\mathrm{\Theta }}_{\text{eq}}^{*}-{\mathrm{\Theta }}_{\text{n/s}}^{*}(t\left)\right)\cos \left({\displaystyle \frac{\mathit{\pi }\mathit{\phi }}{L_{\mathit{\phi }}}}\right),\text{(B2)}& {\displaystyle }\begin{array}{rl}{S}^{*}\left(\mathit{\phi }\right)=& S_{\text{n/s}}^{*}+(S_{\text{eq}}^{*}-S_{\text{n/s}}^{*})\left(\mathrm{1}+\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }\mathit{\phi }}{L_{\mathit{\phi }}}}\right)\right)/\mathrm{2}\\ & -\mathrm{1.25}{e}^{-{\mathit{\phi }}^{\mathrm{2}}/{\mathrm{7.5}}^{\mathrm{2}}},\end{array}\end{array}$$

where the subscript (…)_n/s denotes the restoration value at the northern or southern boundary, depending on whether φ is located in the Northern or Southern Hemisphere, respectively; L_φ=140° is the approximate meridional extent of the domain in degrees. The time dependency of Θ^⋆ reflects an atmospheric response to the seasonal cycle of solar insolation expressed by

$$\begin{array}{}\text{(B3)}& {\displaystyle }{\mathrm{\Theta }}_{\mathrm{n}}^{*}& {\displaystyle }=\mathrm{5}+\mathrm{3}\cos \left(\mathit{\pi }{\displaystyle \frac{d\left(t\right)-\mathrm{201}}{\mathrm{180}}}\right),\text{(B4)}& {\displaystyle }{\mathrm{\Theta }}_{\mathrm{s}}^{*}& {\displaystyle }=-\mathrm{0.5}-\mathrm{0.5}\cos \left(\mathit{\pi }{\displaystyle \frac{d\left(t\right)-\mathrm{201}}{\mathrm{180}}}\right),\end{array}$$

where d denotes the day of the year. For simplification, 1 year of DINO has 12 months with 30 d each. The seasonal cycle has its maximum and minimum on 21 July and 21 January, respectively, lagging 1 month behind the solstices. It only modulates the meridional boundary values of the restoring temperature and leaves the equatorial value constant as ${\mathrm{\Theta }}_{\text{eq}}^{*}=\mathrm{27.0}$ °C. In Eq. (B2) the northernmost, southernmost, and equatorial restoring salinity are all constant and given by $S_{\mathrm{n}}^{*}=\mathrm{35.0}$ g kg⁻¹, $S_{\mathrm{s}}^{*}=\mathrm{35.1}$ g kg⁻¹, and $S_{\text{eq}}^{*}=\mathrm{37.25}$ g kg⁻¹, respectively.

The idealized meridional profile of short-wave radiative heat flux is given by

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}& Q_{\text{sr}}(t,\mathit{\phi })=\\ & max\left(\mathrm{230}\cos \left({\displaystyle \frac{\mathit{\pi }}{\mathrm{180}}}\left[\mathit{\phi }-\mathrm{23.5}\cos \left(\mathit{\pi }{\displaystyle \frac{d\left(t\right)-\mathrm{171}}{\mathrm{180}}}\right)\right]\right),\mathrm{0}\right).\end{array}\end{array}$$

Appendix C: Grid

Here, we provide details on the domain discretization. NEMO employs a staggered Arakawa C-grid (Arakawa and Lamb, 1977) and, for clarity, all grid points mentioned in this section are defined at the T-point of that grid.

C1 Horizontal discretization

DINO is computed on a sphere and discretized by a horizontally isotropic Mercator grid (see Fig. C1). Therefore, the grid spacing in the latitudinal direction decreases towards the poles. We denote the grid indices in zonal and meridional directions as $i\in [\mathrm{1},I]$ and $j\in [\mathrm{1},J]$, respectively. With a horizontal grid spacing of Δλ, this defines the mesh as

$$\begin{array}{}\text{(C1)}& {\displaystyle }\mathit{\lambda }\left(i\right)& {\displaystyle }={\mathit{\lambda }}_{\mathrm{0}}+\mathrm{\Delta }\mathit{\lambda }\cdot i,\text{(C2)}& {\displaystyle }\mathit{\phi }\left(j\right)& {\displaystyle }={\displaystyle \frac{\mathrm{180}}{\mathit{\pi }}}\cdot \arcsin \left(\tanh \left(\mathrm{\Delta }\mathit{\lambda }{\displaystyle \frac{\mathit{\pi }}{\mathrm{180}}}\cdot j\right)\right),\end{array}$$

where I and J are chosen to span a domain of 50° from eastern to western boundary and approximately 70° from Equator to both northern and southern boundary (not exactly as this cannot be guaranteed by Eq. C2). The longitudinal position of the domain λ₀ is arbitrary, but chosen as 50° W to fit the Atlantic Ocean.

Figure C1 The ocean mesh of DINO at 2° horizontal resolution.

C2 Vertical discretization

In this study we chose K=36 full step, quasi-Eulerian vertical levels. The depth at the vertical level midpoints (see Fig. C2) $k\in [\mathrm{1.5},K+\mathrm{0.5}]$ is given by

$$\begin{array}{}\text{(C3)}& z=a_{\mathrm{2}}+a_{\mathrm{1}}k+a_{\mathrm{0}}{a}_{\text{cr}}\ln \left(\cosh \left({\displaystyle \frac{k-k_{\text{th}}}{a_{\text{cr}}}}\right)\right),\end{array}$$

where

$${\displaystyle }\begin{array}{rl}{a}_{\mathrm{0}}=& \left(\mathrm{\Delta }{z}_{\text{min}}-{\displaystyle \frac{H}{K-\mathrm{1}}}\right)/\left(\tanh \left({\displaystyle \frac{\mathrm{1}-k_{\text{th}}}{a_{\text{cr}}}}\right)\right.\\ & \left.-a_{\text{cr}}\cdot {\displaystyle \frac{\ln \left(\cosh \left(\frac{K-k_{\text{th}}}{a_{\text{cr}}}\right)\right)-\ln \left(\cosh \left(\frac{\mathrm{1}-k_{\text{th}}}{a_{\text{cr}}}\right)\right)}{K-\mathrm{1}}}\right),\\ a_{\mathrm{1}}=& \mathrm{\Delta }{z}_{\text{min}}-a_{\mathrm{0}}\cdot \tanh \left({\displaystyle \frac{\mathrm{1}-k_{\text{th}}}{a_{\text{cr}}}}\right),\\ a_{\mathrm{2}}=& -a_{\mathrm{1}}-a_{\mathrm{0}}\cdot a_{\text{cr}}\cdot \ln \left(\cosh \left({\displaystyle \frac{\mathrm{1}-k_{\text{th}}}{a_{\text{cr}}}}\right)\right),\end{array}$$

and H=4000 m is the bottom depth, Δz_min=10 m the minimum level thickness, k_th=35 the index of the inflection point, and a_cr=10.5 the slope of the tanh.

Figure C2 Adapted from Fig. 5 of Madec and Imbard (1996). The solid black line shows the depth for each vertical level, corresponding to the left y axis. The dashed line shows the respective layer thickness indicated by the right y axis.

Appendix D: Initial Conditions

The ocean is initialized at rest for all experiments described in this work:

$$\begin{array}{}\text{(D1)}& u(\mathit{\lambda },\mathit{\phi },z)=\mathrm{0},v(\mathit{\lambda },\mathit{\phi },z)=\mathrm{0}.\end{array}$$

The vertical temperature and salinity profiles are defined through a combination of hyperbolic tangent functions, taken from the GYRE configuration (Lévy et al., 2010).

$$\begin{array}{}\text{(D2)}& {\displaystyle }\begin{array}{rl}\mathrm{\Theta }\left(z\right)& =\left[\mathrm{16}-\mathrm{12}\tanh \left({\displaystyle \frac{z-\mathrm{400}}{\mathrm{700}}}\right)\right]{\displaystyle \frac{\mathrm{1}-\tanh \left(\frac{\mathrm{500}-z}{\mathrm{150}}\right)}{\mathrm{2}}}\\ & +\left[\mathrm{15}\left(\mathrm{1}-\tanh \left({\displaystyle \frac{z-\mathrm{50}}{\mathrm{1500}}}\right)\right)-\mathrm{1.4}\tanh \left({\displaystyle \frac{z-\mathrm{100}}{\mathrm{100}}}\right)\right.\\ & +\left.\mathrm{7}\left({\displaystyle \frac{\mathrm{1500}-z}{\mathrm{1500}}}\right)\right]{\displaystyle \frac{\mathrm{1}-\tanh \left(\frac{z-\mathrm{500}}{\mathrm{150}}\right)}{\mathrm{2}}},\end{array}\text{(D3)}& {\displaystyle }\begin{array}{rl}{S}_{\mathrm{A}}\left(z\right)& =\left[\mathrm{36.25}-\mathrm{1.13}\tanh \left({\displaystyle \frac{z-\mathrm{305}}{\mathrm{460}}}\right)\right]\\ & \cdot {\displaystyle \frac{\mathrm{1}-\tanh \left(\frac{\mathrm{500}-z}{\mathrm{150}}\right)}{\mathrm{2}}}+{\displaystyle \frac{\mathrm{1}-\tanh \left(\frac{z-\mathrm{500}}{\mathrm{150}}\right)}{\mathrm{2}}}\mathrm{\dots }\\ & \mathrm{\dots }\left[\mathrm{35.55}+\mathrm{1.25}\left({\displaystyle \frac{\mathrm{5000}-z}{\mathrm{5000}}}\right)\right.\\ & -\mathrm{1.62}\tanh \left({\displaystyle \frac{z-\mathrm{60}}{\mathrm{650}}}\right)+\mathrm{0.2}\tanh \left({\displaystyle \frac{z-\mathrm{35}}{\mathrm{100}}}\right)\\ & \left.+\mathrm{0.2}\tanh \left({\displaystyle \frac{z-\mathrm{1000}}{\mathrm{5000}}}\right)\right].\end{array}\end{array}$$

We introduce a horizontal gradient from the Equator to the poles with a simple linear transition of temperature and salinity towards bottom values of Θ(z) and S_A(z) at the poles:

$$\begin{array}{}\text{(D4)}& {\displaystyle }\stackrel{\mathrm{̃}}{\mathrm{\Theta }}(\mathit{\phi },z)=\left(\left(\mathrm{\Theta }\right(z)-\mathrm{\Theta }{|}_{z=\mathrm{0}}){\displaystyle \frac{{\mathit{\phi }}_{\mathrm{1}}-\left|\mathit{\phi }\right|}{{\mathit{\phi }}_{\mathrm{1}}}}+\mathrm{\Theta }{|}_{z=\mathrm{0}}\right),\text{(D5)}& {\displaystyle }{\stackrel{\mathrm{̃}}{S}}_{\mathrm{A}}(\mathit{\phi },z)=\left(\left(S_{\mathrm{A}}\right(z)-S_{\mathrm{A}}{|}_{z=\mathrm{0}}){\displaystyle \frac{{\mathit{\phi }}_{\mathrm{1}}-\left|\mathit{\phi }\right|}{{\mathit{\phi }}_{\mathrm{1}}}}+S_{\mathrm{A}}{|}_{z=\mathrm{0}}\right).\end{array}$$

Appendix E: Spin-up strategy

In Fig. 3 we present a schematic of the spin-up and experiment design employed throughout this study. Here we want to provide insights into how well equilibrated each experiment is with respect to the slowest adjusting processes, e.g. dense water formation. Hence, Fig. E1 shows the water mass proportions of the spin-up and each experiment with respect to our idealized definition of AABW (σ₂>36 kg m⁻³), NADW ($\mathrm{35}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}\le {\mathit{\sigma }}_{\mathrm{2}}\le \mathrm{36}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$), and surface water (σ₂<35 kg m⁻³). We find that the fraction of surface water spins up quickly and remains constant for all experiments, as they are subject to identical surface buoyancy forcing. During the spin-up (panel a), the densest water, AABW, is reduced at the expense of NADW. After approximately 1000 years, the model reaches an equilibrated stratification that remains stable throughout the R1 experiment (panel b). When the horizontal resolution is increased, the eddy-induced tracer stirring intensifies, which leads to a flattening of isopycnals and an increase of AABW at the expense of NADW. This can be seen clearly in the meridional sections of stratification shown in Fig. 6. This formation of dense water is slow. We find that the 400 years of the R4 experiment is not enough to fully adjust the state to the change in resolution (panel d). The 30 years of simulation used to compare the R4 and R16 experiment is far from enough for the states to equilibrate (see panels c and d). In fact, the water mass volumes at depth hardly change during these experiments. Consequently, we compare changes in the mean state only after 400 years of the R1 and R4 experiments (see Sect. 5), and focus on metrics related to the fast-adjusting KE fields for the shorter experiments (see Sect. 6).

Figure E1 Time series of the water masses partition based on an idealized definition of Antarctic Bottom Water (σ₂>36 kg m⁻³), North Atlantic Deep Water ($\mathrm{35}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}\le {\mathit{\sigma }}_{\mathrm{2}}\le \mathrm{36}$ kg m⁻³) and surface water (σ₂<35 kg m⁻³). Values are given as percentages of the total water volume. The 3000 years of spin-up at 1° horizontal resolution (panel a) is continued as experiment R1 for 400 years (panel b). R4 (panel c) and R16 (panel e) are initialized from the same stratification field as R1, indicated by the thick black line. The R4 experiment is continued for 400 years (panel d). Hatched areas indicate the period of data collection for each experiment, corresponding to the gray boxes in Fig. 3.