# Vertical Exchange Induced by Mixed Layer Instabilities

Yangcheng Luo, Jörn Callies

► To cite this version:

Yangcheng Luo, Jörn Callies. Vertical Exchange Induced by Mixed Layer Instabilities. Journal of Physical Oceanography, 2023, 53, pp.2701-2716. 10.1175/JPO-D-23-0059.1. insu-04851310

HAL Id: insu-04851310

https://insu.hal.science/insu-04851310v1

Submitted on 20 Dec 2024

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
# Vertical Exchange Induced by Mixed Layer Instabilities

YANGCHENG LUO<sup>a,b</sup> AND JÖRN CALLIES<sup>a</sup>

(Manuscript received 4 April 2023, in final form 26 September 2023, accepted 28 September 2023)

**ABSTRACT:** Submesoscale turbulence in the upper ocean consists of fronts, filaments, and vortices that have horizontal scales on the order of 100 m to 10 km. High-resolution numerical simulations have suggested that submesoscale turbulence is associated with strong vertical motion that could substantially enhance the vertical exchange between the thermocline and mixed layer, which may have an impact on marine ecosystems and climate. Theoretical, numerical, and observational work indicates that submesoscale turbulence is energized primarily by baroclinic instability in the mixed layer, which is most vigorous in winter. This study demonstrates how such mixed layer baroclinic instabilities induce vertical exchange by drawing filaments of thermocline water into the mixed layer. A scaling law is proposed for the dependence of the exchange on environmental parameters. Linear stability analysis and nonlinear simulations indicate that the exchange, quantified by how much thermocline water is entrained into the mixed layer, is proportional to the mixed layer depth, is inversely proportional to the Richardson number of the thermocline, and increases with increasing Richardson number of the mixed layer. The results imply that the tracer exchange between the thermocline and mixed layer is more efficient when the mixed layer is thicker, when the mixed layer stratification is stronger, when the lateral buoyancy gradient is stronger, and when the thermocline stratification is weaker. The scaling suggests vigorous exchange between the permanent thermocline and deep mixed layers in winter, especially in mode water formation regions.

**SIGNIFICANCE STATEMENT:** This study examines how instabilities in the surface layer of the ocean bring interior water up from below. This interior-surface exchange can be important for dissolved gases such as carbon dioxide and oxygen as well as nutrients fueling biological growth in the surface ocean. A scaling law is proposed for the dependence of the exchange on environmental parameters. The results of this study imply that the exchange is particularly strong if the well-mixed surface layer is thick, lateral density gradients are strong (such as at fronts), and the stratification below the surface layer is weak. These theoretical findings can be implemented in boundary layer parameterization schemes in global ocean models and improve our understanding of the marine ecosystem and how the ocean mediates climate change.

**KEYWORDS:** Instability; Ocean dynamics; Transport; Oceanic mixed layer; Thermocline; Ocean models

## 1. Introduction

Submesoscale turbulence, encompassing fronts, filaments, and vortices with horizontal scales typically ranging from 100 m to 10 km, is an important part of the upper ocean's circulation (Thomas et al. 2008; McWilliams 2016; Taylor and Thompson 2023). Numerical modeling has shown that submesoscale dynamics can locally generate vertical motion up to $100 \text{ m day}^{-1}$—one order of magnitude stronger than that generated by mesoscale eddies (Mahadevan and Tandon 2006; Thomas et al. 2008; Capet et al. 2008). Large vertical velocities at the base of the mixed layer can lead to vertical transport of tracers and particles, such as nutrients, dissolved inorganic carbon, oxygen, and particulate organic carbon, between the surface and interior ocean. The vertical exchange can have significant climatological and biological effects (Mahadevan and Archer 2000; Lévy et al. 2001; Thomas et al. 2008; Lévy et al. 2012a,b; Mahadevan 2014; Omand et al. 2015; Mahadevan 2016; Balwada et al. 2018; Lévy et al. 2018). On one hand, strong downwelling takes heat and particulate organic carbon into the interior ocean (Omand et al. 2015; Balwada et al. 2018).

and the resulting acceleration in air-sea exchange has important implications in the context of global warming. Such subduction of surface water can also export phytoplankton from the surface ocean into the interior, thus reducing primary productivity (Lévy et al. 2012a; Mahadevan 2014). On the other hand, it has been proposed that enhanced upwelling speeds up primary production by bringing nutrient-replete interior water into the euphotic surface layer, which overcomes surface oligotrophic conditions and leads to phytoplankton growth (Mahadevan and Archer 2000; Lévy et al. 2001, 2012a; Mahadevan 2016), although some studies have cast doubt on the importance of submesoscale motion in this process (Ascani et al. 2013; Lévy et al. 2018).

Submesoscale turbulence can be energized by a number of mechanisms, such as mixed layer instabilities, mesoscale strain-induced frontogenesis, or topographic wakes (e.g., Thomas et al. 2008; McWilliams 2016; Taylor and Thompson 2023). This paper focuses on the energization of balanced submesoscale motion by mixed layer instability (MLI; Boccaletti et al. 2007; Callies et al. 2016), which invigorates submesoscale turbulence in deep winter mixed layers (Mensa et al. 2013; Sasaki et al. 2014; Callies et al. 2015). Like baroclinic instability in the thermocline that energizes mesoscale eddies, MLI energizes submesoscale eddies by tapping into potential energy stored in the mixed layer in the presence of a lateral buoyancy gradient.

<footnote><sup>a</sup> Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, California</footnote>

<footnote><sup>b</sup> LMD/IPSL, Sorbonne Université, ENS, PSL, École Polytechnique, Institut Polytechnique de Paris, CNRS, Paris, France</footnote>

<footnote>*Corresponding author:* Yangcheng Luo, yangcheng.luo@lmd.ipsl.fr</footnote>
According to QG theory (Eady 1949), which can be used to roughly characterize MLI (Boccaletti et al. 2007; Callies et al. 2016), the most unstable horizontal scale lies around the deformation radius of the mixed layer, $N_m h_m/f$, where $N_m$ and $h_m$ are the buoyancy frequency and thickness of the mixed layer, respectively, and $f$ is the Coriolis parameter. Under Richardson number $Ri \equiv N^2/\Lambda^2 \sim 1$ conditions (where $\Lambda$ is the vertical shear of the horizontal mean flow), the instability scale is larger than predicted by QG theory by a factor $\sqrt{(1+Ri)/Ri}$ (Stone 1966). Under typical conditions of mixed layer stratification and thickness, the deformation radius of the mixed layer falls in the submesoscale range. QG theory further suggests that the growth rate of MLI scales with the Eady growth rate of the mixed layer, $f\Lambda/N_m$. This scaling is also slightly changed by non-QG effects under $Ri \sim 1$ conditions, modifying the growth rate by a factor $\sqrt{Ri/(1+Ri)}$ (Stone 1966). Weak stratification in the mixed layer leads to large growth rates, which is consistent with the time scales of submesoscale turbulence—on the order of a few days (Callies et al. 2020). Due to variations in the mixed layer depth and available potential energy, the strength of MLI undergoes a seasonal cycle, consistent with the observed seasonality in submesoscale turbulence (Callies et al. 2015; Qiu et al. 2017; Soares et al. 2022). The observed vertical structure of submesoscale energy, which is enhanced throughout the mixed layer and decays below its base, is also consistent with energization by MLI (Callies et al. 2015).

Submesoscale turbulence powered by MLI plays two related roles in the upper ocean: restratifying the mixed layer and potentially speeding up the exchange between the thermocline and mixed layer. The former has received significant attention (e.g., Boccaletti et al. 2007; Mahadevan et al. 2010) and a parameterization for coarse models has been proposed (Fox-Kemper et al. 2008; Fox-Kemper and Ferrari 2008; Fox-Kemper et al. 2011). In the parameterization, a scaling law for an overturning streamfunction that restratifies the mixed layer is presented. The streamfunction is proportional to $\Lambda h_m^2$ (Fox-Kemper et al. 2008), although it is not clear how universal this scaling is (Callies and Ferrari 2018b). The Fox-Kemper et al. (2008) streamfunction also advects tracers, but the parameterization does not produce any exchange between the thermocline and mixed layer because the streamfunction vanishes at the base of the mixed layer. The latter, i.e., MLI's role in vertical exchange, however, has received less attention. A drastic increase in vertical velocities is typically seen in numerical simulations when submesoscales are resolved (e.g., Capet et al. 2008), with a vertical structure as expected from MLI (Callies et al. 2016). But it remains unclear whether and how this increase in vertical velocity enhances vertical exchange.

As MLI is affected by various environmental conditions, such as the stratification of the mixed layer and thermocline, the horizontal buoyancy gradient, and the vertical shear of the mean flow, it is expected that these environmental parameters also influence the vertical exchange caused by MLI-induced submesoscale motion. In this study, we first demonstrate how MLI induces vertical exchange between the thermocline and mixed layer, and then explore the parameter dependence of

the exchange. The physics revealed in this work can potentially be added to state-of-the-art boundary layer parameterization schemes in global ocean models (e.g., Large et al. 1994; Fox-Kemper et al. 2008; Bachman et al. 2017; Van Roekel et al. 2018).

This article consists of two parts, linear stability analysis and nonlinear simulations, which will be discussed in the next two sections, respectively. Linear stability analysis captures the onset of MLI and demonstrates how it induces vertical motion and deforms the interface between the mixed layer and thermocline. Nonlinear simulations show how the deformation of the layer interface eventually forms filaments that transport thermocline water into the mixed layer, causing tracer exchange (Figs. 1 and 6).

## 2. Linear stability analysis

### a. Model formulation

We formulate the simplest model for MLI in a controlled setting, a two-layer Eady model (Fig. 2; Eady 1949). The upper layer represents the mixed layer, with depth $h_m$ and buoyancy frequency $N_m$. The lower layer represents the thermocline, with depth $h_t$ and buoyancy frequency $N_t$ ($N_t > N_m$).

For simplicity, we impose the same horizontal buoyancy gradient (only in the $y$ direction) for the two layers, and we assume that the buoyancy is continuous at a flat interface between the two layers. With $B$ representing the mean buoyancy, we assume $\partial B/\partial y = -f\Lambda$. In the $v-z$ plane, the sloping isopycnals are steep in the mixed layer and relatively flat in the thermocline. Conforming to the thermal wind relation, a mean vertical shear $\Lambda$ exits in both layers. We assume that the sea surface is located at $z = h_m/2$ and apply the rigid-lid approximation. These mean flow conditions in the linear stability analysis will also be the initial conditions in the nonlinear simulations presented in the next section. It is worth noting that a thick thermocline layer can isolate the MLI from the mesoscale instability (Callies et al. 2016), but for the sake of computational feasibility, we add a flat, rigid bottom at the base of a relatively thin thermocline layer. In the linear stability analysis, we set $h_t = 5h_m$. In the nonlinear simulations, we use $h_t = h_m$. We only present cases where MLI remains well isolated.

QG theory offers an approximate understanding of MLI (Callies et al. 2016). Due to the discontinuity of stratification across the interface in the presence of a horizontal buoyancy gradient, a sheet of quasigeostrophic potential vorticity (QGPV) is present at the layer interface, with the horizontal QGPV gradient in the $-y$ direction. This PV gradient gives rise to Rossby waves propagating in the $+x$ direction. The rigid surface acts as a layer with infinitely strong stratification. Due to the “effective” stratification jump across the surface, another QGPV sheet is present at the surface with the horizontal QGPV gradient in the $+y$ direction. This QGPV gradient gives rise to Rossby waves propagating in the $-x$ direction. Due to the vertical shear in the mean flow, the Rossby waves at the surface are Doppler shifted toward the $+x$ direction, and the Rossby waves at the interface are Doppler shifted toward the $-x$ direction.
<figure><img src="image_2.png" /><figcaption>FIG. 1. Thermocline water with a high tracer concentration and high potential vorticity (PV) is entrained into the mixed layer. Shown are snapshots of the vertical cross sections (in the x-z plane) of (a) the binary tracer concentration c and (b) the PV q. The initial interface between the mixed layer and thermocline is located at $z = -0.5h_m$. Data are taken from the nonlinear simulation with mixed layer Richardson number $Ri_m = 1$ and thermocline Richardson number $Ri_t = 100$ (simulation time $t = 89N_m/f\Lambda$). The spatial coordinates and PV are nondimensionalized, and their dimensions are shown in the parentheses. Details are described in section 3.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 2. Schematic of the idealized model of mixed layer instability. A deformable interface separates the mixed layer and thermocline. MLI can be isolated by making the thermocline layer infinitely deep, but here a solid bottom is assumed at the base of the thermocline. The black arrows show the mean flow in the x direction. The contours in the y-z plane represent isopycnals. Data are taken from one of the nonlinear simulations in this study ($Ri_m = 1$, $Ri_t = 100$, simulation time $t = 60N_m/f\Lambda$).</figcaption></figure>

direction. This allows the phase speeds of the two sets of Rossby waves to be similar, such that they can phase lock and mutually amplify, giving rise to the instability (see, e.g., Vallis 2017).

Rigorously speaking, this QG reasoning only applies to circumstances where the Rossby number $Ro = U/fL = 1/\sqrt{Ri} \ll 1$, where $U = \Lambda h_m$ is the characteristic speed of the horizontal flow, and $L$ is the characteristic horizontal length scale. In the ocean, the mixed layer is typically observed to be weakly stratified and have an $Ri_m \le 1$ (Thomas et al. 2013). The observed properties originate from two competing processes. On one hand, mixing induced by wind stress, negative surface buoyancy forcing, and cross-front Ekman transport induced by alongfront wind stress can potentially reduce the stratification of the mixed layer down to zero or even negative values. On the other hand, fast-growing gravitational and symmetric instabilities can restore $Ri_m$ back to 1 (Haine and Marshall 1998; Thomas et al. 2013). With $Ri_m$ restored to ~1, baroclinic instability becomes the fastest growing instability and further increases $Ri_m$ by restratifying the mixed layer (Stone 1966; Stone et al. 1969; Stone 1970; Haine and Marshall 1998; Boccaletti et al. 2007). To be able to capture the non-QG mixed layer dynamics with small $Ri$ during the early stage of the development of the mixed layer baroclinic instability, we will use the hydrostatic Boussinesq equations.

Linearized around the mean state described above, these equations read:

$$ \frac{\partial u}{\partial t} + \Lambda z \frac{\partial u}{\partial x} + w \Lambda - f v = - \frac{1}{\rho_0} \frac{\partial p}{\partial x}, \quad (1) $$

$$ \frac{\partial v}{\partial t} + \Lambda z \frac{\partial v}{\partial x} + f u = - \frac{1}{\rho_0} \frac{\partial p}{\partial y}, \quad (2) $$

$$ \frac{1}{\rho_0} \frac{\partial p}{\partial z} = b, \quad (3) $$

$$ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0, \quad (4) $$

$$ \frac{\partial b}{\partial t} + \Lambda z \frac{\partial b}{\partial x} - f \Lambda v + N^2 w = 0. \quad (5) $$

Here, $u$, $v$, and $w$ represent perturbation velocities in the $x$, $y$, and $z$ directions, respectively; $p$ and $b$ are perturbation pressure and buoyancy, respectively. The full buoyancy is $-g(\rho - \rho_0)/\rho_0 = N^2 z - f \Lambda y + b$, where $\rho$ is the density and $\rho_0$ a constant reference density. We ignore the compressibility of seawater. The buoyancy frequency $N$ is equal to $N_m$ in the mixed layer and $N_t$ in the thermocline. For the boundary conditions, we assume zero vertical velocities at the surface and at the base of the thermocline. At the layer interface, $w$ and $p$ are continuous. For simplicity, the hydrostatic approximation is adopted in the vertical momentum equation. Nonhydrostatic effects are insignificant for baroclinic instability unless the vertical shear $\Lambda$ is far smaller than the Coriolis parameter $f$ (Stone 1971; Mahadevan 2006; Callies and Ferrari 2018a).

We nondimensionalize the equations to reduce the number of model parameters and thus simplify the parameter space. We follow QG scaling in picking the deformation radius of the mixed layer $N_m h_m / f$ as the horizontal length scale and the mixed layer thickness $h_m$ as the vertical length scale. For the horizontal perturbation velocities, we pick $\Lambda h_m$, the difference in the mean flow between the surface and base of the mixed layer. We then choose the scales for all other variables based on dimensional relationships in geostrophic and hydrostatic balance: we pick $\Lambda h_m f / N_m$ for the perturbation vertical velocity, $\rho_0 \Lambda h_m^2 N_m$ for the perturbation pressure, $\Lambda h_m N_m$ for the perturbation buoyancy, and $N_m / f \Lambda$ for time. The nondimensionalization is summarized in Table 1. The resulting dimensionless equations are
<table><caption>TABLE 1. Scales used for the nondimensionalization of the equations of motion. These choices follow QG scaling based on mixed layer quantities (denoted with subscript *m*).</caption><thead><tr><th>Variable name(s)</th><th>Symbol(s)</th><th>Scale</th></tr></thead><tbody><tr><td>Horizontal coordinates</td><td>x, y</td><td>N<sub>m</sub>h<sub>m</sub>/f</td></tr><tr><td>Vertical coordinate</td><td>z</td><td>h<sub>m</sub></td></tr><tr><td>Time</td><td>t</td><td>N<sub>m</sub>/fΛ</td></tr><tr><td>Horizontal velocities</td><td>u, v</td><td>Λh<sub>m</sub></td></tr><tr><td>Vertical velocity</td><td>w</td><td>fΛh<sub>m</sub>/N<sub>m</sub></td></tr><tr><td>Pressure perturbation</td><td>p</td><td>ρ<sub>0</sub>Λh<sub>m</sub><sup>2</sup>N<sub>m</sub></td></tr><tr><td>Buoyancy perturbation</td><td>b</td><td>Λh<sub>m</sub>N<sub>m</sub></td></tr></tbody></table>

$$ \frac{1}{\sqrt{Ri_m}} \left( \frac{\partial u'}{\partial t'} + z' \frac{\partial u'}{\partial x'} + w' \right) - v' = - \frac{\partial p'}{\partial x'}, \quad (6) $$

$$ \frac{1}{\sqrt{Ri_m}} \left( \frac{\partial v'}{\partial t'} + z' \frac{\partial v'}{\partial x'} \right) + u' = - \frac{\partial p'}{\partial y'}, \quad (7) $$

$$ \frac{\partial p'}{\partial z'} = b', \quad (8) $$

$$ \frac{\partial u'}{\partial x'} + \frac{\partial v'}{\partial y'} + \frac{\partial w'}{\partial z'} = 0, \quad (9) $$

$$ \frac{1}{\sqrt{Ri_m}} \left( \frac{\partial b'}{\partial t'} + z' \frac{\partial b'}{\partial x'} - v' \right) + Sw' = 0, \quad (10) $$

where primes denote dimensionless variables and $S = 1$ in the mixed layer and $S = Ri_t/Ri_m$ in the thermocline. The boundary condition is $w' = 0$ at $z' = 1/2$ and at $z' = -1/2 - h_t/h_m$.

The dimensionless equations show that the problem is determined by three dimensionless numbers—the Richardson number of the mixed layer $Ri_m$, the Richardson number of the thermocline $Ri_t$, and the ratio of the thermocline to mixed layer thickness $\gamma = h_t/h_m$. The equations above are transformed into an eigenvalue problem for a complex frequency $\omega$ by plugging in a wave ansatz $u' = \mathcal{R}\hat{u}(z)\exp[i(kx + ly - \omega t)]$, etc. for $u'$, $v'$, $w'$, $p'$, and $b'$, where $k$ and $l$ are the horizontal wavenumbers in the $x$ and $y$ directions, respectively, and $\mathcal{R}$ denotes taking the real part. We then use Dedalus, an open-source code that solves differential equations using spectral methods (Burns et al. 2020), to solve the discretized equations and obtain a numerical solution for the dispersion relation $\omega = \omega(k, l; Ri_m, Ri_t, \gamma)$. The imaginary part of $\omega$, denoted as $\sigma = \mathcal{I}\omega$, is the exponential growth rate of the instability. In addition to the dispersion relation, the vertical structures of the growing modes, namely, $\hat{u}(z)$, $\hat{v}(z)$, $\hat{w}(z)$, $\hat{p}(z)$, and $\hat{b}(z)$, can also be obtained. These inform where MLI-induced motion is the strongest and how energy is distributed vertically.

Experimentation shows that $\gamma = 5$ is sufficient for isolating MLI from the influence of the rigid bottom of the thermocline layer, so we use $\gamma = 5$ in all the linear stability analysis calculations. We resolve the mixed layer, located between $z = 0.5h_m$ and $z = -0.5h_m$, and the thermocline layer, located between $z = -0.5h_m$ and $z = -5.5h_m$, with 128 Chebyshev modes each (Boyd 2001; Burns et al. 2020). In both layers, the vertical

spacing is therefore uneven—the spacing is finer near the surface, near the layer interface, and near the bottom of the thermocline. It is coarser in the interior of each layer.

## b. Results

As shown by Stone (1966), at a given $k$, the largest growth rate takes place at $l = 0$, so we consider that case only. The growth rate curve $\sigma(k)$ (Fig. 3a) is similar to that in Callies et al. (2016). The peak at smaller wavenumbers corresponds to the thermocline baroclinic instability, produced by the QGPV gradient reversal between the surface and the base of the thermocline layer. The maximum growth rate takes place near the Rossby deformation radius of the mixed layer, ranging from $1.1N_mh_m/f$ to $1.6N_mh_m/f$, depending on $Ri$ of the two layers. The maximum growth rate falls between $0.23f\Lambda/N_m$ and $0.31f\Lambda/N_m$. As typically $N_t h_t \gg N_m h_m$, the two growth rate peaks are well separated from each other in terms of spatial scales. Also, as typically $N_t \gg N_m$, MLI outcompetes the thermocline instability in terms of the maximum growth rate.

Figure 3 also shows typical spatial structures associated with the fastest-growing MLI modes in the $x-z$ plane. The flow is largely confined to the mixed layer, from $z = -0.5h_m$ to $z = 0.5h_m$, and the modes tilt in the $-x$ direction. The largest amplitudes of $u$ and $v$ are located at the surface ($z = 0.5h_m$) and at the layer interface ($z = -0.5h_m$; Figs. 3b,c). In contrast, the largest amplitude of $w$ takes place near the midplane of the mixed layer ($z = 0$; Fig. 3d). Importantly, $u$, $v$, and $w$ all extend into the upper part of the thermocline. In particular, $w$ at the layer interface can cause a growing vertical displacement of the interface, eventually resulting in vertical exchange between the thermocline and mixed layer. The dependence of the vertical structure of $w$ on $Ri_t$ is illustrated in Fig. 4. At the layer interface ($z = -0.5h_m$), the amplitude of $w$ decreases as $Ri_t$ increases, indicating that large thermocline Richardson number suppresses vertical motion at the layer interface.

The spatial structure of the buoyancy perturbation $b$ is different from that of $u$, $v$, and $w$ (Fig. 3e). The largest amplitude takes place beneath the layer interface. In the thermocline, the amplitude decays with depth, and the modes tilt toward the $-x$ direction. Within the mixed layer, the amplitude of $b$ is small, and its largest values take place at the surface and just above the interface. The modes in the mixed layer tilt toward the $+x$ direction. There is an abrupt phase shift in $b$ across the interface. For typical Richardson numbers of the mixed layer and thermocline, the phase shift is about $180^\circ$, due to the fact that the MLI-induced motion in the $y-z$ plane is inclined with a slope between that of mean isopycnals in the mixed layer and thermocline. The value of the phase shift depends on $Ri_m$ and $Ri_t$. Note that the discontinuity in $b$ across the interface does not imply a discontinuity in the full buoyancy; in fact, considering the background, the full buoyancy $N^2z - f\Lambda y + b$ is continuous across the interface.

The modal amplitudes shown in Figs. 3 and 4 are normalized such that the horizontal-mean column-integrated eddy energy is $\rho_0 \Lambda^2 h_m^3$. Such a normalization is needed to allow for a comparison of modal amplitudes at different
<figure><img src="image_4.png" /><figcaption>FIG. 3. Growth rates and structure of most unstable mode from the linear stability analysis. (a) The growth rate curve for $Ri_m = 1$ and $Ri_t = 100$. The unstable modes to the right of the MLI peak are the ageostrophic modes described in Stone (1970). (b)–(e) Spatial structures of the most unstable mode for $u$, $v$, $w$, and $b$, respectively. The horizontal axes show the distance in the $x$ direction multiplied by the most unstable wavenumber $k$. The vertical axes show depth. The mixed layer is located between $z = \pm 0.5h_m$. The thermocline reaches to $z = -5.5h_m$ but only the upper part is shown. The quantity against which a variable is nondimensionalized is shown in parentheses. The magnitudes of the variables are scaled such that the horizontally averaged energy in the water column is equal to $\rho_0 \Lambda^2 h_m^3$.</figcaption></figure>

parameters because the linear analysis itself leaves the amplitude of the perturbations undetermined. The eddy energy is defined as the sum of perturbation kinetic energy $K = \rho_0(u^2 + v^2)/2$ and eddy potential energy $P = \rho_0 b^2/2N^2$, where $N$ is evaluated with $N_m$ or $N_t$ accordingly in the two layers. The normalization is motivated physically: the energy associated with the perturbations is scaled to be comparable to the energy associated with the mean state in the mixed layer. This roughly marks the transition point from the linear phase of the instability development to its nonlinear phase. The relationship between the energy normalization and the transition point will be discussed further in section 3.

Figure 5 shows the dependence of the linear model on $Ri_m$ and $Ri_t$. The ranges of the Richardson numbers discussed here, $1 \le Ri_m \le 3000$, $60 \le Ri_t \le 10^5$, cover typical conditions in the real upper ocean. If $Ri_m$ and $Ri_t$ are too close to each

other—i.e., when the stratification jump across the interface between the mixed layer and thermocline is too small—MLI modes extend downward far below the base of the mixed layer and touch the bottom of the thermocline layer. To avoid such interaction with the artificial bottom at the base of the thermocline layer, we only present cases with $Ri_t \ge 30Ri_m$.

Figure 5a shows the dependence of the wavenumber $k_{\text{max}}$ of the most unstable mode on $Ri_m$ and $Ri_t$. In the middle right part of the figure, $Ri_m \sim 100$, and $Ri_t$ approaches $10^5$, so the mixed layer approaches the QG limit, and the base of the mixed layer behaves approximately like a rigid bottom. The most unstable wavenumber takes its maximum value of $k_{\text{max}} \approx 1.6f/N_m h_m$ in the middle right part of the figure, consistent with the QG result (Eady 1949). As $Ri_m$ decreases, the mixed layer moves away from the QG limit, and $k_{\text{max}}$ decreases toward $1.1f/N_m h_m$, consistent with Stone's (1966) formula for a single layer:
<figure><img src="image_5.png" /><figcaption>FIG. 4. Large thermocline Richardson number suppresses the vertical velocity at the interface between the mixed layer and thermocline. Shown are the vertical structures of the magnitude of the MLI-induced $w$. The interface is located at $z = -0.5h_m$. Data are taken from three linear stability calculations with fixed mixed layer Richardson number $Ri_m$ and varying thermocline Richardson number $Ri_t$. The magnitudes of the variables are scaled such that the horizontally averaged energy in the water column is equal to $\rho_0 \Lambda^2 h_m^3$.</figcaption></figure>

$$k_{\max} \approx \sqrt{\frac{5}{21} \frac{Ri_m}{1 + Ri_m} \frac{f}{N_m h_m}} \quad (11)$$

This behavior is observed for a broad range of $Ri_t$ and indicates that non-QG effects in general result in a larger spatial scale of the instability than predicted by QG dynamics. When $Ri_t$ approaches $Ri_m$, $k_{\max}$ also decreases due to the modification of the MLI modes by the thermocline as the MLI modes penetrate deeper into it.

The maximum growth rate $\sigma_{\max}$ depends primarily on $Ri_m$ (Fig. 5b). In the upper half of the figure, where $Ri_m \gg 1$, the result from the QG Eady model is reproduced: $\sigma_{\max} \approx 0.3 f \Lambda / N_m$. As $Ri_m$ decreases, the mixed layer moves away from the QG limit, and the maximum growth rate decreases, again consistent with the result in Stone (1966):

$$\sigma_{\max} \approx \sqrt{\frac{5}{541} \frac{Ri_m}{1 + Ri_m} \frac{f \Lambda}{N_m}} \quad (12)$$

This indicates that non-QG dynamics produce slower growth of baroclinic instability compared to the QG prediction. Note that this is a secondary effect on top of the change in the growth rate dimension $f \Lambda / N_m = f / \sqrt{Ri_m}$, which increases with decreasing $Ri_m$ at fixed $f$. Also note that the trends in $k_{\max}$ and $\sigma_{\max}$ are independent of the energy normalization.

With the column-integrated perturbation energy normalized to $\rho_0 \Lambda^2 h_m^3$, the nondimensional vertical velocity at the layer interface has a strong dependence on $Ri_t$ and a weak dependence on $Ri_m$ (Fig. 5c). Increasing $Ri_t$ from $10^2$ to $10^5$

reduces the $w$ amplitude at the interface from $0.1 \Lambda h_m f / N_m$ to $0.003 \Lambda h_m f / N_m$. The reduction of $w$ at the base of the mixed layer is unrelated to the amplitude of $w$ in the interior of the mixed layer, whose change is small (Fig. 4). Taking the scale of $w$ (which is $\Lambda h_m f / N_m = h_m f / \sqrt{Ri_m}$) into consideration, the amplitude of the dimensional $w$ decreases with increasing $Ri_m$ if $h_m$ and $f$ are held fixed.

Figure 5d provides another perspective for how a large $Ri_t$ suppresses vertical motion at the base of the mixed layer, showing the variance of the vertical displacement of the layer interface $\eta$. The magnitude of the interface deformation is a more straightforward indicator of the vertical exchange across the base of the mixed layer after the perturbation grows to finite amplitude than the magnitude of $w$. The amplification of the interface deformation is caused by the correlation between $w$ and $\eta$ rather than $w$ alone. The growth rate of the variance of the vertical displacement of the interface is

$$\frac{d\eta^2}{dt} = 2\overline{\eta w}, \quad (13)$$

where overbars denote a horizontal average. If $\eta$ and $w$ are in phase, the deformation of the interface grows. The part of $w$ that is out of phase with $\eta$ does not lead to an increase in the deformation of the interface and therefore does not contribute to eventual exchange between the two layers.

Like $w$ at the layer interface, the nondimensionalized (against $h_m^2$) variance of the interface displacement depends strongly on $Ri_t$ and weakly on $Ri_m$ (Fig. 5d). Increasing $Ri_t$ from $10^2$ to $10^5$ decreases $\overline{\eta^2}$ from $0.01 h_m^2$ to $5 \times 10^{-6} h_m^2$. But unlike that of $w$, the scale of $\overline{\eta^2}$ does not contain $Ri_m$, so the dimensional $\overline{\eta^2}$ has negligible $Ri_m$ dependence (with the column-integrated perturbation energy normalized to $\rho_0 \Lambda^2 h_m^3$). This is because the scale of time $N_m / f \Lambda = \sqrt{Ri_m} / f$ cancels out the $1 / \sqrt{Ri_m}$ in the scale of $w$.

It is worth emphasizing that the magnitudes of $w$ and $\overline{\eta^2}$ rely on the energy normalization, as stated above, which we choose to target the transition point from the linear to the nonlinear phase of the instability development. It is also worth noting that the $\sqrt{Ri_t}$ dependence of the interface displacement (Fig. 5d) should be understood as a qualitative measure of exchange between the thermocline and mixed layer because nonlinearity in the later stage of the instability development stretches the sinusoidal waves into filaments of thermocline water, which are subsequently entrained into the mixed layer (Fig. 1). The linear stability analysis only captures the onset of this filament formation.

The dimension of $\overline{\eta^2}$ and the parameter dependence of the nondimensional $\overline{\eta^2}$ nevertheless inform where and when strong exchange between the thermocline and mixed layer can be expected in the real ocean. The dimension $h_m^2$ implies that deep winter mixed layers, as observed in mode water formation regions in the Southern Ocean and western parts of subtropical gyres (e.g., de Boyer Montégut et al. 2004; Dong et al. 2008), are favorable for strong vertical exchange. The dependence of $\overline{\eta^2}$ on $Ri_t = N_t^2 / \Lambda^2$ implies that a weakly stratified thermocline, as observed again in mode water formation regions, and strong horizontal buoyancy gradient (i.e., strong mesoscale eddies), as observed
<figure><img src="image_6.png" /><figcaption>FIG. 5. Dependence of the most unstable linear mode on the Richardson numbers of the mixed layer and thermocline. The four panels show the dependence on the Richardson numbers of the mixed layer and thermocline $Ri_m$ and $Ri_t$ of (a) the most unstable wavenumber $k_{max}$, (b) the associated maximum growth rate $\sigma_{max}$, (c) the magnitude of $w$ at the layer interface, and (d) the variance of the vertical displacement of the interface. The magnitudes of the variables are scaled such that the horizontally averaged energy in the water column is equal to $\rho_0 \Lambda^2 h_m^3$.</figcaption></figure>

in the Southern Ocean and near the western boundary currents, are favorable for strong vertical exchange across the base of the mixed layer.

### 3. Nonlinear simulations

#### a. Model formulation

To characterize the parameter dependence of the vertical exchange in the nonlinear phase of the instability, we pose an initial-value problem, where the system freely evolves from an idealized initial state. The model configuration is similar to that used in the linear stability analysis—only $\gamma$ is reduced from 5 to 1 in order to concentrate more resolution in the mixed layer. Linear stability analysis shows that $\gamma = 1$ is sufficiently large that MLI modes in all of the nonlinear simulations presented here are well isolated. The simulation domain extends from $x = 0$ to $12N_m h_m/f$, from $y = 0$ to $12N_m h_m/f$, and from $z = -1.5h_m$ to $0.5h_m$. The domain is doubly periodic in the horizontal directions. In the initial state, the upper half of the domain (from $z = -0.5h_m$ to $0.5h_m$) represents the mixed layer with a small Richardson number $Ri_m$, and the lower half (from $z = -1.5h_m$ to $-0.5h_m$) represents the thermocline with a large $Ri_t$. We now use one Chebyshev basis for the full domain instead of the compound basis used in the linear

stability analysis because the location of the mixed layer base now changes over time.

The nonlinear terms in the equations of motion are restored (Taylor and Ferrari 2010; Callies and Ferrari 2018a):

$$ \frac{\partial u}{\partial t} + (\Lambda z + u) \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + w \left( \Lambda + \frac{\partial u}{\partial z} \right) - f v = - \frac{1}{\rho_o} \frac{\partial p}{\partial x} + \mathcal{F}(u), \quad (14) $$

$$ \frac{\partial v}{\partial t} + (\Lambda z + u) \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} + w \frac{\partial v}{\partial z} + f u = - \frac{1}{\rho_o} \frac{\partial p}{\partial y} + \mathcal{F}(v), \quad (15) $$

$$ \frac{1}{\rho_o} \frac{\partial p}{\partial z} = b, \quad (16) $$

$$ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0, \quad (17) $$

$$ \frac{\partial b}{\partial t} + (\Lambda z + u) \frac{\partial b}{\partial x} + v \left( -f\Lambda + \frac{\partial b}{\partial y} \right) + w \frac{\partial b}{\partial z} = \mathcal{F}(b), \quad (18) $$
where $\mathcal{F}(u)$, $\mathcal{F}(v)$, and $\mathcal{F}(b)$ represent terms that are used to smooth out gridscale gradients and stabilize the numerical simulations. We apply the same operator

$$ \mathcal{F} = -\mu \left( \frac{\partial^4}{\partial x^4} + 2 \frac{\partial^4}{\partial x^2 \partial y^2} + \frac{\partial^4}{\partial y^4} \right) + \nu \frac{\partial^2}{\partial z^2} \quad (19) $$

for momentum and buoyancy, using hyperdiffusion in the horizontal directions and regular diffusion in the vertical direction.

Note that the decomposition of buoyancy into background buoyancy $B = -\Lambda f y$ and perturbation buoyancy $b$ in the nonlinear simulations is different from that employed in the linear stability analysis, where $B = N^2 z - f \Lambda y$. Now, the initial $N^2 z$ is part of the “perturbation” field $b$, which allows for significant changes in the stratification as the instability evolves.

In addition to the five dynamical variables, we solve for a scalar field to trace water that is initially in the thermocline. The tracer concentration $c$ is advected by the background plus perturbation flow and diffused using the same operator as applied for momentum and buoyancy:

$$ \frac{\partial c}{\partial t} + (\Lambda z + u) \frac{\partial c}{\partial x} + v \frac{\partial c}{\partial y} + w \frac{\partial c}{\partial z} = \mathcal{F}(c). \quad (20) $$

For boundary conditions, we apply $w = 0$ at the top of the mixed layer and at the bottom of the thermocline. No perturbation surface or bottom stress is assumed: $\partial u/\partial z = \partial v/\partial z = 0$. This implies a weak zonal stress $v \Lambda$ that maintains the background flow against viscous dissipation. Similarly, we apply no-flux conditions on buoyancy and the tracer: $\partial b/\partial z = 0$, $\partial c/\partial z = 0$.

The initial state of the model consists of the background flow $U = \Lambda z$ with the perturbation flow $u = v = w = 0$. On top of a background lateral buoyancy gradient $\partial B/\partial y = -f \Lambda$, the vertical gradient of the perturbation buoyancy $b$ is initialized using a tanh function in order to mitigate the Gibbs phenomenon at the stratification jump:

$$ \frac{\partial b}{\partial z} = \frac{N_m^2}{2} \left[ \left( 1 - \frac{R_i}{R_i_m} \right) \tanh \frac{z + 0.5h_m}{z_0} + \left( 1 + \frac{R_i}{R_i_m} \right) \right], \quad (21) $$

where $z_0 = 0.02h_m$ is an arbitrarily chosen parameter that sets the transition scale between the mixed layer and thermocline. Away from the layer interface, $\partial b/\partial z \simeq N_m^2$ in the mixed layer, and $\partial b/\partial z \simeq N_t^2$ in the thermocline. The initial buoyancy is determined by vertically integrating (21). Similarly, we used a tanh function to initialize the tracer concentration $c$:

$$ c = \frac{1}{2} \left( 1 - \tanh \frac{z + 0.5h_m}{z_0} \right), \quad (22) $$

so that $c$ is set to 1 (dimensionless) in the thermocline and 0 in the mixed layer. We will from now on refer to $c$ as the “binary tracer.”

The nonlinear equations are nondimensionalized as in the linear stability analysis (Table 1). The viscosities $\mu$ and $\nu$ are nondimensionalized as follows: $\mu = \Lambda N_m^3 h_m^4 / f^3 \text{Re}$, $\nu = \Lambda h_m^2 f / N_m \text{Re}$,

<table><caption>TABLE 2. Richardson numbers used in the nonlinear simulations. The Richardson number is defined as the ratio of the buoyancy frequency squared to the vertical shear of the horizontal flow squared, $R_i \equiv N^2/\Lambda^2$. $R_{im}$ is the mixed layer Richardson number, and $R_i$ is the thermocline Richardson number.</caption><thead><tr><th>Ri<sub>m</sub></th><th>Ri<sub>t</sub></th></tr></thead><tbody><tr><td>1</td><td>10<sup>2</sup></td></tr><tr><td>1</td><td>10<sup>3</sup></td></tr><tr><td>1</td><td>10<sup>4</sup></td></tr><tr><td>10</td><td>10<sup>3</sup></td></tr><tr><td>10</td><td>10<sup>4</sup></td></tr><tr><td>100</td><td>10<sup>4</sup></td></tr></tbody></table>

where Re is the Reynolds number. Re should be as large as possible to minimize the effect of these terms. We set Re = $10^4$, the largest possible Re that does not incur checkerboard instability at the resolution used. The dimensionless $F'$ then takes the form

$$ \mathcal{F}' = - \frac{1}{\sqrt{R_i} m} \frac{1}{\text{Re}} \left( \frac{\partial^4}{\partial x'^4} + 2 \frac{\partial^4}{\partial x'^2 \partial y'^2} + \frac{\partial^4}{\partial y'^4} \right) + \frac{1}{\sqrt{R_i} m} \frac{1}{\text{Re}} \frac{\partial^2}{\partial z'^2}. \quad (23) $$

We use Dedalus to solve the initial value problem with 192 Fourier modes in the $x$ and $y$ directions and 256 Chebyshev modes in the $z$ direction. Translated into the physical space, these are equivalent to 192 evenly spaced grid points in the $x$ and $y$ directions and 256 unevenly spaced grid points in the $z$ direction. The linear stability analysis shows that the wavelength of mixed layer instability in the $x$ direction ranges from 3.9 to $5.7 N_m h_m / f$. Therefore, with 192 Fourier modes resolving $12 N_m h_m / f$ in the horizontal directions, each wavelength is resolved by 62–91 grid points. The vertical grid points are denser near the surface and near the bottom of the thermocline layer and are sparser in the center of the domain, where the base of the mixed layer is located in the initial state. This resolution was sufficient to close the PV budget (see below) for all simulations.

We add white noise to the initial perturbation buoyancy field $b$ to seed the instability. The noise at all grid points follows a Gaussian distribution with a mean value of zero and a standard deviation of $10^{-3}$. All simulations are run to $t = 100 N_m / f \Lambda$.

We simulate six combinations of $R_i$ and $R_t$ (Table 2). $R_{im}$ and $R_i$ not only determine the initial state but also the magnitude of deviation from QG dynamics. We keep $R_i \gg R_{im}$ to avoid MLI modes penetrating deep into the thermocline layer and interacting with the artificial bottom.

## b. Results

Every nonlinear simulation undergoes four stages of evolution (Fig. 6). Out of the initial white noise emerges the fastest exponentially growing mode predicted by the linear stability analysis in section 2, which is characterized by parallel stripes oriented in the $y$ direction ($l = 0$) with alternate positive and
<figure><img src="image_7.png" /><figcaption>FIG. 6. Snapshots from the nonlinear evolution of (top) the perturbation buoyancy b at the midplane of the mixed layer ($z = 0$) and (bottom) the binary tracer concentration at $y = 0$. The panels show the four different stages of the evolution: (a),(e) the linear phase at $t = 29N_m/f\Lambda$, (b),(f) the transition point between the linear phase and the nonlinear phase at $t = 48N_m/f\Lambda$, (c),(g) the early nonlinear phase at $t = 55N_m/f\Lambda$, and (d),(h) the late nonlinear phase when the inverse cascade of energy has taken place at $t = 85N_m/f\Lambda$. Data are taken from a nonlinear simulation with $Ri_m = 10$ and $Ri_t = 10^3$.</figcaption></figure>

negative buoyancy anomalies in the x direction (Fig. 6a). In this linear stage, vertical velocity is very small and the associated vertical transport of thermocline water is invisible (Fig. 6e). The buoyancy anomalies grow exponentially until the advection of perturbations by perturbation flow becomes important, at which time a transition from the linear phase to the nonlinear phase occurs (Figs. 6b,f). Note the correspondence between the patterns in Figs. 6b and 6f. The subsequent nonlinear phase is characterized by small cyclonic and anticyclonic eddies that transport momentum, buoyancy, and the binary tracer (Figs. 6c,g). Gradually, energy is transferred to larger scales, with vortices merging and growing in size. By the end of each simulation, there is a cyclonic vortex centered around a negative buoyancy anomaly at the midplane of the mixed layer (Figs. 6d,h), and the horizontally averaged buoyancy anomaly at this level has decreased below zero (the mean value in the initial state).

After the nonlinear phase starts around $t = 50N_m/f\Lambda$, vertical displacement of the base of the mixed layer becomes finite amplitude, leading to entrainment of thermocline water into the mixed layer. This is evident from the binary tracer and hydrostatic PV fields (Fig. 1). Hydrostatic PV is defined as

$$q = \left( f + \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} \right) \frac{\partial b}{\partial z} - \frac{\partial v}{\partial z} \frac{\partial b}{\partial x} + \left( \Lambda + \frac{\partial u}{\partial z} \right) \left( -f\Lambda + \frac{\partial b}{\partial y} \right). \quad (24)$$

In the initial state, $q = f\partial b/\partial z - f\Lambda^2 = fN^2 - f\Lambda^2 = f\Lambda^2(Ri - 1)$. Therefore, if we pick $f\Lambda^2$ as the dimension for PV, the nondimensional value in the initial state is $q' = Ri - 1$. Thermocline water has a high PV, and mixed layer water has a low PV. Since PV is approximately conserved on water parcels (it is exactly conserved if the flow is inviscid), high PV serves as a tracer of thermocline water, as does high binary tracer concentration. Figure 1a shows an x-z cross section of the binary tracer field, where water from
<figure><img src="image_8.png" /><figcaption>FIG. 7. The amount of thermocline water entrained into the mixed layer is inversely proportional to the initial PV difference between the thermocline and mixed layer. (a) Time series of the average binary tracer concentration in the interior of the mixed layer (at $z = 0.1h_m$). (b) Times series of the product of the average binary tracer concentration in the interior of the mixed layer (at $z = 0.1h_m$) and the initial PV difference $\Delta q' = R_i - R_m$ between the thermocline and mixed layer. We present the level of $z = 0.1h_m$ instead of the midplane of the mixed layer in the initial state ($z = 0$) because the base of the mixed layer is elevated by diffusion (see section 4).</figcaption></figure>

the shallow thermocline with high tracer concentration forms filaments that extend upward into the mixed layer. The cross section of PV shows a similar pattern to that of the binary tracer in the lower part of the mixed layer, where high-PV thermocline water forms filaments and is entrained upward (Fig. 1b), similar to the PV intrusions into the low-PV intermediate layer caused by vortex rollup in Garner et al. (1992). As the low-buoyancy and high-PV thermocline water is gradually entrained into the mixed layer, at the midplane of the mixed layer ($z = 0$), the average buoyancy decreases, and vorticity has an increasingly positive bias (Fig. 6d).

Unlike the lower part of the mixed layer, where PV and tracer filaments largely mirror one another, the upper part of the mixed layer shows differences between the two. The binary tracer concentration near the sea surface keeps a close-to-zero value because there is no tracer input at the surface (Fig. 1a). By contrast, high-PV blobs and filaments frequently appear near the surface and are subsequently transported downward into the interior of the mixed layer. This occurs because PV can be created or destroyed in the presence of viscosity. The upward PV flux can be decomposed into three terms: an advection term

$$J_a = wq, \qquad (25)$$

a diffusion term

$$J_d = -\left(f + \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}\right) \mathcal{F}(b), \qquad (26)$$

and a viscosity term

$$J_v = \frac{\partial b}{\partial x} \mathcal{F}(v) - \left( -f\Lambda + \frac{\partial b}{\partial y} \right) \mathcal{F}(u). \qquad (27)$$

At the boundaries, the advection term is zero because $w = 0$, but the diffusion and the viscosity terms can be nonzero, leading to a nonzero PV flux into or out of the domain. As expected, the instability produces PV input at the surface (e.g., Nakamura and Held 1989). Garner et al. (1992) pointed out that this PV input and entrainment into the mixed layer from

the surface is analogous to that at the base of the thermocline. The rigid surface acts like an overlying layer with infinite stratification.

Figure 7a shows the horizontally averaged binary tracer concentration in the interior of the mixed layer during the late linear phase and the nonlinear phase of the simulations. At fixed $R_{im}$, increasing $R_i$ suppresses the vertical exchange between the thermocline and mixed layer. This is consistent with the linear stability analysis (section 2), where we found a strong dependence of the variance of the interface displacement on $R_i$ (Fig. 5d).

Garner et al.'s (1992) analogy between a rigid boundary and an interface with an infinitely stratified layer inspires a scaling law for the dependence of the tracer injection into the mixed layer and thermocline stratification. We assume that for realistic conditions, $R_{im} \ll R_i$, such that the mixed layer can be thought of as sandwiched between two boundaries: the rigid boundary at the surface and the quasi-rigid interface with the thermocline. As shown in section 2, the instability in the mixed layer is determined largely by $R_{im}$. If this independence of thermocline properties carries over to the nonlinear phase, the dimensionless hydrostatic PV $q'$ in the mixed layer should only be a function of $R_{im}$, nondimensional depth $z'$, and nondimensional time $t'$:

$$q' = q'(R_{im}, z', t'). \qquad (28)$$

We denote the column integrated PV in the mixed layer as $Q$, and its nondimensional form as $Q' = Q/f\Lambda^2 h_m$. Then, as the vertical integration excludes the dependence on $z'$, $Q'$ is only a function of $R_{im}$ and $t'$:

$$Q' = Q'(R_{im}, t'). \qquad (29)$$

Because PV can be regarded as a conserved quantity carried by fluid parcels in the limit of large Re, and the two boundaries of the mixed layer are approximately symmetric, the above argument implies that the PV increase in the mixed layer is half sourced from the surface and half sourced from the interface with the thermocline, the latter via high-PV thermocline water
peeling off the thermocline and getting entrained into the mixed layer in filaments (Fig. 1b).

Let us denote the thickness of the thermocline layer that is pealed off and entrained into the mixed layer, the “entrainment depth,” as D and its nondimensional form as D' = D/hm. The contribution to the increase of mixed layer PV, half of the total PV increase, can then be written as the product of the PV difference between the two layers and the rate of change of the entrainment depth:

$$ \frac{1}{2} \frac{dQ'}{dt'} = (Ri_t - Ri_m) \frac{dD'}{dt'}, \quad (30) $$

where Ri_t - Ri_m are the nondimensional PV difference between the thermocline and mixed layer.

The left-hand side of the equation depends on Ri_m only, so the right-hand side also only depends on Ri_t. That is to say, dD'/dt' can be written as

$$ \frac{dD'}{dt'} = \frac{\varphi(Ri_m, t')}{Ri_t - Ri_m}, \quad (31) $$

where φ is a dimensionless function of Ri_m and t'. This means that the entrainment depth D' also takes the form

$$ D' = \frac{\Phi(Ri_m, t')}{Ri_t - Ri_m} \approx \frac{\Phi(Ri_m, t')}{Ri_t}, \quad (32) $$

where Φ is the integral of φ over t' and is also a dimensionless function, and we made use of Ri_t ≫ Ri_m.

Figure 7b shows the time series of the horizontally averaged binary tracer concentration, which represents the amount of vertical exchange between the thermocline and mixed layer, multiplied by the initial Ri_t - Ri_m. We can observe that for the same Ri_m, the curves with different Ri_t approximately collapse, thus supporting the proposed scaling law D' ~ Ri_t⁻¹. This is consistent with the results in the linear stability analysis, which predicts a decrease of the nondimensional surface displacement with increasing Ri_t (Fig. 5d).

Back in dimensional form,

$$ D = h_m D' \approx \frac{h_m}{Ri_t} \Phi(Ri_m, \frac{f\Lambda}{N_m} t). \quad (33) $$

Using the entrainment depth as a metric for the vertical exchange between the thermocline and mixed layer, we conclude that the vertical exchange is proportional to the thickness of the mixed layer and inversely proportional to the Richardson number of the thermocline.

Figure 7b also shows the dependence of the product c̄(Ri_t - Ri_m) on Ri_m. The product increases as Ri_m increases, which is at first sight inconsistent with the finding in the linear stability analysis that the variance of the vertical displacement of the mixed layer base is largely independent of Ri_m (Fig. 5d). This discrepancy is due to the fact that normalizing the horizontally averaged column-integrated perturbation energy to ρ₀Λ²h_m³ in the linear stability analysis does not precisely capture the energy at the moment when the instability transitions from the

<figure><img src="image_9.png" /><figcaption>FIG. 8. Time series of the horizontal-mean vertically integrated perturbation kinetic energy in the nonlinear simulations. The energy levels off after reaching levels that depend primarily on the Richardson number of the mixed layer Ri_m.</figcaption></figure>

linear to the nonlinear phase. Figure 8 shows the time series of the horizontally averaged column-integrated perturbation kinetic energy (EKE) in the nonlinear simulations, nondimensionalized against ρ₀Λ²h_m³. The kinetic energy first exponentially increases with time and then levels off at the transition point between the linear and the nonlinear phases. With the same Ri_m, the kinetic energy in the nonlinear phase increases with Ri_m. Since the fraction of kinetic energy in total perturbation energy is similar in different cases (ranging from 0.25 to 0.37), the total energy in the nonlinear phase also increases with Ri_m. This means that when normalizing the horizontally averaged column-integrated perturbation energy to ρ₀Λ²h_m³ in the linear stability analysis (the horizontal dashed line in Fig. 8), the amplitude of all the variables (including w) is overestimated (with respect to the “real” amplitude at the transition point as shown in the nonlinear simulations) in the small-Ri_m cases and underestimated in the large-Ri_m cases. Consequently, the positive correlation between the magnitude of vertical motion and Ri_m, as demonstrated in the nonlinear simulations, appears to be absent in Figs. 5c and 5d.

## 4. Discussion

In the nonlinear simulations, there is a marked asymmetry between the upwelling of thermocline water into the mixed layer and the subduction of mixed layer water into the thermocline (Fig. 1). It appears that MLI primarily peels filaments of thermocline water off the underlying layer and incorporates them into the mixed layer. The phenomenology and dynamics are similar to stratospheric intrusions into the troposphere (e.g., Shapiro 1980). Subduction of mixed layer water into the thermocline, as observed in Omand et al. (2015) or described in Freilich and Mahadevan (2021), is not prominent in our model and may be induced by processes that are not represented here, such as the interaction between mesoscale and submesoscale features or atmospherically forced unbalanced instabilities.
<figure><img src="image_10.png" /><figcaption>FIG. 9. In the nonlinear phase of the instability development, advective fluxes dominate in the vertical transport of the binary tracer and PV above the “effective” base of the mixed layer that is elevated by diffusion. Shown are the horizontally averaged upward fluxes of the binary tracer and PV within the mixed layer as functions of time and depth. Each row shows one simulation, with the Richardson numbers of the simulation indicated on the left. The left two columns show the two terms of the upward tracer flux, with the two plots in the same row sharing the same color bar. The right three columns show the three terms of the upward PV flux, with the three plots in the same row sharing the same color bar. The colors show nondimensional fluxes, with their dimensions shown at the top of the figure.</figcaption></figure>

It should be noted that our calculations capture only the tran-
sient evolution of unforced MLI. In the real ocean, mixed layer
turbulence is forced by winds and buoyancy loss at the surface,
which can both counter the restratification by MLI (e.g., Callies
and Ferrari 2018a) and modify the dynamics of the instability
(e.g., Young 1994; Crowe and Taylor 2019). This mixed layer
turbulence, of course, also entrains thermocline fluid when deep-
ening the mixed layer—although it cannot reach beneath the
maximum wintertime mixed layer depth. Given the short time
scale of MLI compared to the seasonal cycle, MLI can still in-
duce substantial exchange at times when the mixed layer is not
deepening, especially in late winter, when the mixed layer depth
is at its maximum and instabilities are strong.

In the nonlinear simulations, an artificial hyperviscosity
that damps variations on small spatial scales is indispensable

for stabilizing the model. The magnitude of the hyperviscos-
ity, measured by the reciprocal of the Reynolds number, must
be chosen such that the artificial diffusion of momentum,
buoyancy, and tracer is small compared to the advection while
keeping the model numerically stable. With a higher model
resolution, smaller-scale flow can be resolved, and the hyper-
viscosity required for numerical stability is smaller. We have
chosen the largest possible Reynolds number with which the
PV budget in the mixed layer, an indicator of numerical stability,
is closed at the current model resolution. We deem the em-
ployed Re = 10⁴ sufficient and expect the findings to be robust
under further increases of the resolution and Reynolds number.

That said, the artificial diffusive fluxes do elevate the base
of the mixed layer and thereby reduce its thickness (Fig. 9).
The first column in Fig. 9 shows the nondimensionalized
<figure><img src="image_11.png" /><figcaption>FIG. 10. In the nonlinear phase, the tracer and PV are homogenized in the mixed layer. Shown are the horizontally averaged binary tracer concentration (blue curves) and hydrostatic PV (orange curves) at $t = 70N_m/fA$. (a)–(f) Vertical profiles in both the initial mixed layer and thermocline. (g)–(l) Vertical profiles in the initial mixed layer only. The Richardson numbers are shown above each panel. Both the binary tracer concentration and PV are normalized such that the initial values in the mixed layer are equal to 0 and the initial values in the thermocline are equal to 1. Both the tracer and PV are roughly homogeneous between $z = -0.1h_m$ and $z = 0.3h_m$. That is why $z = 0.1h_m$ is chosen to represent the whole mixed layer in terms of the binary tracer concentration in Fig. 7.</figcaption></figure>

adveptive tracer flux $\overline{wc}$, the second column shows the nondimensionalized diffusive tracer flux $-v\partial\overline{c}/\partial z$, the third column shows the nondimensionalized advective PV flux $J_a = \overline{wq}$, the fourth column shows the nondimensionalized diffusive PV flux $J_d$, and the rightmost column shows the nondimensionalized

viscous PV flux $J_v$. In the linear phase, diffusive tracer and PV fluxes dominate, elevating the base of the mixed layer to the same depth. After the transition to the nonlinear phase, advective fluxes dominate in the mixed layer above the elevated base and also in the uppermost part of the thermocline below the
elevated mixed layer base. With the same $Ri_m$, a smaller $Ri_t$ permits stronger vertical motion and deeper penetration of vertical velocity into the (diffusion-elevated) thermocline, so the dominance of the advective fluxes over the diffusive and viscous fluxes extends to a greater depth (Fig. 9). This is consistent with the findings in the linear stability analysis (Fig. 4). Figure 10 shows that in the nonlinear phase, the tracer and PV are both nearly homogeneous in the vertical direction in the mixed layer, so we choose the level $z = 0.1h_m$, which is slightly above the initial center of the mixed layer at $z = 0$, to represent the mixed layer in terms of tracer concentration.

## 5. Conclusions

To investigate how mixed layer instabilities produce vertical exchange between the thermocline and mixed layer and how the exchange depends on environmental parameters, we perform a linear stability analysis and nonlinear simulations of an idealized configuration. We find that the instability induces vertical flow reaching into the upper thermocline, deforms the interface between the mixed layer and thermocline, and produces filaments of thermocline water spreading into and blending with the mixed layer water, similar to how stratospheric air is folded into the troposphere (e.g., Shapiro 1980).

The vertical velocity associated with mixed layer instability reaches deeper into the thermocline if the mixed layer is thicker or the thermocline stratification is weaker, and the vertical velocity at the base of the mixed layer is stronger if the thermocline stratification is weaker. A scaling law for the vertical exchange between the thermocline and mixed layer, characterized by how much thermocline water is entrained into the mixed layer, suggests that the exchange is proportional to the mixed layer depth and inversely proportional to the thermocline Richardson number. The vertical exchange also increases with an increasing mixed layer Richardson number. Our results imply that the tracer exchange between the thermocline and mixed layer is more efficient when the mixed layer is deeper, the lateral buoyancy gradient is larger, the thermocline stratification is weaker, and the mixed layer stratification is stronger. This suggests a vigorous exchange between the permanent thermocline and deep mixed layers in winter, especially in mode water formation regions, where mixed layers are deep and the underlying stratification is weak.

Existing parameterizations of MLI-induced restratification (Fox-Kemper et al. 2008) should be augmented to take the tracer exchange into account. Our findings suggest that an entrainment of thermocline water into the mixed layer should be included. The scaling law proposed here could be a starting point for parameterizing the dependence of such an entrainment velocity on environmental parameters.

It should be emphasized that we here investigated the exchange between the thermocline and mixed layer produced by an unforced, idealized, and isolated mixed layer instability. There are other processes that can exchange fluid between these two layers, such as turbulence produced by gravitational and symmetric instabilities or the deepening of the mixed layer caused by wind-forced turbulence. The relative importance of different processes to the exchange depends on environmental conditions and should be investigated in future work.

*Acknowledgments.* This material is based upon work supported by the National Science Foundation under Grant OCE-1924354.

*Data availability statement.* The numerical model simulations upon which this study is based are too large to archive or to transfer. Instead, we provide all the information needed to replicate the simulations; we used Dedalus as the numerical tool. Dedalus is available online at https://dedalus-project.org/. The equations, initial and boundary conditions are described in this article. Example Python scripts for using Dedalus to solve an eigenvalue problem and an initial value problem in this paper along with a Dedalus configuration file can be accessed at https://doi.org/10.22002/fbq4s-vt646 in the CaltechDATA repository.

## REFERENCES

Ascani, F., K. J. Richards, E. Firing, S. Grant, K. S. Johnson, Y. Jia, R. Lukas, and D. M. Karl, 2013: Physical and biological controls of nitrate concentrations in the upper subtropical North Pacific Ocean. *Deep-Sea Res. II*, **93**, 119–134, https://doi.org/10.1016/j.dsr2.2013.01.034.

Bachman, S. D., B. Fox-Kemper, J. R. Taylor, and L. N. Thomas, 2017: Parameterization of frontal symmetric instabilities. I: Theory for resolved fronts. *Ocean Modell.*, **109**, 72–95, https://doi.org/10.1016/j.ocemod.2016.12.003.

Balwada, D., K. S. Smith, and R. Abernathey, 2018: Submesoscale vertical velocities enhance tracer subduction in an idealized Antarctic circumpolar current. *Geophys. Res. Lett.*, **45**, 9790–9802, https://doi.org/10.1029/2018GL079244.

Boccaletti, G., R. Ferrari, and B. Fox-Kemper, 2007: Mixed layer instabilities and restratification. *J. Phys. Oceanogr.*, **37**, 2228–2250, https://doi.org/10.1175/JPO3101.1.

Boyd, J. P., 2001: *Chebyshev and Fourier Spectral Methods*. 2nd ed. Dover Publications, 665 pp.

Burns, K. J., G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown, 2020: Dedalus: A flexible framework for numerical simulations with spectral methods. *Phys. Rev. Res.*, **2**, 023068, https://doi.org/10.1103/PhysRevResearch.2.023068.

Callies, J., and R. Ferrari, 2018a: Baroclinic instability in the presence of convection. *J. Phys. Oceanogr.*, **48**, 45–60, https://doi.org/10.1175/JPO-D-17-0028.1.

—, and —, 2018b: Note on the rate of restratification in the baroclinic spindown of fronts. *J. Phys. Oceanogr.*, **48**, 1543–1553, https://doi.org/10.1175/JPO-D-17-0175.1.

—, —, J. M. Klymak, and J. Gula, 2015: Seasonality in submesoscale turbulence. *Nat. Commun.*, **6**, 6862, https://doi.org/10.1038/ncomms7862.

—, G. Flierl, R. Ferrari, and B. Fox-Kemper, 2016: The role of mixed-layer instabilities in submesoscale turbulence. *J. Fluid Mech.*, **788**, 5–41, https://doi.org/10.1017/jfm.2015.700.

—, R. Barkan, and A. Naveira Garabato, 2020: Time scales of submesoscale flow inferred from a mooring array. *J. Phys. Oceanogr.*, **50**, 1065–1086, https://doi.org/10.1175/JPO-D-19-0254.1.

Capet, X., J. C. McWilliams, M. J. Molemaker, and A. F. Shchepetkin, 2008: Mesoscale to submesoscale transition in the California current system. Part II: Frontal processes. *J. Phys. Oceanogr.*, **38**, 44–64, https://doi.org/10.1175/2007JPO3672.1.
Crowe, M. N., and J. R. Taylor, 2019: Baroclinic instability with a simple model for vertical mixing. *J. Phys. Oceanogr.*, **49**, 3273–3300, https://doi.org/10.1175/JPO-D-18-0270.1.

de Boyer Montégut, C., G. Madec, A. S. Fischer, A. Lazar, and D. Iudicone, 2004: Mixed layer depth over the global ocean: An examination of profile data and a profile-based climatology. *J. Geophys. Res.*, **109**, C12003, https://doi.org/10.1029/2004JC002378.

Dong, S., J. Sprintall, S. T. Gille, and L. Talley, 2008: Southern Ocean mixed-layer depth from Argo float profiles. *J. Geophys. Res.*, **113**, C06013, https://doi.org/10.1029/2006JC004051.

Eady, E. T., 1949: Long waves and cyclone waves. *Tellus*, **1** (3), 33–52, https://doi.org/10.3402/tellusa.v1i3.8507.

Fox-Kemper, B., and R. Ferrari, 2008: Parameterization of mixed layer eddies. Part II: Prognosis and impact. *J. Phys. Oceanogr.*, **38**, 1166–1179, https://doi.org/10.1175/2007JPO3788.1.

—, —, and R. Hallberg, 2008: Parameterization of mixed layer eddies. Part I: Theory and diagnosis. *J. Phys. Oceanogr.*, **38**, 1145–1165, https://doi.org/10.1175/2007JPO3792.1.

—, and Coauthors, 2011: Parameterization of mixed layer eddies. III: Implementation and impact in global ocean climate simulations. *Ocean Modell.*, **39**, 61–78, https://doi.org/10.1016/j.ocemod.2010.09.002.

Freilich, M., and A. Mahadevan, 2021: Coherent pathways for subduction from the surface mixed layer at ocean fronts. *J. Geophys. Res. Oceans*, **126**, e2020JC017042, https://doi.org/10.1029/2020JC017042.

Garner, S. T., N. Nakamura, and I. M. Held, 1992: Nonlinear equilibration of two-dimensional Eady waves: A new perspective. *J. Atmos. Sci.*, **49**, 1984–1996, https://doi.org/10.1175/1520-0469(1992)049<1984:NEOTDE>2.0.CO;2.

Haine, T. W. N., and J. Marshall, 1998: Gravitational, symmetric, and baroclinic instability of the ocean mixed layer. *J. Phys. Oceanogr.*, **28**, 634–658, https://doi.org/10.1175/1520-0485(1998)028<0634:GSABIO>2.0.CO;2.

Large, W. G., J. C. McWilliams, and S. C. Doney, 1994: Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization. *Rev. Geophys.*, **32**, 363–403, https://doi.org/10.1029/94RG01872.

Lévy, M., P. Klein, and A.-M. Tréguier, 2001: Impact of submesoscale physics on production and subduction of phytoplankton in an oligotrophic regime. *J. Mar. Res.*, **59**, 535–565, https://doi.org/10.1357/002224001762842181.

—, R. Ferrari, P. J. S. Franks, A. P. Martin, and P. Rivière, 2012a: Bringing physics to life at the submesoscale. *Geophys. Res. Lett.*, **39**, L14602, https://doi.org/10.1029/2012GL052756.

—, D. Iovino, L. Resplandy, P. Klein, G. Madec, A.-M. Tréguier, S. Masson, and K. Takahashi, 2012b: Large-scale impacts of submesoscale dynamics on phytoplankton: Local and remote effects. *Ocean Modell.*, **43–44**, 77–93, https://doi.org/10.1016/j.ocemod.2011.12.003.

—, P. J. S. Franks, and K. S. Smith, 2018: The role of submesoscale currents in structuring marine ecosystems. *Nat. Commun.*, **9**, 4758, https://doi.org/10.1038/s41467-018-07059-3.

Mahadevan, A., 2006: Modeling vertical motion at ocean fronts: Are nonhydrostatic effects relevant at submesoscales? *Ocean Modell.*, **14**, 222–240, https://doi.org/10.1016/j.ocemod.2006.05.005.

—, 2014: Eddy effects on biogeochemistry. *Nature*, **506**, 168–169, https://doi.org/10.1038/nature13048.

—, 2016: The impact of submesoscale physics on primary productivity of plankton. *Annu. Rev. Mar. Sci.*, **8**, 161–184, https://doi.org/10.1146/annurev-marine-010814-015912.

—, and D. Archer, 2000: Modeling the impact of fronts and mesoscale circulation on the nutrient supply and biogeochemistry of the upper ocean. *J. Geophys. Res.*, **105**, 1209–1225, https://doi.org/10.1029/1999JC900216.

—, and A. Tandon, 2006: An analysis of mechanisms for submesoscale vertical motion at ocean fronts. *Ocean Modell.*, **14**, 241–256, https://doi.org/10.1016/j.ocemod.2006.05.006.

—, —, and R. Ferrari, 2010: Rapid changes in mixed layer stratification driven by submesoscale instabilities and winds. *J. Geophys. Res.*, **115**, C03017, https://doi.org/10.1029/2008JC005203.

McWilliams, J. C., 2016: Submesoscale currents in the ocean. *Proc. Roy. Soc., 472A*, 20160117, https://doi.org/10.1098/rspa.2016.0117.

Mensa, J. A., Z. Garraffo, A. Griffa, T. M. Özgökmen, A. Haza, and M. Veneziani, 2013: Seasonality of the submesoscale dynamics in the Gulf Stream region. *Ocean Dyn.*, **63**, 923–941, https://doi.org/10.1007/s10236-013-0633-1.

Nakamura, N., and I. M. Held, 1989: Nonlinear equilibration of two-dimensional Eady waves. *J. Atmos. Sci.*, **46**, 3055–3064, https://doi.org/10.1175/1520-0469(1989)046<3055:NEOTDE>2.0.CO;2.

Omand, M. M., E. A. D'Asaro, C. M. Lee, M. J. Perry, N. Briggs, I. Cetini, and A. Mahadevan, 2015: Eddy-driven subduction exports particulate organic carbon from the spring bloom. *Science*, **348**, 222–225, https://doi.org/10.1126/science.1260062.

Qiu, B., T. Nakano, S. Chen, and P. Klein, 2017: Submesoscale transition from geostrophic flows to internal waves in the northwestern Pacific upper ocean. *Nat. Commun.*, **8**, 14055, https://doi.org/10.1038/ncomms14055.

Sasaki, H., P. Klein, B. Qiu, and Y. Sasai, 2014: Impact of oceanic-scale interactions on the seasonal modulation of ocean dynamics by the atmosphere. *Nat. Commun.*, **5**, 5636, https://doi.org/10.1038/ncomms6636.

Shapiro, M. A., 1980: Turbulent mixing within tropopause folds as a mechanism for the exchange of chemical constituents between the stratosphere and troposphere. *J. Atmos. Sci.*, **37**, 994–1004, https://doi.org/10.1175/1520-0469(1980)037<0994:TMWTFA>2.0.CO;2.

Soares, S. M., S. T. Gille, T. K. Chereskin, E. Firing, J. Hummon, and C. B. Rocha, 2022: Transition from balanced to unbalanced motion in the eastern tropical Pacific. *J. Phys. Oceanogr.*, **52**, 1775–1795, https://doi.org/10.1175/JPO-D-21-0139.1.

Stone, P. H., 1966: On non-geostrophic baroclinic stability. *J. Atmos. Sci.*, **23**, 390–400, https://doi.org/10.1175/1520-0469(1966)023<0390:ONGBS>2.0.CO;2.

—, 1970: On non-geostrophic baroclinic stability: Part II. *J. Atmos. Sci.*, **27**, 721–726, https://doi.org/10.1175/1520-0469(1970)027<0721:ONGBSP>2.0.CO;2.

—, 1971: Baroclinic stability under non-hydrostatic conditions. *J. Fluid Mech.*, **45**, 659–671, https://doi.org/10.1017/S0022112071000260.

—, S. Hess, R. Hadlock, and P. Ray, 1969: Preliminary results of experiments with symmetric baroclinic instabilities. *J. Atmos. Sci.*, **26**, 991–996, https://doi.org/10.1175/1520-0469(1969)026<0991:PROEWS>2.0.CO;2.

Taylor, J. R., and R. Ferrari, 2010: Buoyancy and wind-driven convection at mixed layer density fronts. *J. Phys. Oceanogr.*, **40**, 1222–1242, https://doi.org/10.1175/2010JPO4365.1.
—, and A. F. Thompson, 2023: Submesoscale dynamics in the upper ocean. *Annu. Rev. Fluid Mech.*, **55**, 103–127, https://doi.org/10.1146/annurev-fluid-031422-095147.

Thomas, L. N., A. Tandon, and A. Mahadevan, 2008: Submeso-scale processes and dynamics. *Ocean Modeling in an Eddying Regime, Geophys. Monogr.*, Vol. 177, Amer. Geophys. Union, 17–38, https://doi.org/10.1029/177GM04.

—, J. R. Taylor, R. Ferrari, and T. M. Joyce, 2013: Symmetric Instability in the Gulf Stream. *Deep-Sea Res. II*, **91**, 96–110, https://doi.org/10.1016/j.dsr2.2013.02.025.

Vallis, G. K., 2017: *Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation*. 2nd ed. Cambridge University Press, 964 pp.

Van Roekel, L., and Coauthors, 2018: The KPP boundary layer scheme for the ocean: Revisiting its formulation and benchmarking one-dimensional simulations relative to LES. *J. Adv. Model. Earth Syst.*, **10**, 2647–2685, https://doi.org/10.1029/2018MS001336.

Young, W. R., 1994: The subinertial mixed layer approximation. *J. Phys. Oceanogr.*, **24**, 1812–1826, https://doi.org/10.1175/1520-0485(1994)024<1812:TSMLA>2.0.CO;2.