NeverWorld2: an idealized model hierarchy to investigate ocean mesoscale eddies across resolutions

We describe an idealized primitive-equation model for studying mesoscale turbulence and leverage a hierarchy of grid resolutions to make eddy-resolving calculations on the finest grids more affordable. The model has intermediate complexity, incorporating basin-scale geometry with idealized Atlantic and Southern oceans and with non-uniform ocean depth to allow for mesoscale eddy interactions with topography. The model is perfectly adiabatic and spans the Equator and thus fills a gap between quasi-geostrophic models, which cannot span two hemispheres, and idealized general circulation models, which generally include diabatic processes and buoyancy forcing. We show that the model solution is approaching convergence in mean kinetic energy for the ocean mesoscale processes of interest and has a rich range of dynamics with circulation features that emerge only due to resolving mesoscale turbulence.

1 Introduction

Mesoscale eddies have a profound impact on the transport of properties in the ocean. They affect the currents, stratification, ocean dynamic sea level variability, and uptake of physical and biogeochemical tracers. Eddies thus play an important role in regulating climate on regional and global scales and on timescales of weeks to centuries. Mesoscale eddies form on spatial scales near the baroclinic Rossby deformation radius (Smith and Vallis, 2002; Arbic and Flierl, 2004; Thompson and Young, 2007; Hallberg, 2013). The deformation radius varies regionally between 10–100 km horizontally (Chelton et al., 1998). These scales are too small to be resolved globally in routinely used structured-grid ocean climate simulations and, therefore, must be parameterized.

The most notable scheme of mesoscale eddy transport implemented in climate models is the Gent McWilliams (GM) parameterization, which mimics the effects of baroclinic instability by flattening isopycnals and acting as a net sink of available potential energy (Gent and McWilliams, 1990; Gent et al., 1995). Motivated by the effect on available potential energy, eddy kinetic energy parameterizations have been developed (Cessi, 2008; Eden and Greatbatch, 2008; Marshall and Adcroft, 2010) to keep track of the mechanical energy in idealized simulations (Marshall et al., 2012) and/or scale the GM parameter in global climate simulations (Adcroft et al., 2019). The specific goal of this paper is to present a test for such energy-based parameterizations, although mesoscale parameterizations based on other approaches can also be tested in the same framework.

As the horizontal grid spacing of climate models is refined, such that the grid box size becomes comparable to the deformation scale, a regime commonly referred to as the “gray zone” is reached. A gray zone is present in virtually all eddying simulations with large meridional extent and continental slopes (Hallberg, 2013). In this regime, some eddies are being partially resolved, but the resolution does not allow for their effects on the large-scale current and stratification to be fully accounted for. In particular, the inverse kinetic energy cascade (or backscatter) and the barotropization of the flow remain too weak in both idealized (Jansen and Held, 2014) and global models (Kjellsson and Zanna, 2017). Recent mesoscale parameterizations focus on these two aspects with novel momentum closures (Bachman, 2019; Jansen et al., 2019); the model hierarchy described here is designed to evaluate and contrast such approaches in an affordable way.

The majority of these mesoscale parameterizations have been developed independently, using different dynamical assumptions (e.g., quasi-geostrophic dynamics or primitive equations) and different idealized configurations with limited spatial extent (e.g., double gyre or channel). This approach has led to a lack of coherent and robust analysis on the effect of eddies and their parameterizations on the ocean dynamics.

Here, we present an idealized model to capture the essence of mesoscale eddy dynamics at varying horizontal grid resolutions, investigate the effect of mesoscale eddies on the large-scale dynamics, and provide a framework for testing and evaluating eddy parameterizations. The model allows for a clean and extensive analysis of the dynamics and energetics of the flow as a function of horizontal resolution, which is often limited in primitive-equation and diabatic global models due to computational resources (Hewitt et al., 2020; McClean et al., 2011).

We introduce a model configuration – referred to as NeverWorld2 (NW2), which is an extension of the Southern Hemisphere-only NeverWorld configuration presented in Khani et al. (2019) and Jansen et al. (2019). NW2 is a stacked shallow-water model configuration with idealized geometry comprising a single cross-equatorial basin and a re-entrant channel in the Southern Hemisphere. The configuration is broadly similar to that of Wolfe and Cessi (2009), except NW2 is strictly adiabatic, on a spherical grid, and forced only by winds. The global volume of water in each density layer is set by the initial conditions, with the dynamics determining the spatial distribution of stratification which can adjust locally. A hierarchy of horizontal grid resolutions allows us to consider mesoscale eddies in distinct dynamical regimes, e.g., Southern Ocean-like dynamics, mid-latitude gyres, and equatorial flows. This hierarchy also encompasses coarse (unresolved), gray (partially resolved), and mesoscale eddying (fully resolved) flow regimes in portions of the domain controlled by the selected model resolution.

We discuss the model equations, the convergence of the simulations as a function of resolution, and the energetics of the flow. This paper is meant to be an introduction to the configuration and the datasets for use by the community to understand, test, and evaluate mesoscale dynamics and novel closures.

Figure 1 NeverWorld2: (a) bathymetry (depth in meters); (b) cross-section of bathymetry on the Equator; (c) stratification, $\mathrm{2}{g}_{k+\mathrm{1}/\mathrm{2}}^{\prime }/(h_k+h_{k+\mathrm{1}})$, at initial depths of each interface; (d) latitudinal profile of zonal wind stress forcing (in pascals), constructed from piecewise cubics with nodes indicated by crosses.

Figure 2 (a) Meridional and (b) zonal resolution parameters, for the $\mathrm{1}/\mathrm{32}{}^{\circ }$ configuration, shown as the deformation radius L_d divided by the model meridional Δy or zonal Δx extent, respectively. The deformation scale is computed from the spun-up state. The aspect ratio between meridional and zonal extent is close to unity for most of the domain except for the very highest latitudes of the Southern Ocean, where the deformation radius is least well resolved.

2 Model description

2.1 Model equations

We consider an adiabatic and hydrostatic fluid system with a single thermodynamic constituent, simplified to a linear equation of state. This system can be approximated by N stacked layers of piecewise constant density. The equations of motion for the primary prognostic model variables, which are the zonal flow (u), meridional flow (v), and layer thickness (h), are written in vector-invariant form:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\partial }_t{u}_k-q_k{v}_k{h}_k+{\partial }_x{K}_k+{\partial }_x{M}_k& {\displaystyle }={\mathcal{F}}_k^x,\text{(2)}& {\displaystyle }{\partial }_t{v}_k+q_k{u}_k{h}_k+{\partial }_y{K}_k+{\partial }_y{M}_k& {\displaystyle }={\mathcal{F}}_k^y,\text{(3)}& {\displaystyle }{\partial }_t{h}_k+{\partial }_x\left(u_k{h}_k\right)+{\partial }_y\left(v_k{h}_k\right)& {\displaystyle }=\mathrm{0},\end{array}$$

with the vertical stress divergence and horizontal friction given by

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathcal{F}}_k^x& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_o{h}_k}}\left({\mathit{\tau }}_{k-\mathrm{1}/\mathrm{2}}^x-{\mathit{\tau }}_{k+\mathrm{1}/\mathrm{2}}^x\right)-\mathrm{\nabla }\cdot {\mathit{\nu }}_{\mathrm{4}}\mathrm{\nabla }\left({\mathrm{\nabla }}^{\mathrm{2}}{u}_k\right)\text{(5)}& {\displaystyle }{\mathcal{F}}_k^y& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_o{h}_k}}\left({\mathit{\tau }}_{k-\mathrm{1}/\mathrm{2}}^y-{\mathit{\tau }}_{k+\mathrm{1}/\mathrm{2}}^y\right)-\mathrm{\nabla }\cdot {\mathit{\nu }}_{\mathrm{4}}\mathrm{\nabla }\left({\mathrm{\nabla }}^{\mathrm{2}}{v}_k\right),\end{array}$$

where a subscript k indicates the vertical layer number with k=1 the topmost and k=N the bottommost. We use the short-hand ∂_t, ∂_x, and ∂_y for partial derivatives in time and in zonal and meridional directions, respectively. ∇⋅ is the horizontal divergence, ∇ the horizontal gradient, and ${\mathrm{\nabla }}^{\mathrm{2}}=\mathrm{\nabla }\cdot \mathrm{\nabla }$ the horizontal Laplacian. The MOM6 code (Adcroft et al., 2019) that discretizes these equations makes use of horizontal orthogonal curvilinear coordinates, though here we write the more concise Cartesian coordinate notation for brevity.

Other dynamic quantities that are derived from the primary variables include

• the interface positions, ${\mathit{\eta }}_{k-\mathrm{1}/\mathrm{2}}=-D+{\sum }_{l=k}^N{h}_l$, indicated with half-integer labels;

• the potential vorticity, $q_k=\frac{\mathrm{1}}{h_k}\left(f+{\partial }_x{v}_k-{\partial }_y{u}_k\right)$;

• the kinetic energy per mass, $K_k=\frac{\mathrm{1}}{\mathrm{2}}\left(u_k^{\mathrm{2}}+v_k^{\mathrm{2}}\right)$;

• the Montgomery potential, $M_k={\sum }_{l=\mathrm{1}}^k{g}_{l-\mathrm{1}/\mathrm{2}}^{\prime }{\mathit{\eta }}_{l-\mathrm{1}/\mathrm{2}}$;

• the dynamic lateral viscosity, ${\mathit{\nu }}_{\mathrm{4}}=C_{\mathrm{4}}\frac{{\mathrm{\Delta }}^{\mathrm{4}}}{\mathrm{8}{\mathit{\pi }}^{\mathrm{2}}}\sqrt{{\left({\partial }_x{u}_k-{\partial }_y{v}_k\right)}^{\mathrm{2}}+{\left({\partial }_y{u}_k+{\partial }_x{v}_k\right)}^{\mathrm{2}}},$ where Δ⁴ is the fourth power of the grid-spacing which follows a particular discretization, as proposed in the Appendix of Griffies and Hallberg (2000) and is different than simply using the square of the cell area;

• the vertical stress, ${\mathit{\tau }}_{k-\mathrm{1}/\mathrm{2}}=-A_{\mathrm{v}}\frac{{\mathit{\rho }}_o}{h_{k-\mathrm{1}/\mathrm{2}}}\left({\mathit{u}}_{k-\mathrm{1}}-{\mathit{u}}_k\right)$;

• the bottom stress, ${\mathit{\tau }}_{N+\mathrm{1}/\mathrm{2}}=-C_{\mathrm{d}}{\mathit{\rho }}_o\left|{\mathit{u}}_{\mathrm{B}}\right|{\mathit{u}}_N$, which uses a quadratic drag law and where u_B is the flow averaged over the bottommost 10 m.

The surface wind stress, ${\mathit{\tau }}_{\mathrm{1}/\mathrm{2}}$, is prescribed, fixed in time, and distributed over the top 5 m. The remaining parameters are the reduced gravity of each layer, $g_{k-\mathrm{1}/\mathrm{2}}^{\prime }$; a reference density, ρ_o=1000 kg m⁻³; the Coriolis parameter, f=2Ωsin ϕ (with $\mathrm{\Omega }=\mathrm{7.2921}\times {\mathrm{10}}^{-\mathrm{5}}$ s⁻¹, and ϕ is latitude); the background kinematic vertical viscosity, $A_{\mathrm{v}}=\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{4}}$ m² s⁻¹; a dimensionless bottom drag coefficient, C_d=0.003; and the bottom depth, $z=-D(x,y)$. We have chosen to use a biharmonic dissipation operator with a dimensionless Smagorinsky coefficient of C₄=0.2, which is larger than the recommended range suggested by Griffies and Hallberg (2000). This is to ensure sufficient dissipation in the absence of any other parameterizations of lateral friction. Other scale-aware alternatives (Bachman et al., 2017) arrived at comparable results.

We provide analysis of the energetics here and in subsequent papers (Loose et al., 2022b; Yankovsky et al., 2022). To facilitate such analysis, it is useful to write out the energy budget equations. The kinetic energy (KE) in layer k is given by KE_k=h_hK_k. To obtain the KE equation, we add u_kh_k× Eq. (1), v_kh_k× Eq. (2), and K_k× Eq. (3), which gives

$$\begin{array}{}\text{(6)}& {\partial }_t\left(K_k{h}_k\right)+\mathrm{\nabla }\cdot \left(K_k{\mathit{u}}_k{h}_k\right)+{\mathit{u}}_k{h}_k\cdot \mathrm{\nabla }{M}_k={\mathit{u}}_k{h}_k\cdot {\mathcal{F}}_k+C_k.\end{array}$$

The contribution from the Coriolis terms, $C_k=q_k{h}_k^{\mathrm{2}}(u_k{v}_k-v_k{u}_k)$, should be zero in Eq. (6). However, the numerical model uses a C-grid staggering of variables so that locally the numerical Coriolis terms do not drop out, thus affecting KE. We use the Arakawa and Hsu (1990) discretization of Coriolis terms, which, when integrated over the whole domain, conserves total KE for horizontally non-divergent flow. The term u_kh_k⋅∇M_k is the conversion between potential energy and KE.

The potential energy (PE) at interface $k-\frac{\mathrm{1}}{\mathrm{2}}$ is given by ${\text{PE}}_{k-\mathrm{1}/\mathrm{2}}=\frac{\mathrm{1}}{\mathrm{2}}{g}_{k-\mathrm{1}/\mathrm{2}}^{\prime }{\mathit{\eta }}_{k-\mathrm{1}/\mathrm{2}}^{\mathrm{2}}$. The corresponding PE equation is obtained by vertically summing equation (Eq. 3) from the bottom up to interface $k-\frac{\mathrm{1}}{\mathrm{2}}$ and then multiplying by $g_{k-\mathrm{1}/\mathrm{2}}^{\prime }{\mathit{\eta }}_{k-\mathrm{1}/\mathrm{2}}$:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& {\partial }_t\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{g}_{k-\mathrm{1}/\mathrm{2}}^{\prime }{\mathit{\eta }}_{k-\mathrm{1}/\mathrm{2}}^{\mathrm{2}}\right)+\mathrm{\nabla }\cdot \left(g_{k-\mathrm{1}/\mathrm{2}}^{\prime }{\mathit{\eta }}_{k-\mathrm{1}/\mathrm{2}}\sum _{l=k}^N{\mathit{u}}_l{h}_l\right)\\ & =g_{k-\mathrm{1}/\mathrm{2}}^{\prime }\left(\sum _{l=k}^N{\mathit{u}}_l{h}_l\right)\cdot \mathrm{\nabla }{\mathit{\eta }}_{k-\mathrm{1}/\mathrm{2}}.\end{array}\end{array}$$

When summed in the vertical, the right-hand side of Eq. (7) and the PE conversion term in Eq. (6) are equal.

The available potential energy (APE) is the domain-integrated PE minus the domain-integrated PE of the resting state. The resting state and interface positions used at initialization are given by ${\mathit{\eta }}_{k-\mathrm{1}/\mathrm{2}}(t=\mathrm{0})=max\left(z_{k-\mathrm{1}/\mathrm{2}}^{\mathrm{0}},-D\right)$, where $z_{k-\mathrm{1}/\mathrm{2}}^{\mathrm{0}}$ is a constant nominal position for each interface. In this adiabatic model, changes in the domain-integrated PE are exactly the changes in APE, with no approximations or ambiguity.

2.2 Configuration

The NW2 configuration is set up as follows. The domain is a sector of a sphere with angular width 60^∘ and with a single basin and a re-entrant channel in the Southern Hemisphere. The basin is bounded by solid coasts at 70^∘ N and 70^∘ S. Not extending to the poles avoids infinitesimal spherical coordinate cells as the meridians converge. We use a regular spherical grid rather than a Mercator grid so that the placement of boundaries (extents of the grid) is exactly the same for all resolutions. The regular spherical grid means the cells are distorted with cell-wise average aspect ratio $\mathrm{\Delta }y/\mathrm{\Delta }x\sim \mathrm{1.4}$, exceeding 2 only poleward of 60^∘ N and 60^∘ S. The full bathymetry is shown in Fig. 1a, with a cross-section along the Equator in Fig. 1b. The depth of the continental shelf is 200 m, and it has a nominal width of 2.5^∘. A cubic profile, of width $\mathrm{1}/\mathrm{8}$ of the shelf, connects the shelf to the beach (which has zero depth out to $\mathrm{1}/\mathrm{8}$ of the shelf width). Another cubic profile of width 2.5^∘ connects the shelf to the abyss with nominal depth of 4000 m. A large abyssal ridge of height 2000 m runs north–south down the middle of the basin with a cubic profile and radius of 20^∘. The re-entrant channel spans 60–40^∘ S, and a semi-circular ridge of height 2000 m, radius of 10^∘, and thickness of 2^∘ is centered on the channel opening in the west to block the deep flow through the channel. These deep ridges are idealizations of the Mid-Atlantic Ridge and the Scotia Arc that acts as a sill across the Drake Passage to remove momentum via form drag.

The wind stress is strictly zonal and fixed in time (Fig. 1b) and is an idealization of the mean zonal wind profile (see for example Fig. A1 of Chaudhuri et al., 2013). We construct the wind stress from piecewise cubic functions that interpolate between the values 0, 0.2, −0.1, −0.02, −0.1, 0.1, and 0 Pa at latitudes −70, −45, −15, 0, 15, 45, and 70^∘, respectively. Each interpolation node has zero derivative so that both the wind stress and the curl of wind stress are zero at the north and south boundaries.

The initial stratification and vertical resolution are intimately linked (Fig. 1c). We use 15 layers and choose nominal thicknesses of 25, 50, 100, 125, 150, 175, 200, 225, 250, 300, 350, 400, 500, 550, and 600 m. The adiabatic conditions relieve the model of resolving surface boundary layer processes, but finer resolution near the surface is preserved to accurately capture surface intensification of mesoscale energy (Smith and Vallis, 2002). Actual initial thicknesses are the shallower of this nominal profile and whatever is clipped by topography to yield flat interfaces in the interior. The reduced gravity at each interface has values of 10, 0.0021, 0.0039, 0.0054, 0.0058, 0.0058, 0.0057, 0.0053, 0.0048, 0.0042, 0.0037, 0.0031, 0.0024, 0.0017, and 0.0011 m s⁻². The first value corresponds to the gravitational acceleration at the surface. The implied density profile is approximately exponential with a depth scale of 1000 m (approximate due to rounding input parameter values to two significant digits).

The ratio between the first baroclinic deformation radius (L_D) to the meridional and zonal grid spacing is shown in Fig. 2a and b, respectively. Assuming that to resolve eddies, $L_{\mathrm{D}}/\mathrm{\Delta }x\ge \mathrm{2}$ (Hallberg, 2013), Fig. 2 highlights that the $\mathrm{1}/\mathrm{32}{}^{\circ }$ horizontal resolution simulation explicitly resolves mesoscale eddies, except on the shelves and at very high latitudes of the Southern Ocean.

3 Spinup

The model is initialized from rest at $\mathrm{1}/\mathrm{4}{}^{\circ }$ horizontal grid spacing, with initial conditions described in Sect. 2.2 and depicted in Fig. 1c. The $\mathrm{1}/\mathrm{4}{}^{\circ }$ simulation is run to a quasi-steady state reached by about 3×10⁴ d, in which the total kinetic energy is no longer drifting. Next, the layer thicknesses are interpolated to the $\mathrm{1}/\mathrm{8}{}^{\circ }$ horizontal grid, and the simulation is run again to a quasi-steady state for a few thousand days. Note that for expediency, velocities and transports are not interpolated but rather reset to zero. This step introduces a mild shock, but the model quickly spins up mechanically, as seen in the recovery of kinetic energy levels in Fig. 3. This procedure is repeated for $\mathrm{1}/\mathrm{16}$ and $\mathrm{1}/\mathrm{32}{}^{\circ }$ horizontal grids until the simulation has reached convergence (see Sect. 5 for further description).

Figure 3 Time series of total energy (a), kinetic energy (b), and available potential energy (c) as a function of time for all horizontal grids considered during spinup and equilibration.

Figure 4 Time series of ACC transport, defined as the total (zonal) transport across longitude 0^∘ E for each of the model horizontal grids during spinup and equilibration. Thin lines are 5 d means (where available), and thicker lines are 100 d means. The inset shows the entire time series, which starts from the $\mathrm{1}/\mathrm{4}{}^{\circ }$ configuration at rest.

Figure 5 Time series of the difference in the barotropic transport stream function at 10 ^∘ E between 30 and 50^∘ N for each of the model horizontal grids during spinup and equilibrium. Thin lines are 5 d means (where available), and thicker lines are 100 d means. The inset shows the entire time series, which starts from the $\mathrm{1}/\mathrm{4}{}^{\circ }$ configuration at rest.

Figure 6 Time series of equatorial thermocline depth, defined here as the depth of the fourth layer, at longitude 15^∘ E for each of the model horizontal grids during spinup and equilibration. Thin lines are 5 d means (where available), and thicker lines are 100 d means.

As the model horizontal grid is refined, the total kinetic energy (KE) increases at each transition to finer spacing (Fig. 3); this behavior is expected since the finer dynamical modes contain more kinetic energy. In addition, the available potential energy (APE) decreases at each transition to finer grids possibly because eddies are numerically better resolved and become more efficient at extracting APE. The total energy of the system drops at each transition in grid spacing, since APE dominates the total energy reservoir, but the drop is less than for APE due to the increase in KE.

The circumpolar transport spins up in ∼10⁴ d (Fig. 4) and fluctuates with no discernible drift for the remaining 2×10⁴ d at $\mathrm{1}/\mathrm{4}{}^{\circ }$. There is a small reduction in circumpolar transport at each refinement in grid spacing, just as for APE. Since the baroclinic component of circumpolar transport is related to APE, arguments about the improved efficiency of eddies with refined grid spacing are relevant (Marshall and Radko, 2003).

The difference in barotropic transport stream function between sub-tropical and subpolar regions spins up in ∼ 3000–4000 d in the $\mathrm{1}/\mathrm{4}{}^{\circ }$ simulation and fluctuates without discernible drift for the remaining 2.5×10⁴ d (Fig. 5). This adjustment time is consistent with basin travel times for mid-latitude baroclinic Rossby waves (∼2 cm s⁻¹). The transition to finer grid spacing does lead to an increase in the variability. However, the time mean of this transport index does not change significantly with grid spacing, as expected via Sverdrup theory.

The equatorial thermocline depth at 15^∘ E adjusts on multiple timescales with large oscillations at the start of the $\mathrm{1}/\mathrm{4}{}^{\circ }$ simulation that damp out over ∼ 4000 d (Fig. 6). A statistical equilibrium depth is reached on the order of ∼10⁴ d. There is an adjustment at each transition in grid spacing, but it is smaller than the dynamical noise.

During the spinup, the model exhibits multiple timescales in the diagnostics described. The cost of the spinup to day 4.2×10⁴ (covering $\mathrm{1}/\mathrm{4}$, $\mathrm{1}/\mathrm{8}$, and $\mathrm{1}/\mathrm{16}{}^{\circ }$) is ∼ 17 % of the cost of the 2800 d segment at $\mathrm{1}/\mathrm{32}{}^{\circ }$ resolution. The initialization procedure of consecutively allowing adjustment and then interpolating to finer grids is approximately 75× cheaper than spinning up the model solution entirely at $\mathrm{1}/\mathrm{32}{}^{\circ }$.

4 Mean circulation and mesoscale turbulence

We illustrate the behavior of the model in terms of key properties of the flow, with a focus on transport and energy reservoirs. We compare the high-resolution NW2 with a $\mathrm{1}/\mathrm{32}{}^{\circ }$ horizontal grid, in which mesoscale eddies are explicitly resolved in the majority of the domain, to lower-resolution NW2 with horizontal grids of $\mathrm{1}/\mathrm{4}$, $\mathrm{1}/\mathrm{8}$, and $\mathrm{1}/\mathrm{16}{}^{\circ }$. No mesoscale eddy closures are used in any of the simulations beyond the Smagorinsky friction.

Figure 7 Time-mean barotropic transport stream function (Sv) averaged over the last 500 d of each simulation. Transport contours are shown every 10 Sv.

Figure 8 (a) Mean hydrography across the southern channel (Drake Passage) for $\mathrm{1}/\mathrm{4}$ and $\mathrm{1}/\mathrm{32}{}^{\circ }$ models, depicted by time-mean position of interfaces (isopycnals). Time-mean zonal flow (contour interval of 0.05 m s⁻¹), with interface positions for models at (b) $\mathrm{1}/\mathrm{4}$ and (c) $\mathrm{1}/\mathrm{32}{}^{\circ }$ horizontal grids, respectively. All results are averages over 500 d at longitude 0^∘ E, with eastward flow red and westward flow blue.

Figure 9 Snapshots of depth-integrated kinetic energy, computed as $\frac{\mathrm{1}}{\mathrm{2}}{\sum }_{k=\mathrm{1}}^N{h}_k(u_k^{\mathrm{2}}+v_k^{\mathrm{2}})$.

Figure 10 Depth-integrated kinetic energy as in Fig. 9, but averaged over 500 d.

Most of the large-scale features of the finest-resolution configuration are recognizable in the coarsest-resolution solution. The time-mean circulation is represented by subtropical gyres in both hemispheres, a subpolar gyre in the Northern Hemisphere (Fig. 7), and a series of circumpolar jets in the Southern Hemisphere reminiscent of the Antarctic Circumpolar Current (ACC) (Figs. 8 and 9). As the grid spacing is refined, the strength of the western boundary currents and their extensions seems to increase (Fig. 7), while the circumpolar jet (Fig. 4) decreases slightly (proportionally).

The thermocline depth along the Equator also barely changes with resolution (Fig. 6). Overall, the character of the large-scale circulation is relatively invariant with resolution even though the metrics of various features converge with finer resolution.

Some details of the circulation that differ across resolution include, for example, major differences in upper-ocean stratification in the Southern Ocean between the coarsest and finest resolutions. The interfaces in the fine-resolution simulation are less steep, and the interface outcrops move southward as a result of re-stratification by eddies (Fig. 8a). As the grid spacing is refined, the vertical extent of jets changes (Fig. 8). The time-mean zonal velocity shows either a change in the number of distinct jets or a migration of mean jet position. Most notable is the appearance of a deep westward circulation in the channel, below the depth of the blocking topography (the Scotia Arc downstream of the passage).

Snapshots of depth-integrated KE reveal the eddying behavior of the simulations (Fig. 9). For the coarsest grid ($\mathrm{1}/\mathrm{4}{}^{\circ }$), the flow permits mesoscale eddies, in particular at low latitudes, where the deformation scale is largest. As the grid spacing is refined, the first deformation radius is better resolved, and mesoscale eddies become more ubiquitous and widespread. Note that a casual glance does not readily distinguish the $\mathrm{1}/\mathrm{16}$ and $\mathrm{1}/\mathrm{32}{}^{\circ }$ models.

The time mean of depth-integrated total KE (Fig. 10) becomes spatially smoother at finer grid spacing as eddies become increasingly ubiquitous. The same is true for sea surface height variance (not shown), a frequently used indicator of eddy activity.

5 Convergence

Defining convergence for a turbulent cascade is challenging when using a dynamic viscosity because finer grid spacing will permit ever more variability on finer scales. For the purposes of mesoscale eddies, we are concerned with convergence as manifested by invariance of the large-scale properties as resolution changes, thus implying the upscale transfer of kinetic energy is not changing so that details of the small scales are not affecting the large-scale solution.

Figure 11 Equilibrium total energy, kinetic energy (KE), and available potential energy (APE) as a function of horizontal resolution ($\mathrm{1}/N$), averaged over the last 500 d of each simulation. Note the different scale used for KE (b). The dotted lines are a fit of the form $A+B(N-N_o{)}^p$, and the asymptotic limit, A, is indicated by the dashed line.

Figure 12 (a) Definition of our band-pass filter to extract the mesoscale range (green shading). The small-scale end (solid red line) of the mesoscale range is defined as max(0.5 R_d,Δx), where Δx is the longitudinal grid spacing in the 1/4^∘ degree run. The large-scale end (blue line) of the mesoscale range is defined as min(5 R_d, 500 km), where R_d is the first baroclinic deformation radius diagnosed from the NeverWorld2 simulations; (b) 500 d averaged, zonally and depth-integrated mesoscale kinetic energy for the four simulations in our model hierarchy. The mesoscale kinetic energy is computed via the band-pass filter defined on the left.

Figure 13 The 100 d averaged zonal (one-dimensional) spectra of surface EKE for all resolutions at three latitudes (top to bottom): 50^∘ N (weak mean flow region), 15^∘ N (mid-latitude gyre), and 50^∘ S (ACC). We employ the XRFT Python package (Uchida et al., 2021) for the calculations of spectra. Spectra are computed from the surface meridional eddy velocity fields, defined as deviations from a 100 d averaged meridional velocity. Meridional velocities are less influenced by temporal fluctuations in the large-scale zonal currents and allow a better visualization of mesoscale eddy influences. The 2.5^∘ of data are cut off from the eastern and western boundaries before computing the spectra, and linear detrending and a Hann smoothing window are applied. The wavenumber k_D corresponding to the deformation radius for the $\mathrm{1}/\mathrm{32}{}^{\circ }$ simulation is shown as a dotted black line, and a −3 spectral slope is shown in gray. Note that the high-wavenumber spectral tails are sensitive to smoothing choice and boundary effects.

The total mechanical energy of the system (APE + KE reservoirs) is dominated by the APE. APE is an integral metric of the system and the reservoir of energy that generates mesoscale eddies. In equilibrium, there is a balance of APE generation by winds (pumping/heaving the isopycnals) and the conversion of APE to mesoscale eddy energy and the subsequent damping of mean and mesoscale energy by bottom drag, Smagorinsky friction, and vertical viscosity. The equilibrium levels of energy reservoirs shown in Fig. 11 indicate a diminishing APE transition for each change in grid spacing, with an empirical fit suggesting convergence towards a value that is about 3 % lower than the value in the $\mathrm{1}/\mathrm{32}{}^{\circ }$ simulation. We suggest that this behavior implies we are approaching convergence for resolving the mean APE-to-eddy energy pathway.

The kinetic energy increases less than linearly as grid spacing is refined (middle panel of Figs. 3 and 11); there is a 5-fold increase in the kinetic energy reservoir between the $\mathrm{1}/\mathrm{4}$ and the $\mathrm{1}/\mathrm{32}{}^{\circ }$ simulations. A power law fit indicates convergence but not until kinetic energy is about another factor of 2 higher than in the $\mathrm{1}/\mathrm{32}{}^{\circ }$ simulation.

To test convergence of the mesoscale, Fig. 12b shows the 500 d averaged, zonally and depth-integrated mesoscale kinetic energy. The mesoscale kinetic energy is computed with a band-pass filter as $\frac{\mathrm{1}}{\mathrm{2}}{\sum }_{n=\mathrm{1}}^N\left({\stackrel{\mathrm{̃}}{h}}_n({\stackrel{\mathrm{̃}}{u}}_n^{\mathrm{2}}+{\stackrel{\mathrm{̃}}{v}}_n^{\mathrm{2}})-{\overline{h}}_n({\overline{u}}_n^{\mathrm{2}}+{\overline{v}}_n^{\mathrm{2}})\right)$. Here, $\stackrel{\mathrm{̃}}{\cdot }$ and $\overline{\cdot }$ denote low-pass filters, with filter scales defined by the red and blue solid lines in Fig. 12a, respectively. For the respective filter operations, we use the Python package gcm-filters (Loose et al., 2022a) with a Gaussian filter shape.

As the horizontal grid spacing is refined from $\mathrm{1}/\mathrm{4}$ to $\mathrm{1}/\mathrm{32}{}^{\circ }$, we observe steadily growing mesoscale eddy activity in all regions of the domain (Fig. 12b), but the increase from $\mathrm{1}/\mathrm{16}$ to $\mathrm{1}/\mathrm{32}{}^{\circ }$ is smaller compared to the increase from $\mathrm{1}/\mathrm{4}$ to $\mathrm{1}/\mathrm{8}{}^{\circ }$ and from $\mathrm{1}/\mathrm{8}$ to $\mathrm{1}/\mathrm{16}{}^{\circ }$. An exception is the energetic recirculation region north of Drake Passage near 38^∘S, where the mesoscale kinetic energy increases considerably as we change the grid spacing from $\mathrm{1}/\mathrm{16}$ to $\mathrm{1}/\mathrm{32}{}^{\circ }$. We speculate that this increase is due to large-scale standing meanders that develop north of Drake Passage due to more finely resolved topography (Kong and Jansen, 2020; Barthel et al., 2017).

The zonal wavenumber spectra of EKE at various latitudes (Fig. 13) show that at all latitudes, there is a clear gain in kinetic energy between the $\mathrm{1}/\mathrm{4}$ and the $\mathrm{1}/\mathrm{32}{}^{\circ }$ simulations, at all wavenumbers. At higher resolutions, a pronounced peak in the EKE at scales at or slightly above the deformation scale develops. The peak is the result of an increasingly resolved mesoscale eddy-driven inverse KE cascade. The notable exception is the ACC region, where there is less sensitivity to resolution. As hinted by the snapshots of KE, convergence at low latitude is achieved faster than at high latitudes. The spectra generally suggest the large scales are more similar between $\mathrm{1}/\mathrm{16}$ and $\mathrm{1}/\mathrm{32}{}^{\circ }$, although full convergence has not yet been obtained.

6 Conclusions

In this paper, we introduce an idealized ocean model configuration, NeverWorld2 (NW2), that resolves mesoscale eddy dynamics in a basin-scale context. NW2 is a stacked shallow-water model using the MOM6 dynamical core (Adcroft et al., 2019), configured with a single cross-equatorial basin and a re-entering channel in the Southern Hemisphere with idealized geometry. Because NW2 is strictly adiabatic, with no parameterizations in the vertical direction, the timescales of adjustment are controlled entirely by dynamics rather than far slower thermodynamic processes. The paper serves as an introduction to the model and grid resolution hierarchy and to the datasets for use by the community.

For the purposes of analyzing the role of mesoscale eddies and deriving and testing parameterizations, we provide evidence that the finest grid spacing shown, $\mathrm{1}/\mathrm{32}{}^{\circ }$, is practically converged. The large-scale metrics of APE and gyre transport, the latitudinal analysis of KE, and the spectral analysis of EKE all suggest convergence of the largest scales when comparing solutions with the $\mathrm{1}/\mathrm{16}$ and $\mathrm{1}/\mathrm{32}{}^{\circ }$ grid spacings.

We have used the Smagorinsky dynamic biharmonic viscosity of Griffies and Hallberg (2000) for the momentum closure in these simulations, including for the finest grid spacing that defines our converged “truth”. Previous work has shown a sensitivity of the details of the forward and inverse turbulent cascade to the form of dissipation (e.g., Smith and Vallis, 2002; Arbic and Flierl, 2004; Thompson and Young, 2007; Arbic et al., 2012; Pietarila Graham and Ringler, 2013; Soufflet et al., 2016; Pearson et al., 2017; Bachman et al., 2017; Treguier et al., 1997). In refining grid spacing here, we seek to converge on resolving the mechanism of interaction between the mesoscale eddies and the large-scale circulation, by showing it to be diminishingly dependent on the details of dissipation near the bottom of the forward enstrophy cascade; this dissipation scale is more and more separated from the eddy production scale. This assumption could be tested by trying alternative closures and evaluating what tuning is necessary to give the same upscale energy flux. We could have used one of several closures proposed as scale-aware schemes to use in the eddy-permitting regime (e.g., Anstey and Zanna, 2017; Bachman et al., 2017). However, using them in the baseline configuration described here would hinder a fair evaluation of those schemes, which we plan to do in the near future. Nevertheless, we acknowledge that the absolute magnitude of eddy energy, the span of the inverse cascade, and other metrics depend somewhat on the choice of viscous closure.

The model spatial-resolution hierarchy, from coarse to eddy-rich, allows for a clean and extensive analysis of the dynamics and energetics of the flow as a function of horizontal grid spacing, which will be upcoming. The coarser grid configurations shown here can serve as a test bed for the evaluation of scale-aware mesoscale eddy parameterizations in the “gray” zone of eddy-permitting resolution. Even coarser, non-eddying configurations, not shown here because they need an eddy parameterization to look sensible, will be used to evaluate parameterizations of subgrid mesoscale eddies.

One virtue of the NW2 model is that the adiabatic limit isolates mesoscale eddies from other processes. However, there are strong interactions between mesoscale eddies and the surface mixed layer which have been recognized since early evaluations of mesoscale parameterizations (Danabasoglu et al., 1994). The dynamical core and algorithm used are the same as for a the full primitive-equation general circulation model (Adcroft et al., 2019; Griffies et al., 2020), so adding diabatic processes is relatively straightforward. The NW2 model can be modified and developed further to explore the large-scale overturning circulation (e.g., Wolfe and Cessi, 2009) or the connection with other parameterizations, for example, mixed-layer and sub-mesoscale parameterizations.