ISOM 1.0: a fully mesoscale-resolving idealized Southern Ocean model and the diversity of multiscale eddy interactions

We describe an idealized Southern Ocean model (ISOM 1.0) that contains simplified iconic topographic features in the Southern Ocean and conduct a fully mesoscale-resolving simulation with the horizontal resolution of 2 km, based on the Massachusetts Institute of Technology general circulation model. The model obtains a fully developed and vigorous mesoscale eddying field with a k⁻³ eddy kinetic energy spectrum and captures the topographic effect on stratification and large-scale flow. To make a natural introduction of large eddy simulation (LES) methods to ocean mesoscale parameterization, we propose the concept of mesoscale ocean direct numerical simulation (MODNS). A qualified MODNS dataset should resolve the first baroclinic deformation radius and ensure that the affected scales by the dissipation schemes are sufficiently smaller than the radius. Such datasets can serve as the benchmark for a priori and a posteriori tests of LES schemes or mesoscale ocean large eddy simulation (MOLES) methods in ocean general circulation models. The 2 km simulation can meet the requirement of MODNS and also capture submesoscale effects. Therefore, its output can be a type of MODNS and provide reliable data support for relevant a priori and a posteriori tests. We demonstrate the diversity of multiscale eddy interactions, validate the crucial role of mesoscale-related strain in submesoscale processes, and uncover the bridge effect of submesoscale processes between mesoscale entities and in the eddy–jet interaction. In addition, we use the model to conduct multipassive tracer experiments and reveal guidelines for the initial settings of passive tracers to delay the homogenization process and ensure the mutual independence of tracers over a long period.

1 Introduction

Oceanic mesoscale processes have motions with spatial scales of O(10 km)–O(100 km) or near the baroclinic Rossby deformation radius, including quasi-geostrophic eddies and meandering jets (Chelton et al., 1998; Hallberg, 2013; Siedler et al., 2013; Thompson and Naveira-Garabato, 2014; Youngs et al., 2017). These processes encompass more than 80 % of the oceanic kinetic energy and play a crucial role in material transport, heat transport, momentum budget, and air–sea interactions. They also modulate climate variability across multiple timescales and engage in drastic multiscale interactions with both large-scale and submesoscale processes (Stammer, 1998; Ferrari and Wunsch, 2009; Zhai et al., 2010; Chelton et al., 2011; Dong et al., 2014; Ma et al., 2016; Vallis, 2017; Busecke and Abernathey, 2019; Schubert et al., 2020; Taylor and Thompson, 2023; Yang et al., 2024). Fully resolving oceanic processes at the mesoscale require ocean general circulation models (OGCMs) with kilometer-scale horizontal resolution (Marques et al., 2022). Such models require massive computational and storage resources for long-term integration or large-ensemble experiments. Therefore, we still need parameterizations that capture the collective effects of the unresolved parts related to oceanic mesoscale processes in lower-resolution OGCMs.

The classic works on parameterizing oceanic mesoscale processes include the isoneutral diffusion scheme by Redi (1982) and the revolutionary Gent–McWilliams (GM) scheme (Gent and McWilliams, 1990; Gent et al., 1995) that represents the effect of eddy-induced adiabatic advection and of generating a net sink of available potential energy. These two schemes are widely used in coarse-resolution OGCMs and can be represented by a flux–gradient relationship with an asymmetric transport tensor (Griffies et al., 1998). Since its inception, scholars have made considerable advancements in the field based on the GM–Redi framework. Additional specific constraints or properties, such as the stratification state (Visbeck et al., 1997), anisotropy (Smith and Gent, 2004), geometrical information (Mak et al., 2018; Wei et al., 2024a), and energetic constraints (Cessi, 2007; Eden and Greatbatch, 2008; Jansen and Held, 2014; Mak et al., 2018; Bachman, 2019; Jansen et al., 2019; Yankovsky et al., 2024), have been embedded in the scheme to produce spatiotemporal variations in the transport coefficients.

In addition to the traditional theory-driven schemes described above, other studies have revealed the potential of large eddy simulation (LES) methods in oceanic mesoscale parameterization (Fox-Kemper and Menemenlis, 2008; Graham and Ringler, 2013; Pearson et al., 2017; Khani et al., 2019; Khani and Dawson, 2023; Xie et al., 2023; Perezhogin and Glazunov, 2023). The application is sometimes called the mesoscale ocean large eddy simulation (MOLES) (Fox-Kemper and Menemenlis, 2008; Graham and Ringler, 2013). First, one advantage of LES is in addressing the Reynolds averaging issue (Khani and Dawson, 2023; Xie et al., 2023; Perezhogin and Glazunov, 2023). Many parameterizations and related diagnostics originating from the GM–Redi framework are based on Reynolds averaging, which may simplify the derivation. However, the Reynolds averaging method inherently suppresses cross-scale interactions near the grid scale, leading to a loss of local information, and its mathematical properties are not fully satisfied by the grid discretization of numerical models (Leonard, 1975; Germano et al., 1991; Germano, 1992; Pope, 2001; Xie et al., 2023). Using Reynolds averaging as a grid discretization approximation in very coarse resolution OGCMs might not cause significant issues. However, as the horizontal resolution of OGCMs increases, mesoscale or even submesoscale dynamics with multiscale interactions enter the model grid-scale regime, and local features should be considered when parameterizing subgrid-scale effects. Within the LES framework, subgrid-scale stress and flux terms that include local interactions can be fully expressed, thereby improving the simulation results. Second, the LES framework can also explicitly involve the stationary eddying effect as a supplement to traditional schemes that mainly focus on the transient eddying effect related to instabilities (Khani and Dawson, 2023; Xie et al., 2023). In addition, constructing parameterization schemes that combine LES with machine learning has become a frontier field in developing OGCMs (Bolton and Zanna, 2019; Zanna and Bolton, 2020; Guillaumin and Zanna, 2021; Frezat et al., 2022; Srinivasan et al., 2024).

There are two types of tests for examining the performance of LES models (i.e., parameterization schemes in oceanography): a priori and a posteriori tests (Meneveau, 1994; Moser et al., 2021). In an a priori test, direct numerical simulation (DNS) for a specific flow is required first. Then, we perform the scale separation of DNS data through coarse-graining methods (e.g., spatial filtering). We regard the filtered field as an approximation of the coarser-resolution model output, and we directly diagnose the “true” subgrid-scale terms following their definition. Finally, we reconstruct the subgrid-scale terms by substituting the filtered field into the LES model. By investigating the performance of the reconstructed and the “true” subgrid-scale terms under given metrics (e.g., spatial correlations and the energy transfer rate), we can find out the properties of the LES model or parameterization scheme. In an a posteriori test, we embed the given LES scheme into a lower-resolution numerical model, run simulations, and test parameter sensitivity if needed. At this point, the macroscopic features of the corresponding flow from DNS data serve as the benchmark for evaluating the lower-resolution simulation results when applying the LES model. For both a priori and a posteriori tests, DNS data for the studied flow are essential. Therefore, when introducing LES methods to ocean mesoscale parameterization, it is necessary to generate DNS datasets for mesoscale processes, thereby facilitating the systematic work of developing, testing, and implementing any LES schemes into OGCMs.

In the literature of computational fluid dynamics, DNS requires the numerical model resolution to be at least close to the Kolmogorov scale (Moin and Mahesh, 1998; Pope, 2001; Kaneda and Ishihara, 2006; Alfonsi, 2011). It is the scale at which molecular viscosity becomes important (Pope, 2001; Vallis, 2017). However, the classical definition of DNS is inapplicable in the context of implementing LES methods into OGCMs. Carrying out simulations with a resolution close to the Kolmogorov scale in OGCMs is not feasible in the foreseeable future, and it is also unnecessary to adopt such a high resolution for merely simulating mesoscale motions. For oceanic mesoscale flow, the dynamically indicative scale is the first baroclinic deformation radius. The radius is in the range of approximately 10–40 km in mid-latitude oceans (e.g., the Southern Ocean) (Chelton et al., 1998; LaCasce and Groeskamp, 2020), necessitating a horizontal resolution of at least $\mathrm{1}/\mathrm{30}$° to resolve it explicitly (Hallberg, 2013; Marques et al., 2022). In addition, OGCMs often adopt dissipative schemes near the grid scale to ensure numerical stability. If the scale at which the dissipative scheme plays a significant role cannot be well separated from the first baroclinic deformation radius, then the intrinsic mesoscale dynamics would be contaminated artificially. Some works (e.g., Graham and Ringler, 2013; Radko and Kamenkovich, 2017) use the term DNS in the context of ocean mesoscale dynamics, which might lead to misunderstandings. Therefore, we hereby explicitly propose the concept of mesoscale ocean direct numerical simulation (MODNS). A qualified MODNS dataset not only requires the model grid to explicitly resolve the first baroclinic deformation radius but also demands that the affected scales by the dissipative scheme employed are significantly smaller than the radius, making it the benchmark for a priori and a posteriori tests of LES schemes (or more specifically, MOLES methods) in OGCMs.

To highlight the oceanic mesoscale dynamics in the simulation while reducing computational and storage costs, we develop an idealized Southern Ocean model (ISOM 1.0) and conduct fully eddy-resolving experiments to generate a type of MODNS dataset. Though it is fundamental for the simulation to conform to the realistic Southern Ocean in terms of basic dynamical features (i.e., quantitatively consistent with the observed Antarctic Circumpolar Current (ACC) transport value to provide a reasonable background flow for eddying processes), we emphasize that the focus of the simulation should be on controlling the dynamics of the idealized model rather than on precise comparisons with observations or realistic model results. We hope that the model can describe processes most closely associated with the mesoscale in the Southern Ocean, including mesoscale motions (mesoscale eddies and meandering jets), large-scale background processes (stratification and eastward transport similar to that in the realistic Southern Ocean), eddy–eddy interactions, eddy–jet interactions, large-scale topographic effects, and mesoscale–submesoscale interactions. ISOM enables us to achieve a type of MODNS in an idealized Southern Ocean with topography, thereby providing reliable supporting data for the design, testing, and application of any potential LES-related mesoscale parameterization schemes and the theoretical exploration of the dynamics.

We introduce the design philosophy and implementation of ISOM 1.0. We verify that the oceanic mesoscale regime is fully resolved and barely contaminated by the dissipative schemes in the highest-resolution simulation. Therefore, the simulation can serve as a type of MODNS dataset with vigorous mesoscale activities. We also offer vivid examples of multiscale eddy interactions, especially the eddy–jet interaction, to intuitively demonstrate the capability of MODNS to capture all the mesoscale-related processes. In addition, we conduct multipassive tracer experiments and explore the principles for setting the initial field of passive tracers to offer technical references for relevant works.

2 Model description

2.1 Model equations and configurations

We establish ISOM 1.0 using the Massachusetts Institute of Technology general circulation model (MITgcm; Marshall et al., 1997). We refer to and improve upon the case of the Southern Ocean Reentrant Channel Example in the MITgcm manual (Adcroft et al., 2021) that is closest to our needs, as well as similar idealized works (e.g., Abernathey et al., 2011; Bischoff and Thompson, 2014), ultimately achieving ISOM 1.0 with topography and intermediate complexity. Similar to these works, we consider a hydrostatic, incompressible Boussinesq fluid on the β plane, with an implicitly linearized free surface and a linearized equation of state (only potential temperature, no salinity). We employ the Cartesian coordinate, and then the governing equations (without specifying concrete parameterizations) are as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{D}u}{\mathrm{D}t}}& -fv+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_c}}{\displaystyle \frac{\partial p^{\prime }}{\partial x}}+{\mathrm{\nabla }}_{\mathrm{h}}\cdot (-A_{\mathrm{h}}{\mathrm{\nabla }}_{\mathrm{h}}u)+{\displaystyle \frac{\partial }{\partial z}}\left(-A_{\mathrm{v}}{\displaystyle \frac{\partial u}{\partial z}}\right)\\ & ={\mathcal{F}}_u,\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{D}v}{\mathrm{D}t}}& +fu+{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_c}}{\displaystyle \frac{\partial p^{\prime }}{\partial y}}+{\mathrm{\nabla }}_{\mathrm{h}}\cdot (-A_{\mathrm{h}}{\mathrm{\nabla }}_{\mathrm{h}}v)+{\displaystyle \frac{\partial }{\partial z}}\left(-A_{\mathrm{v}}{\displaystyle \frac{\partial v}{\partial z}}\right)\\ & ={\mathcal{F}}_{\mathrm{v}},\end{array}\text{(3)}& {\displaystyle \frac{\partial \mathit{\eta }}{\partial t}}+{\mathrm{\nabla }}_{\mathrm{h}}\cdot \left(H\widehat{u}\right)=\mathrm{0},\text{(4)}& {\displaystyle \frac{\mathrm{D}\mathit{\theta }}{\mathrm{D}t}}+{\mathrm{\nabla }}_{\mathrm{h}}\cdot (-{\mathit{\kappa }}_{\mathrm{h}}{\mathrm{\nabla }}_{\mathrm{h}}\mathit{\theta })+{\displaystyle \frac{\partial }{\partial z}}\left(-{\mathit{\kappa }}_{\mathrm{v}}{\displaystyle \frac{\partial \mathit{\theta }}{\partial z}}\right)={\mathcal{F}}_{\mathit{\theta }},\text{(5)}& {\displaystyle }{p}^{\prime }=g{\mathit{\rho }}_c\mathit{\eta }+\underset{z}{\overset{\mathrm{0}}{\int }}g{\mathit{\rho }}^{\prime }\mathrm{d}z.\end{array}$$

Here, u and v are the x and y components of the velocity vector, η is the free surface height, θ is the potential temperature, and p^′ is the pressure field. For more details, please refer to the MITgcm manual.

The forcing terms of the horizontal momentum equation ℱ_u and ℱ_v include the steady zonal surface wind stress and the quadratic bottom drag. The wind stress is set as follows:

$$\begin{array}{}\text{(6)}& {\mathit{\tau }}_{\mathrm{s}}\left(y\right)={\mathit{\tau }}_{\mathrm{0}}\sin (\mathit{\pi }y/L_y),\end{array}$$

where τ₀=0.2 N m⁻². The dimensionless coefficient of quadratic bottom drag C_d is 0.01.

The forcing term of the potential temperature equation ℱ_θ includes relaxation to a prescribed surface temperature profile and a sponge layer at the northern side of the domain. The specified sea surface temperature profile increases linearly from 0 °C in the south to 20 °C in the north, with a relaxation timescale of 30 d (except in the sponge layer region). The sponge layer, or a three-dimensional subdomain where restoring boundary conditions are applied, is confined within 160 km at the northern side of the domain, and the potential temperature is relaxed to the following profile:

$$\begin{array}{}\text{(7)}& {\mathit{\theta }}_p(y,z)=\left[{\mathit{\theta }}_{\mathrm{s}}\right(y)-{\mathit{\theta }}_{\mathrm{b}}]\left(e^{z/h_{\mathrm{0}}}-e^{-H/h_{\mathrm{0}}}\right)/\left(\mathrm{1}-e^{-H/h_{\mathrm{0}}}\right).\end{array}$$

θ_s is the prescribed surface temperature profile that varies linearly and meridionally, and θ_b is the bottom temperature set to 0 °C. The depth of the domain H is 4000 m, and the scaling height h₀ is 1200 m (note that the z-coordinate origin is at the surface, and z in the equation is minus). The setup can represent the stratification well on the northern side of the ACC (Abernathey et al., 2011). The relaxation timescale within 80 km at the northernmost part of the sponge layer is 7 d, and from 80 to 160 km it is 14 d.

We use the combination of the horizontal Laplacian viscosity A_h and biharmonic viscosity A₄ for the closure of the fourth term in horizontal momentum equations. We set the background vertical Laplacian viscosity $A_{\mathrm{v}}=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$ m s⁻² for the fifth term in horizontal momentum equations. We also set the background vertical Laplacian diffusivity κ_v and use no horizontal diffusivity κ_h in the potential temperature equation. In addition, we use the nonlocal K-profile parameterization (KPP) scheme (Large et al., 1994) for vertical mixing to generate the mixed layer. Tables 1 and 2 offer the setting values for all the relevant parameters.

Table 1 Basic parameters of the idealized Southern Ocean simulation.

Table 2 Parameters of simulations with different horizontal resolution at their statistical steady state.

2.2 Model bathymetry

The computational domain is a 18 000 km × 3000 km channel with a prescribed topography and zonal periodic boundary conditions (Fig. 1). The domain depth is 4000 m, with 75 vertical levels, and the vertical grid spacing increases from 2 m at the surface to 125 m at the bottom. The channel mimics the hemispheric Southern Ocean that spans from the west of the Drake Passage to the east of the Kerguelen Plateau. This region has a highly complex topography and contains most of the iconic bathymetric features of the Southern Ocean. These features exert profound impacts on the holistic Southern Ocean flow. Using these iconic features enables the idealized model to preserve the complicated topographic effects on oceanic mesoscale processes. Although we could fill the entire domain with realistic topography, it would cause the bathymetric features to be too close to each other and excessively suppress the development of flow. Moreover, reducing computational and storage costs is another reason for simulating only half rather than the entire Southern Ocean.

Figure 1 The bathymetry of (a) the Southern Ocean State Estimate (SOSE; Mazloff et al., 2010a; Verdy and Mazloff, 2017a) and (b) ISOM 1.0. The snapshots of the surface kinetic energy of (c) SOSE and (d) 2 km resolution ISOM simulation.

We adopt five types of topography in the domain (Fig. 1a, b). They are described in detail from left to right as follows.

1. The first type is an idealized Drake Passage that is a passage from x=800 km to x=3200 km. It smoothly narrows from the sides to form a 1600 km × 600 km rectangular subchannel. The four corners of the passage are $\mathrm{1}/\mathrm{4}$ circular arcs with a radius of 400 km. The bottom topography within the passage rises from −4000 to −2500 m within the range of x=800 km to x=2600 km and descends from −2500 m to −4000 m within the range of x=2600 km to x=3200 km.

   The design helps the idealized model qualitatively reproduce the flow characteristics near the realistic Drake Passage despite its high simplification (Fig. 1c, d). In our early design, we tested the passage with a flat bottom or the semicircular ridge similar to Marques et al. (2022) but found that they somehow suppress eddy activity near the passage.

   We also tested the passage's depth and width in the early spinup stage. Deepening the passage (e.g., changing the highest point to −3000 m) and widening the passage (e.g., increasing the north–south span to 750 km) can significantly increase the ACC transport value to over 200 Sv. Different observations and simulations (e.g., Cunningham et al., 2003; Park et al., 2009; Donohue et al., 2016; Xu et al., 2020; Artana et al., 2021a) yield varying values for ACC transport, but generally, they fall within the range of 130 to 180 Sv, which is also consistent with our calculation using Southern Ocean State Estimate (SOSE). To keep the ACC transport of ISOM within the reasonable range, we decided to make the design described above so that the value in the early spinup stage is around 170 Sv. As the simulation continued, the ACC transport slowly decreased and reached an equilibrium state of around 145 Sv, which is still within the reasonable range after 30 model years of integration. The time-averaged ACC transport can remain at the same level after horizontal resolution increases (Figs. 2a and 3a).

   In addition, we also found in the early design (Fig. S3 and Tables S1 and S2 in the Supplement) that when the depth of the passage is too shallow (e.g., −1000 m), ACC transport is highly sensitive to the bottom drag coefficient. Reducing (increasing) the coefficient can significantly increase (decrease) ACC transport. In our current design (i.e., −2500 m), ACC transport is no longer sensitive to the bottom drag coefficient.

2. The second type is a tilted mid-ocean ridge (ridge 1) with zonal and meridional slopes downstream of the idealized Drake Passage with a zonal width w=3500 km and height h₁=800 m. The ridge is confined within a parallelogram with vertices at $(x,y)=(\mathrm{3500},\mathrm{3000}\mathrm{km})$, (7000, 3000 km), (4500, 0 km), and (8000, 0 km). In the zonal cross-section, the ridge is composed of piecewise cubic functions with height h expressed as follows:

   $$\begin{array}{}\text{(8)}& h\left(x\right)=\left\{\begin{array}{ll}-\frac{\mathrm{16}{h}_{\mathrm{1}}}{w^{\mathrm{3}}}{x}^{\mathrm{3}}-\frac{\mathrm{12}{h}_{\mathrm{1}}}{w^{\mathrm{2}}}{x}^{\mathrm{2}}+h_{\mathrm{1}},& -\frac{w}{\mathrm{2}}<x<\mathrm{0},\\ \frac{\mathrm{16}{h}_{\mathrm{1}}}{w^{\mathrm{3}}}{x}^{\mathrm{3}}-\frac{\mathrm{12}{h}_{\mathrm{1}}}{w^{\mathrm{2}}}{x}^{\mathrm{2}}+h_{\mathrm{1}},& \mathrm{0}<x<\frac{w}{\mathrm{2}}.\end{array}\right.\end{array}$$

   The ACC flow intensity in the realistic Atlantic sector is weaker than in the vicinity of the Drake Passage, the Agulhas region, and the Kerguelen Plateau region. To maintain this feature in ISOM, we set a low mid-ocean ridge with an undulation (800 m) weaker than the realistic ridge in the South Atlantic and our early design (1500 m; Fig. S3 in the Supplement). We explain that the idealized mid-ocean ridge is a continuous and complete topography in the north–south orientation, which has a more concentrated stimulation effect on the activity of eddies and jets compared to the realistic ridge with fractures and breakpoints. The low mid-ocean ridge avoids its topographic effect from excessively enhancing eddy activities, thereby avoiding intense downstream flushing and improving the results of the retroflection feature and eddy-shedding process in the Agulhas region.

3. The third type is an idealized African continent. We notice the asymmetry of southern Africa and define the profile with the quadratic function $y\left(x\right)=\frac{\mathrm{3}}{\mathrm{1000}}(x-x_{\mathrm{0}}{)}^{\mathrm{2}}+y_{\mathrm{0}}$ from x=8500 km to x=9000 km and $y\left(x\right)=\frac{\mathrm{3}}{\mathrm{4000}}(x-x_{\mathrm{0}}{)}^{\mathrm{2}}+y_{\mathrm{0}}$ from x=9000 km to x=10 000 km and the vertex (x₀, y₀) = (9000 km, 2250 km). Compared to the early version with a symmetric design (Fig. S3), the introduction of asymmetry enhances the simulation of Agulhas Retroflection feature.

4. The fourth type is another mid-ocean ridge (ridge 2) with opposite tilt direction and a zonal width w=1500 km and height h₁=1000 m. The ridge is confined within a parallelogram with vertices at $(x,y)=(\mathrm{11}\mathrm{750},\mathrm{3000}\mathrm{km})$, (13 250, 3000 km), (8750, 0 km), and (10 250, 0 km). In the zonal cross-section, the expression of height h is the same as Eq. (8).

5. The fifth type is an idealized Kerguelen Plateau. We set a large-scale elliptical Gaussian plateau centered on $(x_m,y_m)=(\mathrm{14}\mathrm{500},\mathrm{1150}\mathrm{km})$, and its height h is expressed as follows:

   $$\begin{array}{}\text{(9)}& h(x,y)=h_{\mathrm{2}}\exp \left[-{\displaystyle \frac{(x-x_m{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_x^{\mathrm{2}}}}-{\displaystyle \frac{(y-y_m{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_y^{\mathrm{2}}}}\right].\end{array}$$

   σ_x is taken as 350 km, σ_y is 250 km, and h₂ is 4500 m. We further limit h₂ to be no higher than 3800 m to form an elliptical flat plateau underwater.

Figure 2 (a) The monthly averaged ACC transport under three wind stress profiles. Panels (b–d) are zonally averaged zonal velocity (shading) and zonally averaged potential temperature (contours) outside the idealized Drake passage under three wind stress profiles. Panels (e)–(g) are as (b)–(d) but inside the idealized Drake passage.

Though the model bathymetry is highly idealized, ISOM qualitatively reproduces the major flow characteristics of the Southern Ocean (Fig. 1c, d), such as intense eddy activity, enhanced currents downstream of large-scale topographic features, and the eastward extension of the current from southern Africa, and successfully simulates energetic mesoscale and submesoscale phenomena and the vivid multiscale eddy–eddy and eddy–jet interactions (see Sect. 3). By observing the differences between the simulated results and the realistic Southern Ocean and referring to other relevant works, we can gain experience and lessons for potential optimization of the model. We share the following insights and hope they can be helpful to any researchers who conduct similar work in the future.

1. Other idealized models of similar complexity, such as NeverWorld in Khani et al. (2019) and NeverWorld-2 in Marques et al. (2022), often apply a steep slope transition from the landmass to the domain depth to mimic the continental shelf. We find in tests that this leads to the formation of intense currents driven by the topographic β effect near the side boundaries. Unfortunately, these currents, especially near the southern boundary, interfere with the simulation of other processes in ISOM. Marques et al. (2022) mentioned their application of significant lateral dissipation. It might be a way to weaken the overly abundant boundary currents. However, the goal of our simulation is to obtain the MODNS dataset. To avoid directly contaminating mesoscale dynamics, we are not inclined to use excessive viscous dissipation in fully eddy-resolving simulations. Therefore, we abandoned the design of continental shelves. Nevertheless, after we completed the simulation, we found that the result of the nearshore Agulhas Current was weaker than in SOSE. Thus, we recommend that future researchers add a moderate topographic slope around the African continent if necessary.

2. The complicated bathymetry near the South American continent has a decisive influence on the surrounding flow, especially the presence of the Falkland Islands (Islas Malvinas) shelf, which leads to a narrow and extensive boundary current (Artana et al., 2021b). ISOM 1.0 prioritizes simulating the Drake Passage. If researchers conduct similar works in the future, we suggest optimizing the topographical expression along the South American continent.

3. Previous studies (e.g., Speich et al., 2006; Lutjeharms and Van Ballegooyen, 1984; Lutjeharms, 2007) have revealed that the topography near the Agulhas Retroflection region, such as the Agulhas Bank (the continental shelf extending southwest from the African continent) and the Agulhas Plateau (a large-scale seamount offshore of southeastern Africa), plays a crucial role in controlling the flow state of eddies and jets. Altering topographic features, such as the slope of the continental shelf and the degree of topographic undulations, substantially influence the flow path, eddy-shedding process, and cross-basin transport. ISOM 1.0 drastically simplifies the bathymetric environment near South Africa. We recommend adding the topographic details around the continent if future works want to improve the simulation of the retroflection characteristics and the eddy-shedding process of the Agulhas rings.

4. Using a sufficiently wide north–south domain span allows for a more reasonable arrangement of the relative positions of topographic features. In the early design (Fig. S3), we used a north–south span of 2400 km, which could not realize a correct relative position between the African continent and the Drake Passage (i.e., if we set the African continent far north, the flow would be directly and undesirably affected by the sponge layer), thus distorting the holistic pattern of flow. The current design, with a north–south span of 3000 km, offers us more space to optimize the arrangement of iconic bathymetric features, dramatically improving the resemblance of ISOM results to the realistic situation. In addition, a sufficiently long east–west span is also conducive to the full development of flow. Researchers may enhance the resemblance by further enlarging the domain and optimizing the topographic arrangement based on this article. However, one has to consider the compromise and balance between the extra computational costs and experimental aims. For example, our early design uses a domain of 14 400 km × 2400 km with 40 vertical levels (Fig. S3), while the current design uses a domain of 18 000 km × 3000 km with 75 levels. The overall computational load increased by nearly 3 times, which exerted a significant burden during the simulation.

2.3 Sensitivity of wind stress profile

In this section, we discuss the sensitivity of the flow and stratification to wind stress settings. From Eq. (6), wind stress consists of two parts: the amplitude τ₀ and the profile form. Many previous studies (e.g., Abernathey and Cessi, 2014; Balwada et al., 2018; Khani et al., 2019; Marques et al., 2022) set τ₀=0.2 N m⁻² consistent with the observation (Chaudhuri et al., 2013; Abernathey et al., 2011), and the effect of amplitude τ₀ has been extensively studied (Tansley and Marshall, 2001; Bischoff and Thompson, 2014; Youngs et al., 2017). Thus, we focus on the effect of the wind stress profile. Idealized simulations often employ the profile related to the sinusoid function. For example, Abernathey et al. (2011) used the standard sinusoid form, while Balwada et al. (2018) adopted the quadratic form in the sinusoid. We tested half-power (${\mathit{\tau }}_{\mathrm{0}}{\sin }^{\mathrm{1}/\mathrm{2}}[\mathit{\pi }y/L_y]$), standard (${\mathit{\tau }}_{\mathrm{0}}\sin [\mathit{\pi }y/L_y]$), and quadratic (${\mathit{\tau }}_{\mathrm{0}}{\sin }^{\mathrm{2}}[\mathit{\pi }y/L_y]$) forms in the sinusoid because the wind stress profile setting for the idealized simulations is unlikely to exceed this range. Otherwise, it would deviate significantly from the observed situation. We integrated 4 model years under three wind stress profiles, starting from the last moment of the 51st model year of the 8 km simulation.

Figure 2a shows the time series of ACC transport under three wind stress profiles. The response of ACC transport to the wind stress profile is quick (compared to the whole spinup process). In less than 1 model year, the ACC transport under the half power of sinusoidal wind stress forcing is prominently stronger than others and basically reaches a new statistical steady state. The higher the power exponent of the sinusoidal function, the smaller the overall ACC transport is. The last 2-year time-averaged ACC transport is 134.6, 141.4, and 159.5 Sv for quadratic, standard, and half-power forms, respectively. The fundamental reason for this result is that the higher the power exponent of the sinusoidal wind stress profile, the sharper it decays to both sides (north and south). When it reaches the latitude (y coordinate) where the Drake Passage is, the magnitude of the wind stress with the high-power form is smaller than that of the low-power form.

We further demonstrate the time-averaged and zonally averaged potential temperature and zonal velocity under three wind stress profiles (Fig. 2b–g). Since the northern boundary geostrophic flow of the Drake Passage is very energetic with a steep isothermal slope, if we zonally average the entire channel, a discontinuity appears on the figure. But in reality, the fields are continuous with large meridional gradients. Therefore, to avoid the misleading visual illusion and to show in detail the response of flow and stratification to the wind stress profile, we divide the entire channel into two parts for zonal average. One is the rectangular subchannel within the Drake Passage with $x\in [\mathrm{1200},\mathrm{2800}\mathrm{km}]$ and $y\in [\mathrm{400},\mathrm{600}\mathrm{km}]$. The other is the rest of the domain.

Outside the Drake Passage (Fig. 2b, c, d), as the power exponent of the sinusoidal wind stress increases, the zonal flow structure becomes more compact. An intense jet core condenses around the domain center under the quadratic form of forcing (Fig. 2d). In the passage (Fig. 2e, f, g), the northern boundary flow is more vigorous than the southern boundary flow. As the power exponent of the sinusoidal wind stress decreases, the boundary flows become more energetic and penetrate more deeply (Fig. 2e).

In summary, if one wants to study oceanic processes on a concentrated zonal jet like Balwada et al. (2018), a high-power sinusoidal wind stress profile would be a good choice. If one wants to generate a holistically stronger ACC (especially within the Drake Passage), we recommend the profile with a low-power sinusoidal function. We choose the standard form in this article to maintain consistency with previous relevant studies and control the ACC transport value in the reasonable range.

Figure 3 (a) The time series of ACC transport spans 1500 model days. The results are monthly averages for the first 3 model years, 3 d averages for the fourth model year, and daily snapshots for the last 60 model days. The blue line for the 2 km simulation starts from the fourth model year. The red line for the 4 km simulation starts from the 13th model year. The black line for the 8 km simulation starts from the 48th model year. The ACC transport in the 8 km simulation reached a statistically steady state in the 48th model year. After refining the horizontal resolution, the ACC transport can also maintain a statistically steady state close to the result of the 8 km simulation. Panels (b)–(d) are for the Rossby number (see Sect. 3.3 for detailed definition) in the idealized Agulhas region for simulations with different horizontal resolutions. Panels (e)–(g) are as (b)–(d) but for the Gaussian plateau region.

2.4 Implementation of MODNS

The fundamental requirement for MODNS is that the horizontal resolution can explicitly resolve the first baroclinic deformation radius. Under the model configuration, it is basically greater than 15 km. If we follow the experience that the model discretization can adequately represent processes exceeding 5 times the grid spacing, then the horizontal resolution required for MODNS should reach at least 3 km. Additionally, Marques et al. (2022) found that when the horizontal resolution reaches $\mathrm{1}/\mathrm{32}$° in NeverWorld-2, the mesoscale model performance converges, meaning that the resolution is sufficient to fully resolve mesoscale phenomena in their model. Based on these considerations, we conduct a simulation with an even finer horizontal resolution of 2 km to achieve a certain type of MODNS.

Directly running the 2 km simulation (e.g., spinup from rest) is costly. Therefore, our spinup strategy is as follows:

1. Integrate an 8 km simulation from rest for 51 model years to reach the quasi-equilibrium state.

2. Interpolate the final output of the 8 km simulation as the initial field for the 4 km simulation and integrate for 16 model years to reach the corresponding quasi-equilibrium state.

3. Interpolate the final output of the 4 km simulation as the initial value field for the 2 km simulation and integrate for several model years.

We employ the third-order direct space–time (DST) flux limiter scheme (the MITgcm advection scheme code is 33) for all simulations. Since the 2 km simulation has already been submesoscale-permitting, we retain the 3 d average output from the seventh model year and daily instantaneous fields from the first 2 months in the eighth model year in the 2 km simulation to demonstrate the model performance for eddying processes of different scales, including mesoscale and submesoscale.

Figure 3 displays the ACC transport series and the Rossby number (the vertical component of relative vorticity normalized by local Coriolis parameter; see Sect. 3.3 for detailed definition) snapshots at three horizontal resolutions. The ACC transport reached quasi-steady state after integrating the 8 km simulation into 48 model years (Fig. 3a). After improving the horizontal resolution, in the last few years of the 4 and 2 km simulations, the ACC transport stayed in equilibrium with no substantial change compared to the results at 8 km. The magnitude of the Rossby number considerably increases with finer horizontal resolution (Fig. 3b–g). This indicates that the model can explicitly describe more submesoscale phenomena, especially those related to meandering jets, filamentation, and multiscale interactions. We will discuss these in detail in Sect. 3.3.

Figure 4 The meridional cross-sections of 3-year-averaged potential temperature (contours) and zonal velocity (shading) of the 2 km simulation. Panels (a–c) are for the western ridge from upstream to downstream, (d–f) for the Gaussian plateau, and (g–i) for the eastern ridge.

3 Results

In this section, we present the stratification, jets, kinetic energy and enstrophy spectra, and the eddying field from the 2 km simulation to examine the model performance for processes relevant to oceanic mesoscale dynamics.

3.1 Stratification and jets

We show the time-averaged temperature and zonal velocity component on the meridional cross-sections from upstream to downstream of the two mid-ocean ridges and Gaussian plateau (Fig. 4) and demonstrate that the model can capture the impact of large-scale topography on the stratification and zonal flow state.

Upstream of ridge 1 (the western ridge), the zonal flow exhibits a multi-branched state, accompanied by gentle isotherm slopes (Fig. 4a). The flow drifts to the north on the ridge and gets intensified (Fig. 4b). Then, the flow drifts to the south further downstream and becomes more compact (Fig. 4c) with vigorous eddy activities (Fig. 1d) and potentially high eddy diffusivity (Abernathey and Cessi, 2014), forming a distinct jet core corresponding to steeper isotherm slopes under the geostrophic constraint.

As for the flow in the Gaussian plateau region, multiple branches upstream of the topography exist (Fig. 4d), with the jet core at y=1900 km. When the flow passes the plateau (Fig. 4e), it is steered along the northern side and forms a single energetic jet core with a steep isotherm slope. Downstream of the plateau (Fig. 4f), the zonal flow scatters into several branches again, and the jet core moves further south to y=1300 km, where the eddy kinetic energy (Fig. 1d) and mixing induced by stationary eddy would be drastic (Bischoff and Thompson, 2014; Xie et al., 2023).

Ridge 2 (the eastern ridge) has a more complicated bathymetric environment. Its upstream section (Fig. 4g) passes through the Agulhas Retroflection zone. Thus, the northern eastward jet reflects the Agulhas Retroflection current. When crossing the ridge (Fig. 4h), the strong topographic slope remarkably strengthens the jets. Due to its narrow spatial scale, the jets do not completely merge when crossing topography and downstream (Fig. 4i).

The above phenomena are qualitatively consistent with previous studies that used two-layer quasi-geostrophic models with bottom slopes or topography (Tansley and Marshall, 2001; Thompson, 2010; Chen et al., 2015; Khatri and Berloff, 2018), idealized mid-ocean ridge experiments (Abernathey and Cessi, 2014; Youngs et al., 2017), idealized Gaussian plateau experiments (Bischoff and Thompson, 2014), realistic Southern Ocean topography simulations (Thompson and Naveira-Garabato, 2014), Southern Ocean reanalysis data (Lu and Speer, 2010; Abernathey and Cessi, 2014), laboratory experiments (Rhines, 2006), and observational data (Orsi et al., 1995; Thompson and Sallée, 2012; Chapman et al., 2020). This indicates that ISOM 1.0 is capable of describing the large-scale background processes that are closely associated with mesoscale phenomena in the Southern Ocean.

Figure 5 The averaged spectral density function of surface EKE and enstrophy (relative vorticity variance) for simulations with different resolutions. Panels (a)–(c) show the surface EKE spectrum using the 3 d averaged output of 1 model year. Panels (d)–(f) show surface EKE spectrum using daily snapshots of 60 model days. Panels (g)–(i) show the surface enstrophy spectrum using daily snapshots of 60 model days. Lines with different colors represent different sampling positions or reference slopes. The eddy velocity is defined by subtracting the time and zonal average of a 1024 km zonal segment of meridional velocity in each position.

3.2 EKE and enstrophy spectrum

To better extract eddy signals, we take 1024 km zonal segments at given locations and compute the eddy kinetic energy (EKE) spectrum from the meridional component of eddy velocity, similar to Marques et al. (2022). We define the eddy velocity by subtracting the temporal and zonal segment-mean velocity. Figure 5a–c show the EKE spectrum of 3 d averaged model output from 1 model year with three horizontal resolutions. The high-frequency processes (e.g., submesoscale) are thus filtered out, and the output highlights the model performance in simulating mesoscale processes. In regions with active mesoscale eddies, including the Drake Passage, downstream of the two mid-ocean ridges, the Agulhas region, and downstream of the Gaussian plateau, the results consistently exhibit a −3 spectral slope for most of the spectrum of all simulations. As the resolution improves, the range of −3 spectral slope expands. For the 2 km simulation (Fig. 5c), the range of the −3 spectral slope (>8 km) completely covers the deformation radius scale (>15 km), which also means that the dissipative schemes barely contaminate the modeled mesoscale process. This result suggests that the model can fully describe oceanic mesoscale motions that can be deemed quasi-geostrophic turbulence and theoretically possesses a k⁻³ (k is wavenumber) kinetic energy spectrum (Charney, 1971; Fu and Morrow, 2013; Vallis, 2017).

Figure 5d–f display the EKE spectrum computed from the daily instantaneous fields of 60 model days of simulations with three resolutions. The 8 and 4 km simulations maintain a spectral slope −3 related to mesoscale processes over most scale ranges (Fig. 5d, e). The 2 km simulation (Fig. 5f) shows a slight uplift at scales below 40 km, resulting in a spectral slope weaker than −3. The spatial distribution and statistics of sea surface variables (Fig. 3 and Tables 3 and 4) confirm that the 2 km simulation models remarkably more submesoscale signals compared to the lower-resolution simulations. Thus, the flattening of the slope is likely related to the burst of submesoscale phenomena. However, classical submesoscale theories, such as the surface quasi-geostrophic dynamics (Blumen, 1978; Held et al., 1995; Lapeyre, 2017) and the surface quasi-geostrophic dynamics with ageostrophic advection (Boyd, 1992; Callies and Ferrari, 2013), predict a $-\mathrm{5}/\mathrm{3}$ and −2 kinetic energy spectrum, respectively. The flattening degree of the EKE spectral slope in the 2 km simulation has not reached a slope close to −2, which means that the current ISOM cannot sufficiently describe submesoscale processes.

In addition, in our early design (Tables S1 and S2 in the Supplement), ISOM did not use the KPP scheme, which resulted in underdeveloped mixed layers and highly intensified surface flow. Interestingly, it obtains the surface EKE spectrum that exhibits distinct submesoscale features at a scale of 8–30 km with a $-\mathrm{5}/\mathrm{3}$ spectral slope (Fig. S5 in the Supplement) consistent with the surface quasi-geostrophic theory. Incorporating the KPP scheme promotes vertical mixing in the upper layer and weakens the possibly unrealistic submesoscale intensification phenomenon in the surface layer. Therefore, if researchers hope to improve the simulation of submesoscale based on the 2 km resolution ISOM in the future, we suggest introducing effective submesoscale parameterization schemes, such as that of Fox-Kemper and Ferrari (2008). Otherwise, one needs to adopt a finer horizontal resolution of at least 1 km, which can also effectively enhance the model performance of the submesoscale motions, especially the related instabilities (Wei et al., 2024b). Nevertheless, we emphasize again that since the goal of this work is the oceanic mesoscale part, the presented results are sufficient to demonstrate that the model can fully describe mesoscale processes and generate a type of MODNS dataset.

Similarly, we take 1024 km zonal segments at given locations and compute the surface enstrophy (or relative vorticity variance) spectrum Ens(k) (Fig. 5g–i). It shows a spectral slope of −1 for all simulations on large scales, and the dissipation effect removes enstrophy on small scales. The result is qualitatively consistent with Chassignet and Xu (2017) (their Fig. 23), who studied the high-resolution model on the Gulf Stream. The scale range with a −1 spectral slope expands with finer horizontal resolution. In the 2 km simulation (Fig. 5i), the range of −1 spectral slope can extend to nearly 10 km, which fully covers the deformation radius scale. Further examination of the potential enstrophy conversion term between the large-scale flow and eddies in the enstrophy budget shows a holistically forward potential enstrophy cascade (Fig. S9). The result further confirms that the 2 km resolution ISOM can generate a type of MODNS dataset.

Table 3 Statistics (maximum values and time-averaged spatial coverage) of the surface Rossby number for three resolutions in the Agulhas (Fig. 6) and Gaussian plateau (Fig. 7) regions for snapshots of 60 model days.

Table 4 Same as Table 3 but for the normalized strain rate.

3.3 Examples of multiscale eddy interaction

In this section, we examine the performance of the 2 km ISOM in simulating the evolution of eddying processes. We choose the Agulhas and Gaussian plateau regions to show snapshots of the sea surface variables. We especially discuss the interaction between the stationary eddy and the meandering jet in the Gaussian plateau region.

Figures 6 and 7 show sea surface snapshots of the 2 km ISOM in the Agulhas and Gaussian plateau regions, respectively. The variables include sea surface temperature (SST), free sea surface height (SSH), kinetic energy, the normalized relative vorticity, the magnitude of horizontal temperature gradient, and the normalized strain rate. The normalized relative vorticity ($\mathit{\zeta }/f$) is the vertical component of relative vorticity ($\mathit{\zeta }=v_x-u_y$) divided by the local Coriolis parameter ($f=f_{\mathrm{0}}+\mathit{\beta }y$). The Rossby number (Ro) is often defined as the absolute normalized relative vorticity ($|\mathit{\zeta }/f|$), and Ro∼O(1) refers to active submesoscale processes (Thomas et al., 2008; Schubert et al., 2020). Since both forms can describe the richness of submesoscale activities and the version without the absolute value can also reflect the sign of vorticity, we refer to the normalized relative vorticity as the Rossby number in the article for convenience. The normalized strain rate is defined as $\sqrt{{\left(u_x-v_y\right)}^{\mathrm{2}}+{\left(u_y+v_x\right)}^{\mathrm{2}}}/\left|f\right|$.

Tables 3 and 4 present the statistics of the Rossby number and normalized strain rate under three horizontal resolutions in the Agulhas and Gaussian plateau regions. As the horizontal resolution improves, both variables undergo remarkable enhancement. Compared with the 8 km simulation, the maximum values of Ro and normalized strain rate in the 2 km simulation increased by 4.1 times in the Agulhas region and by 5.6 times and 6.1 times in the Gaussian region. In addition, the time-averaged spatial coverage of large Ro or strain rate (>0.5) increases by tens of times from 8 to 2 km simulation. The spatial coverage of large values in the Gaussian plateau region is lower than that in the Agulhas region at all resolutions (because of the higher latitude with a larger $\left|f\right|$ in the denominator) but has more prominent growth with resolution. These results indicate that although the current ISOM cannot sufficiently resolve submesoscale processes according to the spectral analysis in Sect. 3.2, there is still a significant improvement trend with finer resolution.

The following sections will further demonstrate the powerful capability of the 2 km resolution ISOM 1.0 to express mesoscale-related processes and thus become a qualified mesoscale ocean DNS dataset.

Figure 6 Surface snapshot of day 23 in the eighth model year of the 2 km simulation in the Agulhas region. (a) Potential temperature, (b) sea surface height, (c) kinetic energy, (d) Rossby number (normalized relative vorticity), (e) the magnitude of the horizontal temperature gradient, and (f) normalized strain rate.

3.3.1 Agulhas region

Though the Agulhas region (Fig. 6) is experiencing complicated processes, including mesoscale eddies, meandering jets, submesoscale processes, and their multiscale interactions, our model successfully simulates the Agulhas Retroflection (linked to the topographic steered large mesoscale eddies) and the eddy-shedding process of generating Agulhas rings (e.g., at x=8400 km, y=2500 km). These characteristics are consistent with Schubert et al. (2020), who used a regional model with $\mathrm{1}/\mathrm{60}$° horizontal resolution. The warm and cold water masses have extensive meridional displacement and drastically communicate with each other, accompanied by multiple mesoscale eddies and vigorous jets (Fig. 6a, c). The existence of the topography makes the southward intrusion of warm water on the eastern side more prominent, which is consistent with the realistic situation. The warm eddies bear a large size, and their flow matches the SSH field (Fig. 6b, c) but with a loose internal structure and moderate Ro in the core (Fig. 6d). The cold eddies have small scales but compact cores that often reach $\left|\mathit{Ro}\right|\sim O\left(\mathrm{1}\right)$. The temperature field between mesoscale eddies is severely stretched and deformed (Fig. 6e, f) to generate multiple fronts and filaments (Gula et al., 2014; McWilliams, 2016).

Interestingly, our simulation captures the interaction among two large warm eddies and a small cold eddy adjacent to the land at x=9100 km and y=2100 km. These two warm eddies exhibit complete and independent cores in the SSH field (Fig. 6a, b) and can be regarded as typical geostrophic eddies. Albeit squeezed by the large eddies, the small cold eddy keeps a compact internal structure with large Ro in the core (Fig. 6d) and locally enhanced flow (Fig. 6c). Strong, narrow temperature fronts exist surrounding the warm eddies (Fig. 6e), with a width of less than 40 km and the most vigorous surface kinetic energy in the region (Fig. 6c). The edges of mesoscale eddies have an intense strain rate (Fig. 6f) and large Rossby numbers (Fig. 6d). The spatial pattern shows a hint of submesoscale-related instabilities with alternating signs and comma-like structures. However, a horizontal resolution of 2 km (as well as without any specialized submesoscale parameterization) does not seem sufficient to clearly resolve all the submesoscale instabilities. Even with the $\mathrm{1}/\mathrm{60}$° horizontal resolution adopted by Schubert et al. (2020), the development of submesoscale instabilities is far from complete, compared to Gula et al. (2014) and McWilliams (2016) using regional models with 100 m level horizontal resolutions. Recently, a set of global simulations (Wei et al., 2024b) based on the LASG/IAP Climate system Ocean Model of version 3 (LICOM3) showed that when the horizontal resolution reaches 0.01°, the model's capability is explosively enhanced to express the submesoscale instabilities and the “submesoscale soup” (McWilliams, 2016) between mesoscale entities. The submesoscale level of the 2 km simulation in the article basically conforms to the 0.02° LICOM3 and the $\mathrm{1}/\mathrm{48}$° LLC4320 MITgcm simulation (Su et al., 2018).

Figure 7 Same as Fig. 6 but for the Gaussian plateau region.

3.3.2 Gaussian plateau region

Figure 7 shows sea surface snapshots of the 2 km ISOM in the Gaussian plateau region. We can see the topographic imprint, the meandering jet, and the large eddies downstream. The meandering jet is reinforced and steered towards the southeast by the topographic slope and generates a highly chaotic region with energetic eddy activities and potentially elevated mixing. The warm water masses are carried southward by jets and eddies, and the cold water is drastically transported towards the plateau from the south (Fig. 7a), which is consistent with the realistic situation near the Kerguelen Plateau. The SSH field (Fig. 7b) implies the existence of energetic geostrophic currents (Fig. 7c) but is too smooth to reflect the diversity of multiscale eddy interactions that are vividly shown in the Ro field (Fig. 7d). As the jet shoots downstream, intensive strain emerges and drives rich submesoscale filaments along the track. However, the spatial pattern of the strain rate and temperature gradient near the topographically reinforced jet differs from that around large mesoscale eddies. Multiple hotspots and chains with locally excited strain rates and temperature gradients interweave the holistic network-like pattern within the meandering jet (Fig. 7e, f). Rings or elongated lines distinct from the surrounding environment exist around the mesoscale eddy in the Agulhas and Gaussian plateau regions.

Our model focuses on oceanic mesoscale processes (eddies and meandering jets). The mesoscale eddies can be categorized as transient eddies with shorter life cycles and stationary eddies that are long-term and often related to topography (Bischoff and Thompson, 2014; Lu et al., 2016; Khani et al., 2019; Xie et al., 2023). Both of them substantially contribute to the eddy transport process. All types of mesoscale processes are shown vividly in Fig. 7, especially the stationary eddy downstream of the plateau (at x=15 500 km, y=900 km). We will zoom in on this area in Sect. 3.3.3 and demonstrate the evolution details of the flow to explore the complex multiscale interactions between stationary eddies, transient eddies, meandering jets, and submesoscale phenomena.

Figure 8 The evolution of stationary eddies, transient eddies, the meandering jet, and their interactions from the 2 km simulation in a zoomed-in region downstream of the Gaussian plateau.

3.3.3 Eddy–jet interaction

We now zoom in on a region east of the plateau in Fig. 7 and demonstrate the eddy–jet interaction using the model output in the 2 km simulation. In the period displayed, the jet is in a stage of drastic meridional movement, making it easy for eddy formation and intimate eddy–jet interactions.

On day 37, there exists an energetic jet along the temperature front with large meridional displacement, an individual warm mesoscale eddy with a diameter of over 350 km on the southern side, and a cold eddy on the northern side of the jet (Fig. 8a, f, k). From the Rossby number and normalized strain rate (Fig. 8p, u), the jet had a network-like pattern mentioned in Sect. 3.3.2 upstream of $(x,y)=(\mathrm{15}\mathrm{900},\mathrm{1000}\mathrm{km})$, where submesoscale signals are locally excited. However, downstream of this point, the pattern becomes a much simpler curved pattern similar to a mesoscale eddy, and another mesoscale eddy soon forms here later (Fig. 8f–j). Concerning the large eddy southwest to the jet, it has a prominent core on the SSH field (Fig. 8f) with vigorous circulation (Fig. 8k) that tightly locks in a bulk of warm water mass (Fig. 8a). It forms a ring with a high strain rate and Ro in the periphery (Fig. 8p, u), but it is relatively hollow inside. Though this eddy may sometimes couple with the meandering jet during its lifespan, it can survive for over 3 months with an almost unchanged geographical location (Fig. 8k–o and Fig. S1), making it a typical stationary eddy. The cold eddy north of the jet is a transient eddy prone to environmental flow due to its limited spatial scale and circulation intensity. It struggles to maintain its structure and move westward (Fig. 8p–t).

On day 40, the jet stream cannot maintain such a large meridional displacement and starts to shift to the zonal mode (Fig. 8l). A warm eddy begins to separate from the jet stream (Fig. 8b, g). The connection part becomes narrow and filamentary in the SST field and owns enhanced Ro and strain rate fields (Fig. 8q, v). The stationary eddy upstream is still intact and almost becomes a perfect ring.

On day 49, the jet stream finishes throwing out the eddy on the right and moves close to the stationary eddy on the left (Fig. 8c, h, m). They trigger a drastic strain field and intense submesoscale signals in between but with undermined fields for the rest of the ring (Fig. 8r, w).

On day 55, an interesting phenomenon happens. The stationary eddy is temporarily open, and the meandering jet instantly generates a remarkable meridional displacement through the eddy circulation (Fig. 8n), shooting the network-like structure of Ro and strain rate field towards the west side of the eddy (Fig. 8s, x). However, the state is unstable and unsustainable because the submesoscale network can maintain the linkage between the upstream and downstream parts of the jet, thereby helping to stitch and rebuild the zonal geostrophic current and push the stationary eddy disconnection with the jet stream several days later (Fig. 8j, t, y).

Once the western stationary eddy forms and detaches from the jet stream, its geographical location barely changes for over 3 months (Fig. S1 in the Supplement). It can maintain a complete eddy-like pattern most of the time. Even though the jet stream sometimes borrows the eddy circulation, the main body of the stationary eddy does not suffer a devastating consequence. The eastern warm eddy is also generated from the excessive meridional variation of the jet and has a broader range of motion. Further downstream, as the jet stream bends and breaks, another warm eddy with the motion even more wide-ranging is also easy to form. These warm eddies and northern smaller cold eddies can form a spectacular asymmetric eddy street with a submesoscale network on the upstream part of the jet and strain-related submesoscale activities around and between mesoscale eddies (e.g., snapshot shown in Fig. S2).

To sum up, the 2 km resolution ISOM simulation demonstrates reliable capability in expressing mesoscale-related processes and thus is qualified to become a type of MODNS. We obtain qualitatively consistent features with the realistic situation and comprehensively illustrate the diversity of mesoscale processes and multiscale interactions. Previous works (e.g., Capó et al., 2021; Gula et al., 2022) reveal that the mesoscale-related strain is a crucial mechanism for submesoscale motions. We validate this point and intuitively uncover that submesoscale motions can bridge mesoscale entities and affect the eddy–jet interaction. Therefore, an appropriate representation of submesoscale effects on the mesoscale (i.e., MOLES schemes) in eddy-rich or eddy-permitting OGCMs (e.g., 0.1°) is crucial for improving the simulations and is cost-effective compared to using kilometer-scale resolutions or even finer horizontal resolutions. ISOM 1.0 provides a platform for testing such parameterizations. The MODNS data based on the 2 km simulation can also serve as a benchmark for a prior and a posterior tests, thereby guiding modelers to employ suitable schemes.

4 Multipassive tracer tests

We use ISOM 1.0 to conduct multipassive tracer tests in this section. Multiple passive tracers can be leveraged as samples to construct an overdetermined linear system of equations to estimate the eddy transport tensor (Bachman et al., 2015, 2020). A two-dimensional diagnosis (i.e., on the neutral surface) requires at least two nonparallel samples, and a three-dimensional diagnosis (i.e., on the z coordinate) needs at least three nonparallel samples (Xie et al., 2023). We conduct online tests of passive tracers using the 8 km simulation of the 52nd to 55th model years (the flow is the same for all tracers). We uncover the properties of passive tracers in ISOM 1.0 and propose guidelines for selecting passive tracer combinations, thereby providing technical references for works that employ relevant methods for eddy transport diagnosis and parameterization design.

All passive tracers in this article obey the advection equation without diffusion and source–sink terms (i.e., no restoration to a prescribed profile):

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\partial C_i}{\partial t}}+\mathbf{u}\cdot \mathrm{\nabla }{C}_i=\mathrm{0}.\end{array}$$

Figure 9 The initial fields for the recommended passive tracer combination.

Figure 10 Snapshot of day 360 in the multipassive tracer experiment.

Figure 11 Snapshot of day 1440 in the multipassive tracer experiment.

The only difference among tracers lies in their initial fields. To ensure that the tracers are mutually independent (i.e., have very low spatial correlation), we recommend the initialization of a combination of four passive tracers as follows (Fig. 9):

$$\begin{array}{}\text{(11)}& {\displaystyle }{C}_{\mathrm{1}}(x,y,z,t_{\mathrm{0}})=y/L_y+{\mathit{\sigma }}_{\mathrm{1}}(x,y,z)\text{(12)}& {\displaystyle }{C}_{\mathrm{2}}(x,y,z,t_{\mathrm{0}})=\sin (\mathit{\pi }y/L_y)+{\mathit{\sigma }}_{\mathrm{2}}(x,y,z)\text{(13)}& {\displaystyle }{C}_{\mathrm{3}}(x,y,z,t_{\mathrm{0}})=\sin (\mathit{\pi }x/L_x)+{\mathit{\sigma }}_{\mathrm{3}}(x,y,z)\text{(14)}& {\displaystyle }{C}_{\mathrm{4}}(x,y,z,t_{\mathrm{0}})=\left|\sin \right(\mathrm{2}\mathit{\pi }x/L_x+\mathit{\pi }/\mathrm{4}\left)\right|+{\mathit{\sigma }}_{\mathrm{4}}(x,y,z).\end{array}$$

The random terms σ_i ($i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}$) follow a uniform distribution between 0 and 0.1. We further limit the initial fields within the range of 0 to 1. The setup ensures that the absolute Pearson correlation coefficients between tracers are far less than 0.1, and the initial fields can be considered mutually independent. Their spatial standard deviations are at the same level (approximately 0.3). This allows each tracer to contribute almost equally when solving the overdetermined linear system related to the transport tensor.

After long-term stirring by the flow, the homogenization of individual passive tracers occurs, and the spatial pattern among tracers becomes correlated (Figs. 10 and 11). This leads to the disappearance of the tracer gradient and the local alignment of the eddy flux vectors among tracers. Both effects can cause the failure of the multipassive tracer method for diagnosing transport tensors based on the flux–gradient relationship. There are two approaches to addressing these issues. The first approach is to add a restoration term to relax tracers to prescribed profiles. When diagnosing the transport tensor, one must address the problems caused by the restoration term (Bachman et al., 2015, 2020; Haigh and Berloff, 2021). The second approach is to continue to adopt passive stirring but release tracers several times (Wei and Wang, 2021). We prefer the latter approach and conduct validation tests. The key is to explore the timescale of homogenization and the correlation of the passive tracer combination. As long as the time interval of the release is shorter than the timescale at which unacceptable homogenization and correlation occur, the tracer output is suitable for subsequent diagnosis.

Figure 12 The temporal evolution of the standard deviation of tracers (solid line) and the absolute spatial correlation between tracers (dashed line) at different depths.

We test various initial tracer fields to explore how to delay the homogenization process and ensure the mutual independence of the tracers. We share noteworthy guidelines as follows.

1. If the flow is quasi-two-dimensional (e.g., large-scale oceanic motion), the initial tracer field should vary horizontally. One should avoid setting more than one initial field that changes only in the z direction, such as $C\left(z\right)=\left|z\right|/H$ and $C\left(z\right)=\sin \left(\mathit{\pi }\right|z|/H)$. Although these two expressions appear to be independent, they are highly correlated in the horizontal plane.

2. If the domain has periodic boundaries, discontinuities in the initial fields at the corresponding boundaries should be avoided. For example, a passive tracer with an initial field of $C\left(x\right)=x/L_x$ does not seem to cause significant numerical issues under high-order advection schemes, and there is no guarantee that such an initial setting would not be problematic when using lower-order schemes or spectral methods.

3. A crucial principle is that the initial tracer fields should be dominated by structures with large spatial scales. Otherwise, the tracer fields would be rapidly mixed by flow stirring. An extreme example is completely random initialization, which fills the tracer field with small-scale noise and can be homogenized within days. Another example is the initial field with two meridional half-period harmonic waves, $C\left(y\right)=\left|\sin \right(\mathrm{2}\mathit{\pi }y/L_y\left)\right|$. Though its spatial scale is far larger than that of random noise, it can be significantly mixed over a few months. Among the four tracers recommended in the paper, tracer C₄, which has an initial field with two zonal half-period harmonic waves, homogenizes the fastest in the upper layer and reduces its spatial standard deviation by about 40 % in the first model year (the solid purple line in Fig. 12). Tracer C₂, with an initial field of a meridional half-period wave, homogenizes the fastest in the deep level and reduces its spatial standard deviation by about 50 % within 2 model years (the solid red line in Fig. 12). The homogenization rates of tracers C₁ and C₃, which have larger spatial scales, are significantly lower (the solid blue and green lines in Fig. 12).

By further analyzing the temporal evolution of the absolute Pearson correlation between the four passive tracers (dashed lines in Fig. 12), we find the following:

1. The spatial correlation between tracers generally increases over time (though the evolution is not vigorously monotonic). This indicates that their spatial patterns tend to be similar under the continuous stirring of the same flow.

2. The evolution of correlations at different depths varies. Deeper locations generally have a slower increase in absolute correlation (though the monotonicity is not vigorously guaranteed). The increase rate in the upper mixed layer is much faster than in the deeper levels.

3. For convenience, we deem the absolute spatial correlation of less than 0.2 as low correlation. Within 360 model days, at least four tracer pairs for all levels and three tracers for levels under 200 m can maintain low correlation, and the spatial distribution of each tracer has not yet visually shown homogenization (Fig. 10). By the 720th model day, at least four tracer pairs can maintain low correlation for levels under the upper mixed layer.

Therefore, under the condition of ISOM 1.0, controlling the duration of tracer release to less than 2 model years enables homogenization and correlation to be limited to an acceptable level so that the multiple tracer combination can serve as samples for accurately estimating the transport tensor.

5 Conclusions

In this paper, we have introduced an idealized Southern Ocean model (ISOM 1.0) that contains a simplified version of iconic topographic features in the Southern Ocean. We have conducted a fully eddy-resolving (2 km) simulation. The prominent feature of the model is the successful simulation of a fully developed and vigorous mesoscale eddying field. We reproduce the EKE spectrum of k⁻³ predicted by geostrophic turbulence theory. In addition, the simulated ACC transport and geographical distribution of eddy activities are qualitatively consistent with the realistic situation. The model can also describe the topographic effect on stratification and large-scale flow and show reliable capability of simulating mesoscale processes.

To facilitate a smooth introduction of LES methods to ocean mesoscale parameterization, we propose the concept of MODNS. Its model grid should explicitly resolve the first baroclinic deformation radius, and the scales where dissipative schemes play a significant role are distant from the mesoscale dynamical regime, making it the benchmark for a priori and a posteriori tests of LES models or MOLES schemes in OGCMs. The 2 km idealized simulation satisfies the demands for MODNS and captures the submesoscale effects on the mesoscale entities. Therefore, it can serve as a type of MODNS and offer reliable data support for conducting relevant a priori and a posteriori tests of mesoscale-related parameterization schemes.

We demonstrate the diversity and high complexity of multiscale eddy interactions related to mesoscale processes by examining the evolution of mesoscale eddies, submesoscale phenomena, and meandering jets. We validate that mesoscale-related strain plays a crucial role in submesoscale processes. We also uncover that submesoscale motions can exert a bridge effect between mesoscale entities and affect eddy–jet interactions. Therefore, we expect an appropriate expression of the collective submesoscale effects on the mesoscale part (e.g., MOLES schemes) in eddy-rich or eddy-permitting OGCMs to benefit the simulation of mesoscale variability.

We also use the idealized model to conduct multipassive tracer experiments. We reveal some guidelines for the initialization settings of passive tracers. We discover a combination of four passive tracers that can delay the homogenization process and ensure the mutual independence of tracers for a long time. With this combination and in ISOM 1.0, controlling the duration of each experiment of tracer release to less than 2 model years can ensure that the spatial standard deviation of the tracers and the correlations among the tracers are limited to acceptable levels. This allows the results from multiple release experiments of the tracer combination to form a qualified sample for solving the eddy transport tensor based on the flux–gradient relationship. The tracer initialization guidelines can offer important references for similar ocean simulations.

The global ocean is a sophisticated system incorporating multiple spatiotemporal scales, and the large-scale oceanic dynamics vary across different regions. The idealized model in this paper provides a type of MODNS of the Southern Ocean with simplified topography and intermediate complexity. Similarly, one may design other idealized models, such as a double-gyre basin or a cross-equatorial basin model with bottom topography and multiple vertical levels (e.g., Khani et al., 2019; Marques et al., 2022), to involve all mesoscale-related processes and generate corresponding MODNS datasets. The LES schemes and their matching parameters might differ if the background flow changes. Therefore, we should evaluate LES-related methods across various idealized and realistic models.

ISOM 1.0 uses the simple Laplacian and biharmonic dissipation schemes to handle the closure of the momentum equations. Nevertheless, in the case of the kilometer-scale high resolution, employing a localized scheme (e.g., Leith, 1996) or a coupled anisotropic horizontal and vertical dissipation scheme (e.g., Khani and Waite, 2020) might potentially optimize the model performance near the grid-scale range (e.g., submesoscale). In addition, the advection schemes can also influence the oceanic numerical simulations (Uchida et al., 2022; Thiry et al., 2024). Testing the impact of advection schemes on kilometer-scale high-resolution oceanic simulations is a challenging task that requires thorough planning and joint effort of the field.

Finally, we do not use any effective submesoscale-oriented parameterization in the simulation, and the model resolution is insufficient to explicitly resolve some crucial submesoscale processes (e.g., the symmetric instability) in the upper ocean. Thus, using MODNS data to design pure submesoscale parameterization schemes is inappropriate. Improving the horizontal resolution and optimizing the simulation of submesoscale processes in the mixed layer to generate the submesoscale ocean DNS is also an orientation for future work.