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.


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 ∇ h ) of the background density field,ρ (measured by the horizontal analogue to the buoyancy frequency, M 2 ≡ g/ρ 0 |∇ hρ |, with g the acceleration due to gravity, and ρ 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/ρ 0 ∇ hρ ×ẑ (whereẑ is the vertical unit vector). These baroclinic torques tend to flatten isopycnals, but may be counterbalanced by a Coriolis torque, f ∂ zūg (where f is the Coriolis parameter), arising from a vertical shear in the geostrophic velocity,ū 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) (defined with the velocity, u, and buoyancy, b ≡ −gρ/ρ 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 where ω z is the vertical component of the relative vorticity and M 2 ≡ ∂ xb (as above) is the horizontal analogue to the buoyancy frequency, N 2 ≡ ∂ zb . 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 ω z + f > 0 whereas 'inertial instability' occurs when N 2 > 0 and ω z + f < 0. Therefore SI only occurs when (ω 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 ≡ N 2 f 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  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) andlater 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 2009Thomas & 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, Γ ≡ 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, Γ .
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.

Linear stability analysis
Perhaps the simplest model of a front, the Eady model was first introduced by Eady (1949) and later used by Stone (1966Stone ( , 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.
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 Here, ν is the kinematic viscosity and κ 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, Γ = M 2 /f 2 ≈ [1, 100], which covers a wide variety of ocean fronts. Although very strong fronts with Γ > 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 (Γ > 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 tof = 2Ω cos φ, where Ω is the angular velocity and φ 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 Consistent with the non-dimensional parameters (2.1a-d) introduced above, the dimensionless ( * ) variables here are The dimensionless pressure head acceleration, ∇ * Π * , absorbs the hydrostatic pressure gradient, and is eliminated when writing (2.2a) with the along-front streamfunction. We choosex to be the across-front direction (parallel to ∇ hb ). The background basic state (denoted by an overbar) used for linearisation and the initial condition for the numerical simulations isv 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).

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 (ŷ) direction (Stone 1966), we consider linear perturbations that vary only in x and z: 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 ω) of the form where the eigenfunction,χ(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 ) = ∇ × ψŷ, becomes where D ≡ d/dz 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ψ at z = 0, 1 on this sixth-order ODE areψ 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, 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 ⎡ which has eigensolutions of the form that match the boundaries if λ 1 − λ 2 = 2π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 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, σ , is the imaginary part of ω. 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 2π. 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).
We can now consider how the fastest growing mode of (2.10) varies with Γ and Ri, as shown in figure 2(b). For Ri = 0 the energy growth rate relative to f increases nearly 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 (θ m , dot-dashed line) and isopycnals (θ b , dotted lines), which for N 2 /M 2 = 0, θ b = 90 • . This unstratified case has modes nearly equally spaced between the isopycnals and absolute momentum surfaces for large Γ , but with increasingly horizontal isopycnals the SI modes grow more along these isopycnals. While the angle of the contour ψ(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 θ 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, σ . 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 τ c /2.
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 (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 • 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 (m =v g + Γ −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 (θ m = tan −1 Γ −1 ) for small front strength. This permits a larger buoyancy production of energy (B = w b ), as shown in figure 3(b), while the geostrophic shear production (P g = − v w ∂ zvg ) is the dominant energy source in the rest of the parameter space. Here and throughout the rest of this paper, · indicates a volume average over the entire domain, and primes represent local departures from the horizontally averaged fields denoted by·.

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, σ SI , is equal to the KHI growth rate, σ KH . Of course, σ KH is a monotonically increasing function of the shear, and thereby of U SI which exponentially grows at a rate σ 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 σ KH,max (U SI ) with a pseudo-spectral solver until finding the critical SI mode amplitude, U c . While the most unstable SI wavevector, |k SI |, increases as the mode number (n) and Re increase, the scaling for U c appears to be dominated by σ SI and so remains largely unchanged.
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, (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, σ KH ∝ U c |k SI |) with the SI growth rate in the limit of large Γ , σ SI ∝ Γ −1/2 . We see that this simple scaling argument fails for small Γ where the growth rate in these weak fronts is slow compared with f .

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 Γ = {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, τ 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 τ 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.

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, τ c = σ −1 SI 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,ψ, found by iteratively solving for λ 1 and λ 2 in (2.11), we determined the full structure of these modes:û,v,ŵ andb 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 |û| 2 + |ŵ| 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), (4.1) The first two terms on the right-hand side represent the shear production, P, converting energy from the mean flow into TKE. Specifically, the along-front contribution (P y ) is split into a geostrophic shear production term, P g , energised by the thermal wind shear, and an ageostrophic part. The other potential source of TKE comes from buoyancy production, 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 τ c is shown in figure 5(b). As expected for SI, the contribution from P g exceeds B, except for small Γ . 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, 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 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, τ c . We do this by introducing the energy production ratio, as plotted in figure 7, where B is the buoyancy flux, P g is the geostrophic shear production and λ 1 and λ 2 describe the vertical mode structure (2.11) and depend on Ri and Γ (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 θ b . Non-hydrostatic effects flatten the SI modes particularly for small Γ . This permits 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 Γ ) 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.
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 Γ 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). In the 'slantwise convection' regime, where B > 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 τ 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 C e = Γ −1 w b MLI = 0.06 − 0.08. (4.5) In comparison, a typical SI buoyancy flux at τ c for the slantwise convective regime (specifically at Γ = 1 and Ri = 0) is Note that the buoyancy flux increases in time during the growing phase of SI and MLI. The fact that w b SI at t = τ c is smaller than w b 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 (Γ = 1). For stronger fronts with weak vertical stratification (i.e. large Γ 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.

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.

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 thex andŷ directions. Using continuity and geostrophic balance, the horizontal ageostrophic momentum equations are 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), where we have also assumed exponential growth in time, ∝ exp(σ 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, τ , then cancels with the exponential evaluated at τ in the numerator. Similarly for the terms in the y-momentum equation (5.1b), the ratio is 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 1. We found this to be the case for Γ 1 and Ri 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 ∂ tva ≈ −∂ z v w .

Loss of geostrophic balance
This dominant balance with the ∂ z 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, we can estimate the time scale required to mix the thermal wind shear for SI momentum fluxes evaluated at τ c . This value is plotted for each Γ 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 Γ 10. Tandon & Garrett (1994) showed that in the limit of instantaneous mixing (here for Γ 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 ∂ zv | t=0 = s∂ zvg . (5.6) 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 These inertial oscillations draw closed circular orbits and have a linear structure in z The amplitudes of these inertial shear oscillations are dimensionally (1 − s)M 2 /f .

Inertial oscillation amplitude
The reduction in thermal wind shear before τ 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 τ c 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 τ 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.

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 Γ = 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 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 (∂ xv = 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 Γ 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.
in the inviscid limit by Colin de Verdière (2012) and for unbounded modes by Zeitlin (2018).
One consequence of the traditional approximation we used in the analysis thus far is that the dynamics is independent of the front orientation. We therefore only specified that x points across the front (parallel to ∇ hb ). However, the non-traditional terms break this horizontal isotropy, and so we must specify the angle, ϑ, of the background buoyancy gradient relative to north. We still take x to be across front (i.e. |∇ hb | = ∂ xb = M 2 ), but we now orient the entire front (andx) at an angle ϑ from north.
Including the northward horizontal component of Earth's rotation, the Boussinesq momentum equation (2.2a) becomes The importance of these non-traditional terms is measured by which accounts for both the latitude (φ) and the orientation (ϑ) of the across front (x-axis) relative to north. α is the 'symmetric' component off in the along-front (ŷ) direction, and drops out upon writing (A1) with the streamfunction. So while the front orientation and latitude are both important when considering non-traditional effects, these can be reduced into the single parameter γ . It becomes apparent now that the traditional approximation (γ → 0) used to simplify the analysis in § 2 holds better at mid to high latitudes and for fronts with a nearly east/west lateral density gradient. Additionally, the importance of this horizontal component is diminished in the large shear regime of the strong fronts we considered, where the vorticity from the thermal wind shear (M 2 /f ) greatly exceedsf (i.e. when γ /Γ 1).
A.2. Non-traditional results at φ = 45 • To demonstrate the effects of the non-traditional terms on the main results in this paper, we present a selection of these results for φ = 45 • , and for ϑ = 0 • and 180 • . We find that while the horizontal component of Earth's rotation quantitatively influences the SI growth and transport properties, it does not qualitatively change the observed trends and our conclusions.
The non-traditional terms impact the stability of SI by changing the contours of absolute momentum,m (Recall x is still the across-front coordinate, but now the entire front has been oriented ϑ from north.) This means that for the range Γ < γ , the front is stable to SI. This is written equivalently as a sub-critical Richardson number, Ri c = 1 − γ /Γ . Of course it should be emphasised that at φ = 45 • , γ = 1 only if the high buoyancy side is further north (ϑ = 0 • ). In the opposite orientation (when the buoyancy gradient points south) then γ = −1. Thus non-traditional effects can either increase or decrease the region of instability (in Ri-Γ space) and consequently influence the growth rate. This is apparent in figure 9(a), where for strong yet unstratified fronts, the non-traditional effects have a uniform influence of increasing (decreasing) the growth rate by ∼ 25 % when the  (1), as shown in figure 7(b) in the traditional approximation. The black solid line is the sub-critical Richardson number, Ri c = 1 − γ /Γ , which is no longer equal to 1. The white line separates regions of parameter space where SI modes are aligned closer to isopycnals (inside) from regions (outside) where they are more along absolute momentum surfaces. Comparing these two contour plots with figure 7(b) shows similarly distinct regions that could be characterised as 'slantwise convection' separated from the 'slantwise inertial instability'.
buoyancy gradient is north (south). This effect is much less pronounced with even a weak stratification of Ri = 0.25, in agreement with Colin de Verdière (2012). By changing the contours of absolute momentum (A3), the non-traditional Coriolis terms also influence the angle of the SI modes. As shown in figure 9(b), the SI mode angle becomes steeper with increasing γ , and tends to align more with isopycnals as the tilted rotation vector steepens the absolute momentum contours.
926 A6-19 We finally consider how the energy source for SI changes with varying γ . This is best seen by the generalised production ratio (B13) which is plotted for γ = ±1 in figure 10. Compared with figure 7(b) under the traditional approximation, we note similarly distinct regions of slantwise convection (black) and slantwise inertial instability (red). These regions are largely unchanged for large Γ (where the strong thermal wind shear means that all rotation is less important), and the slantwise convective character still persists near Ri = 1 compared with the slantwise inertial instability region for small Ri. Still, the energy source for weaker fronts appears to be influenced more by the non-traditional effects. optimisation problem to find the minimum U SI that satisfies σ KH (k KH ; U SI ) = σ SI , and plot this U c (Γ ) as a solid line in figure 5(a). This value for U c can then be used to calculate the various transport and energetic quantities for SI in Appendix B.