# Exploring Baroclinic Instability of the Computational Kind (BICK) in Numerical Simulations of the Ocean

Claire Ménesguen, Nicolas Ducousso, Clement Vic, Sylvie Le Gentil

► To cite this version:

Claire Ménesguen, Nicolas Ducousso, Clement Vic, Sylvie Le Gentil. Exploring Baroclinic Instability of the Computational Kind (BICK) in Numerical Simulations of the Ocean. Journal of Advances in Modeling Earth Systems, 2025, 17 (4), e2024MS004600 (23p.). 10.1029/2024MS004600. hal-05105673

HAL Id: hal-05105673

https://hal.science/hal-05105673v1

Submitted on 11 Jun 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

10.1029/2024MS004600

**Key Points:**

* Ocean models that generally use the Lorenz vertical grid staggering can suffer from the spurious numerical instability called Baroclinic Instability of the Computational Kind (BICK)

* The nonlinear evolution of BICK produces spurious small scale eddies and reduces frontal sharpness

* In order to avoid the development of BICK, we recommend to design the grid such that $δx/δz > 2N/f$

**Correspondence to:**

C. Ménesguen,
claire.menesguen@ifremer.fr

**Citation:**

Ménesguen, C., Ducoussou, N., Vic, C., & Le Gentil, S. (2025). Exploring baroclinic instability of the computational kind (BICK) in numerical simulations of the ocean. *Journal of Advances in Modeling Earth Systems, 17*, e2024MS004600.
https://doi.org/10.1029/2024MS004600

Received 26 JUL 2024

Accepted 11 MAR 2025

© 2025 The Author(s). Journal of Advances in Modeling Earth Systems published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

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

# Exploring Baroclinic Instability of the Computational Kind (BICK) in Numerical Simulations of the Ocean

Claire Ménesguen¹,² ID, Nicolas Ducoussou³, Clément Vic¹ ID, and Sylvie Le Gentil¹

¹Laboratoire d'Océanographie Physique et Spatiale, UMR6523, CNRS, IRD, Ifremer, Université de Bretagne Occidentale, Brest, France, ²ODYSSEY, Rennes, France, ³Service Hydrographique et Océanographique de la Marine, Brest, France

**Abstract** Primitive-equation models are essential tools for studying ocean dynamics and their ever-increasing resolution uncovers ever finer scales. At mesoscales and submesoscales, baroclinic instability is one of the main drivers of turbulence, but spurious numerical instabilities can also arise, leading to nonphysical dynamics. This study investigates a spurious instability termed Baroclinic Instability of Computational Kind (BICK), discovered in Arakawa and Moorthi (1988, https://doi.org/10.1175/1520-0469(1988)045<1688:BIIVDS>2.0.CO;2) and Bell and White (2017, https://doi.org/10.1016/j.ocemod.2017.08.001), through idealized configurations using a vertical (Modified) Lorenz grid. Here, we explore the growth of BICK within quasi-geostrophic (QG) and hydrostatic primitive-equation (HPE) frameworks for different setups: the canonical Eady configuration, stratification-modified Eady configurations, and a surface-intensified jet configuration. Our results confirm that the emergence of BICK is specific to the vertical staggering of the (Modified) Lorenz grids. Its growth is consistent with linear QG theory, and BICK is confined near the surface and bottom boundaries. In HPE simulations, the nonlinear evolution of BICK generates small-scale spurious eddies and reduces frontal sharpness. Increasing the number of levels reduces BICK's horizontal scale down to below the model's effective resolution. We illustrate this property using regional HPE simulations with a varying number of levels. BICK is found to significantly affect the vertically under-resolved simulations by introducing small-scale noise from both the bottom and surface boundaries. Our recommendation is to keep the ratio between the model horizontal ($δx$) and vertical ($δz$) resolution greater than $2N/f$, where $N$ is the Brunt-Väisälä frequency and $f$ the Coriolis parameter, to minimize the impact of BICK on the dynamics.

**Plain Language Summary** Numerical simulations of the ocean circulation are routinely used to investigate regional dynamics. In recent years, increases in their resolution have allowed the community to explore new ranges of fine-scale dynamics. However, these new regimes of dynamics come with new numerical challenges inherent to the increase in resolution. In addition to physical instabilities, numerical instabilities artificially introduce spurious fine-scale dynamics. Here, we investigate such an instability called Baroclinic Instability of the Computational Kind (BICK), identified in previous studies. We study BICK in different configurations, from idealized to realistic setups. We show that BICK is initially excited close to the bottom or surface boundaries, especially where the density changes rapidly over small scales along the boundaries. We illustrate the effects of this instability on a regional simulation of the Mozambique Channel circulation. Specifically, we show that, contrary to common practices, the vertical resolution of the grid has to be refined hand-in-hand with the horizontal resolution in order to tone down BICK. By identifying the conditions that trigger BICK, we aim to unravel the physical dynamics and numerical artifacts of small-scale ocean simulations.

## 1. Context and Motivation

Primitive-equation models are powerful tools for studying regional ocean dynamics. They have been widely used in recent decades. Their resolution has increased and the physics they represent has become more realistic (e.g., state of the art is reviewed in Fox-Kemper et al., 2019). A ubiquitous feature of ocean dynamics is baroclinic instability (BI, e.g. Chapter 6 in Vallis, 2006). BI is a major source of balanced turbulence, acting through a wide range of spatial and temporal scales, from basin scales to mesoscales and submesoscales (Capet et al., 2016; Hochet et al., 2015; Smith, 2007). It was furthermore shown to be active in surface (Boccaletti et al., 2007) and bottom (Wenegrat et al., 2018) boundary layers. In addition to physical instabilities, numerical (nonphysical) instabilities can arise spontaneously in simulations (e.g., Barham et al., 2018; Ducoussou et al., 2017; Hallberg, 2005). These instabilities lead to spurious dynamics with consequences that remain difficult to quantify on, for example, energy balance or tracer transport and mixing.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Surface relative vorticity snapshots after 105 days of simulations within the Mozambique Channel with horizontal resolutions of 500 m, and with three distinct vertical discretizations, stretched to increase the resolution near the surface: (a) 30 levels, (b) 90 levels, and (c) 360 levels (with vertical resolution varying from 4 cm at surface to 60 m at depth). Spurious noise is mainly observed for the lowest vertical resolutions.</figcaption></figure>

The motivation for this work is as follows. Using a regional simulation of the Mozambique Channel circulation based on the primitive-equation Coastal and Regional Ocean Community model (CROCO, Auclair et al. (2024), detailed in Section 3), we detected a potentially spurious instability developing in the surface layer. This is illustrated in Figure 1, which shows snapshots of relative vorticity in simulations that only differ by the number of vertical levels. One can see that at low vertical resolution (30 levels), thin stripes of vorticity patterns ubiquitously populate the field. Increasing the resolution to 90 levels and further to 360 levels dramatically reduces these noisy patterns. We hypothesized this numerical artifact to be related to Baroclinic Instability of the Computational Kind (BICK), as first named by Arakawa and Moorthi (1988). They showed that BICK can develop when the Lorenz vertical staggering of variables is used (hereafter L-grid, as illustrated in Figure 2). This staggering is such that the densities are not at the same points as the vertical velocities but between them. Analysis of Arakawa and Moorthi (1988) has shown that discretization of the Quasi-Geostrophic (QG) equations based on this vertical staggering not only gives rise to physical BI of the Eady or Green type, but also allows unstable spurious modes to develop on small horizontal scales and at the upper and lower boundaries. Furthermore, their study indicates that the discretization of the same equations using the Charney-Phillips vertical staggering (hereafter CP-grid), where the density is collocated with the vertical velocity, prevents the growth of short-wave BICK. The same authors demonstrated that the BICK occurrence is not limited to the discretization of QG equations but also extends to the discretization of primitive equations. This has prompted some to prefer the Charney-Phillips staggering when designing atmospheric general circulation models (e.g., Arakawa & Konor, 1996; Konor & Arakawa, 1997). This approach seems to have been little adopted for the design of quasi-Eulerian vertical coordinate ocean models, for which Lorenz staggering is most often used. This is the case with CROCO, which is therefore likely to let BICK grow. An important property of BICK shown by Arakawa and Moorthi (1988) is that its characteristic horizontal scale shifts to small scales with increasing vertical resolution, until it may no longer be supported by the finite grid size of the model. The sensitivity of the noise to vertical resolution in our Mozambique simulations (Figure 1) follows this rule and leads us to believe that it is BICK.

Bell and White (2017) developed an analytical framework for the QG Eady problem that provides growth rates and scales for BICK. They developed their analysis on a Modified-Lorenz grid (hereafter ML-grid, see Figure 2). ML-grid discretization is a natural choice for hydrostatic balance and its properties are very similar to those of the L-grid (Bell, 2003).

This paper aims to identify the conditions that lead to BICK and to provide recommendations on numerical choices to mitigate it. We revisit and expand on previous work related to BICK, with a special attention toward more realistic applications. The structure of the paper is as follows. In Section 2, we analyze linear instabilities within the QG and Hydrostatic Primitive Equation (HPE) frameworks through various case studies. These include the classical

Eady configuration, and modified Eady configurations with variable surface stratification to better represent oceanic conditions. We then extend the analysis to a surface-intensified jet configuration. In Section 3, we investigate the nonlinear evolution of BICK within the aforementioned idealized configurations through HPE simulations performed with CROCO. Finally, in Section 4, we summarize our findings and discuss their implications.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The staggering of variables used by the Charney-Phillips grid (left), the original (middle) and modified (right) Lorenz grids for QG (blue) and HPE (black) variables.</figcaption></figure>

## 2. Linear Instabilities

In this section, we perform a linear stability analysis of four idealized configurations that are representative of open ocean currents. The first case is the well-known Eady configuration (e.g., Chapter 6 in Vallis, 2006). The next two cases are slightly modified from the Eady configuration to account for modified stratification near the surface. The final case is a surface intensified jet on an *f* plane. Figure 3 shows the different configurations.

### 2.1. 1D Linearized Inviscid Equations

#### 2.1.1. QG Equations

We consider the QG equations, linearized around a base state comprising a zonal flow $U(z)$, with no lateral variation, and a temperature field $T^b(z) + T(y, z)$, where $T^b(z)$ represents the background temperature profile and $T(y, z)$ is the component in geostrophic balance with the zonal flow. We assume a fluctuation eigenmode solution with its streamfunction of the form $\phi(z)\exp(i(k_x x + k_y y) + \sigma t)$ and the QG potential vorticity $q(z)\exp(i(k_x x + k_y y) + \sigma t)$. The real part of $\sigma$, $\Re(\sigma)$, is the growth rate and the opposite of its imaginary part, $-\Im(\sigma/k_x)$, is the phase speed in the direction of increasing $x$.

$$ (\sigma + ik_x U) q + ik_x \left[ -f^2 \partial_z \left( \frac{U_z}{N^2(z)} \right) \right] \phi = 0 \quad (1a) $$

$$ q = \left[ -k_x^2 - k_y^2 + f^2 \partial_z \left( \frac{1}{N^2(z)} \partial_z \right) \right] \phi, \quad (1b) $$

with

$$ N^2 = -\frac{g}{\rho_0} \rho_z, \quad (2a) $$

$$ \rho = \rho_0 [1 - \alpha_T (T - T_{ref})], \quad (2b) $$
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The four different idealized base states studied in section. Colors are for the zonal velocity and continuous lines are isopycnals. The left side of each plot shows the vertical profile of the Brunt-Väisälä frequency and generalized QG-potential vorticity gradient (cf. Equation 4) in the center of the domain. Insets show zooms of the generalized QG-potential vorticity gradient down to 1,500 m. In the case of the surface intensified jet (d), the y-dependent base state is used as the initial conditions of the numerical integrations (Section 3.5), while the linear stability analysis (Section 2.4) is computed at the center of the domain, without considering the y-dependence.</figcaption></figure>

g the gravity, $\rho_0$ the density in the Boussinesq approximation, and $\alpha_T$ the thermal expansion.

The set of equations leads to matrix eigenvalue problems, which we solve using the eigensolver provided by the SciPy package in Python. The boundary conditions are derived by requiring $w = 0$ at the upper and lower boundaries. Two vertical discretisations are implemented: the CP-grid, as described by Arakawa and Moorti (1988), and the ML-grid, as detailed in Bell and White (2017). Discretisations are written in Appendix A1 and A2.

### 2.1.2. Hydrostatic Primitive Equations

We have considered HPE linearized around the base state to compute the full set of eigenmodes beyond the QG framework for the two vertical discretisations, L-grid, as described in the appendix of Molemaker et al. (2005) and CP-grid which is a modified version of the L-grid code following Figure 2. Discretisations are written in Appendix A3 and A4.

The base state is the same as in the previous section with $U(z), P(z), T^h(z) + T(y, z)$ in geostrophic balance. We assume a fluctuation eigenmode solution form of $\{u(z), v(z), w(z), p(z), \theta(z)\} \exp(i(k_x x + k_y y) + \sigma t)$. The linearized inviscid HPE equations are:
$$ \sigma u + ik_x U u + U_y v + U_z w - fv + \frac{1}{\rho_0} ik_x p = 0 \quad (3a) $$

$$ \sigma v + ik_x U v + fu + \frac{1}{\rho_0} ik_y p = 0 \quad (3b) $$

$$ \frac{1}{\rho_0} p_z - g \alpha_T \theta = 0 \quad (3c) $$

$$ ik_x u + ik_y v + w_z = 0 \quad (3d) $$

$$ \sigma\theta + T_y v + ik_x U\theta + (T_z^b + T_z) w = 0 \quad (3e) $$

The HPE matrix eigenvalue problem is also solved with the eigensolver of the SciPy package in Python. For inviscid equations, the boundary condition is reduced to the no-normal flow condition to the surface such that $w = 0$.

## 2.2. The Eady Case

The base state is a zonal flow with speed $U(z)$ that depends linearly on the vertical coordinate $z$, and a stratification that is the sum of a background linear profile and a field that depends on the meridional coordinate $y$ and that is in geostrophic balance with the zonal flow (Figure 3a). In other words, the background temperature and zonal velocity are written as $T^b(z) = \frac{T_0}{H}(z - z_{bottom})$, and $U(z) = \frac{\Delta U}{H}(z - z_{mid})$. In our study, the parameters are: the Brunt-Väisälä frequency $N = 1.3 \times 10^{-3}$ s$^{-1}$, the Coriolis parameter $f = 4 \times 10^{-4}$ s$^{-1}$, the difference between top and bottom velocity amplitudes $\Delta U = 1$ ms$^{-1}$, the total depth $H = 4000$ m, the gravity $g = 9.81$ s$^{-1}$, $\rho_0 = 1025$ kgm$^{-3}$, $\alpha_T = 2.7 \times 10^{-4}$°C$^{-1}$ and $T_0 = \frac{N^2 H}{g \alpha_T} = 2.5$°C.

We checked that the set up was unstable regarding the generalized potential vorticity gradient Charney-Stern-Pedlosky criterion (e.g., Chapter 6 in Vallis, 2006). It requires a change of sign of the meridional gradient of the QG potential vorticity,

$$ Q_y = f \partial_z \frac{b_y}{N^2} + \frac{f^2}{N^2} \frac{dU}{dz} (\delta_{upper} - \delta_{lower}), \quad (4) $$

with $b = -\frac{\rho g}{\rho_0}$ (see Equation 6.1 in Arakawa & Moorthi, 1988). Figure 3a shows the change of sign between the surface and the bottom for the classic BI in the Eady case.

### 2.2.1. Analytical BICK QG Eigen Modes

On the L-grid, in a QG framework, Arakawa and Moorthi (1988) showed that the Eady flow develops BICK modes, while Bell and White (2017) further formulated on the ML-grid an analytical approximation to compute their characteristics. As the approximate solutions reproduce accurately the numerical solutions, we present a few results from this approximation to give an initial overview of the problem. For a perturbation streamfunction on an integer level $k$, written as $\psi'_k = \phi_k \exp[ik_x(x - ct)] \sin(k_y y)$, the growth rate of the perturbation is $k_x \Re(c)$ and its phase speed is $\Re(c)$. The BICK modes can be characterized by the position of their critical level between $k = j$ and $k = j + 1$, in which the mean flow $U_{j+1/2}$ is closest to the growing mode's phase speed $\Re(c)$.

Since the approximation in Bell and White (2017) was developed for BICK, where the most unstable mode occurs when $j$ is the closest to the boundaries, they focus on critical layers near the boundaries. Consequently, they select small values of $j$, and because the Eady problem is vertically symmetric, the results apply equally well to the opposite boundary.

Bell and White (2017) consider constant $N^2$ and $\delta_z$, and make two approximations: (a) they neglect the potential vorticity gradient of the mean flow and the variation of $U$ within each vertical layer other than the critical layer; (b) they neglect the component of the solution that grows exponentially for $k > j + 1$ in the equations for
$k < j + 1$. Using (a) the equation for the streamfunction of the perturbation at all half-levels other than $k = j + 1/2$ is given by

$$ \left(1 - \frac{m^2}{2}\right) \left(\phi_{k+3/2} + \phi_{k-1/2}\right) - (2 + m^2)\phi_{k+1/2} = 0, \qquad (5a) $$

$$ m^2 = \frac{N^2 \delta z^2 (k_x^2 + k_y^2)}{2f^2}, \qquad (5b) $$

where $\delta z$ is the vertical grid spacing. Bell and White (2017) give the solution of this equation as

$$ \phi_k = a s^{-(k-j-1/2)} + b s^{(k-j-1/2)}, \qquad (6) $$

where $a$ and $b$ are arbitrary constants (determined by the boundary conditions) and

$$ s^2 - 2\tau s + 1 = 0, \quad s < 1, \qquad (7a) $$

$$ \tau = \frac{2 + m^2}{2 - m^2}. \qquad (7b) $$

From there, they obtain an equation for the complex phase speed $c$ in the case of an arbitrary vertical shear (cf. their Equation 24) and also one corresponding to the Eady case of constant vertical shear which is considered here. A non-dimensional complex number $c_j$ is introduced such that

$$ c_j = \frac{c - U_{j+1/2}}{\delta U}, \qquad (8) $$

with $\delta U = U_{k+1} - U_k = \Delta U/K$, $K$ being the number of levels and $\Delta U = U_{K+1/2} - U_{1/2}$. Hence, Bell and White (2017) end up with the growth rate of the perturbation as the solution of a simple quadratic equation:

$$ \mu_2 c_j^2 + \mu_1 c_j + \mu_0 = 0, \qquad (9) $$

with

$$ \mu_2 = 16(1 - s), \qquad (10a) $$

$$ \mu_1 = a_{\mu} + b_{\mu}, \qquad (10b) $$

$$ a_{\mu} = -4(1 + s)^2 + 16\left(j - \frac{1}{2}\right)(1 - s), \qquad (10c) $$

$$ b_{\mu} = 2s^{2j-1}(1-s)^3, \qquad (10d) $$

$$ \mu_0 = \frac{1}{2} s^{2(j-1)} (1-s)^2 \left[ (1+s)^2 + \left(j - \frac{1}{2}\right) 4s(1-s) \right]. \qquad (10e) $$

In the case of BI and BICK in the Eady case, the most unstable eigenmode is found for $k_y = 0$ (not shown here). As a first approximation, we keep this assumption in our four case studies. For perturbations with no lateral variation ($k_y = 0$), the growth rate $k_x \Im(c)$ is given by:

$$ k_x \Im(c) = m \Im(c_j) \frac{\sqrt{2f\Delta U}}{NH}. \qquad (11) $$

Figure 4a shows approximate growth rates computed from Equation 9 for our Eady configuration. Growth rates are calculated for modes with critical levels in layers 1 to 5. The most unstable mode is the one whose critical level
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Growth rate of the BICK instability according to Bell and White (2017)'s short-wave approximation on a Lorenz grid in a QG regime. (a): computation for 16 levels, $\Delta U = 1 \cdot m \cdot s^{-1}$ and $N = 1.3 \times 10^{-3} \cdot s^{-1}$. Vertical dashed dot lines are aligned with the highest growth rates resulting from Bell and White (2017) calculations for the first 3 critical levels with $m_1 = 1.1$, $m_2 = 0.37$ and $m_3 = 0.24$. Dashed purple and black vertical lines are aligned with Nyquist frequencies for high and coarse resolution simulations respectively, discussed in Section 3.3. Dotted purple and black vertical lines are aligned with effective resolutions for the same simulations. (b–d) show the influence of the number of levels, velocity shear amplitude and stratification, respectively, on the predicted growth rates.</figcaption></figure>

is closest to the boundary ($j = 1$). However, BICK instabilities have positive growth rates for several wavelength ranges. As we move away from the boundary (increasing $j$), the modes of BICK that are likely to be destabilized have lower growth rates, and affect larger and larger scales.

In the QG framework, the analytical development of Bell and White (2017) gives us some insight into how the setup of our configuration can affect the BICK properties. Figure 4b shows that the wavelengths of the most unstable spurious modes become smaller as $\delta z$ decreases (increasing number of levels). Figure 4c shows that the larger the shear, the faster the spurious modes grow. Finally, Figure 4d illustrates that as $N$ decreases, BICK occurs at smaller scales with higher growth rates. Note that this is exactly the opposite for $f$, since $N$ and $f$ always appear in the calculation as the ratio $N/f$.

### 2.2.2. QG and HPE Eigen Modes

In Figure 5a, the top panel shows the growth rates of unstable modes of BI and BICK solved from the 1D linearized QG equations. For both vertical discretisations, the BI mode is found at scales predicted by theory $L_{Eady} = (2\pi NH/f)/1.61 \sim 50$ km (cf. Gill, 1982). Considering the ML-grid, the BICK modes are in very close agreement with the Bell and White (2017) calculations. Considering the CP-grid, no BICK mode appears. The ML-grid supports two types of BICK modes: one with a positive phase speed and one with a negative phase speed
<figure><img src="image_7.png" /><figcaption>**Figure 5.** BI and BICK eigenmodes solved from linearized QG and HPE 1D equations for 16 levels. (a) Growth rates from the QG equations in the top panel and the HPE equation in the bottom panel. Blue dots are for the L- or ML-grid and red and coral crosses are for the CP-grid. Black dashed curves are the analytical QG approximations from Bell and White (2017) of the BICK growth rates. Blue vertical dashed lines are aligned with their maxima. Gray vertical dashed line is the expected highest growth rate of the BI in a QG framework. (b) Corresponding phase speeds. Gray horizontal dashed lines are aligned with discrete values of *U* initialized at integer levels. (c) Vertical profiles of the HPE eigenmodes for the BI and the first three BICK eigenmodes, respectively. The BICK eigenmode with the highest growth rate is confined very close to the boundary.</figcaption></figure>

(Figure 5b), corresponding to the BICK modes near the upper boundary and near the lower boundary, respectively. The critical levels, defined as the levels where the phase speed of the BICK modes is equal to the base flow, are at half-integer levels for the ML-grid, as shown in Bell and White (2017). In Figure 5a, the lower panel shows the growth rates of unstable modes of BI and BICK solved from the 1D linearized HPE. The growth rates of the BICK modes in the ML-grid are still in good agreement with the QG predictions of Bell and White (2017). We evaluate dimensionless numbers as the Rossby number $Ro = U_0/fL$ and the Burger number $Bu = (NH/fL)^2$, with *L* a characteristic length. This characteristic length is 50 km for the BI mode, and gives $Ro = 0.02$ and $Bu = 0.07$. For the first BICK mode, the characteristic length is 3.3 km, giving $Ro = 0.4$ and $Bu = 15.5$. Yet,
QG approximation is relevant for $Ro \ll 1$ and $Bu = O(1)$. In both cases, we are not too far from QG conditions and HPE formulation does not change QG results. Considering the CP-grid, some modes appear between those of the L-grid, with very small growth rates and at critical levels at $j$ integer levels (Figure 5b, lower panel). However, their growth rates are small compared to the growth rates of BI modes, hence they are not likely to develop in a simulation.

Figure 5c shows in the first panel the vertical structure of the canonical Eady BI modes with a smooth structure over the whole column. Conversely, other panels show the BICK modes relatively confined to a boundary (to the upper boundary in the figure, those confined to the lower boundary are not shown, but are symmetric). The most unstable BICK mode is the one most tightly confined to a boundary.

## 2.3. Two Modified Eady Flows

### 2.3.1. With a Surface-Intensified Stratification

We set up a modified Eady configuration where the background stratification is increased near the surface to mimic springtime restratification (Figure 3b). The configuration is the same as the previous Eady configuration, except that the background stratification is changed to: $T^b(z) = \frac{T_0}{H}(z - z_{bottom}) + 3.1T_0 \exp\left(-\left(1 - 4\left(\frac{z - z_{surf}}{H}\right)\right)^2\right)$. The vertical variation of $Q_y$ is smoother at the surface than in the classic Eady case (Equation 3b). This prevents the growth of BICK at the surface.

Results of the stability analysis are shown in Figures 6a and 6b. Essentially, BICK modes with critical levels near the surface do not grow or grow slowly, while BICK modes with critical levels near the bottom and BI modes are very similar to those obtained in the reference Eady case. In other words, the increased stratification mutes the growth of BICK modes. Actually, $N$ appears in the denominator of the growth rate Equation 11. It implies that as $N$ increases, the growth rate decreases (Figure 4d). Moreover, as demonstrated in the results of Bell and White (2017), BICK modes are triggered by two adjacent levels and are confined to a few levels. Therefore, it is reasonable to consider the local $N$ in the growth rate rather than the depth average as it is the case in the Eady problem. That is why changes in stratification near a boundary have a significant impact on BICK modes.

### 2.3.2. With a Surface Low Stratification

In contrast, we set up a modified Eady configuration with a stratification that decreases at the surface to account for destratification effects (Figure 3c). The configuration is the same as the previous Eady configuration, except that the background stratification is modified as: $T^b(z) = \frac{T_0}{H}(z - z_{bottom}) - 0.31T_0 \exp\left(-\left(1 - 4\left(\frac{z - z_{surf}}{H}\right)\right)^2\right)$. There are two scales over which $Q_y$ changes sign (Figure 3c). We retrieve the large scale BI, and there is an additional change of sign at small vertical scale that corresponds to surface mixed layer instability (Boccaletti et al., 2007).

Results of the linear stability analysis are shown in Figures 6c and 6d. CP-grid has now unstable modes at small scales. Their growth rates have shifted maxima compared to those of the L-grid, but the growth rates of both grids increase drastically at small scales, converging to a same curve when increasing the number of levels (see light blue dots and orange crosses in Figure 6c). This small scale instability is not spurious because it is due to the change of sign of $Q_y$ close to the surface, caused by the low subsurface stratification. It shares similarities with mixed layer instabilities (Boccaletti et al., 2007).

## 2.4. The Surface Intensified Jet Case

We set up a last idealized configuration that features a more realistic water column dynamics than the simple Eady problem. This consists in a surface intensified zonal flow with no lateral shear $U(z) = \Delta U \exp\left(\frac{zH}{H_{jet}}\right)$ with a constant background stratification $T^b(z) = z - z_{bottom}$. This set up, even if no $y$-dependence is taken into account, is relevant to interpret the result of the HPE simulations in the next section (Figure 3d). We have set $N = 4$. $\times 10^{-3} \text{ s}^{-1}$, $f = 1$. $\times 10^{-4} \text{ s}^{-1}$, $\Delta U = 0.6 \text{ m}\cdot\text{s}^{-1}$, $H = 3000 \text{ m}$, $H_{jet} = 1000 \text{ m}$.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** BI and BICK eigenmodes solved from linearized 1D HPE with 16 levels in the cases where: (a, b) the stratification increases at the surface, (c, d) the stratification decreases at surface, and (e, f) the velocity is a surface intensified jet. Blue and light blue dots are for the L-grid. Red and coral crosses are for the CP-grid. In panels (c–f), the light blue dots are the computation for 100 levels on the L-grid, while the orange crosses are the computation for 100 levels on the CP-grid.</figcaption></figure>

In this case (Figures 6e and 6f), BI develops around 250 km. On the L-grid, the largest growth rate of BICK ($\lambda_x$ = 30 km) is more than twice as large as the BI growth rate and the most unstable wavelengths are still accurately estimated by Bell and White (2017)'s theory. The CP-grid also supports a subsidiary maximum at small scales. However, as the number of levels increases, this subsidiary maximum is transformed into multiple subsidiary maxima and shifted to smaller scales with decreasing growth rates. This is not the case for the L-grid, where the highest growth rate is shifted to smaller scales as the number of levels is increased, but is maintained at the same magnitude. The L-grid supports BICK, while the CP-grid supports subsidiary maxima due to the discretization, as in Bell and White (1988), where spurious subsidiary maxima are identified with the following characteristics:“as the number of levels increases, (a) the number of subsidiary maxima increases and (b) the magnitude of each subsidiary maxima decreases, so that the growth rates approach the monotonic variation obtained with a very large number of levels.”
### 3. Hydrostatic Primitive-Equation Simulations

We now analyze the free evolution of the different idealized base states in the primitive-equation ocean model CROCO.

#### 3.1. Ocean Model Configuration

CROCO (Auclair et al., 2024) solves the Hydrostatic Primitive Equations on a horizontal C-grid and a L-grid in the vertical (see Figure 2). The domain is a flat-bottom re-entrant channel in the zonal direction. It has a free surface, a free-slip bottom and free-slip walls in the north and south. The Coriolis parameter is set constant over the domain (f-plane). The advection scheme in the horizontal is a 5th-order upstream biased scheme (UP5) for both momentum and tracers, and a fourth-order centered parabolic spline reconstruction (SPLINES) in the vertical. We do not use explicit dissipation and diffusion in the horizontal and vertical.

For the Eady configuration, the number of points on the horizontal grid is 70 × 35 (corresponding to a grid spacing of 1.4 km in both directions) and 560 × 280 (corresponding to a grid spacing of 180 m in both directions) for 16 vertical levels (corresponding to a constant vertical spacing of 250 m). For the surface intensified jet configuration, the number of points is 400 × 400 (corresponding to a grid spacing of 2 km in both directions), for 16, 32, 54 and 200 levels (corresponding to a constant vertical spacing of 187.5 m, 93.7 m, 55.5 m and 15.0 m). In this study case, the vertical grid spacing is set constant using standard CROCO parameters (Auclair et al., 2024): $\theta_s = 0$, $\theta_b = 0$, $h_c = 10^6$.

Given a grid spacing, we consider two different length scales. First, the Nyquist length scale is the minimum wavelength that can theoretically be represented in the model without aliasing. It is twice the grid spacing. Second, the effective resolution is defined as the cutoff length scale below which numerical dissipation progressively damps the dynamics (Ménesguen et al., 2018; Soufflet et al., 2016). Given the horizontal advective schemes (a 5th order upstream scheme), the horizontal effective resolution is 8 times the grid spacing (Ménesguen et al., 2018). In the vertical, effective resolution is less documented.

#### 3.2. The Base States

In all cases, the initial flow is in geostrophic balance, such that the surface height $\eta$ and the temperature $T$ are defined as:

$$ 
\begin{aligned}
\eta(x, y) &= \frac{-f}{g} \int_{0}^{y} U_{surf} dy, \\
T(x, y, z) &= T^{b}(z) - \frac{f}{g \alpha_{T}} \int_{0}^{y} \partial_{z} U \ dy,
\end{aligned} 
\qquad (12) 
$$

with

$$ 
\begin{aligned}
x \in [-L_x/2, L_x/2], \quad y \in [-L_y/2, L_y/2], \quad z \in [-H, 0], \\
T^b(z) = T_{surf} + \frac{T_{surf} - T_{bottom}}{H} z + \tau(z).
\end{aligned} 
$$

$T^b$ is the background temperature for the three stratifications discussed in the study. It has a surface expression $\tau(z)$, which is defined as $\tau(z) = 10 \exp(-(1 - 4z/H)^2)$ in the surface-intensified stratification case, $\tau(z) = 0$ in the constant stratification case and $\tau(z) = -\exp(-(1 - 4z/H)^2)$ in the surface low stratification case. We use a linear equation of state $\rho = \rho_0[1 - \alpha_T(T - T_{ref})]$. Parameters are defined in Table 1.

In the Eady configuration,

$$ U_{Eady} = U_{surf} \left( \frac{2z - H}{H} \right) (1 + 5 \times 10^{-3} \times R), $$

with $R$ a random white noise between 0 and 1.
<table><caption>**Table 1**
Parameters for the Different Configurations</caption><thead><tr><th>General</th><th>Eady configuration</th><th>Surface-intensified jet</th></tr></thead><tbody><tr><td>g = 9.81 m.s<sup>-2</sup></td><td>U<sub>surf</sub> = 0.5 m.s<sup>-1</sup></td><td>U<sub>surf</sub> = 0.6 m.s<sup>-1</sup></td></tr><tr><td>α<sub>T</sub> = 2.8.10<sup>-4</sup> °C<sup>-1</sup></td><td>H = 4000 m</td><td>H = 3000 m</td></tr><tr><td>ρ<sub>0</sub> = 1025 kg.m<sup>-3</sup></td><td>L<sub>x</sub> = 100 km</td><td>L<sub>x</sub> = 800 km</td></tr><tr><td>T<sub>ref</sub> = 25°C</td><td>L<sub>y</sub> = 50 km</td><td>L<sub>y</sub> = 800 km</td></tr><tr><td></td><td>f = 4.10<sup>-4</sup> s<sup>-1</sup></td><td>f = 1.10<sup>-4</sup> s<sup>-1</sup></td></tr><tr><td></td><td>T<sub>surf</sub> = 26°C</td><td>T<sub>surf</sub> = 23°C</td></tr><tr><td></td><td>T<sub>bottom</sub> = 23°C</td><td>T<sub>bottom</sub> = 4°C</td></tr><tr><td></td><td></td><td>H<sub>jet</sub> = 1000 m</td></tr><tr><td></td><td></td><td>L<sub>jet</sub> = 40 km</td></tr></tbody></table>

In the surface intensified jet configuration,

$$U_{jet} = U_{surf} \exp\left(-\left(\frac{y - y_0}{L_{jet}}\right)^2\right) \exp\left(\frac{z}{H_{jet}}\right) (1 + 5 \times 10^{-3} \times R),$$

with $y_0$, centered in the domain.

### 3.3. Eady Configuration

In Section 2.2.2, we built upon Arakawa and Moorthi (1988)'s and Bell and White (2017)'s studies to show that BICK modes grow in the HPE Eady configuration if the horizontal resolution is fine enough to resolve their horizontal wavelengths. We showed that the QG framework provides a good approximation of BICK properties in our case studies.

The theory predicts that the most unstable Eady mode in our configuration is $L_{Eady} \sim 50$ km, that is, half the domain in the x direction. With the

70 × 35 × 16 grid (Figure 7), the Nyquist length scale is about 3 km and the effective resolution is about 12 km (see Section 3.1 for definitions). Consequently, the model configuration cannot accurately resolve the fastest growing modes of BICK ($\lambda_x \sim 3.6$ km) because it falls within the dissipative range of the model (cf. Figure 4a). Indeed, BICK is absent from the dynamics at time scales relevant to its growth (Figure 7 for the dynamics and Figure 5a for the growth rate), and the overall pattern aligns well with expectations for the Eady configuration, featuring the development of deep modes and the emergence of a coherent structure. With the 560 × 280 × 16 grid, the Nyquist length scale is 360 m and the effective horizontal resolution is about 1.4 km. The fastest growing mode of BICK is resolved and Figures 8a and 8c show additional small structures blurring the BI signal near the

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Fields after 10 days (a and c) and 20 days (b and d) for a 70 × 35 × 16 grid in a canonical Eady configuration. As predicted by theory, only baroclinic instability is triggered.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Fields after 10 days (a and c) and 20 days (b and d) for a 560 × 280 × 16 grid in a canonical Eady configuration. With an effective horizontal resolution of 1.4 km, BICK can develop alongside BI.</figcaption></figure>

upper and lower boundaries. This is consistent with the findings of Arakawa and Moorthi (1988), showing that the expression of BICK modes is highly dependent on the horizontal resolution of a model.

In Figure 9 we show energy wavenumber spectra versus time at the level nearest to the surface and at middepth level. BI and BICK modes are identified by their theoretical length scale. In the two grids the energy grows at the scale of the BI mode (about 50 km) with a growth rate corresponding to $\sigma = 1.35 \text{ days}^{-1}$, which is close to the growth rate predicted by linear theory (about $2.1 \text{ days}^{-1}$). On the 560 × 280 × 16 grid, the fastest growing mode of BICK with a scale of 3.6 km grows rapidly with a growth rate corresponding to $\sigma = 2 \text{ days}^{-1}$, which is the growth rate predicted by linear theory. On the 70 × 35 × 16 grid, the fastest growing mode of BICK is damped by the dissipation of the model. However, the scale of the second fastest growing mode of BICK (11 km) is at the limit of the effective resolution (about 12 km), and the energy at the surface shows the growth of this second mode, with a growth rate corresponding to $\sigma = 0.6 \text{ days}^{-1}$, which is also well predicted by the linear theory. This second BICK mode is far less energetic than the BI mode and is visually absent in physical fields (Figure 7).

In summary, we have shown that linear QG theory is a good approximation for predicting BICK properties in the HPE framework for the Eady configuration. However, in the 3D simulations, the Nyquist scale of the configuration sets the lower bound on the horizontal length scale of the mode that can grow. The lower bounds sets by the Nyquist scale is potentially shifted to slightly larger scales by the effective horizontal resolution of the model.

### 3.4. Eady Modified Configurations With Variable Stratification

Starting from the Eady configuration, we slightly modified it by increasing and decreasing the background stratification near the surface (cf. Section 2). In the first case, Figure 10 shows that, as expected from linear studies, BICK is significantly reduced at the surface as compared to the bottom layer that features typical BICK modes. In the second case, as the stratification decreases near the surface, BICK develops very rapidly as compared to the bottom layer (Figure 11). The small-scale instability observed at the surface at a wavelength of
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Eady case's power spectra versus time of the total energy field $E_{tot} = 0.5(u'^2 + v'^2) + 0.5 \frac{g^2}{\rho_0^2 N^2} \rho'^2$, where ' is the perturbation relative to the initial state. The top row is for the 70 × 35 × 16 grid and the bottom line is for the 560 × 280 × 16 grid. (a, d) are spectra of the surface field where BICK modes are enhanced. (b, e) are mid-depth field spectra where BI is present and BICK modes have reduced amplitude. The black dashed lines indicate the length scale of the BI unstable mode. Blue dashed-dot lines indicate the length scales of the three most unstable BICK modes. The white dotted lines indicate the horizontal effective resolution of the configurations. (c, f) are spectra along the BI and BICK length scales. Fields are saved every 8 hr in the 12 first days.</figcaption></figure>

about 6 km is similar to the second unstable mode of BICK (Figure 6c). The most unstable linear BICK mode is at the limit of the effective resolution of this configuration and may be damped or superimposed on the second.

### 3.5. Surface Intensified Jet Configuration

The jet case simulations presented in this section confirm the results of the Eady cases, and allow us to develop an intuition on the effect of varying the vertical resolution. Figure 12 shows the relative vorticity at the surface, around the jet position, for the four different vertical resolutions used in the study. Following Section 2.2.1, the length scale of the BICK fastest growing mode is given by $\lambda_x = 2\pi/k_x \sim 4N\delta z/f$. Consequently, as we observe in Figure 12, in the two lower vertical resolutions where BICK is resolved, the horizontal length scale of BICK is smaller in the 32-level configuration than in the 16-level configuration. The 54-level configuration should feature BICK instability with a characteristic length scale of 9.3 km. This scale lies in the range where the dynamics are progressively dissipated by the model between the effective resolution and the Nyquist cutoff. BICK mode is actually strongly damped. The last configuration with 200 levels sets BICK's length scale to be less than the Nyquist cutoff and prevents BICK from developing in the simulation.

Figure 13 shows surface relative vorticity fields in the four configurations at later stages, when the nonlinear dynamics have fully developed. The coarser vertical resolution has developed many small-scale eddies that we suspect to be the nonlinear evolution of BICK. Note that these simulations do not have a surface mixed layer, so the typical small scales eddies in Figure 13a cannot be attributed to mixed layer instabilities (Boccaletti et al., 2007). The finer vertical resolution shows sharper fronts and fewer small scale eddies (Figure 13d). It is possible that BICK is acting to destroy sharp fronts and create small scales, potentially increasing dissipation. In the different configurations, a single change in the number of levels produces very different nonlinear
<figure><img src="image_12.png" /><figcaption>**Figure 10.** In a configuration with a surface intensified stratification, fields after 10 days for a 560 × 280 × 16 grid. (a) At surface, only BI is triggered. (b) Bottom fields are similar to what observed in the canonical Eady configuration—Figure 8a. (c) Velocity vertical section at the same time.</figcaption></figure>

developments in turbulence. Note that although only the linear phase around the initial state of the jet is documented in this paper, we suspect that the growth of BICK in the coarser vertical resolution configuration generates more small scales throughout the run that have an impact at later stages than in the finer vertical resolution configuration.

## 4. Conclusion and Discussion

This study is a follow-up to previous studies, which demonstrated the existence of BICK on the original or modified Lorenz grid in an Eady configuration (Arakawa & Moorthi, 1988; Bell & White, 2017). We highlighted the growth of BICK in different idealized configurations, combining linear analyses and numerical simulations with CROCO. The strongest BICK is characterized by critical levels close to the surface and bottom boundaries and its growth depends on the local Coriolis parameter, the local velocity shear and the local stratification. It has been shown that flattened isopycnals at the surface reduce or even inhibit BICK growth near the surface (Figures 6a and 10), while steeper isopycnals increase the instability (Figures 6c and 11). Similar observations can be made at the bottom, for example, in the case of a slope. When the slope of a topography follows the bottom isopycnals, BICK is greatly reduced, whereas a topography with a slope opposite to that of the isopycnals increases instability (not shown). The direct consequence of BICK is the growth of small scales near the boundaries (top or bottom). The long-term effects are substantial, distorting the life cycle of physical instabilities. As Figure 13 shows, where BICK is likely to grow, a smoothing of the fronts occurs. This feature is reminiscent of the impact of other numerical instabilities, such as the Symmetric Instability of the Computational Kind (SICK) (Ducoussou et al., 2017).

A pertinent question is how to configure numerical models to prevent the development of BICK instabilities. Our findings, consistent with previous studies, indicate that the horizontal length scale of BICK decreases with reduced vertical cell thickness. A straightforward way to prevent BICK in simulations is to use a fine vertical grid
<figure><img src="image_13.png" /><figcaption>**Figure 11.** In a configuration with a surface decreased stratification, temperature (a, b) and velocity fields after 5 days for a 560 × 280 × 16 grid. Small scale features are rapidly present at surface (a, c).</figcaption></figure>

spacing near boundaries to dampen the dynamics at BICK's horizontal length scale. In practice, modelers often aim to resolve horizontal processes at specific length scales, which may include BICK if the vertical grid spacing allows its growth. Based on the Eady flow model, the maximum growth rate is expected at a horizontal length scale of $\lambda_x \sim 4N\delta z/f$. Scales smaller than $\delta x_{eff}$ are progressively dissipated numerically, and scales below $2\delta x$ (Nyquist scale) remain unresolved. Consequently, $\delta x/\delta z > 2N/f$ is required to safely prevent BICK from being resolved at the grid scale. However, the numerical horizontal dissipation also acts to dampen the growth of scales between the Nyquist scale and the effective resolution scale. Consequently, depending on the numerical choices of the model, the threshold can be slightly crossed and still prevent BICK growth. In other words, the vertical resolution should be chosen such that $\delta z < \delta x \times f/2N$ to keep BICK safely unresolved but can be slightly crossed, depending on the numerical dissipation.

This concept is illustrated using realistic simulations in the Mozambique Channel with varying numbers of vertical levels. The vertical grid is stretched to increase the resolution near the surface (the vertical resolution is stretched from 4 cm at the surface to 60 m at depth, keeping 10 m resolution at 800 m depth, with code parameters $\theta_s = 6$, $\theta_b = 0$ and $h_c = 10$ m, Auclair et al. (2024)). Consistent with the idealized simulations, we use the UP5 advection scheme for both momentum and tracers in the horizontal, and the SPLINE scheme in the vertical, without explicit horizontal dissipation and diffusion, but adding a vertical turbulent scheme in the form of the k-omega equations (Umlauf & Burchard, 2003). Figure 1 demonstrates the impact of vertical resolution on surface dynamics. Figure 14 further reveals that although BICK is inherent to boundaries, and in particular to the bottom boundary in these simulations, we see that their nonlinear evolution contaminates the ocean interior with grid-scale noise-like patterns, mainly in cases with lower vertical resolution. In cases with higher vertical resolution, these scales are confined to regions with steep topography. To compare small-scale formation with the proposed criterion, we have added contours of the quantity $2N\delta z/f\delta x$, which should be less than 1 to safely avoid BICK. This condition is not met in the 30-level configuration. In the 360-level configuration, the resolution appears to be safe from BICK above approximately 100 m depth. However, the ratio deviates from the ideal value
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Relative vorticity over the Coriolis parameter (*f*), depicted at the surface after 20 days in simulations of surface intensified jets. The simulations employ a horizontal resolution of *δx* = 2 km, with progressively finer resolutions (from a to d subpanels). Subfigures provide a close-up view of the mid-latitude region within the domain.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Surface-relative vorticity over the Coriolis parameter (*f*) during the final stage of simulations spanning 200 days, showing increasing vertical resolutions (from a to d subpanels).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Snapshots of relative vorticity (at constant sigma-levels) along a given latitude in the middle of our domain after 105 days of simulations conducted within the Mozambique Channel. The horizontal resolution is 500 m and the vertical discretizations uses: (a) 30 levels, (b) 90 levels, and (c) 360 levels. Black dashed lines are few isopycnals. Green contour denotes $(2N\delta z)/(f\delta x) = 1$ and warm colors contours are when this ratio equals 5, 10, 15, 35, and 50. In the 360-levels configuration, the ratio is less than 1 above approximatively 100 m and is more than 5 at few spot close to the topography were noise is observed.</figcaption></figure>

of 1, particularly near steep topographies, coinciding with observed small-scale features. This supports the hypothesis that BICK forms when the criterion is not satisfied. For models using sigma coordinates, such as CROCO, a compromise must be found to choose the number of levels and their distribution: thin enough bottom
<figure><img src="image_17.png" /><figcaption>**Figure 15.** *fδx/2N* ratio at the ocean bottom. It represents the maximum cell height to avoid BICK near topography. Maps are computed at the bottom of two simulations: LLC4320 (a) and GIGATL3 (b). LLC4320 is a global simulation of the MITgcm with a spatial resolution of about 2.3 km. The date shown in panel (a) is 14 November 2012. GIGATL3 is a simulation of the Atlantic performed by CROCO with a spatial resolution of about 3 km. The date shown in panel (b) is 20 September 2008.</figcaption></figure>

cells would avoid an over-representation of BICK; but bottom cells have also to be thick enough to prevent pressure gradient errors (Shchepetkin & McWilliams, 2003).

In the design of primitive-equation model configurations, it is now standard to set up a horizontal grid spacing of a few kilometres or less, and down to a few hundreds of meters (e.g., Chassignet & Xu, 2017; Chassignet & Xu, 2021; Rocha et al., 2016; A. L. Stewart et al., 2018; Uchida et al., 2022). Nonetheless, the vertical grid spacing often remains an overlooked parameter, despite its importance for resolving both physical processes (Nelson et al., 2020; K. Stewart et al., 2017) and impeding numerical instabilities such as BICK. Our results highlight the importance of carefully setting the vertical grid as the horizontal resolution increases. We propose a criterion to design the grid aspect ratio in order to avoid BICK in such simulations. Figure 15 shows examples of the *fδx/2N* ratio for two simulations: LLC4320 (Rocha et al., 2016; A. L. Stewart et al., 2018) and GIGATL3 (Gula et al., 2021). LLC4320 is a global simulation of the MITgcm with a spatial resolution of about 2.3 km. GIGATL3 is a CROCO simulation of the Atlantic with a spatial resolution of about 3 km. The ratio corresponds to a minimum *δz* required at the bottom of the ocean to avoid BICK. We should also note that BICK grows faster when a strong velocity shear is present. That is why special care should be taken in the grid design near western boundaries for example, as BICK could feed on the deep currents that feature intense shear. With a resolution of *δx* = 3km, the criterion is not restrictive in the abyssal plains, as a minimum *δz* of 500 m is required. Over the Mid-Atlantic Ridge, the criterion specifies a minimum *δz* of 200 m. The requirement becomes stricter in the South Atlantic, where *δz* values below 50 m are needed at the bottom, and even more so on the continental shelves and western margins, where *δz* at the bottom can drop below 20 m in some areas. Note that these threshold values must be divided by three if the horizontal grid spacing is divided by three to reach ~1 km. Thus, as the horizontal resolution increases, the threshold becomes more restrictive.

# Appendix A: Discretized Equations

## A1. Linearized QG Equations on a CP-Grid

Following Arakawa and Moorthi (1988), by combining the simplest vertical discretisations of the relative vorticity equation centered at full levels and the density equation centered at half levels, we obtain the following discrete equation of the QG potential vorticity centered at full levels
$$ \delta z_1 \frac{D_{g1}}{Dt} (\nabla_h^2 \psi_1 + f) + \frac{D_{g1}}{Dt} \left( \frac{f_0^2}{N_{3/2}^2} \frac{\psi_2 - \psi_1}{\delta z_{3/2}} \right) = 0, \quad (A1a) $$

$$ \delta z_k \frac{D_{gk}}{Dt} (\nabla_h^2 \psi_k + f) + \frac{D_{gk}}{Dt} \left( \frac{f_0^2}{N_{k+1/2}^2} \frac{\psi_{k+1} - \psi_k}{\delta z_{k+1/2}} - \frac{f_0^2}{N_{k-1/2}^2} \frac{\psi_k - \psi_{k-1}}{\delta z_{k-1/2}} \right) = 0, \quad \text{for } 2 \le k \le K-1, \quad (A1b) $$

$$ \delta z_K \frac{D_{gK}}{Dt} (\nabla_h^2 \psi_K + f) + \frac{D_{gK}}{Dt} \left( -\frac{f_0^2}{N_{K-1/2}^2} \frac{\psi_K - \psi_{K-1}}{\delta z_{K-1/2}} \right) = 0, \quad (A1c) $$

where

$$ \frac{D_{gk}}{Dt} \phi = \partial_t \phi + \mathbf{u}_{gk} \cdot \nabla_h \phi = \partial_t \phi + J(\psi_k, \phi) \quad (A2) $$

is the rate of change of any quantity $\phi$ following the quasi-geostrophic flow. Linearizing these equations around a zonal flow $U$ that depends only on the vertical, and looking for a solution in the form $\psi_k(x, y, t) = \phi_k e^{i(k_x x + k_y y) + \sigma t}$, we obtain the following eigenvalue problem

$$ (U_1 k_x - i\sigma) \left( \delta z_1 (-k_x^2 - k_y^2) \phi_1 + \frac{f_0^2}{N_{3/2}^2} \frac{\phi_2 - \phi_1}{\delta z_{3/2}} \right) + k_x \phi_1 \left( \delta z_1 \beta - \frac{f_0^2}{N_{3/2}^2} \frac{U_2 - U_1}{\delta z_{3/2}} \right) = 0 \quad (A3a) $$

$$ (U_k k_x - i\sigma) \left( \delta z_k (-k_x^2 - k_y^2) \phi_k + \frac{f_0^2}{N_{k+1/2}^2} \frac{\phi_{k+1} - \phi_k}{\delta z_{k+1/2}} - \frac{f_0^2}{N_{k-1/2}^2} \frac{\phi_k - \phi_{k-1}}{\delta z_{k-1/2}} \right) + k_x \phi_k \left( \delta z_k \beta - \frac{f_0^2}{N_{k+1/2}^2} \frac{U_{k+1} - U_k}{\delta z_{k+1/2}} + \frac{f_0^2}{N_{k-1/2}^2} \frac{U_k - U_{k-1}}{\delta z_{k-1/2}} \right) = 0 \quad (A3b) $$

for $2 \le k \le K-1$,

$$ (U_K k_x - i\sigma) \left( \delta z_K (-k_x^2 - k_y^2) \phi_K - \frac{f_0^2}{N_{K-1/2}^2} \frac{\phi_K - \phi_{K-1}}{\delta z_{K-1/2}} \right) + k_x \phi_K \left( \delta z_K \beta + \frac{f_0^2}{N_{K-1/2}^2} \frac{U_K - U_{K-1}}{\delta z_{K-1/2}} \right) = 0. \quad (A3c) $$

Note that $\beta = 0$ in the idealized cases studied here.

## A2. Linearized QG Equations on a ML-Grid

Following Bell (2003), Bell and White (2017), by combining the simplest discretisations of the relative vorticity and density equations both centered at full levels, we obtain the following discrete equation of the QG potential vorticity centered at half levels

$$ \frac{\delta z_1}{2} \frac{D_{g1}}{Dt} (\nabla_h^2 \psi_1 + f) + \frac{D_{g1}}{Dt} \left( \frac{f_0^2}{N_1^2} \frac{\psi_{3/2} - \psi_{1/2}}{\delta z_1} \right) = 0, \quad (A4a) $$

$$ \frac{\delta z_k}{2} \frac{D_{gk}}{Dt} (\nabla_h^2 \psi_k + f) + \frac{\delta z_{k+1}}{2} \frac{D_{gk+1}}{Dt} (\nabla_h^2 \psi_{k+1} + f) + \frac{D_{gk+1/2}}{Dt} \left( \frac{f_0^2}{N_{k+1}^2} \frac{\psi_{k+3/2} - \psi_{k+1/2}}{\delta z_{k+1}} - \frac{f_0^2}{N_k^2} \frac{\psi_{k+1/2} - \psi_{k-1/2}}{\delta z_k} \right) = 0 \quad (A4b) $$

for $1 \le k \le K-1$
$$ \frac{\delta z_K}{2} \frac{D_{gK}}{Dt} (\nabla_h^2 \psi_K + f) + \frac{D_{gK}}{Dt} \left( -\frac{f_0^2}{N_K^2} \frac{\psi_{K+1/2} - \psi_{K-1/2}}{\delta z_K} \right) = 0, \quad (A4c) $$

where the streamfunction at full levels is chosen to be expressed as

$$ \psi_k = \frac{1}{2}(\psi_{k-1/2} + \psi_{k+1/2}). \quad (A5) $$

Linearizing these equations around a zonal flow U that depends only on the vertical, and looking for a solution in the form $\psi_k(x, y, t) = \phi_k e^{i(k_x x + k_y y) + \sigma t}$, we obtain the following eigenvalue problem

$$ \begin{aligned} & (U_1 k_x - i\sigma) \left( \frac{\delta z_1}{2} (-k_x^2 - k_y^2)\phi_1 + \frac{f_0^2}{N_1^2} \frac{\phi_{3/2} - \phi_{1/2}}{\delta z_1} \right) \\ & \qquad + k_x \phi_1 \left( \frac{\delta z_1}{2} \beta - \frac{f_0^2}{N_1^2} \frac{U_{3/2} - U_{1/2}}{\delta z_1} \right) = 0, \end{aligned} \quad (A6a) $$

$$ \begin{aligned} & (U_k k_x - i\sigma) \left( \frac{\delta z_k}{2} (-k_x^2 - k_y^2)\phi_k \right) + (U_{k+1} k_x - i\sigma) \left( \frac{\delta z_{k+1}}{2} (-k_x^2 - k_y^2)\phi_{k+1} \right) \\ & \qquad + (U_{k+1/2} k_x - i\sigma) \left( \frac{f_0^2}{N_{k+1}^2} \frac{\phi_{k+3/2} - \phi_{k+1/2}}{\delta z_{k+1}} - \frac{f_0^2}{N_k^2} \frac{\phi_{k+1/2} - \phi_{k-1/2}}{\delta z_k} \right) \\ & \qquad + k_x \left( \frac{\delta z_k}{2} \phi_k + \frac{\delta z_{k+1}}{2} \phi_{k+1} \right) \beta \\ & \qquad + k_x \phi_{k+1/2} \left( - \frac{f_0^2}{N_{k+1}^2} \frac{U_{k+3/2} - U_{k+1/2}}{\delta z_{k+1}} + \frac{f_0^2}{N_k^2} \frac{U_{k+1/2} - U_{k-1/2}}{\delta z_k} \right) = 0 \end{aligned} \quad (A6b) $$

$$ \begin{aligned} & (U_K k_x - i\sigma) \left( \frac{\delta z_K}{2} (-k_x^2 - k_y^2)\phi_K - \frac{f_0^2}{N_K^2} \frac{\phi_{K+1/2} - \phi_{K-1/2}}{\delta z_K} \right) \\ & \qquad + k_x \phi_K \left( \frac{\delta z_K}{2} \beta + \frac{f_0^2}{N_K^2} \frac{U_{K+1/2} - U_{K-1/2}}{\delta z_K} \right) = 0. \end{aligned} \quad (A6c) $$

We take U at half levels as known, evaluated analytically.

### A3. Linearized HPE Equations on a CP-Grid

$$ \begin{aligned} & \sigma u_k + i k_x U_k u_k + (U_y)_k v_k + \frac{1}{2} \left( (U_z)_{k+1/2} w_{k+1/2} + (U_z)_{k-1/2} w_{k-1/2} \right) \\ & \qquad - f v_k + \frac{1}{\rho_0} i k_x p_k = 0 \end{aligned} \quad (A7a) $$

$$ \sigma v_k + i k_x U_k v_k + f u_k + \frac{1}{\rho_0} i k_y p_k = 0 \quad (A7b) $$

$$ \frac{1}{\rho_0} \frac{p_{k+1} - p_k}{\delta z} - g \alpha_T \theta_{k+1/2} = 0 \quad (A7c) $$

$$ i k_x u_k + i k_y v_k + \frac{w_{k+1/2} - w_{k-1/2}}{\delta z} = 0 \quad (A7d) $$
$$ \sigma\theta_{k+1/2} + (T_y)_{k+1/2} \frac{1}{2}(v_{k+1} + v_k) + ik_x U_{k+1/2}\theta_{k+1/2} + (T_z^b + T_z)_{k+1/2}w_{k+1/2} = 0 \quad (A7e) $$

All equations are defined for $k \in [1, K]$, with $w_{1/2}$ is set to 0 in (Equation A7a) when $k = 1$ and with the exception of (Equations A7c and A7e) which are defined for $k \in [1, K - 1]$. Then, supplemental equations for boundary conditions are added: $w_{K+1/2} = 0$ and $v_K + 1 = v_K$.

## A4. Linearized HPE Equations on a L-Grid

$$ \sigma u_k + ik_x U_k u_k + (U_y)_k v_k + \frac{1}{2}((U_z)_{k+1/2} w_{k+1/2} + (U_z)_{k-1/2} w_{k-1/2}) -fv_k + \frac{1}{\rho_0}ik_x p_k = 0 \quad (A8a) $$

$$ \sigma v_k + ik_x U_k v_k + fu_k + \frac{1}{\rho_0} ik_y p_k = 0 \quad (A8b) $$

$$ \frac{1}{\rho_0} \frac{p_{k+1} - p_k}{\delta z} - g\alpha_T \frac{1}{2}(\theta_{k+1} + \theta_k) = 0 \quad (A8c) $$

$$ ik_x u_k + ik_y v_k + \frac{w_{k+1/2} - w_{k-1/2}}{\delta z} = 0 \quad (A8d) $$

$$ \sigma\theta_k + (T_y)_k v_k + ik_x U_k \theta_k + \frac{1}{2}((T_z^b + T_z)_{k+1/2} w_{k+1/2} + (T_z^b + T_z)_{k-1/2} w_{k-1/2}) = 0 \quad (A8e) $$

All equations are defined for $k \in [1, K]$, with $w_{1/2}$ is set to 0 in (Equation A8a) when $k = 1$ and with the exception of (Equation A8c), which is defined for $k \in [1, K - 1]$. Then, a supplemental equation for a boundary condition is added: $w_{K+1/2} = 0$.

## Data Availability Statement

The CROCO (Coastal and Regional Ocean Community) model code is available at Auclair et al. (2024).

## Acknowledgments

We acknowledge the support of the French LEFE (Les Enveloppes Fluides et l'Environnement), CLIMAGO (Etude du système Climatique et de ses composantes Atmosphère, Glace et Océan) program through project DIDON (DIssipation Des tOurbillons et des oNdes internes dans l'océan). We would like to thank Aurélien Ponte for helpful discussions and help with LLC4320 data processing, and Jonathan Gula and Dimitris Menemenlis for making available some fields from GIGATL3 and LLC4320 simulations. We thank M. J. Bell and two other anonymous reviewers for their constructive reviews, which helped us to consolidate the paper.

## References

Arakawa, A., & Konor, C. S. (1996). Vertical differencing of the primitive equations based on the Charney–Phillips grid in hybrid &sigma;-p vertical coordinates. *Monthly Weather Review*, **124**(3), 511–528.

Arakawa, A., & Moorthi, S. (1988). Baroclinic instability in vertically discrete systems. *Journal of the Atmospheric Sciences*, **45**(11), 1688–1708. https://doi.org/10.1175/1520-0469(1988)045<1688:BIIVDS>2.0.CO;2

Auclair, F., Benshila, R., Bordois, L., Boutet, M., Brémond, M., Caillaud, M., et al. (2024). Coastal and regional ocean community model. *Zenodo*. https://doi.org/10.5281/zenodo.13898339

Barham, W., Bachman, S., & Grooms, I. (2018). Some effects of horizontal discretization on linear baroclinic and symmetric instabilities. *Ocean Modelling*, **125**, 106–116. https://doi.org/10.1016/j.ocemod.2018.03.004

Bell, M. J. (2003). Conservation of potential vorticity on lorenz grids. *Monthly Weather Review*, **131**(7), 1498–1501. https://doi.org/10.1175/1520-0493(2003)131<1498:copvol>2.0.co;2

Bell, M. J., & White, A. A. (1988). Spurious stability and instability in n-level quasi-geostrophic models. *Journal of the Atmospheric Sciences*, **45**(11), 1731–1738. https://doi.org/10.1175/1520-0469(1988)045<1731:ssaii>2.0.co;2

Bell, M. J., & White, A. A. (2017). Analytical approximations to spurious short-wave baroclinic instabilities in ocean models. *Ocean Modelling*, **118**, 31–40. https://doi.org/10.1016/j.ocemod.2017.08.001

Boccaletti, G., Ferrari, R., & Fox-Kemper, B. (2007). Mixed layer instabilities and restratification. *Journal of Physical Oceanography*, **37**(9), 2228–2250. https://doi.org/10.1175/jpo3101.1

Capet, X., Roullet, G., Klein, P., & Maze, G. (2016). Intensification of upper-ocean submesoscale turbulence through charney baroclinic instability. *Journal of Physical Oceanography*, **46**(11), 3365–3384. https://doi.org/10.1175/jpo-d-16-0050.1
Chassignet, E. P., & Xu, X. (2017). Impact of horizontal resolution (1/12 to 1/50) on gulf stream separation, penetration, and variability. *Journal of Physical Oceanography*, 47(8), 1999–2021. https://doi.org/10.1175/jpo-d-17-0031.1

Chassignet, E. P., & Xu, X. (2021). On the importance of high-resolution in large-scale ocean models. *Advances in Atmospheric Sciences*, 38(10), 1621–1634. https://doi.org/10.1007/s00376-021-0385-7

Ducousso, N., Le Sommer, J., Molines, J.-M., & Bell, M. (2017). Impact of the “symmetric instability of the computational kind” at mesoscale-and submesoscale-permitting resolutions. *Ocean Modelling*, 120, 18–26. https://doi.org/10.1016/j.ocemod.2017.10.006

Fox-Kemper, B., Adcroft, A., Böning, C. W., Chassignet, E. P., Curchitser, E., Danabasoglu, G., et al. (2019). Challenges and prospects in ocean circulation models. *Frontiers in Marine Science*, 6, 65. https://doi.org/10.3389/fmars.2019.00065

Gill, A. E. (1982). *Atmosphere-ocean dynamics* (Vol. 30). Academic Press.

Gula, J., Theetten, S., Cambon, G., & Roullet, G. (2021). Description of the gigat simulations. *Zenodo*. https://doi.org/10.5281/zenodo.4948523

Hallberg, R. (2005). A thermobaric instability of Lagrangian vertical coordinate ocean models. *Ocean Modelling*, 8(3), 279–300. https://doi.org/10.1016/j.ocemod.2004.01.001

Hochet, A., Huck, T., & De Verdière, A. C. (2015). Large-scale baroclinic instability of the mean oceanic circulation: A local approach. *Journal of Physical Oceanography*, 45(11), 2738–2754. https://doi.org/10.1175/jpo-d-15-0084.1

Konor, C. S., & Arakawa, A. (1997). Design of an atmospheric model based on a generalized vertical coordinate. *Monthly Weather Review*, 125(7), 1649–1673. https://doi.org/10.1175/1520-0493(1997)125<1649:doaamb>2.0.co;2

Ménesguen, C., Le Gentil, S., Marchesiello, P., & Ducousso, N. (2018). Destabilization of an oceanic meddy-like vortex: Energy transfers and significance of numerical settings. *Journal of Physical Oceanography*, 48(5), 1151–1168. https://doi.org/10.1175/jpo-d-17-0126.1

Molemaker, M. J., McWilliams, J. C., & Yavneh, I. (2005). Baroclinic instability and loss of balance. *Journal of Physical Oceanography*, 35(9), 1505–1517. https://doi.org/10.1175/jpo2770.1

Nelson, A., Arbic, B. K., Menemenlis, D., Peltier, W. R., Alford, M. H., Grisouard, N., & Klymak, J. M. (2020). Improved internal wave spectral continuum in a regional ocean model. *Journal of Geophysical Research: Oceans*, 125(5), e2019JC015974. https://doi.org/10.1029/2019jc015974

Rocha, C. B., Chereskin, T. K., Gille, S. T., & Menemenlis, D. (2016). Mesoscale to submesoscale wavenumber spectra in drake passage. *Journal of Physical Oceanography*, 46(2), 601–620.

Shchepetkin, A. F., & McWilliams, J. C. (2003). A method for computing horizontal pressure-gradient force in an oceanic model with a nonaligned vertical coordinate. *Journal of Geophysical Research*, 108(C3). https://doi.org/10.1029/2001jc001047

Smith, K. S. (2007). The geography of linear baroclinic instability in earth’s oceans. *Journal of Marine Research*, 65(5), 655–683. https://doi.org/10.1357/002224007783649484

Soufflet, Y., Marchesiello, P., Lemarié, F., Jouanno, J., Capet, X., Debreu, L., & Benshila, R. (2016). On effective resolution in ocean models. *Ocean Modelling*, 98, 36–50. https://doi.org/10.1016/j.ocemod.2015.12.004

Stewart, A. L., Klocker, A., & Menemenlis, D. (2018). Circum-antarctic shoreward heat transport derived from an eddy-and tide-resolving simulation. *Geophysical Research Letters*, 45(2), 834–845. https://doi.org/10.1002/2017gl075677

Stewart, K., Hogg, A. M., Griffies, S., Heerdegen, A., Ward, M., Spence, P., & England, M. H. (2017). Vertical resolution of baroclinic modes in global ocean models. *Ocean Modelling*, 113, 50–65. https://doi.org/10.1016/j.ocemod.2017.03.012

Uchida, T., Le Sommer, J., Stern, C., Abernathey, R. P., Holdgraf, C., Albert, A., et al. (2022). Cloud-based framework for inter-comparing submesoscale-permitting realistic ocean models. *Geoscientific Model Development*, 15(14), 5829–5856. https://doi.org/10.5194/gmd-15-5829-2022

Umlauf, L., & Burchard, H. (2003). A generic length-scale equation for geophysical turbulence models. *Journal of Marine Research*, 61(2), 235–265. https://doi.org/10.1357/002224003322005087

Vallis, G. K. (2006). *Atmospheric and oceanic fluid dynamics: Fundamentals and large-scale circulation*. Cambridge University Press.

Wenegrat, J. O., Callies, J., & Thomas, L. N. (2018). Submesoscale baroclinic instability in the bottom boundary layer. *Journal of Physical Oceanography*, 48(11), 2571–2592. https://doi.org/10.1175/jpo-d-17-0264.1