The influence of front strength on the development and equilibration of symmetric instability. Part 1. Growth and saturation

Abstract

Submesoscale fronts with large horizontal buoyancy gradients and $O(1)$ Rossby numbers are common in the upper ocean. These fronts are associated with large vertical transport and are hotspots for biological activity. Submesoscale fronts are susceptible to symmetric instability (SI) – a form of stratified inertial instability which can occur when the potential vorticity is of the opposite sign to the Coriolis parameter. Here, we use a weakly nonlinear stability analysis to study SI in an idealised frontal zone with a uniform horizontal buoyancy gradient in thermal wind balance. We find that the structure and energetics of SI strongly depend on the front strength, defined as the ratio of the horizontal buoyancy gradient to the square of the Coriolis frequency. Vertically bounded non-hydrostatic SI modes can grow by extracting potential or kinetic energy from the balanced front and the relative importance of these energy reservoirs depends on the front strength and vertical stratification. We describe two limiting behaviours as ‘slantwise convection’ and ‘slantwise inertial instability’ where the largest energy source is the buoyancy flux and geostrophic shear production, respectively. The growing linear SI modes eventually break down through a secondary shear instability, and in the process transport considerable geostrophic momentum. The resulting breakdown of thermal wind balance generates vertically sheared inertial oscillations and we estimate the amplitude of these oscillations from the stability analysis. We finally discuss broader implications of these results in the context of current parameterisations of SI.

1. Introduction

The upper ocean is a dynamically active and important region, relevant not only to Earth's climate due to exchanges at the air–sea interface, but to biogeochemical processes. Turbulence acts to vertically homogenise this upper-most layer of the ocean down to typical depths of $10$ to $100$ m, driven by wind stresses, surface waves, heat or salinity fluxes or internal flow instabilities. Dynamics in the mixed layer influences exchanges of heat, momentum, carbon, oxygen and other important biogeochemical tracers with the ocean interior.

Fronts, or regions with large lateral density gradients, are common in the upper ocean. These lateral gradients (denoted $\nabla_h$) of the background density field, $\bar {\rho }$ (measured by the horizontal analogue to the buoyancy frequency, $M^2\equiv g/\rho _0 |\nabla _h \bar {\rho }|$, with $g$ the acceleration due to gravity, and $\rho _0$ a reference density) are often in near-geostrophic balance and may be generated by the frontogenetic strains of mesoscale eddies, by coastal upwelling, intrusions into intermediate waters or river discharges. Additionally, persistent frontal systems in the ocean include western boundary currents (e.g. the Gulf Stream and Kuroshio) and the Antarctic Circumpolar Current.

Horizontal density gradients can drive across-front flow due to baroclinic torques, $g/\rho _0 \nabla _h \bar {\rho } \times \hat {\boldsymbol {z}}$ (where $\hat {\boldsymbol {z}}$ is the vertical unit vector). These baroclinic torques tend to flatten isopycnals, but may be counterbalanced by a Coriolis torque, $f \partial _z \bar{\boldsymbol{u}}_{g}$ (where $f$ is the Coriolis parameter), arising from a vertical shear in the geostrophic velocity, $\bar {\boldsymbol {u}}_g$. This geostrophic balance with the horizontal gradient of hydrostatic pressure arising from the background density field is often called the thermal wind balance. The reservoir of potential energy associated with the horizontal density gradient and kinetic energy associated with the thermal wind is available to energise secondary motions. The dynamics within fronts (if not the entirety of frontal systems), however, is often unresolved in global and regional numerical models. A better understanding of these self-regulating frontal dynamics is therefore crucial to modelling the up-scale influence of unresolved processes.

Fronts are susceptible to a number of linear instabilities which drive submesoscale ($100$ m–$10$ km) motions.

Baroclinic instability releases the potential energy stored in the horizontal density gradient, rather than extracting it from the thermal wind shear (Charney 1947; Stone 1972), and is a major mechanism behind the generation of submesoscale eddies (e.g. Boccaletti, Ferrari & Fox-Kemper 2007; Fox-Kemper, Ferrari & Hallberg 2008; Callies et al. 2016). Symmetric instability (SI) is an ageostrophic instability that can develop in frontal regions when the Ertel potential vorticity (PV)

(1.1)\begin{equation} q \equiv (\,f \hat{\boldsymbol{z}} + \boldsymbol{\nabla} \times \boldsymbol{u}) \boldsymbol{\cdot} \boldsymbol{\nabla} b, \end{equation}

(defined with the velocity, $\boldsymbol {u}$, and buoyancy, $b\equiv -g\rho /\rho _0$) is of the opposite sign to the Coriolis parameter, $f$ (Hoskins 1974). The destabilising contributions of a balanced flow are evident if we decompose the PV into a vortical and baroclinic component, respectively

(1.2)\begin{equation} q = (\omega_z + f) N^2 - M^4/f, \end{equation}

where $\omega _z$ is the vertical component of the relative vorticity and $M^2 \equiv \partial _x \bar {b}$ (as above) is the horizontal analogue to the buoyancy frequency, $N^2 \equiv \partial _z \bar {b}$. A negative PV does not necessarily imply SI, however. In the absence of a frontal buoyancy gradient (i.e. $M^2 = 0$) ‘gravitational instability’ occurs when $N^2 < 0$ and $\omega _z + f > 0$ whereas ‘inertial instability’ occurs when $N^2 > 0$ and $\omega _z + f < 0$. Therefore SI only occurs when $(\omega _z + f) N^2 > 0$ but $M^4/f$ is sufficiently large so that $fq < 0$. Much of the ocean interior is sufficiently stratified such that $fq > 0$. However, as noted by Thomas et al. (2016), a frictional stress or diabatic flux at the surface and bottom boundaries lead to $fq < 0$ and trigger SI.

In the context of the Eady model with uniform horizontal and vertical buoyancy gradients, Stone (1966) found that symmetric modes, defined as those independent of the along-front direction (i.e. perpendicular to the horizontal buoyancy gradient), grow faster than baroclinic modes (independent of the cross-front direction) for $Ri < 0.95$, where $Ri\equiv N^2f^2/M^4$ is the balanced Richardson number. Stone (1971) considered the non-hydrostatic contributions to symmetric and baroclinic instabilities in the ageostrophic Eady model, showing that the vertical inertia suppresses both baroclinic and symmetric instabilities. Viscous contributions to the bounded non-hydrostatic SI problem were then included by Weber (1980) and approximated by a viscosity acting on a vertically unbounded normal mode. Beyond the Eady model other types of instability are possible. For example, Wang, McWilliams & Ménesguen (2014) describe a variety of instabilities that develop in more general vertically sheared flows and how they relate to symmetric and baroclinic instability in the Eady model.

Recent observational studies have accumulated evidence of SI in the ocean. For example, increased turbulence and dissipation (exceeding that from atmospheric forcing) in regions where $fq < 0$ has been attributed to SI (Thomas et al. (2016) in the Gulf Stream, and D'Asaro et al. (2011) in the Kuroshio). This negative PV is generated by atmospheric forcing – either by upward buoyancy fluxes (for example cooling) (Haine & Marshall 1998; Thomas et al. 2013) or wind stresses (Thomas & Lee 2005) – which can reduce the PV, and sustain ‘SI turbulence’ and mixing stronger than what the forcing alone could generate. Thompson et al. (2016b) and later Yu et al. (2019) have also found evidence for SI in glider and mooring observations of the open ocean away from major frontal systems. Recently, Savelyev et al. (2018) captured aerial images of SI in the North Wall of the Gulf Stream (cf. figure 2), which constitutes the only visual evidence of the structure of SI to date.

Vertically sheared inertial oscillations of the isopycnals can result from the rapid mixing of geostrophic momentum, and were present following the saturation of SI in the simulations of Taylor & Ferrari (2009). Tandon & Garrett (1994) modelled the response of a mixed layer front to impulsive vertical mixing using the inviscid hydrostatic equations. After a mixing event, the front undergoes inertial oscillations and modulates the background stratification about the average steady-state position ($Ri = 1$). Tandon & Garrett (1994) also considered the case when the vertical stratification is perfectly homogenised (for example by a passing storm), but where the geostrophic shear is partially mixed leaving only a fraction, $s$, of the balanced shear profile. We will show that when acting on times short relative to the inertial period, then SI can generate sufficient geostrophic momentum transport needed to prompt adjustment. We quantify this mixing fraction, $(1-s)$, resulting from the effects of SI.

A number of previous numerical process studies of SI have investigated its nonlinear evolution with varying set-ups, but most have focused only on a single value of the non-dimensional horizontal buoyancy gradient (Thomas & Lee 2005; Taylor & Ferrari 2009, 2010; Thomas & Taylor 2010; Stamper & Taylor 2016). Nonetheless, between persistent fronts, transient fronts and mid-ocean fronts, the strength of these horizontal buoyancy gradients span a large range in the ocean (Hoskins & Bretherton 1972; Jinadasa et al. 2016; Thompson et al. 2016a). We therefore vary the front strength, $\varGamma \equiv M^2/f^2$, rather than changing the vertical stratification as measured by $Ri$.

In this paper, we investigate the equilibration of SI-unstable fronts. We focus on the development and saturation of SI in the Eady model configuration to determine the transport by SI and explain how the rate of energy extraction and amplitude of the resulting inertial oscillations vary with frontal strength. To do this, we first extend the non-hydrostatic and bounded linear analysis of Stone (1971) to include viscosity. We represent the vertical viscous terms using the wave-mode approximation of Weber (1980) to find an analytic solution, but further solve the full numerical eigensystem to verify this approximation in the regime of interest. Compared with Stone (1966) and ensuing papers which studied instability of the Eady model in the inviscid, hydrostatic limit, our analysis is no longer a function only of $Ri$, but now also depends on the front strength, $\varGamma$.

These purely linear analyses are unable to determine the finite contribution of SI to the momentum transport, buoyancy fluxes and energetics of the flow. We analyse the weakly nonlinear problem by considering the growth of secondary instabilities on the growing finite-amplitude SI modes. Our analysis formally extends the work by Taylor & Ferrari (2009) who implicitly considered the secondary shear instability of (unbounded) SI modes by applying the Miles–Howard theorem. We are thereby able to compute a critical amplitude beyond which SI transitions to turbulence and calculate the efficiency with which SI mixes geostrophic momentum prior to transition. To our knowledge this is the first calculation of the mixing fraction, $(1-s)$, (as used by Tandon & Garrett (1994) to describe the geostrophic response of a front) associated with SI.

We begin in § 2 by introducing the problem set-up and primary linear stability analysis for SI. In § 2.3, we consider the stability of these growing SI modes to secondary shear instability and find a critical mode amplitude beyond which the front transitions to turbulence. We finally combine these two stability analyses in §§ 4 and 5 to determine the finite-amplitude contributions of SI to the energetics and momentum transport, respectively.

In a companion paper, we explore the nonlinear consequences of these findings beyond the saturation point. We extend the numerical simulations (from § 3 here used for validation) to study the evolution of these fronts following SI. We use the framework of Tandon & Garrett (1994) to shed light on the effects of dissipation and a finite mixing time on the adjustment response and resulting inertial oscillations.

2. Linear stability analysis

Perhaps the simplest model of a front, the Eady model was first introduced by Eady (1949) and later used by Stone (1966, 1970) to study ageostrophic instabilities. As illustrated in figure 1, the Eady model can be viewed as a local idealisation of a submesoscale mixed layer front where the bottom of the mixed layer is replaced with a flat, rigid boundary. Specifically, an incompressible flow in thermal wind balance with uniform horizontal and vertical buoyancy gradients is bounded between two rigid, stress-free horizontal surfaces.

Figure 1. Schematic of a model frontal region showing coloured contours of density varying both across the front and vertically. The across-front stratification is balanced by the thermal wind shear in $\bar {v}_g$, shown on the top face. A local horizontally homogeneous model can be constructed by considering the region within the grey box, where the buoyancy gradient is approximately uniform.

Non-dimensionalising the Eady problem such that the thermal wind shear, $M^2/f$, is unity in units where the vertical domain size, $H = 1$, brings out four dimensionless parameters

(2.1a–d)\begin{equation} \varGamma \equiv \frac{M^2}{f^2}; \quad Re \equiv \frac{H^2 M^2}{f \nu}; \quad Ri \equiv \frac{N^2 f^2}{M^4}; \quad Pr \equiv \frac{\nu}{\kappa}.\end{equation}

Here, $\nu$ is the kinematic viscosity and $\kappa$ is the diffusivity of buoyancy, but we take $Pr = 1$. It should be noted that the Rossby number is not a parameter in this local frontal zone configuration because there is no horizontal length scale.

We consider a range of front strengths, $\varGamma =M^2/f^2 \approx [1,100]$, which covers a wide variety of ocean fronts. Although very strong fronts with $\varGamma >100$ have been observed (e.g. Sarkar et al. 2016), these fronts are typically very narrow and hence our assumption of a uniform horizontal density gradient is expected to break down. The development of SI at very strong ($\varGamma >100$) and narrow fronts will be reserved for future work.

2.1. Governing equations

Here, we invoke the Boussinesq approximation with a linear equation of state. We further assume that the Coriolis parameter, $f$, is constant and neglect the non-traditional Coriolis terms (i.e. those proportional to $\tilde {f} = 2\varOmega \cos \phi$, where $\varOmega$ is the angular velocity and $\phi$ is latitude). This ‘traditional’ approximation is made here for simplicity but is shown in Appendix A to not qualitatively change our conclusions. The resulting non-dimensionalised Boussinesq equations are

(2.2a)$$\begin{gather} \frac{\textrm{D} \boldsymbol{u}^*}{\textrm{D} t^*} ={-}\boldsymbol{\nabla}^*\varPi^* - \frac{1}{\varGamma} \hat{\boldsymbol{z}}\times \boldsymbol{u}^* + \frac{1}{Re} \nabla^{*2}\boldsymbol{u}^* + b^* \hat{\boldsymbol{z}}, \end{gather}$$

(2.2b)$$\begin{gather}\frac{\textrm{D} b^*}{\textrm{D} t^*} = \frac{1}{Re} \nabla^{*2} b^*, \end{gather}$$

(2.2c)$$\begin{gather}0 = \boldsymbol{\nabla}^*\boldsymbol{\cdot} \boldsymbol{u}^*. \end{gather}$$

Consistent with the non-dimensional parameters (2.1a–d) introduced above, the dimensionless ($^*$) variables here are

(2.3a–e)\begin{equation} \boldsymbol{u}^* \equiv \boldsymbol{u}\frac{f}{H M^2}; \quad b^* \equiv b \frac{f^2}{H M^4}; \quad t^* \equiv t \frac{M^2}{f}; \quad \boldsymbol{x}^* \equiv \boldsymbol{x}\frac{1}{H}; \quad \boldsymbol{\nabla}^* \equiv H \boldsymbol{\nabla}. \end{equation}

The dimensionless pressure head acceleration, $\boldsymbol{\nabla}^* \varPi ^*$, absorbs the hydrostatic pressure gradient, and is eliminated when writing (2.2a) with the along-front streamfunction.

We choose $\hat {\boldsymbol {x}}$ to be the across-front direction (parallel to $\nabla _h \bar {b}$). The background basic state (denoted by an overbar) used for linearisation and the initial condition for the numerical simulations is

(2.4)\begin{equation} \left. \begin{array}{c@{}} \bar{v}^* = z^* - 1/2 \\ \bar{b}^* = \varGamma^{{-}1} x^* + Ri z^* \end{array} \right\}, \end{equation}

as shown in the grey shaded region of figure 1. Following the Eady model, we will also use solid horizontal boundaries at $z^* = 0$ and $1$ which are taken to be insulating and stress free (Eady 1949). In what follows we will omit the appended asterisks for notational simplicity. All variables are dimensionless unless the units are explicitly stated (as in some figures).

2.2. Primary instability

We begin by linearising the Boussinesq equations (2.2) about the basic state (2.4) to describe the evolution of small anomalies in buoyancy and momentum, denoted with a prime. Since the most unstable mode of SI is independent of the along-front ($\hat {\boldsymbol {y}}$) direction (Stone 1966), we consider linear perturbations that vary only in $x$ and $z$:

(2.5)\begin{equation} \left. \begin{array}{c@{}} \displaystyle \partial_t u' ={-}\partial_x \varPi' + \dfrac{1}{\varGamma} v' + \dfrac{1}{Re}\nabla^2 u' \\ \displaystyle \partial_t v' + w'\partial_z\bar{v} ={-}\dfrac{1}{\varGamma} u' + \dfrac{1}{Re}\nabla^2 v' \\ \displaystyle \partial_t w' ={-}\partial_z \varPi' + \dfrac{1}{Re}\nabla^2 w' + b'\\ \displaystyle \partial_t b' + u'\partial_x\bar{b} + w'\partial_z\bar{b} = \dfrac{1}{Re}\nabla^2 b'\\ 0 = \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}' \end{array}\right\}. \end{equation}

We transform this set of partial differential equations (PDEs) into a set of ordinary differential equations (ODEs) by further assuming normal mode perturbations autonomous in $x$ (with wavenumber $k_x$) and in time (with frequency $\omega$) of the form

(2.6)\begin{equation} \chi'(x,z,t) = \Re [ \hat{\chi}(z) \, \mathrm{e}^{\mathrm{i}(k_x x - \omega t)} ], \end{equation}

where the eigenfunction, $\hat {\chi }(z)$, must then be chosen to satisfy the relevant boundary conditions. The set (2.5), after substitution and simplification using the streamfunction defined by $(u',w') = \boldsymbol {\nabla }\times \psi \hat {\boldsymbol {y}}$, becomes

(2.7)\begin{equation} \left(\mathrm{i}\omega + \frac{1}{Re}({-}k_x^2 + D^2)\right)^2({-}k_x^2+D^2)\hat{\psi} = \left(-\frac{1}{\varGamma^2}D^2 -\frac{2\mathrm{i}k_x}{\varGamma}D + k_x^2 Ri\right) \hat{\psi}, \end{equation}

where $D \equiv \mathrm {d}/\mathrm {d}z$ for notational ease. Note that this equation is closely related to (14) in Grisouard & Thomas (2016) who formulated the equation in terms of pressure and neglected horizontal diffusion. The boundary conditions for $\hat {\psi }$ at $z = 0, 1$ on this sixth-order ODE are

(2.8)\begin{equation} \left. \begin{array}{c@{}} \hat{\psi} = 0\\ D^2 \hat{\psi} = 0 \end{array} \right\}. \end{equation}

To make this system tractable, we follow the method of Weber (1980) and approximate equation (2.7) as a second-order ODE by writing the vertical diffusion terms as spatially invariant wave modes,

(2.9)\begin{equation} \frac{1}{Re} D^2 \hat{\psi} \approx{-}\frac{1}{Re} k_z^2 \hat{\psi}, \end{equation}

with vertical wavenumber $k_z$. By neglecting the vertical variations in $k_z$, this approximation constrains the SI mode angle to be uniform in $z$. This is a good approximation for large $Re$ and $k_z$, when the effects of diffusion are dominated by the interior of the domain. This does consequently prohibit the boundaries from generating vorticity, but it is found to not influence the selection or stability of SI, which is only energised by the bulk background buoyancy and shear. Equation (2.7) then becomes

(2.10) \begin{align} \left[\left(\mathrm{i}\omega - \frac{k_x^2 + k_z^2}{Re}\right)^2 + \frac{1}{\varGamma^2}\right] D^2\hat{\psi} + \frac{2\mathrm{i}k_x}{\varGamma} D\hat{\psi} -k_x^2\left[\left(\mathrm{i}\omega - \frac{k_x^2+k_z^2}{Re}\right)^2 + Ri\right] \hat{\psi} = 0,\end{align}

which has eigensolutions of the form

(2.11)\begin{equation} \hat{\psi} = \exp(\mathrm{i} \lambda_1 z) - \exp(\mathrm{i} \lambda_2 z),\end{equation}

that match the boundaries if $\lambda _1 - \lambda _2 = 2{\rm \pi} n$, for the chosen eigenmode number, $n$. Equation (2.10) is thus reduced to a quadratic eigenproblem which may be solved by numerical iteration while enforcing the vertical viscous wave-mode approximation that

(2.12)\begin{equation} k_z^2 = \tfrac{1}{2} (\lambda_1^2 + \lambda_2^2). \end{equation}

Complete details of this solution are included in Appendix B.1.

The exact numerical eigensolution to the linear set (2.5) was also computed using a pseudo-spectral eigenvalue solver written in Matlab. The computed solutions to (2.10) give good agreement with this numerical solution, as shown in figure 2(a), where the growth rate, $\sigma$, is the imaginary part of $\omega$. This new solution correctly accounts for both the limiting effects of the vertical boundaries at low wavenumber, and of viscosity at high wavenumber. Accurate in the low wavenumber limit, Stone (1971) determined this inviscid, bounded solution, where the mode growth becomes suppressed as it feels the constraint of the boundaries for $k_x \lesssim 2{\rm \pi}$. In the other limit of unbounded, viscous and hydrostatic motions, Taylor & Ferrari (2009) (and later Bachman & Taylor (2014) for non-hydrostatic motions) found that the most unstable mode has a vanishing wavenumber. The structure of the exact ($n = 1$) viscous, bounded SI mode ($u'$) is shown in the background of figure 4. Due to viscous and non-hydrostatic effects, the modes are no longer parallel to isopycnals as they were in e.g. Stone (1966).

Figure 2. (a) The growth rate for the $n=1$ SI mode in a vertical front with $\varGamma = 10$ and $Re = 10^5$. The real part of $\omega$ for the SI modes are everywhere $0$ except where linearly stable at very small wavenumber. (b) The growth rate of the fastest growing SI mode ($n=1$) and wavenumber at $Re = 10^5$, as a function of $\varGamma$. The vertical front ($Ri = 0$) is shown in black and also for increasing stratification.

We can now consider how the fastest growing mode of (2.10) varies with $\varGamma$ and $Ri$, as shown in figure 2(b). For $Ri=0$ the energy growth rate relative to $f$ increases nearly linearly with front strength. However, for strong fronts stratification significantly reduces the growth rate of the most unstable modes.

In a vertically unbounded domain with an inviscid, hydrostatic dynamics, the maximum release of energy can be achieved by motion aligned with $b$ surfaces, with ${\theta _b = \tan ^{-1}(M^2/N^2)}$ from the horizontal (i.e. $k_x/k_z = M^2/N^2$) (Taylor & Ferrari 2009), effectively precluding any buoyancy flux. However, in a vertically bounded front with weak stratification, the most unstable modes become very inclined to the isopycnals as shown in figure 3(a), and reach nearly $45^\circ$ for $N^2 = 0$. While the angle of the unstable SI modes must still be between the angle of the isopycnals and surfaces of constant absolute momentum ($\bar {m} = \bar {v}_g + \varGamma ^{-1} x$) (dotted and dash-dotted curves in figure 3a), the most unstable modes approach more closely to the angle of the absolute momentum surfaces ($\theta _m = \tan ^{-1}\varGamma ^{-1}$) for small front strength. This permits a larger buoyancy production of energy ($\mathcal {B} = \langle w'b' \rangle$), as shown in figure 3(b), while the geostrophic shear production ($\mathcal {P}_g = -\langle \overline {v'w'} \partial _z \bar {v}_g \rangle$) is the dominant energy source in the rest of the parameter space. Here and throughout the rest of this paper, $\langle \cdot \rangle$ indicates a volume average over the entire domain, and primes represent local departures from the horizontally averaged fields denoted by $\bar {\cdot }$.

Figure 3. (a) The angle of the fastest growing SI mode as measured from horizontal, plotted as a function of front strength, $\varGamma$, and for different background stratifications measured by the inverse isopycnal slope, $N^2/M^2$. The shaded grey region indicates where $fq>0$ and the front is stable to SI. The unstable SI mode inclination must remain between the angle of absolute momentum surfaces ($\theta _m$, dot-dashed line) and isopycnals ($\theta _b$, dotted lines), which for $N^2/M^2 = 0$, $\theta _b = 90^\circ$. This unstratified case has modes nearly equally spaced between the isopycnals and absolute momentum surfaces for large $\varGamma$, but with increasingly horizontal isopycnals the SI modes grow more along these isopycnals. While the angle of the contour $\psi (x,z)=0$ is a weak function of $z$ in the full numerical eigensolution (decreasing by at most $5\,\%$ at the boundaries), the mode angle of the solutions (2.11) are independent of $z$. (See Appendix B.1 for details on the calculation of $\theta$ and the eigenfunctions.) (b) Contribution of the most unstable linear SI mode to the energy budget (4.1) of the vertical front for $Re = 10^5$. Normalised by the kinetic energy, the geostrophic shear production and buoyancy flux are relatable to the growth rate, $\sigma$. As expected with SI, the instability still primarily draws energy from the thermal wind shear into the kinetic energy of the mode through the TKE production term. The grey dotted line indicates the growth rate of baroclinic instability for this choice of parameters (Stone 1966). Symbols correspond to the numerical simulations discussed in § 3, computed as a time average from $t=0$ to $\tau _c/2$.

2.3. Secondary instability

Secondary instability plays a key role in the equilibration of SI. Here, we explore the onset of secondary instability to determine the cumulative effects of SI in the front equilibration energetics and the contribution to mixing down the thermal wind shear.

As described in Taylor & Ferrari (2009), shear associated with the growing SI modes becomes unstable to a secondary Kelvin–Helmholtz instability (KHI) which prompts a transition to turbulence. We identify this critical SI mode amplitude, $U_{SI} = U_c$, at which the SI modes themselves break down as the time when the SI growth rate, $\sigma _{SI}$, is equal to the KHI growth rate, $\sigma _{KH}$. Of course, $\sigma _{KH}$ is a monotonically increasing function of the shear, and thereby of $U_{SI}$ which exponentially grows at a rate $\sigma _{SI}$. We therefore iteratively compute the secondary linear stability of the combined Eady and growing SI mode basic state to determine this critical amplitude that is plotted in figure 5(a).

We formulate the one-dimensional (1-D) linear Kelvin–Helmholtz stability problem using a sinusoidal extension of the structure of the full SI mode (evaluated at the mid-plane) in the rotated coordinates shown in figure 4. As described in Appendix C, this basic state includes the constant vertical and horizontal buoyancy gradients associated with the basic state in the Eady model as well as the buoyancy changes induced by the SI modes. We iteratively compute $\sigma _{KH,max}(U_{SI})$ with a pseudo-spectral solver until finding the critical SI mode amplitude, $U_{c}$. While the most unstable SI wavevector, $|\boldsymbol {k}_{SI}|$, increases as the mode number ($n$) and $Re$ increase, the scaling for $U_{c}$ appears to be dominated by $\sigma _{SI}$ and so remains largely unchanged.

Figure 4. Diagram showing the secondary stability analysis coordinate transformation drawn over the linear SI mode ($u'$). The primary SI basic state is also indicated, with grey isopycnal lines showing the linearly increasing buoyancy from left to right (for $Ri = 0$), as well as the thermal wind vectors into the page which balance the baroclinic torques.

We demonstrate this here just for the unstratified ($Ri = 0$) front, but a general analysis is provided in Appendix C. Figure 5(a) shows the classic KH stability analysis (i.e. ignoring rotation and neglecting the along-mode component of the background stratification) alongside the full solution for $U_{c}$. The dashed line shows the resulting scaling,

(2.13)\begin{equation} U_{c} \propto (\sqrt{\varGamma} |\boldsymbol{k}_{SI}|)^{{-}1}, \end{equation}

(in our same dimensionless units of the velocity associated with thermal wind shear). We obtain this scaling by balancing the KHI growth rate (proportional to the non-dimensional shear in the SI mode, $\sigma _{KH}\propto U_{c} |\boldsymbol {k}_{SI}|$) with the SI growth rate in the limit of large $\varGamma$, $\sigma _{SI} \propto \varGamma ^{-1/2}$. We see that this simple scaling argument fails for small $\varGamma$ where the growth rate in these weak fronts is slow compared with $f$.

Figure 5. (a) The critical amplitude of the most unstable SI mode velocity at which secondary instability begins to dominate, shown in units of the thermal wind. The dotted line shows this critical amplitude when rotation and along-shear stratification (i.e. $x^{\dagger}$ in Appendix C) are neglected in the KHI stability analysis. The dashed line shows the scaling (2.13) achieved by taking the KHI growth rate directly proportional to the shear and matching $U_c$ in the limit of large $\varGamma$. (b) The cumulative kinetic energy (KE) budget contributions from the $n = 1$ linear SI mode of the unstratified front, integrated through $U_{c}$. Coloured symbols show the value derived from the 2-D simulations. Due to weak scale and mode selection, these simulations contain a range of $n$ and $k_x$, yet with increasing front strength the values calculated from the simulations approach the $n = 1$ line shown due to stronger mode selection as the higher modes are damped by viscosity.

3. Numerical simulations

We employed the non-hydrostatic hydrodynamics code, diablo, to verify the conclusions of the preceding linear primary and secondary instability theory as well as the results in the following two sections. diablo solves the fully nonlinear Boussinesq equations (2.2) on an $f$-plane (Taylor 2008). Second-order finite differences in the vertical and a collocated pseudo-spectral method in the horizontal periodic directions are employed, along with a third-order accurate implicit–explicit time-stepping algorithm using Crank–Nicolson and Runge–Kutta with an adaptive step size. Rigid, stress-free and insulating horizontal boundaries are enforced to match the linear analysis in § 2. Following Taylor & Ferrari (2009), the simulations are run in a 2-D ($x$–$z$) domain while retaining all three components of the velocity vector. This choice allows us to focus on the evolution and saturation of the symmetric modes.

While the presented analytical results in this paper are general, we focus these numerical verification experiments on an initially unstratified front ($Ri=0$) with ${Re = 10^5}$. It should be noted that the along-front flow would be susceptible to KHI in a 3-D simulation, but this is not considered for the purpose of this study. Each of the simulations were initialised as a balanced front (2.4) with strength ${\varGamma = \{1, 10, 100\}}$ and white noise was added to the velocity with a (dimensionless) amplitude of $10^{-4}$. The simulations were run through the linear phase until secondary instability breaks down the SI modes at the critical time, $\tau _c$, as shown in the right column of figure 6. At this point, we measure the cumulative effects of SI on the front – the integrated shear production, buoyancy fluxes and momentum transport – and present these values alongside the analytical results of §§ 4 and 5. While we restrict these verification simulations to initially unstratified fronts and do not consider times after $\tau _c$, we extend these simulations in a companion paper to explore the SI-induced re-stratification and geostrophic adjustment of the fronts at later times.

Figure 6. Slices across each front show the along-front vorticity, $\omega _y$, along with buoyancy contours (black lines), for $\varGamma = 1$ (top), $10$ (centre) and $100$ (bottom). Two snapshots are shown, at $t = \tau _c/2$ (left) when the fastest linear SI mode has emerged, and at $t = \tau _c$ (right) when secondary KHI first begins to break the coherent energy of the SI modes into small-scale turbulence. Note that the vorticity is normalised by $M$, which keeps the amplitude similar across the range of $\varGamma$ (consistent with the scaling (2.13)). The vorticity normalised by $f$ can be obtained by multiplying the values shown here by $\varGamma ^{-1/2}$. During the linear growth phase (left panels), the SI modes do not align with the isopycnals, and rather become increasingly flat for larger $\varGamma$, consistent with the results shown in figure 3(a).

4. Energetics of SI

In light of these stability analyses, a natural question is: What impact does the linear SI phase and ensuing turbulence have on the resulting equilibration of the front, and how does it depend on the frontal strength? To answer this, we combine the primary linear instability results of § 2.2 with the details of SI saturation from § 2.3 to determine the cumulative contribution of SI modes up to the critical time, ${\tau _c = \sigma _{SI}^{-1}\log (U_c/U_0)}$, when SI has grown to an amplitude $U_{c}$. This allows us to quantify the energetics of the linear SI modes and their influence on the evolution of the front.

With the complex eigenfunction, $\hat {\psi }$, found by iteratively solving for $\lambda _1$ and $\lambda _2$ in (2.11), we determined the full structure of these modes: $\hat {u}$, $\hat {v}$, $\hat {w}$ and $\hat {b}$ as given in Appendix B.1. With these, we compute the correlations relevant to the transport and energetics of the development of SI. We first must normalise each of the modes by $\sqrt {|\hat {u}|^2 + |\hat {w}|^2}$, and then rewrite them in the normal mode form, (2.6), using the parameter and eigenvalues of (2.10).

We will first consider the contribution of SI to the turbulent kinetic energy (TKE), ${E_K \equiv \frac {1}{2} \langle u_i' u_i'\rangle }$:

(4.1)\begin{equation} \frac{\partial E_K}{\partial t} ={-}\underbrace{\left\langle \overline{u'w'}\frac{\partial\bar{u}}{\partial z}\right\rangle }_{-\mathcal{P}_x} - \underbrace{\left\langle \overline{v'w'}\left(\frac{\partial\bar{v}_a}{\partial z} + \frac{\partial\bar{v}_g}{\partial z}\right)\right\rangle }_{-\mathcal{P}_y} + \underbrace{\left\langle w'b'\right\rangle }_\mathcal{B} - \underbrace{\frac{1}{Re}\left\langle \frac{\partial u'_i}{\partial x_j}\frac{\partial u'_i}{\partial x_j}\right\rangle }_{\varepsilon_t}. \end{equation}

The first two terms on the right-hand side represent the shear production, $\mathcal {P}$, converting energy from the mean flow into TKE. Specifically, the along-front contribution ($\mathcal {P}_y$) is split into a geostrophic shear production term, $\mathcal {P}_g$, energised by the thermal wind shear, and an ageostrophic part. The other potential source of TKE comes from buoyancy production, $\mathcal {B}$, which represents the transfer of energy from PE into TKE. The cumulative generation of TKE by each of these terms in (4.1), integrated from $t = 0$ up to transition at $\tau _c$ is shown in figure 5(b). As expected for SI, the contribution from $\mathcal {P}_g$ exceeds $\mathcal {B}$, except for small $\varGamma$. Interestingly, even for these SI modes that are very flat (i.e. inclined to the isopycnals) in strong fronts, the energetics are still dominated by geostrophic shear production which relies on the vertical velocity to exchange geostrophic momentum. We confirm this result using the numerical simulations described in § 3. Even though the initial white noise and weak mode selection mean that a range of wavenumbers are represented in the simulations, these predictions still remain robust.

Following Haine & Marshall (1998), it is possible to re-frame the SI stability criterion, $fq < 0$, in terms of the energy sources driving growth: the background buoyancy gradient and the geostrophic kinetic energy. First consider fluid parcels that are constrained to move along isopycnals (and thus incur no gravitational penalty). The criterion, $fq < 0$, for instability then becomes the Rayleigh criterion describing inertial instability,

(4.2) \begin{align} f\left(\frac{\partial m}{\partial x}\right)_b < 0.\end{align}

where the subscript here indicates the gradient is measured along isopycnals. Since they are aligned with the isopycnals, these modes do not extract potential energy from the front and instead grow by drawing energy from the thermal wind. Contrast this with the other limiting angle that SI modes can take, when their motion is aligned with surfaces of constant absolute momentum. Now, instability requires that the vertical buoyancy gradient measured along these absolute momentum surfaces is negative

(4.3)\begin{equation} \left(\frac{\partial b}{\partial z}\right)_m < 0. \end{equation}

Therefore, motions that are constrained to follow absolute momentum surfaces can extract potential energy, analogously to ‘upright convection’ (Haine & Marshall 1998). For hydrostatic perturbations in an unbounded domain, the most unstable mode of SI is aligned with the isopycnals and hence grows by extracting kinetic energy from the thermal wind through geostrophic shear production (Stone 1972; Haine & Marshall 1998; Taylor & Ferrari 2009).

As shown previously in figure 3(a), for non-hydrostatic modes in a bounded domain, the most unstable mode of SI is not necessarily aligned with isopycnals and hence these modes can grow through a non-trivial combination of buoyancy production and geostrophic shear production. We can quantify the energetic influences on the most unstable mode of SI using the linear stability analysis up to the critical time, $\tau _c$. We do this by introducing the energy production ratio,

(4.4)\begin{equation} \frac{\mathcal{B}}{\mathcal{B}+\mathcal{P}_g} = \frac{\lambda_1 + \lambda_2}{2 k_x \varGamma}, \end{equation}

as plotted in figure 7, where $\mathcal {B}$ is the buoyancy flux, $\mathcal {P}_g$ is the geostrophic shear production and $\lambda _1$ and $\lambda _2$ describe the vertical mode structure (2.11) and depend on $Ri$ and $\varGamma$ (details of which are given in Appendix B.2). The production ratio suggests the expected character of SI. For strong fronts with weak vertical stratification, SI extracts energy from shear production, and so we refer to this flavour of SI as ‘slantwise inertial instability’. In the inviscid and hydrostatic limits, the linear analysis indicates that energy is always fully derived from geostrophic shear production, with modes aligned perfectly with isopycnals at $\theta _b$. Non-hydrostatic effects flatten the SI modes particularly for small $\varGamma$. This permits buoyancy production to contribute to the energy more than shear production, and so we call this flavour of SI ‘slantwise convection’. Note that this term has sometimes been used synonymously with SI in the literature, although it is not always congruous with the energetics of SI (Haine & Marshall 1998). The boundary-permitted viscous limit in figure 7(b) each for large $\varGamma$ and large $Ri$ also exhibits slantwise convection modes. It is perhaps then surprising that within the white outlined region, indicating where SI modes are more aligned with isopycnals, the instability does not always extract a majority of energy from the shear production (i.e. red shading).

Figure 7. Contours of the production ratio (4.4) distinguish regions where geostrophic shear production dominates ($0$) and regions where buoyancy production dominates ($1$). The white line separates regions of parameter space where SI modes are more aligned with isopycnals, i.e. $|\theta -\theta _b|< |\theta -\theta _m|$ (inside), from the regions (outside) where they are more closely aligned with absolute momentum surfaces. (a) The production ratio plotted in parameter space with $N^2/f^2$ on the $y$-axis, chosen so that the axes are only interdependent on $f$. A black dashed line designates the contour $Ri = 0.25$. Lines of constant isopycnal slope ($M^2/N^2$) are straight lines of slope $1$ in this log–log scale. Strong fronts with weak stratification (equivalently, large isopycnal slope) derive energy primarily from geostrophic shear production. Thus, rapid frontogenesis (moving horizontally to the right), or rapid de-stratification via mixing (moving vertically downwards) will tend the SI modes to slantwise inertial instability. (b) The parameter space is rescaled with $Ri$ on the $y$-axis to emphasise the region near $Ri = 1$ where SI in a balanced front becomes stabilised. Non-hydrostatic effects (for small $\varGamma$) and boundary viscous effects (for large $Ri$) influence the SI modes to derive this portion of energy from the background buoyancy gradient. Non-traditional effects also influence how SI extracts energy, as shown by figure 10 in Appendix A.2.

In the ‘slantwise convection’ regime, where $\mathcal {B}>\mathcal {P}_g$, SI tends to be weak and the total energy production is small. This raises the question of whether it is important to account for the SI-driven buoyancy flux in parameterisations of SI. To provide context, we compare the SI-driven buoyancy flux at $\tau _c$ with the buoyancy flux associated with mixed layer instability (MLI) in the parameterisation from Fox-Kemper et al. (2008). They empirically estimated a constant efficiency factor for the finite-amplitude MLI, which in our non-dimensional variables can be written

(4.5)\begin{equation} C_e = \varGamma^{{-}1}\langle \overline{w'b'}\rangle _{MLI} = 0.06-0.08. \end{equation}

In comparison, a typical SI buoyancy flux at $\tau _c$ for the slantwise convective regime (specifically at $\varGamma = 1$ and $Ri = 0$) is

(4.6)\begin{equation} \varGamma^{{-}1}\langle \overline{w'b'}\rangle_{\mathrm{SI},c} = 0.0074. \end{equation}

Note that the buoyancy flux increases in time during the growing phase of SI and MLI. The fact that $\langle \overline {w'b'}\rangle _{\mathrm {SI}}$ at $t=\tau _c$ is smaller than $\langle \overline {w'b'}\rangle _{MLI}$ highlights the comparatively early saturation of SI through secondary instabilities. Thus even though MLI grows more slowly for this set of parameters (cf. the dotted line in figure 3b), the finite-amplitude buoyancy flux associated with MLI has a significantly larger influence on the rate of re-stratification compared with SI for weak fronts ($\varGamma =1$). For stronger fronts with weak vertical stratification (i.e. large $\varGamma$ and small $Ri$) where the geostrophic shear production is larger than the buoyancy flux, SI can indirectly induce re-stratification by first generating large vertical fluxes of geostrophic momentum. This will be discussed in the next section.

5. Momentum transport by SI

We now consider the effect that SI has on the geostrophic shear and the implications for the subsequent response of the front.

5.1. Dominant momentum balance

We would like to determine the dominant terms in the mean horizontal momentum equations to understand what is driving the evolution of the front. Subtracting off the background geostrophic velocity and buoyancy gradient from the Boussinesq equation (2.2a) gives a horizontally autonomous system allowing us to Reynolds average in the $\hat {\boldsymbol {x}}$ and $\hat {\boldsymbol {y}}$ directions. Using continuity and geostrophic balance, the horizontal ageostrophic momentum equations are

(5.1a)$$\begin{gather} \partial_t \bar{u}_a + \partial_z\overline{u'w'} = \hphantom{-}\varGamma^{{-}1}\bar{v}_a, \end{gather}$$

(5.1b)$$\begin{gather}\partial_t \bar{v}_a + \partial_z\overline{v'w'} ={-}\varGamma^{{-}1}\bar{u}_a. \end{gather}$$

To determine the dominant balance at early times arising from the growing SI modes, we first assume the Coriolis term in (5.1b) is small. With this approximation, we construct a ratio from the terms in (5.1a),

(5.2)\begin{align} \frac{\partial_z \overline{u'w'}}{\varGamma^{{-}1} \bar{v}_a} &\approx{-} \frac{\partial_z \overline{u'w'}}{\displaystyle \varGamma^{{-}1} \int_0^{\tau} \partial_z\overline{v'w'} \,\mathrm{d}t} = 2\varGamma^2\sigma^2\frac{\lambda_1+\lambda_2}{\lambda_1 + \lambda_2 - 2k_x \varGamma} \nonumber\\ & \quad \sim 2 \quad \mathrm{for }\ \varGamma \gg 1 \ \mathrm{and} \ Ri = 0, \end{align}

where we have also assumed exponential growth in time, $\propto \exp (\sigma t)$. We take $U_{SI}$ at $t = 0$ to be infinitesimal so that the lower limit of integration evaluates to $0$. The arbitrary upper limit, $\tau$, then cancels with the exponential evaluated at $\tau$ in the numerator. Similarly for the terms in the $y$-momentum equation (5.1b), the ratio is

(5.3)\begin{align} \frac{\partial_z \overline{v'w'}}{-\varGamma^{{-}1} \bar{u}_a} &\approx \frac{\partial_z \overline{v'w'}}{\displaystyle \varGamma^{{-}1} \int_0^{\tau} \partial_z\overline{u'w'} \, \mathrm{d}t} = 2\frac{\lambda_1 + \lambda_2 - 2k_x \varGamma}{\lambda_1+\lambda_2} \nonumber\\ &\quad \sim 2 \varGamma \quad \mathrm{for} \ \varGamma \gg 1 \ \mathrm{and} \ Ri = 0, \end{align}

where we again use the solution for the eigenfunctions (B7) derived in Appendix B.1 to evaluate these integrals. Each of these expressions in (5.1b) are self-consistent with our assumption to neglect the Coriolis term if both ratios are $\gg 1$. We found this to be the case for $\varGamma \gtrsim 1$ and $Ri \lesssim 0.5$ and so we conclude that the mean ageostrophic $y$-momentum is driven more strongly than the $x$-momentum – i.e. the dominant balance is initially $\partial _t\bar {v}_a \approx -\partial _z\overline {v'w'}$.

5.2. Loss of geostrophic balance

This dominant balance with the $\partial _z\overline {v'w'}$ Reynolds stress term suggests that, at first order, SI can influence the large-scale evolution of the front by rearranging the momentum of the balanced thermal wind. The rate at which this geostrophic shear profile is reduced will give hints as to the type of adjustment that follows SI.

Taking the vertical gradient of the dominant momentum balance,

(5.4)\begin{equation} \frac{\partial }{\partial t} ( \partial_z\bar{v} ) \approx{-}\partial^2_z\overline{v'w'},\end{equation}

we can estimate the time scale required to mix the thermal wind shear

(5.5)\begin{equation} \tau_{mix} = \frac{\partial_z \bar{v}_g}{\langle \partial_z^2\overline{v'w'}\rangle _c}, \end{equation}

for SI momentum fluxes evaluated at $\tau _c$. This value is plotted for each $\varGamma$ in figure 8(a), and details of the calculation are saved for Appendix B.3. If this time scale is long compared with $f$ (as for very weak fronts), then we might expect the front to slowly slump over while remaining quasi-balanced. In contrast, when the vertical fluxes rapidly (relative to $f$) mix down the thermal wind shear before inertial effects can influence the large-scale dynamics, then the response can be viewed as a form of geostrophic adjustment. This is the case for $\varGamma \gtrsim 10$.

Figure 8. (a) The mixing time scale (5.5) required for SI-driven fluxes to entirely destroy the thermal wind shear. Strong fronts with steep isopycnals are rapidly mixed (relative to an inertial period) whereas this time scale becomes longer than an inertial period for small $\varGamma$. (b) The thermal wind shear mixing fraction, ($1-s$), induced by SI (5.9). Colours correspond to the different values of inverse isopycnal slope, $N^2/M^2$, shown in the legend at left. Symbols show the extracted values from the 2-D simulations. Momentum is only rearranged within the domain by the linear modes, and so the domain average is conditioned on positive fluxes.

Tandon & Garrett (1994) showed that in the limit of instantaneous mixing (here for $\varGamma \gg 1$) this resulting geostrophic adjustment of the front results in inertial shear oscillations. They considered the evolution of a mixed layer front when a fraction $(1-s)$ of the vertical shear is removed, such that initially

(5.6)\begin{equation} \partial_z \bar{v}|_{t=0} = s \partial_z \bar{v}_g. \end{equation}

The subsequent horizontally invariant inertial oscillations modulate the background stratification by differentially advecting the lateral buoyancy gradient across the front. Assuming the PV remains constant, the (dimensionless) stratification evolves according to

(5.7)\begin{equation} \partial_z \bar{b}_i = (1-s) (1 - \cos (\varGamma^{{-}1}t) ) + Ri. \end{equation}

These inertial oscillations draw closed circular orbits and have a linear structure in $z$

(5.8a)$$\begin{gather} \bar{u}_i ={-}(1-s) (z - 1/2 ) \sin (\varGamma^{{-}1}t), \end{gather}$$

(5.8b)$$\begin{gather}\bar{v}_i = s (z - 1/2 ) + (1-s) (z - 1/2 ) (1 - \cos (\varGamma^{{-}1}t) ). \end{gather}$$

The amplitudes of these inertial shear oscillations are dimensionally $(1-s)M^2/f$.

5.3. Inertial oscillation amplitude

The reduction in thermal wind shear before $\tau _c$ thus should dictate the amplitude of these inertial oscillations in a front following SI. We can estimate this mixing fraction, $(1-s)$, as introduced in Tandon & Garrett (1994). Again using the vertical derivative of the dominant momentum balance (5.4), we compute the cumulative contribution of the SI modes through to $\tau _c$

(5.9)\begin{equation} 1-s = \int_0^{\tau_c} \langle \partial_z^2 \overline{v'w'} \rangle \, \mathrm{d}t ,\end{equation}

(detailed in Appendix B.3). Note that the term on the right-hand side has been non-dimensionalised by $M^2/f$ (consistent with the dimensionless units used throughout this paper) so that $(1-s)$ is interpreted as a fraction of the thermal wind shear. This mixing fraction is shown in figure 8(b). We see that with increasing front strength the linear SI modes are able to remove a larger fraction of the thermal wind shear before $\tau _c$, setting up larger inertial oscillations. While these results combine the analysis of SI with the theory of Tandon & Garrett (1994), in a companion paper we consider the direct and indirect nonlinear effects of SI on the evolution of these inertial oscillations.

6. Conclusions

SI occurs at density fronts in the ocean and atmosphere when the PV takes the opposite sign to the Coriolis parameter, i.e. $fq < 0$. While previous studies have focused on the effect of Richardson number on SI, here we have explored the dependence of SI on front strength, parameterised by $\varGamma = M^2/f^2$, where $M^2$ is the horizontal buoyancy gradient. To that end, we have analysed an idealised model of a frontal region initially in thermal wind balance with a uniform horizontal buoyancy gradient and a constant background vertical stratification. Although highly idealised, this configuration was motivated by rapid mixing events such as the passage of a storm or an event which vertically mixes the buoyancy profile.

Using a linear stability analysis in a vertically bounded domain with viscous and non-hydrostatic effects, we have shown that SI can grow via two routes: by converting kinetic energy associated with the balanced thermal wind into the growing perturbations, or by extracting potential energy from the front via the buoyancy flux. For strong fronts and where $Ri \lesssim 0.5$, the larger contribution energising the instability comes from geostrophic shear production, but for large $Ri$ and/or weak fronts the buoyancy flux is also important. We have characterised the two limiting behaviours of SI distinguished by the dominant energy source: ‘slantwise convective instability’ extracts energy from the background potential energy via buoyancy production with modes tending along absolute momentum surfaces, while ‘slantwise inertial instability’ is energised by shear production and has more upright modes nearly along isopycnals. This finding provides context to the work by Grisouard (2018) on mixed ‘inertial–symmetric instability’. By varying the Rossby number, they found that while the two limiting instabilities extract energy via shear production, buoyancy fluxes can still be important for the mixed modes. Here, we have focussed on pure SI ($\partial _x \bar {v} = 0$), and found that even in this limit the dominant energy source depends on the details of the front. However, for the parameters where the buoyancy flux is the largest energy source (the ‘slantwise convection’ regime), the SI-driven buoyancy flux is small compared with the mixed layer eddy parameterisation of Fox-Kemper et al. (2008). Nonetheless, at stronger fronts SI can induce rapid re-stratification by first generating large vertical fluxes of geostrophic momentum, as parameterised by Bachman et al. (2017).

By extracting energy from the balanced thermal wind, SI leads to re-stratification, and can induce vertically sheared inertial oscillations depending on the strength of the front. The mixing time scale for SI to homogenise the thermal wind shear decreases with front strength, and is faster than an inertial period for $\varGamma \gtrsim 10$. Thus the response to rapid mixing of the thermal wind shear at strong fronts can be described in terms of geostrophic adjustment. We analysed this behaviour in the context of the model used in Tandon & Garrett (1994) which assumed that the PV was constant throughout the adjustment process. Using the linear stability analysis, we estimated the degree to which SI mixes the thermal wind shear and concluded that SI can generate large amplitude inertial oscillations at strong fronts.

In Part 2 of this series (Wienkers, Thomas & Taylor 2021), we consider the nonlinear consequences of these findings well beyond the saturation point of SI. We continue the numerical simulations presented here to study the long-term evolution of initially unstratified fronts. In particular, we focus on the equilibration of the front and how the details depend on the particular flavour of SI and the front strength.