1. Introduction
At length scales larger than
$\sim \!100 \text{ km}$
, flows in the ocean are close to geostrophic balance – a state of approximate balance between the Coriolis acceleration and pressure gradient – and evolve slowly, over many days. At smaller length scales, beyond those that can currently be resolved in global ocean models, lies the submesoscale (
$\sim\! 0.1{-}10\ \text{ km}$
). Understanding of the impact of submesoscale flows on the ocean is emerging, with observational and numerical evidence that processes at this scale, which is poorly represented in parametrisations, are crucial components of the ocean’s physical and biogeochemical processes (Mahadevan Reference Mahadevan2016; Calil et al. Reference Calil, Suzuki, Baschek and da Silveira2021; Gula et al. Reference Gula, Taylor, Shcherbina and Mahadevan2022; Aravind et al. Reference Aravind, Verma, Sarkar, Freilich, Mahadevan, Haley, Lermusiaux and Allshouse2023; Pham et al. Reference Pham, Verma, Sarkar, Shcherbina and D’Asaro2024; Zhu et al. Reference Zhu, Yang, Li, Chen, Ma, Cai and Wu2024).
Images of ocean properties at submesoscale resolution are often peppered with highly anisotropic features. Fronts – confined lateral gradients in density – and filaments – opposing pairs of fronts – remain balanced by the Coriolis acceleration of an along-front jet. In the open ocean, these features are created primarily at the edges of mesoscale vortices such as those created by separating western boundary currents, such as the Gulf Stream or the Kuroshio current (McWilliams Reference McWilliams2016). As fronts and filaments are drawn out by the strain flow in these regions, horizontal gradients in velocities and tracers become larger, eventually leading to the onset of a new dynamical regime, where geostrophic and hydrostatic balances no longer describe the flow well (Shakespeare Reference Shakespeare2016; Barkan et al. Reference Barkan, Molemaker, Srinivasan, McWilliams and D’Asaro2019; Taylor & Thompson Reference Taylor and Thompson2022).
In addition to the along-front geostrophic jet, an ageostrophic secondary circulation may develop at a front due to a background strain, turbulent mixing, instability or otherwise (Crowe & Taylor Reference Crowe and Taylor2018; Sullivan & McWilliams Reference Sullivan and McWilliams2018; McWilliams Reference McWilliams2021). The resulting vertical velocities provide a route for transport of tracers, momentum and energy between the surface layer and deep water (Mahadevan Reference Mahadevan2016; Klein et al. Reference Klein, Lapeyre, Siegelman, Qiu, Fu, Torres, Su, Menemenlis and Le Gentil2019; Verma, Pham & Sarkar Reference Verma, Pham and Sarkar2019; Freilich & Mahadevan Reference Freilich and Mahadevan2021). In particular, the submesoscale is a link between the larger-scale, horizontally non-divergent flows and smaller-scale turbulence that facilitates the observed forward cascade of energy in the ocean (Srinivasan, Barkan & McWilliams Reference Srinivasan, Barkan and McWilliams2023).
In the absence of turbulence, finite time singularities form in the horizontal gradients when forced by a strain flow (Hoskins & Bretherton Reference Hoskins and Bretherton1972). For a realistic front, some effective horizontal diffusion due to turbulence will arrest the frontogenesis before this singularity is reached, though this process is not well understood (Samelson & Skyllingstad Reference Samelson and Skyllingstad2016; Verma et al. Reference Verma, Pham and Sarkar2019). The turbulence, namely the vertical momentum mixing induced by surface cooling or wind stress, acts on the vertically sheared thermal wind velocity. If only these processes act in the submesoscale mixed layer, then the expected momentum balance is the turbulent thermal wind balance between the Coriolis acceleration of an along-front jet, the horizontal pressure gradient and vertical momentum mixing (Gula, Molemaker & McWilliams Reference Gula, Molemaker and McWilliams2014; Crowe & Taylor Reference Crowe and Taylor2018; Bodner et al. Reference Bodner, Fox-Kemper, Johnson, Van Roekel, McWilliams, Sullivan, Hall and Dong2023).
Of particular interest in this study is symmetric instability (SI). A symmetrically unstable flow is unstable to perturbations that exchange fluid parcels along isopycnals, with no variation along the direction of mean flow (hence symmetric; Stone Reference Stone1966). The instability manifests as highly anisotropic rolls, along the direction of the frontal jet, with major axis aligned with isopycnals. The most unstable across-front wavenumber tends to infinity in the inviscid limit, so the wavenumber of the observed mode is presumed to be restricted by diffusion (Harris, Poulin & Lamb Reference Harris, Poulin and Lamb2022). The SI extracts kinetic energy from the thermal wind velocity shear (Haine & Marshall Reference Haine and Marshall1998; Thomas et al. Reference Thomas, Taylor, Ferrari and Joyce2013; Taylor & Thompson Reference Taylor and Thompson2022).
For a Boussinesq fluid with equal momentum and buoyancy diffusivities, the stability of a flow to symmetric perturbations can be diagnosed by the sign of Ertel potential vorticity
$q$
. Fluid unstable to SI has
$fq = f(f\hat {\boldsymbol{z}} + {\boldsymbol\nabla} \times \boldsymbol{u}) \,{\boldsymbol\cdot}\, {\boldsymbol\nabla} b \lt 0,$
where
$\boldsymbol u = u \hat {\boldsymbol{x}} + v \hat {\boldsymbol{y}} + w \hat {\boldsymbol{z}}$
and
$b$
are velocity and buoyancy fields, respectively, and
$f$
is the Coriolis frequency (Emanuel Reference Emanuel1979). The SI competes with gravitational, inertial, baroclinic and Kelvin–Helmholtz instabilities, and dominates in regions where the Richardson number (
$Ri = N^2 / S^2$
, where
$N$
and
$S$
are buoyancy frequency and vertical shear of horizontal velocity, respectively) of the flow lies in the range
$0.25 \lt Ri \lt 0.95$
(Stone Reference Stone1966; Thomas et al. Reference Thomas, Taylor, Ferrari and Joyce2013). The small size and strong baroclinicity of submesoscale fronts often produce favourable conditions for SI (Thomas et al. Reference Thomas, Taylor, Ferrari and Joyce2013; Gula et al. Reference Gula, Molemaker and McWilliams2014).
The behaviour of fluid susceptible to SI (or close to unstable) is suspected to have important consequences for ocean flows at the submesoscale and, for example, has been studied in the contexts of frontogenesis and frontal arrest (Thomas Reference Thomas2012; Verma et al. Reference Verma, Pham and Sarkar2019), mixing and vertical transport (Brannigan Reference Brannigan2016; Bosse et al. Reference Bosse, Testor, Damien, Estournel, Marsaleix, Mortier, Prieur and Taillandier2021), and near-inertial wave generation (Wienkers et al. Reference Wienkers, Thomas and Taylor2021a ). However, the initial and long-time behaviours of a symmetrically unstable flow differ in form, growth rate and energy source. The form of the linear perturbations that optimally extract energy from an Eady background state differs from the most unstable normal mode during initial evolution, and their growth rate is enhanced during an initial, transient phase of instability. This behaviour persists in fully nonlinear, three-dimensional simulations (Heifetz & Farrell Reference Heifetz and Farrell2007; Zemskova, Passaggia & White Reference Zemskova, Passaggia and White2020; Kimura Reference Kimura2024). Initial motions in an unstable flow seeded with noise do not resemble the symmetric normal modes, and exhibit a faster growth rate across parameter space, a behaviour that is general across hydrodynamic stability (Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993; Zemskova et al. Reference Zemskova, Passaggia and White2020). This non-normal, transient growth will play a particularly important role in this study.
In this paper, we present a set of strongly baroclinic (
$Ri \leqslant 0.2$
) simulations resembling the configuration of Sullivan & McWilliams (Reference Sullivan and McWilliams2018), but with coarser spatial resolution, and without the surface cooling that drives vertical mixing. On the other hand, we run more simulations, allowing for a better parameter space exploration, and we do so for longer periods of time, namely, several inertial periods each time.
As in the aforementioned study, the reference state is unstable; and we pay particular attention to the consequences of instability when initial noise is applied, and how this may obfuscate conclusions that assume a balanced turbulent thermal wind state. Specifically, we share the goal of Sullivan & McWilliams (Reference Sullivan and McWilliams2018) to have realistic levels of noise, but while they attempt to generate realistic turbulent noise by continuously forcing the ocean surface, we simply generate random noise that aims to mirror the turbulent kinetic energy profile of mixed layer convection driven by surface cooling. This also results in a larger amplitude than in many idealised studies, which seed the fluid with very small levels of noise (Zemskova et al. Reference Zemskova, Passaggia and White2020; Wienkers et al. Reference Wienkers, Thomas and Taylor2021b ).
We show that instabilities act quickly to remove kinetic energy and excite near-inertial oscillations within the fronts that briefly strengthen horizontal buoyancy gradients at the surface. We argue that we observe a transient response of the filament as it evolves towards normal-mode SI. Finally, we suggest that the initialisation of both turbulence and an unstable reference state must be done with care, and initial conditions of numerical studies of similar strongly baroclinic features must be designed with initial instability and possible subsequent oscillations in mind, or risk drawing conclusions that may not be generalisable to balanced flows in the real ocean.
2. Numerical simulations
The simulations presented in this paper are of a set of buoyancy filaments in a rotating, Boussinesq fluid, namely,
\begin{align} \frac {\text{D}\boldsymbol u}{\text{D}t} + f \hat{\boldsymbol z} \times \boldsymbol u = -{\boldsymbol\nabla} \varphi + b\hat{\boldsymbol z} + {\boldsymbol\nabla} \,{\boldsymbol\cdot}\, (\nu\, {\boldsymbol\nabla} \boldsymbol{u})+ \boldsymbol {F_u}, \end{align}
\begin{align} \frac {\text{D}b}{\text{D}t} = {\boldsymbol\nabla} \,{\boldsymbol\cdot}\, (\nu\, {\boldsymbol\nabla} b) + F_b \quad \text{and} \quad {\boldsymbol\nabla} \boldsymbol{\,{\boldsymbol\cdot}\, } \boldsymbol u = 0, \end{align}
with
$f = 10^{-4}\ \text{s}^{-1}$
, geopotential
$\varphi$
, and Lagrangian derivative
${\text{D}\phi }/{\text{D}t} = \partial \phi / \partial t + \boldsymbol{u} \,{\boldsymbol\cdot}\, {\boldsymbol\nabla} \phi$
. The terms
$\boldsymbol{F_u}$
and
$F_b$
are velocity and buoyancy forcing terms; we elaborate on these quantities below. The equation set (2.1) is integrated using Oceananigans, a Julia package for simulation of incompressible fluid flows in an oceanic context (Ramadhan et al. Reference Ramadhan, Wagner, Hill, Campin, Churavy, Besard, Souza, Edelman, Ferrari and Marshall2020). The grid size is
$(N_x, N_y, N_z) = (1024, 1024, 128)$
, and the simulation domain is
$(x, y)\in [-5L, 5L)^2$
, where
$L=1000\text{ m}$
is a measure of the initial surface separation between the centres of the two counter-flowing jets, and
$z\in [-2.5H, 0]$
, where
$H=100\text{ m}$
is the initial depth of the mixed layer. The grid spacing is uniform in the horizontal directions. The vertical grid spacing smoothly transitions between two uniform values at the pycnocline, such that approximately
$3/4$
of the cells are in the surface mixed layer (
$z\gt -H$
), while the remainder are in the deep water.
The domain is periodic in both horizontal directions. In the vertical, boundary conditions are no-flux – the filament is not forced by wind stress or surface cooling, namely,
$u_{,z} = v_{,z} = w = b_{,z} = 0$
at
$z=0$
and
$u_{,z} = v_{,z} = w = 0$
and
$b_{,z} = N_0^2$
at
$z=-2.5H$
, where
$\phi _{,i}$
denotes differentiation of
$\phi$
with respect to coordinate
$i$
, and
$N_0$
is the deep-water buoyancy frequency.
To prevent internal wave reflections from the bottom boundary, a sponge layer forcing
$\boldsymbol F_{\boldsymbol u}, F_b$
is included for all fields. The fields are relaxed towards the reference state (
$\boldsymbol u_0, b_0$
) at a rate that decays quadratically away from the bottom boundary, namely,
\begin{align} \boldsymbol F_{\boldsymbol u} = - \sigma (z)\, (\boldsymbol u - \boldsymbol u_0) \quad \text{and} \quad F_b = -\sigma (z)\, (b - b_0), \end{align}
\begin{align} \text{with}\quad \sigma (z) = \frac 12\frac {N_0}{2\unicode{x03C0} }\times \left \{\begin{array}{ll} (2z/H + 4)^2, & z \leqslant -2H, \\ 0, & z \gt -2H. \end{array}\right. \end{align}
The isotropic eddy viscosity
$\nu$
is parametrised as a Smagorinksy–Lilly turbulent closure with Smagorinsky constant
$C=0.16$
(Smagorinsky Reference Smagorinsky1963; Lilly Reference Lilly1967), and we set
$\kappa = \nu$
. The adaptive time-stepping routine aims to achieve maximum Courant–Friedrichs–Levy number 0.5.
The reference state of the simulation is a buoyancy profile
$b_0(x, z)$
, and a jet
$\boldsymbol u_0 = v_0(x, z)\, \hat{\boldsymbol y}$
in thermal wind balance, namely,
\begin{align} v_0 = \frac {1}{f}\int _{-\infty }^z \frac {\partial b_0}{\partial x} \, \text{d}z. \end{align}
The reference buoyancy contains a weakly stratified surface layer above a strongly stratified fluid. The height of the transition between these two stratifications varies, providing the horizontal buoyancy gradient that makes up the filament. The reference filament state is shown in figure 1, and the exact definitions of the functional forms of
$v_0$
and
$b_0$
are given in Appendix A.
Figure 1. (a) The local Rossby number for the reference state with
$Ri_{{min}} = 0.1$
. Black contours are isopycnals. Fluid inside the red, dashed contour has
$fq \lt 0$
. (b) The local Richardson number, as in (a).
The reference state is invariant in the down-front direction
$y$
, and the initial evolution introduces little variation along this direction for the large-scale velocity and buoyancy fields (i.e. the filament does not meander as the short down-front domain constrains the growth rate of baroclinic instabilities). As such, we define the mean state to simply be the down-front average of the full three-dimensional state, with averaging operator
\begin{equation} \overline {\phi }(x, z, t):= \frac {1}{10L}\int _{-5L}^{5L} \phi (x, y, z, t) \,\text{d}y, \end{equation}
and down-front fluctuations
$\phi ^{\prime} = \phi - \overline {\phi }$
. The primary non-dimensional numbers that determine the dynamical regime of an inviscid submesoscale mixed layer flow are the Rossby and Richardson numbers (McWilliams Reference McWilliams2016; Bodner et al. Reference Bodner, Fox-Kemper, Johnson, Van Roekel, McWilliams, Sullivan, Hall and Dong2023). For the base state, the local values of each are, respectively,
\begin{gather} Ro = \frac {1}{f}\frac {\partial v_0}{\partial x}\quad \text{and}\quad Ri = \frac {\partial b_0 / \partial z}{(\partial v_0 / \partial z)^2}. \end{gather}
We parametrise our initial conditions by extremal values of the local Rossby and Richardson numbers as follows.
This study is concerned with the submesoscale ocean, where advection and the Coriolis acceleration are comparable in magnitude:
$Ro \sim\! 1$
. To limit the number of possible instabilities, the reference state has minimum vertical vorticity above
$-f$
, which excludes inertial instability. We use
$Ro_{{min}} = -0.8$
here. Given the prescribed filament shape, this results in a maximum
$Ro_{{max}} = 1.9$
. For typical horizontal velocities in submesoscale fronts, this corresponds to jet widths of order
$1{-}10\ \text{km}$
at mid-latitudes (Siegelman et al. Reference Siegelman, Klein, Thompson, Torres and Menemenlis2020; Calil et al. Reference Calil, Suzuki, Baschek and da Silveira2021).
The submesoscale mixed layer often has
$Ri \lesssim 1$
(Zhu et al. Reference Zhu, Yang, Li, Chen, Ma, Cai and Wu2024). For the filament shape used here, the Richardson number is approximately constant with depth in the mixed layer, and attains a minimum value in the horizontal centre of each front. We present the results of a set of low-
$Ri$
filaments, with a range of minimum Richardson numbers:
$Ri_{{min}} \in \{0.0, 0.1, 0.2\}$
. The remainder of the non-dimensional parameters, as well as the relationship between
$Ro_{{min}}$
,
$Ri_{{min}}$
and dimensional parameters, are given in Appendix A. Figures 1(a) and 1(b) show the local
$Ro$
and
$Ri$
of the
$Ri_{min} = 0.1$
reference state, respectively, and the region with
$fq\lt 0$
at the outer edges of the filament.
Before we impose the reference state, there is a brief initialisation process to ensure that the initial conditions, with noise, satisfy the Boussinesq equations to at least first order, as well as to allow the sub-grid parametrisation to reduce very large gradients caused by the pre-initial random noise before the unstable flow is introduced. We pre-initialise the simulation at
$ft/2\unicode{x03C0} =-1$
with a state identical to the reference state at
$|x|\rightarrow \infty$
. We then add noise to all velocities, and run the simulation until
$t=0$
. This produces fields
$\boldsymbol u_{{init}} = \overline {\boldsymbol{u}}_{{init}} + \boldsymbol{u}_{{init}}^{\prime}$
and
$b_{{init}} = \overline {b}_{{init}} + b_{{init}}^{\prime}$
. The pre-initialised noise is sampled from a Gaussian distribution, and is chosen to have magnitude and variation with depth of the turbulent kinetic energy (TKE) similar to that which would be produced by forcing due to surface cooling of order
$\sim\! 100\text{ W m}^{-2}$
(see Sullivan & McWilliams Reference Sullivan and McWilliams2018). A slice of vertical velocity and the horizontally averaged kinetic energy are shown in figure 2. We explore the impact of the choice of initial noise on the conclusions of this study in Appendix B. The initial state of the filament phase of the simulation is taken to be the reference state plus the fluctuation terms from the initialisation phase, namely,
\begin{align} \boldsymbol{u}(t=0) = v_0\hat {\boldsymbol{y}} + \boldsymbol{u}_{{init}}^{\prime}\quad \text{and}\quad b(t=0) = b_0 + b_{{init}}^{\prime}. \end{align}
Figure 2. (a) A slice of vertical velocity
$w_{{init}}$
in the initial condition, just before the filament is imposed at
$t=0$
. (b) Kinetic energy
$\boldsymbol{u}_{{init}}\boldsymbol{\,{\boldsymbol\cdot}\, }\boldsymbol{u}_{{init}} / 2$
at the same time, horizontally averaged, as a function of depth.
The mean, balanced state provides the energy source for instabilities, with the nature of the instability determining the dominant energy source. We consider contributions to the rate of change of the total (kinetic and potential) energy in the mean state, integrated over the whole domain
$\mathcal{V}$
. For the boundary conditions considered here, the rate of change of mean total energy may be partitioned into turbulent flux terms as
\begin{equation} \frac {\text{d}}{\text{d}t}\int _{\mathcal{V}} \text{d}V \left ( \frac {1}{2}\overline {\boldsymbol{u}} \boldsymbol{\,{\boldsymbol\cdot}\, } \overline {\boldsymbol{u}} - \overline {b}z \right) \approx \int _{\mathcal{V}}\text{d}V \left ( {\overline {\boldsymbol{u}}_{,x} \boldsymbol{\,{\boldsymbol\cdot}\, } \overline {u^{\prime}\boldsymbol{u}^{\prime}}} + {\overline {\boldsymbol{u}}_{,z} \boldsymbol{\,{\boldsymbol\cdot}\, } \overline {w^{\prime}\boldsymbol{u}^{\prime}}} - { \overline {w^{\prime}b^{\prime}}} \right), \end{equation}
wherein the direct impacts of the sponge forcing and sub-grid dissipation on the reference state are assumed to be negligible. Energisation of fluctuations by a laterally or vertically sheared velocity profile is given by the lateral or vertical shear production, respectively
\begin{align} \text{LSP} = -\int _{\mathcal{V}} \text{d}V\; \overline {\boldsymbol u}_{,x} \boldsymbol{\,{\boldsymbol\cdot}\, } \overline {u^{\prime} \boldsymbol u^{\prime}}, \quad \text{VSP} = -\int _{\mathcal{V}} \text{d}V\; \overline {\boldsymbol u}_{,z} \boldsymbol{\,{\boldsymbol\cdot}\, } \overline {w^{\prime} \boldsymbol u^{\prime}}. \end{align}
The VSP approximates the geostrophic shear production (GSP) if the down-front averaging can be assumed to separate the geostrophic and ageostrophic components. The GSP provides the energy source for SI, among others. Lateral shear is known to energise inertial and barotropic instabilities. Mean potential energy is transferred to fluctuations via the (vertical) buoyancy flux, typically associated with baroclinic and gravitational instabilities:
\begin{align} \text{BFLUX} = \int _{\mathcal{V}} \text{d}V\; \overline {w^{\prime}b^{\prime}}. \end{align}
3. Results
3.1. General evolution
For the first 2–3 inertial periods (
$\sim\! 54\ \text{h}$
at mid-latitudes), the filament remains approximately uniform in the down-front direction. Towards the end of our simulations (
$ft/ 2\unicode{x03C0} = 5.4$
, not shown), and especially visible at the surface, the isopycnals begin to meander due to baroclinic instability. The domain is not large enough to support many baroclinic modes; following Stone (Reference Stone1966), the fastest-growing baroclinic mode wavelength
$\lambda$
can be estimated for an Eady model with parameters comparable to the centre of a jet,
\begin{equation} \lambda =\frac {2\unicode{x03C0} v_0}{f}\sqrt {\frac {1+Ri_{min}}{5/2}} \approx 5L, \quad \end{equation}
which is half the down-front length of our domain. The small domain limits the baroclinic growth rate. Nevertheless, we ran a simulation with a longer down-front domain that resolved approximately eight wavelengths of this instability (not shown), in which the meandering indeed becomes noticeable sooner, at
$ft/2\unicode{x03C0} \approx 2$
; and we verified that even then, the baroclinic instability does not alter the conclusions that we present in this study. The rest of this discussion is concerned with the evolution of the filament before the meandering, which can be understood from the perspective of the down-front mean state.
Figure 3. Evolution of the mean state for the simulation with
$Ri_{{min}} = 0$
. Three snapshots are shown, at
$ft/2\unicode{x03C0} =\{0.8, 1.5, 2.2\}$
. The heat map shows
$\overline {v}$
. Solid black contours are mean isopyncals
$\overline {b}$
. Dashed orange (purple) contours are (anti)clockwise streamlines of the down-front mean flow.
The vertically sheared jet does not satisfy the no-flux boundary condition at
$z=0$
for the down-front velocity, and there is an onset of restratifying across-front Ekman flow; however, by
$ft/2\unicode{x03C0} = 0.7$
, this is overshadowed by a strong secondary circulation throughout the mixed layer, which is the focus of this section.
Figure 3 outlines the evolution of the down-front mean state during the first inertial periods. Shortly after the filament is allowed to evolve freely, a strong secondary circulation on the scale of the fronts emerges (figure 3
a). The secondary circulation acts to restratify the fronts. At the surface, the circulation increases front strength: horizontal convergence pushes isopycnals together immediately under the surface. Afterwards, the secondary circulation changes sign (figure 3
b). Hovmöller plots of the down-front-averaged horizontal buoyancy gradient and sub-grid dissipation, averaged over the down-front direction and the top 10 % of the mixed layer, are shown in figure 4. Near-surface isopycnals are drawn together, and horizontal buoyancy gradients reach maximum strength at
$ft/2\unicode{x03C0} \approx 1.3$
. This is associated with the maximum level of sub-grid dissipation. Afterwards, the down-front-averaged near-surface buoyancy gradients weaken before somewhat strengthening again, reaching a second maximum approximately one inertial period after the initial maximum. A sub-grid dissipation signal also recurs at this time (
$ft/2\unicode{x03C0} \approx 2.1$
). The fronts approach a final state (disregarding the later meandering due to baroclinic instability) with surface buoyancy gradients that are stronger than those in the reference state, though weaker than the maxima achieved earlier.
Figure 4. (a) A Hovmöller plot of the horizontal buoyancy gradient of the down-front mean state for the simulation with
$Ri_{{min}} = 0$
, averaged over the top 10 % of the mixed layer. Black lines are contours of the buoyancy averaged over the same region. Red dashed lines indicate the times at which the snapshots in figure 3 are taken. (b) The down-front mean sub-grid dissipation of kinetic energy
$\overline {\varepsilon } = \overline {\nu\, {\boldsymbol\nabla} \boldsymbol{u}:{\boldsymbol\nabla} \boldsymbol{u}}$
, as in (a).
3.2. Initial evolution is a fast instability
In this subsection, we show that the dynamics in the first inertial period contains an instability that extracts energy from the the vertical sheared velocity.
Figure 5. Terms extracting energy from the mean state for each of the simulations integrated over the mixed layer. Each term is normalised by the maximum VSP in the simulation, with
$Ri_{min}=0$
. Marked points in (c) indicate the times at which the snapshots of vertical velocity and contribution to VSP in figures 6 and 7 are taken.
Figure 6. At
$t/2\unicode{x03C0} =0.30$
, comparison of an
$x$
–
$z$
slice at
$y=0$
of (a) vertical velocity
$w$
and (b) contribution to VSP in the left-hand side of the filament for the simulation with
$Ri_{{min}}=0$
. The red dashed line in (a) locates
$z=-50\ \text{m}$
, the depth of the
$x$
–
$y$
slice of vertical velocity displayed in (c). The horizontal slice of
$w$
is coarse-grained, with a Gaussian filter with a half-width of one grid cell for clarity. A video time series of this plot appears in supplementary movie 1.
Figure 5 shows an initial sharp peak in the VSP, greater for the simulations with smaller minimum
$Ri$
, before
$t/2\unicode{x03C0} \lesssim 0.8$
. Other sources of TKE are negligible in comparison.
Figures 6 and 7 show the VSP and horizontal (
$x$
–
$y$
) and vertical (
$x$
–
$z$
) slices of the vertical velocity
$w$
at
$ft/2\unicode{x03C0} =0.30$
and
$0.46$
, respectively, for one of the fronts that comprise the filament. At the earlier time, the contribution of the turbulence to VSP is greater near the surface, and the horizontal slice shows some organisation within the front, though there is significant down-front variation.
Once the instability saturates and gives way to turbulence, there is a large production of TKE via VSP. The saturated regions in figure 6(a,b) above
$z=-50\text{ m}$
are associated with a larger VSP signature. The maximum in VSP moves down over time, indicating that the instability grows in a non-normal fashion, and figure 7(c) shows the horizontal state as the somewhat coherent cells are destroyed by secondary instabilities and the turbulence begins to isotropise in the
$x$
–
$y$
plane. The extraction of energy is fast, and can be seen in figure 5 to grow, peak and subside within the first inertial period. The consequence is that the filament is now out of geostrophic balance, as shown in the next subsection.
Supplementary movie 1 is a video time series of the vertical velocity and mean contribution to VSP as shown in figures 6 and 7 or similar. Appendix B presents changes in behaviour as the large initial noise (figure 2) is lowered. Doing so recovers the symmetric normal modes as expected of the low-
$Ri$
initial conditions. While qualitatively distinct, the effect of the normal-mode SI on the filament is largely the same as the high-noise case, insofar as this study is concerned. That is, even though normal-mode SI grows somewhat more slowly than the non-normal instability that we just described, it remains fast enough to cause a rapid imbalance of the jet, followed by a filament-scale inertial oscillation, similar to the one that we are about to describe.
3.3. The initial instability echoes throughout the front as an inertial oscillation
Following the onset of instability, the evolution of the filament is not balanced. Indeed, the process in the previous subsection happens within less than an inertial period, and in the resulting evolution, the front does not stabilise at some minimum width. We highlight the consequences of the short-lived turbulent energy production using the linearised Sawyer–Eliassen equation for the mean secondary circulation streamfunction
$\psi$
forced by by-products of turbulent vertical fluxes as we measure them in our simulations. That is, we use
\begin{equation} {\boldsymbol\nabla} ^2 {\psi }_{,{tt}} + \underbrace {f^2 {\psi }_{,{zz}} - J(fv_0,{\psi }_{,{z}})+ J(b_0, {\psi }_{,{x}})}_{L} = \underbrace {f {\overline {w^{\prime}v^{\prime}}}_{,{zz}} - {\overline {w^{\prime}b^{\prime}}}_{,{xz}}}_{G}, \end{equation}
in which we have neglected
$O(\psi ^2$
) terms, and where
\begin{align} (\overline {u}, 0, \overline {v}) = (-{\psi }_{,{z}}, 0, {\psi }_{,{x}}) \quad \text{and}\quad J(a, b) = {a}_{,z} {b}_{,{x }} - {b}_{,{ z}} {a}_{,{ x}}. \end{align}
In doing so, we neglect terms due to the sub-grid momentum and buoyancy fluxes, as well as the remainder of the turbulent flux terms. The turbulent velocity flux
$\overline {w^{\prime}v^{\prime}}$
induces a forcing on
$\psi$
. Integrating this equation over the left half of the interior mixed layer (IML) gives a leading-order evolution equation for the circulation around the corresponding front
$C_{{IML}} = -\iint {\boldsymbol\nabla} ^2\psi \,\text{d} x\,\text{d} y$
, namely,
\begin{align}\notag &\frac {\text{d}^2C_{{IML}}}{\text{d}t^2} \\ &\quad = \int _{-0.9H}^{-0.1H}\text{d}z \int ^0_{-L_x/2} \text{d}x \left[f^2 {\psi }_{,{zz}} - J(fv_0,{\psi }_{,{z}})+ J(b_0, {\psi }_{,{x}}) - f\, {\overline {w^{\prime}v^{\prime}}}_{,{zz}} + {\overline {w^{\prime}b^{\prime}}}_{,{xz}} \right]. \end{align}
The frontal circulation appears dominated by decaying near-inertial oscillations. The peak in
$G$
is large, and occurs in less than one inertial period, before significant change in the mean state; this indicates that the subsequent behaviour will be unbalanced and lead to oscillations. To confirm, a comparison of the terms in (3.4) is shown for each of the simulations in figure 8. The circulation shows clear decaying oscillatory behaviour with period
$2\unicode{x03C0} / f$
, with linear terms dominating for small
$Ri_{min}$
. The strength of these oscillations is dependent on the Richardson number of the flow: a smaller
$Ri_{min}$
produces a faster and stronger extraction of kinetic energy, and in turn leads to higher-amplitude oscillations.The secondary circulation is correlated with the surface front strength. The maximum front strength in figure 4(a) emerges at
$ft/2\unicode{x03C0} \approx 1$
, at the same time that isopycnals begin to move apart. This coincides with the moment
${\text{d}^2C}/{\text{d}t^2} =0$
, and the circulation changes sign (figure 8
a).
Figure 8. The circulation around the left half of the interior mixed layer as in (3.4). The black line is the actual evolution of the circulation, i.e. the left-hand side of (3.4). The blue line is the contribution of the linear terms, the orange line is the contribution of the by-products of turbulent flux terms, and the green line is their sum. Here,
$L$
and
$G$
are defined as in (3.2). All time series are smoothed with a Gaussian kernel with standard deviation
$0.1/f$
.
4. Discussion and conclusions
We investigated the evolution of low-Richardson number mixed layer filaments with numerical simulations. Energy is quickly extracted, via vertical shear production (VSP), from the down-front velocity. This puts the filament in an unbalanced state, and results in vertically sheared inertial oscillations that strengthen horizontal buoyancy gradients near the surface, while weakening them near the base of the mixed layer. The oscillations are well described by the linear terms in the Sawyer–Eliassen equation for the down-front streamfunction, following a forcing induced by the vertical turbulent flux of down-front momentum.
The specific nature of the instability is outside the scope of this paper, but will no doubt be of interest to readers. In short, we believe it to be a non-normal, transient response towards symmetric instability (SI) that is sensitive to the initial noise, rather than normal-mode SI. While we leave a full demonstration to future studies, we find a comparison with reduced-noise simulations particularly illuminating. We present these results in Appendix B. The initial noise illustrated in figure 2 has a large magnitude, compared to the mean flow; and its vertical profile is responsible for the apparent downward propagation of the VSP contribution from figures 6 and 7. At low levels of initial noise, the unstable modes that emerge are down-front symmetric (figure 9) before being destroyed by secondary instabilities (not shown). Perturbations emerge at the same rate throughout the mixed layer in a normal-mode fashion, rather than growing fastest at the surface. The low-noise behaviour is similar to the uniform-gradient Eady state studied by Wienkers et al. (Reference Wienkers, Thomas and Taylor2021a ). The authors observe a similar instability-to-inertial-oscillations pathway that begins with a strong SI in 2.5-dimensional direct numerical simulations of finite-width, zero Richardson number frontal zones (i.e. constant lateral buoyancy gradients in domains that are periodic in the cross-front direction). The generation of inertial oscillations and subsequent surface convergence and frontal strengthening is agnostic to the method of extraction of kinetic energy, as long as it is rapid enough.
Figure 9. Vertical slices of vertical velocity of the simulation with smallest pre-initial noise amplitude (
$10^{-8}$
) at (a)
$ft/2\unicode{x03C0} = 0.8$
and (b)
$ft/2\unicode{x03C0} = 1.0$
, the time of maximum
$\text{VSP}_0$
(see figure 10). The red dashed line in (b) locates
$z=-50\text{ m}$
, the depth of the
$x$
–
$y$
slice of vertical velocity displayed in (c). (c) A horizontal slice of vertical velocity at
$z=-50\text{ m}$
at the time of maximum
$\text{VSP}_0$
.
Arobone & Sarkar (Reference Arobone and Sarkar2015), Zemskova et al. (Reference Zemskova, Passaggia and White2020) and Kimura (Reference Kimura2024) investigate the initial and long-time behaviour of a symmetrically unstable flow. The initial and normal growths of linear perturbations in such a flow are different. In particular, the perturbations are initially strongly asymmetric and not aligned with isopycnals, before becoming nearly symmetric for large times (Arobone & Sarkar Reference Arobone and Sarkar2015). The instability at the start of the simulations in this paper is such an initial response of the filament associated with SI. Non-normal perturbations, seeded by the initial noise, do not have time to grow into normal-mode SI before becoming significantly nonlinear and being isotropised by secondary instabilities acting on the growing modes, and by interaction with the remainder of the initial noise, leading to the rapid energy extraction from the mean state via the VSP. The role and applicability of linear instability analyses in the real ocean, where internal waves, inertial waves and convection contribute to noise, which may strongly interact with both normal and transient forms of instability, is not clear and requires attention in future studies.
Submesoscale structures in a realistic, time-varying ocean can show inertial oscillations associated with the diurnal variation of vertical diffusivity as shown by Dauhajre & McWilliams (Reference Dauhajre and McWilliams2018). We demonstrate here that the changes in effective diffusivity need not come from changes in surface forcing, but can arise without forcing, from turbulent kinetic energy production due to submesoscale instabilities acting on a front. The resulting inertial oscillations contribute to surface frontal strengthening, but originate from unstable initial conditions, rather than a general feature of mixed layer frontogenesis in the real ocean, where instabilities, fronts and turbulent mixing co-exist and evolve together. With similarly unstable initial conditions, and the additional presence of a persistent turbulence-driven secondary circulation, we suspect that the inertial oscillations will be superimposed on a monotonic, turbulent thermal wind frontogenesis.
Frontal features created in the prototypical sense, via a confluent strain flow, cross a range of scales before they reach the regime of the filament presented here. We believe that this work indicates there is understanding to be gained by spinning up fronts in situ for idealised numerical studies of turbulence, rather than beginning with what we think a submesoscale front ought to look like – in particular, the transition from
$Ro \ll 1$
to
$Ro \sim\! 1$
frontogenesis, as these scales remain poorly resolved by global ocean models. Of particular interest here, Thomas (Reference Thomas2012) showed that the behaviour of near-inertial oscillations is significantly altered in regions of frontogenesis driven by strain. Near-inertial waves are captured by growing fronts, which quickly rotate the wavevector to align with isopycnals, before kinetic energy extraction via the frontogenetic secondary circulation damps them completely. Therefore, we believe that in large-eddy simulations combining vertical convection and mixed layer strain-driven frontogenesis across scales, near-inertial oscillations of the nature described in this paper will not be present.
Submesoscale fronts and filaments are the sites of various instabilities, the consequences on the global ocean and best parametrisation of which remain an open question. Here, we have shown that instability alone is capable of strengthening surface buoyancy gradients in a mixed layer filament via a multi-stage process in which the initial instability, which is likely a non-normal, transient response to initial noise, acts quickly and puts the fronts out of balance. The resulting adjustment begins with a convergent secondary circulation on the scale of the fronts that strengthens them near the surface. For our experiments, arrest of surface convergence is not exclusively due to turbulence, and is at least in part the restoring phase of the inertial oscillation. We also comment on the consequences for studies of submesoscale flows, and suggest that the most fruitful way to investigate how these structures evolve in the real ocean, from a computational perspective, may require spinning up – rather than imposing – initial conditions, e.g. starting with a larger-scale mixed layer front and applying an advective forcing, to avoid energising unwanted dynamics due to instability or otherwise.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.10348.
Acknowledgements
The authors ackowledge fruitful discussions with W.R. Young, L.N. Thomas, J. Taylor and M.N. Crowe, as well as three anonymous reviewers, for their illuminating critique and suggestions. Computations were performed on the Mist supercomputer at the SciNet HPC Consortium (Loken et al. Reference Loken2010; Ponce et al. Reference Ponce2019). SciNet is funded by Innovation, Science and Economic Development Canada, the Digital Research Alliance of Canada, the Ontario Research Fund: Research Excellence, and the University of Toronto.
Funding
E.A. and N.G. acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) (funding reference number RGPIN-2022-04560).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
Replication code for this paper can be found at https://doi.org/10.5281/zenodo.14852129 (Atkinson Reference Atkinson2025).
Appendix A. Parametrisation of the reference state
The reference state of the filament is assumed to be independent of the down-filament direction and described by a buoyancy profile
\begin{equation} b_0(x, z) = N_0^2 z + (N^2 - N_0^2)\lambda H\; g\!\left (\frac {z - h(x)}{\lambda H}\right), \end{equation}
with a jet
$v_0(x, z)$
in thermal wind balance (
$f\, {\partial v_0}/{\partial z} = {\partial b_0}/{\partial x}$
), with
$v_0\rightarrow 0$
as
$z \rightarrow - \infty$
. The above equations are defined in terms of a filament shape function
\begin{align} h(x) = H\gamma \left (\frac {x}{\ell }\right)\quad \text{with}\quad \gamma (s) = -1 + \frac {\delta }{2} \left [ \text{erf} \left ( s + \frac {1}{2\alpha }\right) - \text{erf} \left ( s - \frac {1}{2\alpha }\right)\right]. \end{align}
and mixed layer transition function
$g(\xi) = \ln (1 + e^{\xi })$
. Here,
$N_0$
and
$N$
are deep-water and mixed layer buoyancy frequencies, respectively,
$f$
is the Coriolis frequency,
$L$
is the distance between front centres,
$H$
is the depth of the mixed layer, and
$\ell$
is the frontal half-width. Remaining (non-dimensional) parameters are described in table 1. We parametrise the state by varying
$N$
and
$N_0$
to achieve prescribed minimum local Rossby and Richardson numbers, given respectively by
\begin{align} Ro_{{min}} = \frac {(N^2 - N_0^2)\beta ^2}{f^2\alpha ^2} \max _s \left (\gamma \gamma _{,s}\right)\quad \text{and}\quad Ri_{{min}} = \frac {f^2 \alpha ^2 N^2}{\beta ^2 (N^2 - N_0^2)^2} \min _s ( \gamma _{,s}^{-2}). \end{align}
Target Rossby and Richardson numbers are given in table 1.
Table 1. Filament parameters.
Appendix B. Transient behaviour and the impact of initial noise
In the main body of this paper, we outline the process leading to inertial oscillations in an unstable front. In this appendix, we will attempt to connect the topic of this study with SI literature. We find that the emergence and propagation of the initial disturbance throughout the filament is dependent on the amplitude and form of the fluctuations in the initial conditions. This is beyond the scope of the main body in this paper, which is concerned with the consequences of unstable initial conditions for frontogenesis studies, but we include this secondary set of simulations and analyses to provide a more thorough exploration of the generation of inertial oscillations, as well as to demonstrate that the conclusions in the main body of the text hold for rather different circumstances.
Due to the choice of pre-initial conditions, the initialisation described in § 2 produces a state with a turbulent kinetic energy (TKE) profile that is approximately linear in
$z$
(from zero at the base of the mixed layer, to a maximum near the thermocline). This serves as an analogue for the TKE of a state driven by surface cooling. This profile is the cause of the downward propagation of the VSP from figure 6 to figure 7. Fluctuations simply start larger at the surface. Altering the noise profile to be constant in depth in the mixed layer, for example, causes the VSP to emerge approximately at the same time throughout the unstable region (not shown).
Here, we present a set of three simulations, each with different values for the maximum (surface) amplitude of the pre-initial noise. The simulation in the main body of this paper has maximum pre-initial noise amplitude a factor
$10^{-2}$
of the maximum down-front velocity in the reference state. In this appendix, we repeat the analysis in the main body of the paper, but for simulations pre-initialised with noise of relative magnitudes
$10^{-4}$
,
$10^{-6}$
,
$10^{-8}$
.
For the simulation with smallest pre-initial noise amplitude, the form of the vertical velocity at early times is shown in figure 9, and a time series of vertical velocity slices and VSP (as in figures 6 and 7) is shown in supplementary movie 2. Characteristic SI modes are recovered. These modes are somewhat tilted with respect to the isopycnals in the
$x$
–
$z$
plane, expected for non-hydrostatic SI (Wienkers et al. Reference Wienkers, Thomas and Taylor2021a
).
The down-front mean defined in the earlier sections does not remove the large across-front wavenumber symmetric modes, and as such is not appropriate or illustrative to describe a coarse-grained state. With this in mind, and for this appendix only, we define the VSPs and LSPs, as well as the vertical buoyancy flux, using the departure of the simulation state from the reference filament state, namely,
\begin{equation} \text{LSP}_0 = -\int _{\mathcal{V}} \text{d}V\, \frac {\partial \boldsymbol{u}_0}{\partial x}\boldsymbol{\,{\boldsymbol\cdot}\, } u (\boldsymbol{u} - \boldsymbol{u}_0), \quad \text{VSP}_0 = -\int _{\mathcal{V}} \text{d}V\, \frac {\partial \boldsymbol{u}_0}{\partial z}\boldsymbol{\,{\boldsymbol\cdot}\, } w (\boldsymbol{u} - \boldsymbol{u}_0), \end{equation}
and
\begin{equation} \text{BFLUX}_0 = \int _{\mathcal{V}} \text{d}V\, w(b - b_0). \end{equation}
Figure 10 shows the shear production and buoyancy flux terms for a set of simulations with decreasing amplitude of initial noise. In addition, it includes the VSP as defined in (2.8). The
$\text{VSP}_0$
is the primary source of energy for the fluctuations, and behaviour across all cases up to the extremum in
$\text{VSP}_0$
is largely identical; though, for the largest two noise amplitudes (
$10^{-2}$
, in figure 5, and
$10^{-4}$
), the
$\text{VSP}_0$
growth and peak occur earlier. As expected, VSP associated with the perturbations vanishes as they become more symmetric at lower noise levels.
Figure 10. As in figure 5, but for a series of simulations with decreasing amplitudes of initial noise relative to the maximum velocity, and using the TKE production definitions in (B1) and (B2). (a,b,c) The simulations are pre-initialised with noise of relative magnitudes
$10^{-4}$
,
$10^{-6}$
and
$10^{-8}$
, respectively. Terms are shown relative to the same scale as in figure 5. Additionally, and for comparison, the VSP calculated according to (2.8) is included as a dashed green line. Marked points in (c) indicate the times at which the snapshots of vertical velocity in figure 9 are taken.
Figure 11. As in figure 8, but for a series of simulations with decreasing amplitudes of initial noise relative to the maximum velocity. (a,b,c) The simulations are pre-initialised with noise of relative magnitudes
$10^{-4}$
,
$10^{-6}$
and
$10^{-8}$
, respectively. All time series are smoothed with a Gaussian kernel with standard deviation
$0.3/f$
.
We continue with these alternative definitions, noting that the perturbations from the reference state will include the down-front-mean ageostrophic circulation, and may be appropriate only at early times. Nevertheless, figure 5 was reproduced for the first set of simulations with the above definitions (not shown) and bears only quantitative differences.
Following the TKE production by
$\text{VSP}_0$
, inertial oscillations are also observed in the secondary circulations, as seen in figure 11. Hovmöller plots in figure 12 show that the near-surface isopycnals follow the inertial oscillations, with a large front strength and enhanced dissipation just after the filament reaches a minimum width (
$ft / 2\unicode{x03C0} \approx 1.7$
).
Figure 12. (a) A Hovmöller plot of the horizontal buoyancy gradient of the down-front mean state for the simulation with smallest pre-initial noise amplitude (
$10^{-8}$
), averaged over the top 10 % of the mixed layer. Black contours are contours of the buoyancy averaged over the same region. (b) The down-front mean sub-grid dissipation of kinetic energy
$\overline {\varepsilon } = \overline {\nu\, {\boldsymbol\nabla} \boldsymbol{u}:{\boldsymbol\nabla} \boldsymbol{u}}$
, as in (a).
Open access
