# Vertical Mixing Can Both Induce and Inhibit Submesoscale Frontogenesis

Daniel Dauhajre, Kaushik Srinivasan, M. Jeroen Molemaker, Jonathan Gula, Delphine Hypolite, James Mcwilliams, Roy Barkan, William Young

► To cite this version:

Daniel Dauhajre, Kaushik Srinivasan, M. Jeroen Molemaker, Jonathan Gula, Delphine Hypolite, et al. Vertical Mixing Can Both Induce and Inhibit Submesoscale Frontogenesis. Journal of Physical Oceanography, 2025, pp.1-57. 10.1175/JPO-D-24-0148.1. hal-05196319

HAL Id: hal-05196319

https://hal.univ-brest.fr/hal-05196319v1

Submitted on 1 Aug 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
1

Vertical Mixing Can Both Induce and Inhibit Submesoscale Frontogenesis

2

Daniel P. Dauhajre,ª,ᵇ Kaushik Srinivasan,ª M. Jeroen Molemaker,ª

3

Jonathan Gula,ᶜ,>d Delphine Hypolite,ª James C. McWilliams,ª Roy Barkan, a,e

4

William R. Young ᵈ

5

ª *Department of Atmospheric and Oceanic Sciences, University of California Los Angeles*

6

ᵇ *Department of Geography, University of California, Santa Barbara*

7

ᶜ *Univ. Brest, CNRS, Ifremer, IRD, INRIA, Laboratoire d'Océanographie Physique et Spatiale (LOPS), IUEM, F29280, Plouzané, France*

8

ᵈ *Institut Universitaire de France (IUF), Paris, France*

9

ᵉ *Department of Geosciences, Tel Aviv University*

10

f Scripps Institution of Oceanography, University of California, San Diego

11

12
ABSTRACT: Past studies separately demonstrate that vertical boundary layer turbulence can either sharpen or weaken submesoscale fronts in the surface mixed layer. These studies invoke competing interpretations that separately focus on the impact of either vertical momentum mixing or vertical buoyancy mixing, where the former can favor sharpening (frontogenesis) by generation of an ageostrophic secondary circulation, while the latter can weaken the front (frontolysis) via diffusion or shear dispersion. No study comprehensively demonstrates vertical mixing induced frontogenesis and frontolysis in a common framework. Here, we develop a unified paradigm for this problem with idealized simulations that explore how a front initially in geostrophic balance responds to a fixed vertical mixing profile. We evolve 2D fronts with the hydrostatic, primitive equations over a range of Ekman ($Ek = 10^{-4} - 10^{-1}$) and Rossby numbers ($Ro = 0.25 - 2$), where $Ek$ quantifies the magnitude of vertical mixing and $Ro$ quantifies the initial frontal strength. We observe vertical momentum mixing induced, nonlinear frontogenesis at large $Ro$ and small $Ek$ and inhibition of frontogenesis via vertical buoyancy diffusion at small $Ro$ and large $Ek$. Symmetric instability can dominate frontogenesis at very small $Ek$; however, the fixed mixing limits interpretation of this regime. Simulations that suppress vertical buoyancy mixing are remarkably frontogenetic, even at large $Ek$, explicitly demonstrating that buoyancy mixing is frontolytic. Application of two scalings to quantify the competition between cross-front buoyancy advection and vertical diffusion identifies practically equivalent controlling parameters ($Ro^2/Ek, Ro/Ek^{1/2}$); these ratios approximately map regime transitions across simulations with equal vertical eddy viscosity and diffusivity.
SIGNIFICANCE STATEMENT: This study reconciles competing views on how turbulent vertical mixing on 0.01 - 1 m scales controls the sharpening or weakening of upper-ocean fronts characterized by horizontal changes in density and velocity over 100 m - 1 km scales. This sharpening or weakening modulates frontal circulation that acts to bring heat upwards; given the pervasiveness of upper-ocean fronts, these local dynamics influence global upper-ocean heat content. Utilizing simulations, we identify a measurable parameter that predicts frontal sharpening or weakening via vertical mixing. This new dynamical framework can better inform the necessary parameterization of these fronts in global climate models. However, future work should interrogate the validity of our simplified model, which unrealistically assumes that the vertical mixing does not evolve.

# 1. Introduction

Submesoscale turbulence (Thomas et al. 2008; McWilliams 2016; Gula et al. 2022; Taylor and Thompson 2023) spawns flow patterns with horizontal scales $O(0.1 – 1 \text{ km})$ and vertical scales $O(1 – 100 \text{ m})$ that frequently populate the surface mixed layer. These submesoscale flow patterns encompass mixed layer eddies (Boccaletti et al. 2007) generated by a form of baroclinic instability and density fronts and filaments (McWilliams 2016). This idealized modeling study focuses on submesoscale fronts, which characteristically exhibit strong, dynamically consequential ageostrophic overturning circulations. These ageostrophic overturning circulations can be triggered and fueled by straining currents – supplied for example, by mixed layer or mesoscale eddies (Zhang et al. 2021) – and/or vertical boundary layer turbulence (Gula et al. 2014; McWilliams et al. 2015; Dauhajre and McWilliams 2018; Barkan et al. 2019). Once activated, an overturning circulation can rapidly sharpen a submesoscale front via amplification of horizontal density and velocity gradients in a process known as frontogenesis (Hoskins 1982; McWilliams 2021). Submesoscale fronts undergoing frontogenesis re-stratify the mixed layer (Taylor and Thompson 2023) while simultaneously serving as conduits for a forward energy cascade (Srinivasan et al. 2023; Contreras et al. 2023). The conditions required for submesoscale frontogenesis – available potential energy in a surface mixed layer with vertical boundary layer turbulence and/or ambient straining or deformation currents – typically arise in both the open-ocean and continental shelves, making this seemingly spontaneous process pervasive and significant.
Like many other types of oceanic fronts with ageostrophic motions (e.g., gravity fronts), submesoscale fronts exhibit a Rossby number ($Ro = V/f_l$) that is $O(1-10)$, often quantified by a horizontal velocity gradient ($V/l$) normalized by the Coriolis frequency ($f$). However, unlike other types of fronts, submesoscale fronts (and filaments) uniquely exhibit geostrophic and ageostrophic velocities of comparable of magnitude (Barkan et al. 2019). The horizontal gradients of these velocities result in distinct surface signatures of large cyclonic relative vorticity ($\zeta/f \gg 1$; where $\zeta = v_x - u_y$) and convergence ($\delta/f \ll -1$; where $\delta = u_x + v_y$). The cyclonic vorticity is related to an along-front jet ($\zeta \approx v_x$, where $x$ is the across-front direction) and the convergence indicative of a preferentially downwelling, ageostrophic cross-frontal circulation ($\delta \approx u_x$), commonly referred to as an ageostrophic secondary circulation (ASC). ASCs exhibit extreme vertical velocity ($w \lesssim -100 \text{ m/day}$), relative to larger scale currents (Farrar et al. 2020; Taylor and Thompson 2023) and can regulate a variety of oceanic processes, including: energetic exchanges (Molemaker et al. 2010; Su et al. 2018; Srinivasan et al. 2023; Contreras et al. 2023), smaller-scale turbulence (Buckingham et al. 2019; Wenegrat et al. 2020; Peng et al. 2021; Chor et al. 2022), larger-scale stratification (Fox-Kemper et al. 2011; Su et al. 2018), air-sea fluxes (Siegelman et al. 2020), biogeochemical cycling (Taylor 2016; Freilich et al. 2022; Damien et al. 2023), ecosystem functioning (Levy et al. 2012; Lévy et al. 2018; Fahlbusch et al. 2024), and pollution dispersal (D'Asaro et al. 2018).

Central to understanding the manner in which submesoscale fronts modulate these processes – and fundamental to designing parameterization of submesoscale material fluxes (Young 1994; Fox-Kemper et al. 2008; Zhang et al. 2023; Bodner et al. 2023; Yang et al. 2024) – are the dynamical frameworks that explain how the ASC, and thus the front, strengthens or weakens over a frontal life-cycle, which typically spans hours-to-days. This life-cycle generically encompasses the triggering of a surface convergent ASC that initiates frontogenesis, followed by the subsequent erosion of horizontal density and velocity gradients (frontolysis) by some arresting mechanism (e.g., instability); other currents or atmospheric forcing can externally influence this sequence. Submesoscale frontogenesis spans hours (Barkan et al. 2019), while fronts maintain, undergo instabilities and dissipate on daily time-scales.

This study re-litigates the role of vertical boundary layer turbulence in submesoscale frontogenesis and frontolysis, motivated by competing interpretations (overviewed in Sec. 1a) of whether
and how vertical mixing sharpens or weakens submesoscale fronts. Here, vertical boundary layer turbulence refers generally to motions smaller than the submesoscale that drive vertical mixing, which is often represented as a vertical eddy viscosity ($v_v$) and diffusivity ($\kappa_v$), and considered to be driven by atmospheric forcing or instabilities that can emerge locally at submesoscale fronts (Thomas et al. 2013; Verma et al. 2019; Yu et al. 2019; Carpenter et al. 2020; Peng et al. 2021; Chor et al. 2022). The lack of consensus on this problem stems from dynamical frameworks that separately invoke either the impact of $v_v$ (in the momentum equation) or $\kappa_v$ (in the buoyancy equation), where $v_v$ can setup a convergent ASC that favors frontogenesis (and re-stratifies the mixed layer), while $\kappa_v$ acts to diffuse the front (and maintains the mixed layer).

While observations shed some light on the coupling between vertical mixing and submesoscale fronts (Nagai et al. 2006; Johnston et al. 2011; Johnson et al. 2020a,b; Carpenter et al. 2020; Swart et al. 2020; Peng et al. 2021), the difficulty in simultaneously measuring spontaneously arising submesoscale fronts as well as smaller-scale turbulence over a range of conditions has left the bulk of mechanistic interpretation to theoretical and numerical treatments. These studies comprise analyses of submesoscale fronts and filaments in realistically configured, primitive equation simulations (Gula et al. 2014; Dauhajre et al. 2017; Wang et al. 2021; Barkan et al. 2019; Srinivasan et al. 2023) as well as more idealized or theoretical approaches that span for example, complex (Sullivan and McWilliams 2017; Verma et al. 2019; Sullivan and McWilliams 2024) and simplified (Crowe and Taylor 2019) large-eddy simulations; 2D semi-geostrophic (Thompson 2000) or primitive equation (McWilliams et al. 2015) models; asymptotic expansions (Young 1994; Crowe and Taylor 2018); and perturbation analysis (Bodner et al. 2019).

A general approach that isolates the role of vertical mixing evolves an initial front (or filament) that is forced only by vertical mixing, which can be prescribed (Thompson 2000; Crowe and Taylor 2018, 2019; McWilliams 2017), parameterized (McWilliams et al. 2015) or partially resolved (Sullivan and McWilliams 2017; Verma et al. 2019; Sullivan and McWilliams 2024). The two primary controlling parameters inherent in this posing are the initial frontal strength, which can be quantified with a Rossby number ($Ro = \zeta/f, \delta/f$) and the vertical mixing intensity, which can be quantified with an Ekman number ($Ek = v_v/fh_{ml}^2$, where $h_{ml}$ is a mixed or turbulent boundary layer depth).
Past studies – many of which employ the above-described approach – separately demonstrate that vertical mixing can either sharpen (Thompson 2000; McWilliams et al. 2015; McWilliams 2017; Sullivan and McWilliams 2017, 2024) or weaken (Young 1994; Crowe and Taylor 2018, 2019; Bodner et al. 2019) fronts. However, these studies sample separate regions of the relevant parameter space (Ek, Ro) (see Table A1), with no individual study demonstrating both vertical mixing induced frontogenesis and frontolysis (for Ro > 1) in a common framework. This has led to seemingly disconnected, competing interpretations for the impact of vertical mixing on submesoscale fronts. Below, we heuristically describe these competing views to motivate our attempt to develop a common paradigm for this problem.

### a. Competing interpretations on the role of vertical mixing

It is important to first note that, regardless of the role of vertical mixing, frontogenesis is a non-linear process involving cross-front advection of buoyancy and/or momentum. Barkan et al. (2019) demonstrate—with a submesoscale-resolving realistic model and asymptotic theory—that once the ASC Rossby number ($Ro = \delta/f$) reaches $O(1)$, the convergence becomes the primary determinant of frontal sharpening; Yu et al. (2024) corroborate the existence of this nonlinear frontogenetic regime in observations. That is, the question of how vertical mixing impacts frontogenesis primarily concerns whether or not vertical mixing *induces* or *inhibits* the transition to this nonlinear frontogenetic stage.

The prevailing paradigm predicts that vertical mixing induces submesoscale frontogenesis. This view stems from considering the impact of vertical momentum mixing (via $v_v$) on thermal wind balance, which results in a linear, three-way balance between rotation, pressure gradient, and vertical momentum mixing. This diagnostic balance dates back to Heaps (1972) and in recent literature is referred to as “turbulent thermal wind (TTW) balance” (Gula et al. 2014; McWilliams et al. 2015; Bachman and Taylor 2016; Crowe and Taylor 2018; Lentz 2022), with Garrett and Loder (1981) providing original theoretical treatment for $Ro < 1$. TTW balance is also contained within the (more generalized) subinertial mixed layer model of Young (1994); it arises as the dominant balance there when the vertical mixing is large.

TTW often successfully predicts a surface convergent, frontogenetically favorable ASC for submesoscale fronts and filaments with characteristic $Ro > 1$ in realistic settings (Gula et al.
2014; Dauhajre et al. 2017; Wang et al. 2021; Barkan et al. 2019). A commonly invoked TTW scaling – which assumes that only the geostrophic velocity is vertically mixed (Garrett and Loder 1981; Thomas and Ferrari 2008; McWilliams 2017) – anticipates stronger convergence for larger $v_v$, given the same horizontal density gradient and $f$ (McWilliams 2017). Thus, TTW balance provides a route for vertical momentum mixing to setup a convergent, high $Ro$ ASC that can initiate a transition to nonlinear frontogenesis (Barkan et al. 2019), with a prevailing expectation that larger $v_v$ (and $Ek$) induces a stronger ASC (larger $Ro \approx \delta/f$). Interestingly, Barkan et al. (2019) show that the TTW balance describes the ‘early-time’ (less than one inertial period) convergence of the submesoscale fronts in a realistic submesoscale resolving ocean model solution, which naturally contains straining induced by background mesoscale and mixed-layer eddies.

The framework of ‘TTW frontogenesis’ often assumes (at least heuristically) that the effect of vertical buoyancy mixing is negligible due to weak stratification in the mixed layer (i.e., $\kappa_v \partial b / \partial z \sim 0$, where b is the buoyancy). This heuristic leads to the (unrealistic) prediction that an ‘approximately’ balanced TTW secondary circulation can sharpen a front until a singularity is reached (see Sec. 6 of McWilliams et al. (2015)). Of course, the weak stratification in the mixed layer is fundamentally due to strong vertical buoyancy mixing. Considering the kinematics of vertical buoyancy mixing acting on a front leads to a perhaps more intuitive, although less invoked, expectation that vertical mixing will weaken a front. This frontolytic view (Young 1994; Crowe and Taylor 2018, 2019) implicates $\kappa_v$ in weakening fronts via vertical diffusion (at large $\kappa_v$, $Ek$) or shear dispersion (at intermediate $\kappa_v$, $Ek$) (Crowe and Taylor 2018), where oscillations of the vertically sheared ASC coupled to vertical diffusion lead to an effective horizontal diffusivity that spreads isopycnals apart (Young and Jones 1991; Young 1994; Wenegrat et al. 2020; Swart et al. 2020). The asymptotic theory (Crowe and Taylor 2018) underpinning this view can be considered a specialized demonstration of dynamics in the subinertial mixed-layer model of Young (1994). The asymptotic analysis in Crowe and Taylor (2018) is limited to $Ro < 1$, not typical of submesoscale fronts; although, a numerical test of the theory (Crowe and Taylor 2019) suggests validity at $Ro = 1$. The interpretation is that the momentum balance satisfies a ‘quasi-steady’ (linear) TTW balance, with the buoyancy evolution (vertical diffusion or shear dispersion) dominating frontal evolution. Importantly, this frontolytic interpretation results from focusing primarily on the long-term ($\gtrsim 10$ inertial periods) solution behavior.
No past studies (summarized in Table A1) account for all of the above-described mechanisms in a common framework, with these heuristics made more complex by the fact that fronts actually modulate the boundary layer turbulence (Verma et al. 2019; Carpenter et al. 2020; Peng et al. 2021; Sullivan and McWilliams 2024; Johnson and Fox-Kemper 2024). An additional factor contributing to the present confusion is that separate interpretations generally focus on different time-scales: (TTW) frontogenesis that lasts ≲ 1 inertial period versus vertical diffusion or shear dispersion frontolysis that emerges after ≳ 1 – 10 inertial periods.

## b. This study

Here, we attempt to reconcile competing views on this problem with a comprehensive, idealized exploration of vertical mixing ($v_v$, $\kappa_v$) impacts on submesoscale frontogenesis. By design, we do not consider the impact of straining, which is well-understood to induce frontogenesis (Hoskins and Bretherton 1972; Shakespeare and Taylor 2013); we note that Bodner et al. (2019) find that vertical mixing generally inhibits strain-induced frontogenesis with a perturbation analysis.

Sec. 2 describes the experimental setup, which poses a simple question: how does a surface layer front initially in geostrophic balance evolve in response to the introduction of a prescribed vertical mixing profile ($v_v(z)$, $\kappa_v(z)$)? We primarily focus on time-scales of ≈ 0.1 – 1 inertial periods, motivated by the *O*(hour) frontogenesis observed in more realistic scenarios (e.g., Barkan et al. (2019)). Scaling analysis (Sec. 2b) identifies three controlling non-dimensional parameters: a Rossby number $Ro$; an Ekman number $Ek = v_0/fh_{ml}^2$, where $h_{ml}$ is a mixed-layer depth, and $v_0$ is a mixed-layer average vertical eddy viscosity; and a turbulent Prandtl number $Pr = v_v/\kappa_v$. We demonstrate that two separately formulated scaling choices for the cross-front velocity–originating in previous studies of submesoscale frontogenesis (Thomas and Ferrari 2008; Barkan et al. 2019)– result in a consistent quantification of the controlling buoyancy dynamics (Sec. 2b). This scaling guides the parameter variations in the numerical experiments (Sec. 2c, Table 1).

Our idealized setup prescribes an initial front that is motivated by realism (analogous to Sullivan and McWilliams (2017, 2024); Verma et al. (2019); McWilliams et al. (2015)), but explicitly isolates the role of vertical mixing by artificially holding $v_v$ and $\kappa_v$ fixed in time (as in Crowe and Taylor (2019)). The primary numerical experiment evolves 2D fronts over a range of initial frontal strengths ($Ro = 0.25 – 2$) and vertical mixing intensities ($Ek = 10^{-4} – 10^{-1}$) with $Pr = 1$. This is
a broader parameter space than previous individual studies that targets *Ro* and *Ek* values typical of extratropical submesoscale fronts (**references??**); we intentionally exclude geostrophic adjustment (no initial flow; *Ro* = 0) and gravity currents (*f* ≈ 0, *Ro* → ∞).

We demonstrate that vertical mixing can both induce and inhibit frontogenesis, with all solutions eventually exhibiting frontolysis, and we map these regimes in the (*Ek*,*Ro*) space (Sec. 3). The scaling arguments in Sec. 2b identify two separately derived, yet practically equivalent, non-dimensional numbers (*Ro*²/*Ek* and *Ro*/*Ek*¹/²) that quantify competition between (frontogenetic) horizontal buoyancy advection and (frontolytic) vertical diffusion. We demonstrate that these parameters well predict the transition from frontogenetic inhibition to frontogenesis (for *Pr* = 1). This mapping of regimes to a single parameter leverages understanding of the distinct roles of *ν*_v and *κ*_v in frontogenesis or frontolysis, which is elucidated with simulations that suppress the vertical buoyancy mixing (*κ*_v = 0, *Pr* = ∞). Sec. 4 details the controlling dynamical balances for frontogenetic and frontogenetically inhibited regimes as well as describing the mechanisms controlling late-time frontolysis. Sec. 5 discusses caveats of the idealization and contextualizes our interpretations relative to previous studies and applicability in more realistic environments. Sec. 6 summarizes the results and interpretations of this study.

## 2. Idealized setup

The basic idealized experiment triggers the evolution of a surface layer density front initially in geostrophic balance with the introduction of a fixed vertical mixing profile (*ν*_v(z), *κ*_v(z)). We evolve the fronts for a range of vertical mixing intensities (*Ek*) and initial frontal strengths (*Ro*) with the hydrostatic, primitive equations in a 2D configuration utilizing the Regional Oceanic Modeling System (ROMS; Shchepetkin and McWilliams (2005)). Here, we define the 2D system (Sec. 2a), identify controlling non-dimensional parameters (Sec. 2b), and detail the ROMS idealized setup and experimental design (Sec. 2c).

### a. 2D system

We take the hydrostatic, primitive equations with *x*, *y* as the across and along-front directions, respectively. For simplicity, we assume along-front uniformity ∂/∂*y* = 0 and a vertically variable
mixing profile ($v_v(z), \kappa_v(z)$). This gives the 2D $(x, z)$ system:

$$D_t u - f v = -\phi_x + \frac{\partial}{\partial z} (v_v u_z), \qquad (1a)$$

237

$$D_t v + f u = \frac{\partial}{\partial z} (v_v v_z), \qquad (1b)$$

238

$$D_t b = \frac{\partial}{\partial z} (\kappa_v b_z), \qquad (1c)$$

239

$$\phi_z = b, \qquad (1d)$$

240

$$u_x + w_z = 0, \qquad (1e)$$

where $D_t = \partial_t + u\partial_x + w\partial_z$ is the material derivative; $u$ is the across-front velocity; $v$ is the along-front velocity; $w$ is the vertical velocity; $\phi = p/\rho_0$ is the pressure normalized by a reference density $\rho_0$; and $b = -g\rho/\rho_0$ is the buoyancy.

The horizontal boundary conditions are periodic in $x$ and the vertical boundary conditions are zero buoyancy flux and zero stress at the free-surface ($z = \eta$) and a zero buoyancy flux and a bottom stress ($\vec{\tau}^b$) at the (flat) bottom ($z = -H$).

244

245

246

$$v_v \frac{\partial \vec{u}_h}{\partial z} = 0, \quad \text{at } z = \eta \qquad (2a)$$

247

$$v_v \frac{\partial \vec{u}_h}{\partial z} = \frac{\vec{\tau}^b}{\rho_0}, \quad \text{at } z = -H \qquad (2b)$$

248

$$\kappa_v \frac{\partial b}{\partial z} = 0, \quad \text{at } z = \eta, -H \qquad (2c)$$

249

250

251

The bottom boundary conditions have little significance on the near-surface behavior that is the focus of this study due to the deep bottom and insignificant near-bottom velocities (and bottom stress).

252

253

254

TTW balance (Sec. 1a) is specified by removing the $D_t$ terms in Eq. 1a-1b. Including acceleration ($\partial u/\partial t, \partial v/\partial t$) in TTW gives the ‘transient’ TTW (or ‘$T^3$W’) balance (Dauhajre and McWilliams 2018; Wenegrat and McPhaden 2016), which is discussed in Sec. 4.
b. *Scalings for cross-front velocity and controlling buoyancy dynamics*

Here, we present two separate scalings for the cross-front velocity (u) which fundamentally drives frontogenesis (amplification of $b_x, u_x, v_x$ in Eq. 1) via cross-front advection of velocity and buoyancy. These scalings originate in previous studies of submesoscale frontogenesis (Thomas and Ferrari 2008; Barkan et al. 2019), with distinct arguments guiding their formulations, as detailed below. The purpose of this section is to demonstrate that both scalings lead to a consistent measure of the competition between cross-front buoyancy advection (which sharpens a front) and vertical buoyancy diffusion (which weakens a front; Sec. 1a). This demonstration motivates our experimental design (Sec. 2c) and anticipates a central result of the study: buoyancy dynamics, as opposed to momentum dynamics, fundamentally control whether frontogenesis or frontolysis occurs in our 2D simulations (assuming $v_v = \kappa_v$).

Both scalings choose an along-front velocity scaling $v \sim V$. Thomas and Ferrari (2008) employ a dynamical argument that the cross-front velocity scales with the Ekman transport associated with the geostrophic stress ($\tau_g = \rho v_0 \partial v_g / \partial z$; where $v_g$ is a geostrophic along-front velocity). They define an Ekman layer thickness $\delta_e = \sqrt{2v_0/f}$ and choose $v_v \sim v_0$ and $\partial v_g / \partial z \sim V/h_{ml}$, which gives the cross-front velocity scaling:

$$u_{\text{TF08}} \sim \frac{\tau_g}{\rho_0 f \delta_e} \sim \frac{v_0 |\partial v_g / \partial z|}{f \delta_e} \sim \frac{\delta_e V}{2h_{ml}} \sim Ek^{1/2}V, \quad (3)$$

where $Ek = v_0/h_{ml}^2 f$ is the Ekman number. This scaling essentially captures the TTW dynamics that set up or sustain an ASC (Sec. 1a).

Barkan et al. (2019) instead make an empirically based assumption of comparable along- and across-front velocities at submesoscale fronts with $Ro \sim O(1)$. This assumption is based on observed velocities in a realistic submesoscale resolving simulation (see Fig. 2 in Barkan et al. (2019)) and leads to the scaling:

$$u_{\text{B19}} \sim RoV, \quad (4)$$

where $Ro = V/fl$ is the Rossby number with cross-front length scale $l$. This scaling attempts to capture the observed nonlinearity during frontogenesis (i.e., after an ASC has formed). We note that any scaling of the form $u \sim Ro^nV$ with $n > 0$ would give comparable along- and across
front velocities at $Ro \sim O(1)$, with $n \approx 1/2 - 1$ the most physically plausible for submesoscale frontogenesis (Barkan et al. 2019; Yu et al. 2024).

We use $u_{B19}$ and $u_{TF08}$ to separately scale the competition in Eq. 2c between the (frontogenetic) cross-front buoyancy advection via the ASC ($ub_x$) and (frontolytic) vertical buoyancy diffusion ($\frac{\partial}{\partial z}[\kappa_v b_z]$). Thompson (2000) and Thomas and Ferrari (2008) also highlight these competing buoyancy dynamics in frontogenesis, albeit in a semi-geostrophic framework (Thompson 2000) and in the context of re-stratification at fronts forced by wind and straining (Thomas and Ferrari 2008).

We choose an across-front length scale $x \sim l$, vertical length scale $z \sim h_{ml}$, and vertical eddy viscosity $\nu_v \sim \nu_0$ and diffusivity $\kappa_v \sim \kappa_0$. Applying $u_{TF08}$ gives:

$$ \frac{u_{TF08} b_x}{\kappa_0 b_{zz}} \sim Pr \frac{Ro}{Ek^{1/2}}, \qquad (5) $$

and applying $u_{B19}$ gives:

$$ \frac{u_{B19} b_x}{\kappa_0 b_{zz}} \sim Pr \frac{Ro^2}{Ek}, \qquad (6) $$

where $Pr = \nu_0/\kappa_0$ is the turbulent Prandlt number.

Assuming $Pr = 1$ gives two non dimensional ratios that equivalently quantify the competition between cross-front buoyancy advection and vertical diffusion: $Ro^2/Ek$ (Eq. 5) and $Ro/Ek^{1/2}$ (Eq. 6). Given that $Ro/Ek^{1/2}$ is the square root of $Ro^2/Ek$, these parameters can be used interchangeably; we somewhat arbitrarily choose $Ro^2/Ek$ for the purposes of notation in the rest of the study. We will demonstrate that the competition measured by $Ro^2/Ek$ well-predicts regime transitions (frontolysis to frontogenesis) for the 2D simulations in this study when $Pr = 1$.

It is worth commenting on the seemingly coincidental equivalence that results from these distinct scaling arguments. First, we note that these scaling choices only result in a consistent scaling of the buoyancy equation; application of $u_{B19}$ and $u_{TF08}$ results in differently scaled momentum equations. A heuristic interpretation is that $u_{TF08}$ captures non-conservative (i.e., TTW) dynamics that can create or sustain a frontal ASC, while $u_{B19}$ assumes a front (and ASC) have already formed, and captures the ASC during advective frontogenesis (Barkan et al. 2019). The two different scaled momentum equations delineate this distinction (see Appendix B); a single, generic scaling would need to account for transitions between the $u_{TF08}$ and $u_{B19}$ regimes during a frontal life-cycle.
TABLE 1. Parameters for sets of idealized simulations. The controlling non-dimensional parameters (Sec. 2b) for each case comprise the Rossby number, based on the geostrophic initial condition ($Ro = \zeta_{init}/f$, where $\zeta = v_x$ at $t = 0$; see Fig. 1a-d); the Ekman number ($Ek = v_v/(h_{ml}^2 f)$); and the Prandtl number ($Pr = v_v/\kappa_v$). We define $Ek$ based on a fixed mixed layer depth $h_{ml} = 70$ m and the mixed layer mean of the prescribed vertical eddy viscosity profile (Fig. 1e) that is constant in the cross-front direction ($x$) and time. In all cases, $f = 10^{-4}$ s$^{-1}$. The simulations comprise three solution sets: (1) a primary set of 16 cases (4 $Ro$ numbers $\times$ 4 $Ek$ numbers) with $Pr = 1$; (2) a secondary set where vertical buoyancy mixing is suppressed ($\kappa_v = 0, Pr = \infty$); and (3) two cases with fixed $Ro^2/Ek$ (= 350.87) corresponding to the $Ro = 2, Ek = 1.14 \times 10^{-2}$ case (see Sec. 3c). Movies S1 and S2 illustrate the ASC and density evolution for primary and buoyancy mixing suppression solution sets, respectively. The listed $Ek$ contain a 1.14 multiplicative factor for the first two solution sets that comes from the vertical average of the non-dimensional $v_v$ profile in Fig. 1e, which is 0.56. This gives $Ek = 0.56v_{max}/fh_{ml}^2 = 1.14v_{max}$; we refer to $Ek$ in the text and subsequent figures without this factor (e.g., $Ek = 10^{-2}$) for brevity.

<table><thead><tr><th>Solution set</th><th>Number of solutions</th><th><i>Ro</i></th><th><i>Ek</i></th><th><i>Pr</i></th></tr></thead><tbody><tr><td>Primary</td><td>16</td><td>[0.25 0.5, 1, 2]</td><td>1.14 &times; [10<sup>-4</sup>, 10<sup>-3</sup>, 10<sup>-2</sup>, 10<sup>-1</sup>]</td><td>1</td></tr><tr><td>Buoyancy mixing suppression</td><td>16</td><td>[0.25 0.5, 1, 2]</td><td>1.14 &times; [10<sup>-4</sup>, 10<sup>-3</sup>, 10<sup>-2</sup>, 10<sup>-1</sup>]</td><td>&infin;</td></tr><tr><td>Fixed <i>Ro</i><sup>2</sup>/<i>Ek</i></td><td>2</td><td>1,4</td><td>2.85 &times; 10<sup>-3</sup>, 4.56 &times; 10<sup>-2</sup></td><td>1</td></tr></tbody></table>

### c. ROMS experimental design

We solve Eq. (1) with the UCLA ROMS code, utilizing an idealized, 2D ($x, z$) configuration (Fig. 1). ROMS solves the primitive equations in a terrain-following coordinate system with an implicitly hyper-diffusive (Lemarie et al. 2012), 3rd-order upstream advection scheme for horizontal advection and a parabolic spline scheme for vertical advection of momentum and tracers. The idealized configuration employs a flat bottom ($H = 500$ m); periodic boundary conditions in the cross-front ($x$) direction; constant horizontal resolution ($\Delta x = 50$ m) over a domain length of 51.2 km; 128 vertical levels with grid-stretching parameters $\theta_s = 6, \theta_b = 2, h_c = 25$; and constant Coriolis frequency $f = 1 \times 10^{-4}$ s$^{-1}$. The buoyancy is defined as a linear function of temperature with $\rho_0 = 1000$ kg/m$^3$ and thermal expansion coefficient $\alpha = 2 \times 10^{-4}$. In practice, the 2D configuration is achieved with 4 grid-points in the along-front direction ($y$), periodic boundary conditions in $y$, and an initial condition that is uniform in $y$.
<figure><img src="image_2.png" /></figure>

FIG. 1. Idealized 2D (x,z) double front initial conditions (a-d) and prescribed vertical mixing profile (e) normalized by a maximum ($v_{max}$) (a-d): initial temperature (contour lines) and along-front velocity (v, colors). The initial condition is in geostrophic balance and there is no initial across-front velocity ($u(t=0)=0$). The initial Rossby number is defined based on the surface relative vorticity normalized by the Coriolis frequency ($Ro = v_x/f$). The Ekman number ($Ek$) is defined based on the vertical average of $v_v(z)$ in the mixed layer and is modulated via $v_{max}$. The isotherms in (a-d) are the same in every panel and range from 23.74 °C (black lines) to 24.66 °C (white lines) with 0.051 °C change between each isotherm. Any (x,z) snapshots that follow show the same isotherms as here; the simulations employ a linear equation of state, dependent only on temperature, so these isotherms can be interpreted as isopycnals. All analyses focus on the eastern front ($x > 0$). Both fronts behave the same and are separated enough so as to not influence each other. Appendix C details formulations for the initial condition and vertical mixing profile.

The experimental setup prescribes an initial 2D buoyancy field (Fig. 1a-d) with a geostrophic, along-front velocity (v) and triggers frontal evolution with the introduction of a vertically variable mixing profile ($v_v(z), \kappa_v(z)$) that is constant in the cross-front direction and time (Fig. 1e). This idealized setup (Fig. 1) sits between the intended realism of LES solutions (Sullivan and McWilliams 2017; Verma et al. 2019; Wenegrat et al. 2020; Sullivan and McWilliams 2024; Johnson and Fox-Kemper 2024) – that partially resolve the evolving boundary layer turbulence – and the extreme idealization of Crowe and Taylor (2018, 2019, 2020) that primarily investigates fronts with no initial stratification and constant mixing in space and time with a free-slip bottom<fnref n="1" />. Sullivan and McWilliams (2017, 2024) and Crowe and Taylor (2018, 2019, 2020) employ an initial condition that satisfies TTW balance; this modifies both $v(t=0)$ and introduces a secondary

<footnote n="1">¹The complete theoretical model of Crowe and Taylor (2018) allows for initial stratification and spatially variable $v_v, \kappa_v$, noting that the numerical successor paper (Crowe and Taylor 2019) mainly focuses on the simplest cases with zero initial stratification and constant viscosity and diffusivity.</footnote>
circulation (u, w) at $t = 0$ relative to our geostrophic initial condition. Based on testing a subset of cases with a TTW-balanced initial condition (not shown) we do not expect the choice of initial condition (TTW or geostrophic) to impact the results; the secondary circulation develops via a transient adjustment with a geostrophic initial condition (detailed in Sec. 4) and we observe comparable (nonlinear) frontogenesis for both TTW and geostrophic initial conditions given the same Ro, Ek.

The initial condition (Fig. 1a-d) defines a surface mixed layer with a weak stratification ($b_z = N^2 \sim 10^{-7} s^{-2}$) and horizontal buoyancy gradient $b_x$ that transitions to a pycnocline and interior stratification ($N^2 \sim 10^{-5} s^{-2}$). The surface mixed layer depth $h_{ml}(x)$ varies between 60 m and 75 m over a frontal width that is modulated to set the magnitude of the initial horizontal buoyancy gradient. The prescribed vertical eddy viscosity and diffusivity profile (Fig. 1e) is zero at the surface (to ensure zero-stress), reaches a maxima ($v_{max}, \kappa_{max}$) in the mixed layer interior, and transitions to zero at $z = -70$ m (an approximate base of the mixed layer). This mixing shape approximately maintains the mixed layer density structure while minimizing entrainment from the pycnocline. As a way to ensure across-front periodicity, the initial condition comprises a double front (Fig. 1a-d), as opposed to subtracting off a background $b_x$ (Crowe and Taylor 2019). Appendix C details the initial condition and mixing profile formulations. Analyses focus on the front on the eastern end of the domain ($x > 0$), noting that the fronts evolve symmetrically and are separated enough so as to not interact with each other.

We evolve fronts for a range of initial frontal strengths ($\max[b_x(t=0)]$) and vertical mixing intensities (defined based on $v_{max}, \kappa_{max}$; Fig. 1e). This parameter variation translates to varying the initial Rossby number ($Ro = \max[\zeta_{init}]/f$, where $\zeta_{init} = v_x$ at $t=0$) and an Ekman number ($Ek = v_0/fh_{ml}^2$), where $v_0$ is the average of $v(z)$ in the upper 70 m. We choose Ro, Ek values associated with our phenomenological target of mid-latitude submesoscale fronts; this focus excludes a geostrophic adjustment problem of a front with no initial flow ($Ro = 0$) or a gravity current regime ($f \approx 0$; $Ro \to \infty$).

We run 34 total simulations (Table 1): two solutions sets of 16 solutions each as well as 2 more simulations with fixed $Ro^2/Ek$ (or equivalently $Ro/Ek^{1/2}$; Sec. 2b). The primary solution set consists of four $Ro(= 0.25, 0.5, 1, 2$; Fig. 1a-d) by four $Ek = (1.14 \times [10^{-4}, 10^{-3}, 10^{-2}, 10^{-1}])^2$ with $Pr = v_v/\kappa_v = 1$; this value of the turbulent Prandtl number is consistent with parameterizations

<footnote n="2">²In subsequent text and figures we generally do not list the 1.14 factor for Ek; see Table 1 caption.</footnote>
of oceanic vertical boundary layer turbulence (Large et al. 1994) and present conceptions of atmospheric boundary layers (Li 2019). The second set of solutions suppresses the buoyancy diffusivity ($\kappa_v = 0, Pr = \infty$) for the same $Ek, Ro$ as the primary solution set. These solutions allow us to distinguish the roles of $v_v$ and $\kappa_v$ in governing frontal evolution (Sec. 3d). The third set of solutions comprises two cases with fixed $Ro^2/Ek$ (see Sec. 2b) corresponding to a case in the primary solution set with $Ro = 2, Ek = 1.14 \times 10^{-2} (Ro^2/Ek = 350.87)$.

All cases are run for $\approx 5$ inertial periods ($T_i = 2\pi/f$) with a time-step $\Delta t = 15$ s. Model output ($u, v, w, \rho, \eta$) is stored as instantaneous snapshots and saved every 15 minutes. We compute diagnostic terms in the momentum equations online with the model time-step, which we use to describe the controlling balances during frontal evolution (Sec. 4).

The model resolution ($\Delta x = 50$ m) is higher or comparable to previous realistic modeling studies of submesoscale fronts (Gula et al. 2014; Barkan et al. 2019; Qu et al. 2022; Srinivasan et al. 2023) and adequately resolves a submesoscale frontal width of $\sim O(100\ \text{m})$. This resolution allows for a computationally efficient, large parameter sweep (Table 1), while noting that some cases exhibit frontal sharpening that temporarily reaches the grid-scale. We also note that the initial condition contains very weak mixed layer stratification ($N^2 \sim 10^{-7}\ \text{s}^{-2}$), strong vertical shear, and negative potential vorticity (Fig. C1). These conditions make some solutions susceptible to symmetric or shear instabilities. Symmetric instability dominates solution behavior at $Ek = 10^{-4}$, where the fixed vertical mixing is too weak to suppress these motions. While such instabilities are known to occur at submesoscale fronts with stratification and shear similar to our initial conditions (Yu et al. 2019; Peng et al. 2021), they are incompletely resolved in the present simulations due to the inability of the vertical mixing to respond to their onset, the hydrostatic assumption, and resolution limitations. Sec. 5 discusses these instabilities and their impact on our interpretations.

### 3. Vertical mixing induces and inhibits frontogenesis

#### a. Illustrative solutions

Here, we overview the evolution of two solutions (Fig. 2-3) to illustrate two typical responses to vertical mixing: frontogenesis (Fig. 2,3 top; $Ek = 10^{-2}, Ro = 2, Pr = 1$) and frontogenetic inhibition (Fig. 2,3 bottom; $Ek = 10^{-1}, Ro = 2, Pr = 1$). The snapshot sequences in Fig. 2 and 3 show the evolution of the density (contour lines) and ageostrophic secondary circulation (ASC)
<figure><img src="image_3.png" /></figure>

FIG. 2. Snapshots of isopycnals (black lines) and ageostrophic secondary circulation (ASC) streamfunction (colors, $\Psi(x, z)$, where $u = \partial\Psi/\partial z$ and $w = -\partial\Psi/\partial x$) for two cases in the primary solution set ($Pr = 1$; see Table 1); $Ek$ and $Ro$ are indicated in the first snapshot (left column) for each solution. These solutions (both with initial $Ro = 2$, Fig. 1d) demonstrate vertical mixing induced frontogenesis (top; $Ek = 10^{-2}$) and vertical mixing inhibition of frontogenesis (bottom; $Ek = 10^{-1}$). The isopycnal contours for both cases are the same as in Fig. 1. Note the surface-intensified sharpening of the front at $Ek = 10^{-2}$ (top; $t = 1 T_i$) compared to the stronger mixing case (bottom; $Ek = 10^{-1}$), despite the presence of non-zero ASC ($\Psi > 0$). In both cases, introduction of the vertical momentum mixing ($v_v$) induces the generation of the ASC (Sec. 4 details this adjustment). Small scale motions at the base of the mixed layer (most clearly visible in the bottom panels) are due to the prescribed vertical mixing tending to zero at $z = -70$ m; these motions do not significantly impact the target ASC behavior in the mixed layer.

streamfunction in Fig. 2 (colors, $\Psi(x, z)$, where $u = \partial\Psi/\partial z$ and $w = -\partial\Psi/\partial x$) and normalized relative vorticity in Fig. 3 ($\zeta/f = v_x/f$). Movie S1 shows the evolution of density and $\Psi$ for all cases in the primary solution set (Table 1).

In both cases in Fig. 2 (with $Ro = 2$, $Pr = 1$ and differing $Ek$) an ASC develops, indicating a generic momentum adjustment to the introduction of vertical momentum mixing. The ASC acts counterclockwise ($\Psi > 0$) in the $(x, z)$ plane in a manner that favors re-stratification of the mixed layer (pushing lighter water over heavier water). While both cases exhibit this ASC, the frontal evolution is different. With stronger mixing ($Ek = 10^{-1}$; Fig. 2 bottom), the ASC weakens and there is no discernible frontal sharpening (or spreading), indicated by the relatively constant $\zeta/f$ (Fig. 3). With weaker mixing ($Ek = 10^{-2}$; Fig. 2 top), the front sharpens, (over approximately 1 inertial period), with this sharpening characterized by amplification of $b_x$ (contour lines; top row Fig. 2 and 3), $\zeta/f$ (colors; Fig. 3 top), and $\Psi$ (colors; Fig. 2 top). We characterize the top row
<figure><img src="image_4.png" /><figcaption>FIG. 3. As in Fig. 2, but showing the relative vorticity $\zeta/f = v_x/f = Ro$ in colors.</figcaption></figure>

of Fig. 2,3 as frontogenesis and the bottom row of Fig. 2,3 as frontogenetic inhibition, noting that frontolysis occurs at later time for a majority of the cases with $Pr = 1$ (see weakening of the ASC in Fig. 2 top right, Movie SI, or Fig. 4). Sec. 4 provides more detailed descriptions of the dynamics controlling these characteristic sequences.

b. Solution regimes

Fig. 4 provides a visualization of solution regimes in the $(Ek, Ro)$ space with $(x, t)$ Hovmöller plots of the surface, cross-front buoyancy gradient that is normalized relative to the maximum $b_x$ in the initial condition:

$$ \tilde{b}_x = \frac{b_x(x, t)}{\max[b_x(x, t = 0)]} - 1. \qquad (7) $$

We observe three regimes in the primary solution set ($Pr = 1$; Fig. 4): ‘typical’ frontogenesis via the ASC (Fig. 2 top, Fig. 4b,c), frontogenetic inhibition or frontolysis (Fig. 2 bottom, e.g., Fig. 4d,h), and frontogenesis via (incompletely resolved) symmetric instability (Fig. 4i). We define frontogenesis in Fig. 4 as a sustained increase of the (surface) buoyancy gradient over $\approx 1$ inertial period. These frontogenetic cases (e.g., Fig. 4b,c) are characterized by a growth and peak in $\tilde{b}_x$ that generally weakens at later time (e.g., Fig. 4c). This late-time weakening can be due to (viscously damped) inertial oscillations that reverse the ASC and weaken $b_x$ (e.g., see $x \approx 12$ km, $t \approx 2 T_i$ in Fig. 4b). These late-stage oscillations of the ASC, detailed further in Sec. 4, occur in many of the solutions and are less damped with weaker mixing (see Movie S1).

In all cases the front moves laterally (to the left in Fig. 4), with this movement most pronounced at larger $Ro$ and smaller $Ek$. Some of the strongly frontogenetic cases also exhibit the formation
<figure><img src="image_5.png" /><figcaption>FIG. 4. Vertical mixing induced frontogenesis and frontolysis across Ekman (Ek; columns) and Rossby numbers (Ro; rows) in the primary solution set (Table 1). Each panel shows a (x,t) Hovmöller of the surface cross-front buoyancy gradient that is normalized relative to maximum buoyancy gradient in the initial condition ($\tilde{b}_x$) (Eq. 7); $\tilde{b}_x$ is dimensionless. All panels share the same colorbar, which is log-scaled for better visualization. The vertical time axis has units of inertial period $T_i = 2\pi/f$. Here, frontogenesis appears as a red 'streak'. Note the stronger frontogenesis at smaller Ek and larger Ro and extremely weak to no frontogenesis at larger Ek and smaller Ro.</figcaption></figure>

of a secondary front (e.g., Fig. 4a, $x \gtrsim 14$ km). However, we do not focus on this behavior, which is likely a byproduct of the fixed (weak) vertical mixing combined with overturning motions associated with incompletely resolved symmetric instability (discussed below and in Sec. 5a).

Fig. 4 illustrates general trends of stronger frontogenesis (darker reds) for decreasing Ek or increasing Ro. In particular, there is a frontogenetic 'sweet spot' (Fig. 4b,c,f,g) at large Ro ($Ro \gtrsim 1$) and intermediate Ek ($Ek = 10^{-3} - 10^{-2}$). Very strong mixing ($Ek = 10^{-1}$; right column,
<figure><img src="image_6.png" /><figcaption>FIG. 5. Surface (left) and vertically averaged (right) metrics of frontogenesis as a function of $Ro^2/Ek$ (see Eq. 6; Sec. 2b) for solutions with $Pr = 1$. The top left legend in (b) indicates $Ek$ (color) and $Ro$ (marker shape) for each case in the primary solution set (first row, Table 1). The green and orange octagons represent cases with fixed $Ro^2/Ek$ (corresponding to the red star case; $Ro = 2$ $Ek = 10^{-2}$), but different $Ro$ and $Ek$ (with $Ro$ indicated in the legend, see bottom row Table 1). The top panels (a,b) measure frontogenesis as the ratio of the maximum buoyancy gradient ($\max[b_x]$) normalized by the initial condition ($\max[\text{init. } b_x]$); values greater than 1 (horizontal dashed line) indicate frontal sharpening. The bottom panels (c,d) measure frontogenesis as the maximum buoyancy frontogenetic tendency rate ($\max[\mathcal{T}_{\text{tot}}]$); $\mathcal{T}_{\text{tot}}$ is defined in Eq. (8)). For both frontogenesis metrics, maximum values are taken over the entire simulation period (5.5 inertial periods) with most cases with exhibiting maximum $b_x$ before $\approx 3$ inertial periods (see Fig. 4). The inset text in (a) summarizes interpretations, detailed in Sec. 3c. Note that solutions with both low $Ek$ and $Ro$ (e.g., grey circle) exhibit frontogenesis due to symmetric instability (see Sec. 5a) that is distinct from the frontogenesis via an ageostrophic secondary circulation induced by vertical momentum mixing (e.g., red star).</figcaption></figure>

Fig. 4) or very weak initial fronts ($Ro = 0.25$; bottom row, Fig. 4) exhibit extremely weak or no frontogenesis for $t \le 3 T_i$.
The weakest mixing cases ($Ek = 10^{-4}$, Fig. 4) also exhibit strong frontogenesis that takes a distinct three-front structure, particularly for smaller $Ro$ (e.g., Fig. 4i). This frontogenesis is due to the onset of (improperly resolved) symmetric instability (SI), which occurs most prominently at $Ek = 10^{-4}$ and small $Ro$, but also appears at $Ek = 10^{-3}$ to varying degrees (e.g., Fig. 4j). We do not over-interpret the low $Ek$ and low $Ro$ solutions dominated by later-time SI (Fig. 4i,j,m,n where SI is not visible in Fig. 4m,n because it occurs at $t \approx 5 T_i$) primarily due to the fixed vertical mixing assumption that limits their fluid dynamical validity. We discuss these instabilities further in Sec. 5a.

All solutions exhibit a transient response (i.e., non-steady $b_x$), indicating the breaking of geostrophic balance by the vertical momentum mixing (or onset of SI). In the typical (non-SI) frontogenetic cases (e.g., Fig. 4b,c) the increase in $b_x$ is driven by the convergent ASC (e.g., Fig. 2, top) that develops in response to this balance-breaking and subsequently amplifies during nonlinear frontogenesis (Sec. 4). Fig. 4 makes apparent the $Ek$ and $Ro$ dependence on the time-scale over which the ASC and frontogenesis develop. As $Ek$ or $Ro$ increase, the initiation of frontogenesis (red streaks in Fig. 4) occurs earlier. However, larger $Ek$ also results in faster and stronger erosion of the front (e.g., late-time presence or absence of dark red moving from left to right in the top row of Fig. 4). This suggests an intrinsic competition between the frontogenetic ASC – triggered (faster) by (stronger) vertical momentum mixing – and the vertically diffusive erosion of the front triggered (faster) by (stronger) vertical buoyancy mixing.

### c. Controlling buoyancy dynamics

If we assume that the primary function of the momentum dynamics in our system is to trigger the ASC (detailed in Sec. 4), we can assume that the buoyancy equation (Eq. 2c) controls whether or not fronts sharpen or weaken in our 2D simulations. Sec. 2b describes two separately derived nondimensional parameters ($Ro^2/Ek \sim Ro/Ek^{1/2}$) that equivalently measure the relevant competition in the buoyancy equation, despite each parameter originating from a different choice for scaling the cross-front velocity (Barkan et al. 2019; Thomas and Ferrari 2008). These parameters measure the competition between the (frontogenetic) cross-front buoyancy advection ($ub_x$) versus the (frontolytic) vertical buoyancy diffusion ($\frac{\partial}{\partial z}[\kappa_v b_z]$); more generally, these ratios scale the nonlinear term in the buoyancy equation ($D_t B$; Eq. 2c) when $Pr = 1$ (see Appendix B).
Fig. 5 demonstrates that these ratios approximately map the transition between frontogenetically inhibited or frontolytic ($Ro^2/Ek \lesssim 50$) and frontogenetic ($Ro^2/Ek \ge 50$) solutions, where we somewhat arbitrarily choose $Ro^2/Ek$ for illustrative purposes. We map these transitions using two metrics of frontogenesis: the maximum buoyancy gradient normalized by the initial condition maxima (max [$b_x$] / max [init. $b_x$]); Fig. 5a,b) or the maximum buoyancy frontogenetic tendency rate (max [$T_{tot}$]; Fig. 5c,d), where $T_{tot}$ represents the change of amplitude of the buoyancy gradient following a fluid parcel (i.e., the rate of frontal sharpening; see Eq. 8 for a complete definition). These metrics are defined relative to the maximum taken over all time. The surface metrics of frontogenesis (Fig. 5a,c) most successfully map regime transitions. The mixed layer average metrics (Fig. 5b,d) more moderately demonstrate predictive utility of $Ro^2/Ek$ noting that frontogenesis is primarily a near-surface process.

In Fig. 5a, the frontal sharpening increases with $Ro^2/Ek \ge 50$ and approximately plateaus for $Ro^2/Ek \gtrsim 10^3$, indicating a grid-scale constraint on $b_x$ that occurs at weakest $Ek$; this trend of increasing frontogenesis for larger $Ro^2/Ek$ also holds when mapping solutions based on $T_{tot}$ (Fig. 5c). Note that the cases with SI-induced frontogenesis (e.g., $Ro = 0.5, Ek = 10^{-4}$; Fig. 4i) also approximately collapse on these curves of surface frontogenesis (Fig. 5a,c). In particular this is due to taking maximum $b_x$ and $T_{tot}$ over the whole simulation period, which captures very late-stage ($t \approx 4-5 T_i$) onset of SI e.g., in the $Ro = 0.25, Ek = 10^{-4}$ solution (grey circle in Fig. 5, Fig. 4m). Again, we do not place too much importance on these cases.

We perform an additional test of $Ro^2/Ek$ by running two cases at fixed $Ro^2/Ek$ ($\sim Ro/Ek^{1/2}$) corresponding to the characteristic frontogenetic case with $Ro = 2, Ek = (1.14\times)10^{-2}$ (red star Fig. 5, Fig. 4c). The green ($Ro = 1, Ek = 2.85 \times 10^{-3}$) and orange ($Ro = 4, Ek = 4.56 \times 10^{-2}$) octagons in Fig. 5 represent these cases with $Ro^2/Ek = 350.87$. Both of these additional cases exhibit frontogenesis analogous to corresponding case in the primary solution set (red star Fig. 5, Fig. 4c), albeit with different time-evolution (not shown). This qualitative agreement in solution behavior further indicates the utility of $Ro^2/Ek$ (or the practically equivalent $Ro/Ek^{1/2}$) as a predictor of vertical mixing induced frontogenesis. While the larger $Ro$ case (orange octagon) appears as an outlier relative to the other two cases (red star and green octagon) when measuring maximum $b_x$ relative to the initial condition (Fig. 5a), all three cases exhibit comparable frontal sharpening rates (Fig. 5c,d). The reduced max [$b_x$] / max [init. $b_x$] for the $Ro = 4$ case (orange octagon Fig.
<figure><img src="image_7.png" /><figcaption>FIG. 6. As in Fig. 4, but for cases with the vertical buoyancy mixing suppressed ($\kappa_v = 0, Pr = \infty$; second row, Table 1) Note the sustained frontogenesis ($\tilde{b}_x \gg 0$) in the stronger mixing cases (e.g., $Ek = 10^{-1}$), unlike the frontogenetic inhibition for these cases in Fig. 4 (where $\kappa_v = \nu_v$). All panels share the same, log-scaled colorbar as in Fig. 4.</figcaption></figure>

5a) relative to the other two cases is likely due to the grid-scale constraint on frontogenesis that artificially inhibits frontal sharpening earlier for the large initial $b_x(\propto Ro)$.

The approximate regime collapse on $Ro^2/Ek$ in Fig. 5 is a primary result of this study and constitutes an attempted unification of the competing views on the role of vertical mixing in submesoscale frontogenesis (Sec. 1a). In Sec. 5b we discuss how application of a scaling in Crowe and Taylor (2018) (where $D_t b \sim Ro/Ek$) partially succeeds in mapping frontogenesis in our solutions. Because of the idealized setup, particularly the inability of the vertical mixing to respond to the frontal evolution (discussed further in Sec. 5), we do not over-interpret the exact magnitudes in Fig. 5 (e.g., $Ro^2/Ek \sim 50$ as a transition point). For example, the shape of the
mixing profile (e.g., constant $v_v$, $\kappa_v$) could (marginally) modify this transition point. Instead, the
utility of this parameter mapping is the insight it provides into the governing dynamics of the
frontal evolution. That is, because $Ro^2/Ek$ (and not e.g., $1/RoEk$, which arises when scaling the
momentum equations with $u_{B19}$; see Appendix B) appears to explain the solution behavior for all
cases with $Pr = 1$, it suggests a strong control on the frontal evolution by the *non-conservative*
buoyancy equation. The next section (Sec. 3d) explicitly demonstrates this control of vertical
buoyancy mixing on solution behavior via suppression of $\kappa_v$.

### d. Suppression of vertical buoyancy mixing

Here, we distinguish the role of $v_v$ and $\kappa_v$ by suppressing the vertical buoyancy mixing ($\kappa_v = 0$, $Pr = \infty$) in a twin set of solutions (Table 1). Fig. 6-7 demonstrate frontal evolution in these solutions with $\kappa_v = 0$ and allow comparison with their counterparts with $\kappa_v = v_v$. The Hovmöller plots in Fig. 6 are analogous to Fig. 4.

Strikingly, frontogenesis occurs for large *Ek* and small *Ro* when $\kappa_v = 0$ (Fig. 6, right column), remembering that these solutions are frontogenetically inhibited when $\kappa_v = v_v$ (Fig. 4, right column). In particular, note the *sustained* amplification of $\tilde{b}_x$ in Fig. 6 for these larger *Ek* solutions (right two columns). This result explicitly demonstrates that the role of vertical *buoyancy* mixing is to suppress frontogenesis, which is driven by the vertical *momentum* mixing induced ASC. Taking the limit of $Pr \to \infty$ in Eq. 5 or 6 offers one explanation for this result. Note that the $Ek = 10^{-4}, 10^{-3}$ cases with buoyancy mixing suppressed (Fig. 6 left two columns) appear similar to their $Pr = 1$ analogs (Fig. 4), suggesting the limited role of $\kappa_v$ in those solutions with $Pr = 1$.

Fig. 7 compares solutions with $v_v = \kappa_v$ ($Pr = 1$; solid lines) and $\kappa_v = 0$ ($Pr = \infty$; dashed lines) at different *Ek* (colors) and *Ro* (rows). We plot the time-series of the maximum surface buoyancy gradient (left column) and maximum surface convergence normalized by *f* (an ASC Rossby number; right column). Again, note the striking, sustained amplification of $\max[b_x]$ for the cases with $\kappa_v = 0$ compared to their $v_v = \kappa_v$ counterparts (Fig. 7 left column). For example, note the difference between the black ($Ek = 10^{-1}$) solid ($\kappa_v = v_v$) and dashed ($\kappa_v = 0$) lines in Fig. 7a,c,e. For $\kappa_v = 0$, the frontogenetic rate (i.e., the slope of $\max[b_x]$) is strongest for largest *Ro* and, interestingly, smaller *Ek*. Fig. 6 and 7 also demonstrate, particularly at large *Ek*, that when $\kappa_v = 0$
<figure><img src="image_8.png" /><figcaption>FIG. 7. Demonstration of the impact of vertical buoyancy mixing ($\kappa_v$) on frontogenesis with comparison of solutions with $Pr = 1$ ($\kappa_v = \nu_v$, solid) and $Pr = \infty$ ($\kappa_v = 0$, dashed). The panels show time-series of (left) the maximum surface buoyancy gradient and (right) the maximum surface convergence normalized by the Coriolis frequency $f$ for all $Ro$ (rows) at $Ek = 10^{-3}$ (red), $Ek = 10^{-2}$ (blue) and $Ek = 10^{-1}$ (black). The time axis has units of inertial period $T_i = 2\pi/f$. The horizontal dashed line on the left panels indicates the initial condition value, which is the same for each $Ro$. Note that all the dashed lines ($\kappa_v = 0$, $Pr = \infty$) exhibit frontogenesis, even at large $Ek$, and a plateau in $b_x$ (a,b) indicates that the front has reached the grid-scale.</figcaption></figure>

frontogenesis continues until reaching the grid-scale; this occurs earlier and is most apparent at larger $Ro$. This again emphasizes how vertical buoyancy mixing inhibits frontogenesis.
All solutions generally exhibit the same initial growth rate of the the convergence (Fig. 7b,d,f,h) at early time ($t \lesssim 1 T_i$). This indicates an inertial control on the initial adjustment (detailed in Sec. 4). The time-series of convergence also show that this initial convergent ASC is strongest for the strongest mixing (largest $Ek$, black lines). However, in the cases with $\kappa_v = 0$ (dashed lines), the smaller $Ek$ exhibit stronger later-time frontogenetic rate. That is, the dashed red line in Fig. 7a is larger in magnitude and slope than the dashed black line. This seemingly counter-intuitive trend in $Ek$ for $\kappa_v = 0$ solutions – stronger frontogenesis for smaller $Ek$, despite stronger initial convergence for larger $Ek$ – indicates a frontolytic role of vertical momentum mixing *at later time*, described in Sec. 4.

# 4. Dynamical balances

Here, we detail controlling dynamical balances from the perspective of buoyancy and momentum evolution. We diagnose terms in various evolutionary equations (defined in Sec. 4a) that collectively demonstrate the mechanisms governing frontal evolution in frontogenetic and frontolytic or frontogenetically inhibited solution regimes (Sec. 3b). We exemplify these mechanisms with detailed analysis of the two characteristic cases in Fig. 2-3 ($Ro = 2, Pr = 1$ and $Ek = 10^{-1}, 10^{-2}$) as well as their $Pr = \infty$ analogs (Fig. 8-11). The latter allows us to further distinguish the roles of $v_v$ and $\kappa_v$.

A generic, mechanistic description of frontal evolution can be summarized as follows (for $Pr = 1$): vertical momentum mixing induces a convergent ASC via a transient, inertial adjustment–i.e., the linear, transient TTW (or $T^3W$) balance (Wenegrat and McPhaden 2016; Dauhajre and McWilliams 2018; Johnson et al. 2020b)–with $Ro, Ek$ modulating the resulting ASC magnitude and structure. Subsequent frontogenesis, if it occurs, is dominated by nonlinear advection of momentum and buoyancy. Early-stage ($t \lesssim 2 T_i$) frontogenetic inhibition or frontolysis (larger $Ek$) results from vertical buoyancy diffusion that competes comparably with horizontal buoyancy advection. Later-stage ($t \gtrsim 1.5-2 T_i$) frontolysis at smaller $Ek$ can result from viscously damped inertial oscillations that reverse the ASC. These ASC oscillations weaken the front in a manner that is qualitatively consistent with shear dispersion (Young and Jones 1991; Crowe and Taylor 2018; Wenegrat et al. 2020; Swart et al. 2020), although this is not explicitly diagnosed here.
<figure><img src="image_9.png" /><figcaption>FIG. 8. Time-series of horizontally and vertically averaged buoyancy frontogenetic tendency terms (Eq. 8; terms indicated in the legend) for four cases with initial $Ro = 2$: $Ek = 10^{-2}$ (a,b), $Ek = 10^{-1}$ (c,d) and $Pr = 1$ (a,c) and $Pr = \infty$ (b,d). Here, $T_{tot} > 0$ (black) indicates frontogenesis. Note the different y-axis ranges for each panel. Time-series are obtained with spatial averaging in a front-following window that tracks the maximum cross-frontal buoyancy gradient; here, the averaging is done 400 m around the maximum cross-front buoyancy gradient and vertically in the upper $\approx 5$ m. Note the strong frontogenesis at $Ek = 10^{-2}$ (a) that is dominated by horizontal advection ($T_u$); this sharpening is inhibited with stronger buoyancy vertical mixing (c).</figcaption></figure>

a. Definitions

The buoyancy frontogenetic tendency equation quantifies the rate of change of the amplitude of the buoyancy gradient ($|\nabla_h b|$) following a fluid parcel (Hoskins 1982). For the 2D system here, $|\nabla_h b| = b_x$ and the frontogenetic tendency equation is written as:

$$ \underbrace{\frac{1}{2} D_t (b_x)^2}_{T_{tot}} = \underbrace{-b_x b_x u_x}_{T_u} \underbrace{-b_x b_z w_x}_{T_w} + \underbrace{b_x \frac{\partial}{\partial x} \left( \frac{\partial}{\partial z} (\kappa_v b_z) \right)}_{T_{vmix}} . \quad (8) $$

We diagnose the frontogenetic ($T_{tot} > 0$) and frontolytic ($T_{tot} < 0$) contributions from horizontal advection ($T_u$), vertical straining ($T_w$), and vertical diffusion ($T_{vmix}$). In Eq. 8 and other balance equations (defined below), we consider horizontal diffusion negligible. It arises in the present
simulations from the implicit hyper-diffusion of the third-order upwind advection scheme and is always frontolytic.

We interpret the controlling momentum dynamics with the cross- and along-front momentum balances, as well as divergence and vorticity equations. We diagnose terms in the ageostrophic, cross-front momentum balance as follows:

$$ \underbrace{\frac{\partial u}{\partial t}}_{\text{Rate}_u} = \underbrace{-uu_x - wu_z}_{\text{Adv}_u} + \underbrace{fv_{ag}}_{\text{Cor}_{ag,u}} + \underbrace{\frac{\partial}{\partial z}(v_v u_z)}_{\text{Vmix}_u}, \qquad (9) $$

where $Cor_{ag,u} = fv_{ag} = f\nu - \phi_x$ is the ageostrophic Coriolis term. We split the material derivative to both isolate the nonlinear contribution ($Adv_u$) and diagnose the importance of the linear diagnostic TTW balance ($Cor_{ag,u} = -Vmix_u$) and transient TTW balance ($T^3W; Rate_u = Cor_{ag,u} - Vmix_u$).

Similarly, we diagnose terms in the along-front momentum balance, which makes clear the role of $v_v$ in breaking geostrophic balance as well as the dynamics controlling the initial generation of the secondary circulation. The along-front momentum balance is defined as follows:

$$ \underbrace{\frac{\partial v}{\partial t}}_{\text{Rate}_v} = \underbrace{-uv_x - wv_z}_{\text{Adv}_v} - \underbrace{fu}_{\text{Cor}_v} + \underbrace{\frac{\partial}{\partial z}(v_v v_z)}_{\text{Vmix}_v}, \qquad (10) $$

where terms are analogous to those in Eq. 9, with TTW balance ($Cor_v = -Vmix_v$)

The divergence equation provides additional, useful perspective on the ASC evolution (with divergence $\delta = u_x$):

$$ \underbrace{D_t \delta}_{\text{Rate}_\delta} = \underbrace{-\delta^2}_{\text{Hadv}_\delta} + \underbrace{f \zeta_{ag}}_{\text{Cor}_{ag,\delta}} + \underbrace{w_x u_z}_{\text{Vadv}_\delta} + \underbrace{\frac{\partial}{\partial x} \left( \frac{\partial}{\partial z} (v_v u_z) \right)}_{\text{Vmix}_\delta}. \qquad (11) $$

$Rate_\delta < 0$ near the surface indicates amplification of a frontogenetic ASC (i.e., intensification of surface convergence), with contributions from horizontal advection ($Hadv_\delta$), the inertial ageostrophic residual ($Cor_{ag,\delta} = f\zeta - \phi_{xx}$), vertical advection ($Vadv_\delta$), and vertical mixing ($Vmix_\delta$).

Finally, diagnosis of terms in the vorticity equation (where $\zeta = v_x$) enables additional insight on the particular role of vertical momentum mixing in breaking geostrophic balance and/or (weakly)
<figure><img src="image_10.png" /><figcaption>FIG. 9. As in Fig. 8 but for the cross-front momentum balance terms (Eq. (9), terms indicated in the legend). The TTW residual (TTW<sub>u</sub> = Cor<sub>ag,u</sub> + Vmix<sub>u</sub>) is shown in dashed grey. Here, the time-series are obtained with spatial averaging 400 m around the maximum cross-front buoyancy gradient and vertically in the upper ≈ 2.6 m. Note the transient, linear adjustment in all cases (Rate<sub>u</sub> ≈ Cor<sub>ag,u</sub> + Vmix<sub>u</sub>) in response to vertical momentum mixing breaking geostrophic balance. At smaller Ek (a), momentum advection (green) dominates frontogenesis, while at larger Ek (c) vertical buoyancy diffusion (Fig. 8c) inhibits nonlinear frontogenesis and TTW balance dominates (approximately equal and opposite red and blue curves).</figcaption></figure>

damping frontogenesis (of velocity gradients) at later time:

$$ \underbrace{D_t \zeta}_{\text{Rate}_\zeta} = \underbrace{-\zeta \delta}_{\text{Hadv}_\zeta} \underbrace{-w_x v_z}_{\text{Vadv}_\zeta} \underbrace{-f \delta}_{\text{Cor}_\zeta} + \underbrace{\frac{\partial}{\partial x} \left( \frac{\partial}{\partial z} (\nu_v v_z) \right)}_{\text{Vmix}_\zeta}, \quad (12) $$

Rate<sub>ζ</sub> > 0 near the surface indicates amplification of the front (increase in cyclonic vorticity) with contributions from vortex stretching (Hadv<sub>ζ</sub>), vortex tilting (Vadv<sub>ζ</sub>), stretching of planetary vorticity (Cor<sub>ζ</sub>), and vertical mixing (Vmix<sub>ζ</sub>).
<figure><img src="image_11.png" /><figcaption>FIG. 10. As in Fig. 9 but for the along-front momentum balance terms (Eq. (10), terms indicated in the legend). The TTW residual ($TTW_v = Cor_v + Vmix_v$) is shown in dashed grey.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 11. As in Fig. 8 but for the divergence balance terms (Eq. (11), terms indicated in the legend). The TTW residual ($TTW_δ = Cor_{ag,δ} + Vmix_δ$) is shown in dashed grey. Here, the time-series are obtained with spatial averaging 400 m around the maximum cross-front buoyancy gradient and vertically in the upper ≈ 5 m.</figcaption></figure>

b. *Frontogenetic and frontolytic balances*

1) BUOYANCY FRONTOGENETIC TENDENCY

The buoyancy frontogenetic tendency diagnostics (Fig. 8) demonstrate that horizontal buoyancy advection ($T_u > 0$) drives frontogenesis (here, illustrated for $Ek = 10^{-2}$, green curve, Fig. 8a).
<figure><img src="image_13.png" /><figcaption>FIG. 12. As in Fig. 8 but for the vorticity balance terms (Eq. (12), terms indicated in the legend). The TTW residual ($TTW_{\zeta} = Cor_{\zeta} + Vmix_{\zeta}$) is shown in dashed grey. The spatial averaging is the same as in Fig. 11.</figcaption></figure>

At stronger mixing ($Ek = 10^{-1}$, Fig. 8b), there is extremely weak advective frontogenesis (note the different y-axis ranges between panels) that competes against the frontolytic impact of vertical mixing ($T_{vmix} < 0$; Fig. 8 orange curve), which is comparable in magnitude to $T_u$.

When buoyancy mixing is suppressed ($Pr = \infty$, Fig. 8b,d), both Ekman number cases exhibit strong frontogenesis driven by horizontal advection ($T_u > 0$). In these cases, the frontolytic impact of vertical straining ($T_w < 0$, purple curves) is more apparent, noting that $T_w < 0$ in a typical, frontogenetic ASC (e.g., Fig. 8a). Note the sustained, increasing frontogenesis in Fig. 8d ($Ek = 10^{-1}, Pr = \infty$), which dominates over the early-time frontogenetic signal that is not visible.

## 2) MOMENTUM, DIVERGENCE, AND VORTICITY

The introduction of vertical eddy viscosity $v_v$ initially acts to mix the along-front geostrophic momentum (via $Vmix_v$ in Eq. 10 and analogously $Vmix_{\zeta}$ in Eq. 12), noting that $u, \delta = 0$ initially. This balance-breaking via $v_v$ is apparent in both the along-front momentum and vorticity balances (red and black curves for $t \sim 0 T_i$ in Fig. 10 and 12), and occurs across all solutions, albeit with differing time-scales dictated by $Ek$.

The breaking of the initial geostrophic balance by the vertical momentum mixing generates an ageostrophic secondary circulation ($u, w, \delta$), as well as an ageostrophic component of the
along-front momentum ($v_{ag}, \zeta_{ag}$). The early-time along-front momentum balance (Fig. 10, Eq. 689) illustrates that the secondary circulation (i.e., $Cor_v = -fu$; blue curve) is balanced by the 690 vertical mixing ($Vmix_v$) and acceleration ($Rate_v$). That is, the adjustment to balance-breaking 691 by $v_v$ drives a secondary circulation via linear, transient turbulent thermal wind ($T^3W$) dynamics 692 ($Rate_v \approx Cor_v + Vmix_v$). This early-time $T^3W$ balance is indicated by the agreement between dashed 693 grey ($TTW_v$, $TTW_u$) and black ($Rate_v$, $Rate_u$) curves in the along- and across-front momentum 694 balances (Fig. 9,10; $t \lesssim 0.5 T_i$). Despite this universal balance-breaking mechanism, Fig. 9-12 695 show a dichotomy in the later-time momentum, divergence, and vorticity balances for solutions 696 with weak (top panels) and strong (bottom panels) vertical mixing. 697

The smaller $Ek$, frontogenetic solutions (e.g., Fig. 9a,b, 10a,b 11a,b, and 12a,b) exhibit a 698 three-stage evolution: 699

1. A linear, transient ($T^3W$) adjustment to balance-breaking via vertical momentum mixing that 700 initially drives the secondary circulation (as described above). This is most clearly seen in the 701 momentum diagnostics (Fig. 9a,b, 10a,b) 702

2. The transition to- and dominance of nonlinear (conservative) frontogenesis via the ASC, 703 which is most strikingly observed in the dominance of horizontal advection in the divergence 704 balance ($Rate_\delta \approx Hadv_\delta$; Fig. 11a,b; Eq. 11). 705

3. Subsequent frontolysis resulting from both the weakening of the ASC ($Cor_v \to 0$ in Fig. 10a 706 and $Rate_\delta \to 0$ in Fig. 11a) and vertical buoyancy diffusion (e.g., Fig. 8a, orange curve). 707

Note that the vorticity balance also demonstrates the dominance of horizontal advection in the 708 strongest, nonlinear frontogenetic stage (Fig. 12a,b; green curve), with the rate of velocity gradient 709 frontogenesis ($Rate_\delta$, $Rate_\zeta$) decreasing after approximately $\approx 0.6 - 0.7 T_i$ (for both $Pr = 1$ 710 and $Pr = \infty$). 711

After frontogenesis in the small $Ek$ solutions, $Rate_u$, $Rate_v$, $Rate_\delta$ and $Rate_\zeta$ move towards zero, 712 interestingly following the nonlinear terms (i.e., agreement between green and black curves after 713 the peak at $t \approx 0.7 - 1.5 T_i$ in Fig. 9a, 11a, 12a). After this decrease of the Rate terms following 714 frontogenesis, the ASC in the small $Ek$ solutions notably oscillates. These oscillations, some of 715 which result in reversals of the ASC (positive to negative $\Psi$), are apparent in Movie S1. The later- 716 time ($t \gtrsim 2 T_i$) momentum balances (not visible in Fig. 9,10) show an approximate $T^3W$ balance 717
during these oscillations of the ASC. We interpret these oscillations as viscously damped inertial
oscillations, which arise more prominently (less damped) with weaker mixing (smaller Ek). Visual
inspection of the Ro = 2, Ek = 10⁻², Pr = 1 case (Fig. 2 top) and other weaker-mixing cases in
Movie S1 illustrates that the weakening (and reversals) of the ASC after the initial frontogenetic
peak can act to spread the previously sharpened front. We interpret this late-stage frontolysis as
analogous to shear dispersion (Young and Jones 1991; Crowe and Taylor 2018; Wenegrat et al.
2020; Swart et al. 2020). However, we note that ASC reversals (and later-stage frontal weakening)
also occur when buoyancy mixing is suppressed (Pr = ∞).

At large Ek (Fig. 9c, 11c), where frontogenesis is inhibited, TTW balance (Cor<sub>ag,u</sub> ≈ −Vmix<sub>u</sub>;
Cor<sub>ag,δ</sub> ≈ −Vmix<sub>δ</sub>) dominates and there is negligible nonlinearity (see red versus blue curves
in both Fig. 9c and 11c). However, there is still a transient adjustment in this case via T³W
(agreement between dashed grey and black at very early time); this indicates a generic balance-
breaking adjustment, regardless of the mixing amplitude. TTW balance is not as apparent in the
along-front momentum or vorticity balances for this case (Fig. 10c, 12c), indicating some transient,
although negligible, evolution of the along-front velocity.

Interestingly, the TTW dominance holds for Pr = ∞ at Ek = 10⁻¹ (Fig. 9d, 11d), noting
the sustained (yet near-constant) amplification of convergence (Rate<sub>δ</sub> < 0; Fig. 11d), vorticity
(Rate<sub>ζ</sub> > 0; Fig. 12d), and buoyancy gradient (T<sub>tot</sub> > 0; Fig.8d) <sup>3</sup>. Additionally, despite the
suppression of vertical buoyancy mixing, the smaller Ek case with Pr = ∞ is more frontogenetic
than the larger Ek case (i.e., T<sub>tot</sub> is larger in Fig. 8b than 8d). The vorticity balance (Fig. 12d)
illustrates that vertical momentum mixing of the vorticity (red curve Fig. 12d) – as well as the TTW
residual (grey curve Fig. 12d) – transition from early-time frontogenetic (amplifying the cyclonic
vorticity; Vmix<sub>ζ</sub> > 0) to later-time frontloytic (eroding the geostrophic vertical shear; Vmix<sub>ζ</sub> < 0).
This highlights a subtle, relatively weak dampening of frontogenesis by vertical momentum mixing
at later time, despite its initial role in inducing the frontogenetic ASC.

<footnote n="3"><sup>3</sup>The buoyancy gradient frontogenesis T<sub>tot</sub> > 0 plateaus similarly to the Rate<sub>δ</sub>, Rate<sub>ζ</sub> terms in Fig. 11d,12d, although this is not visually apparent with the 2 T<sub>i</sub> time-limit in Fig. 8d.</footnote>
## 5. Discussion

### a. Symmetric and gravitational instabilities at low Ek

The frontal initial condition (Fig. 1), while designed to be quasi-realistic in structure and magnitude of buoyancy gradients, has a small Richardson number ($Ri = N^2 f^2 / b_x^2 < 1$) and contains negative potential vorticity (Fig. C1); the negative PV is due to the strong cross-front buoyancy gradient (as opposed to negative absolute vorticity). This negative PV can lead to the onset of (unforced) symmetric instability (SI; Hoskins (1974)), which particularly dominates solution behavior (over motions associated with ASC or shear instability) at $Ek = 10^{-4}$ and $Ro = 0.25, 0.5$. We illustrate the onset of SI in one of these weak mixing solutions ($Ro = 0.5, Ek = 10^{-4}, Pr = 1$; Fig. 13) to exemplify how these motions lead to frontogenesis (e.g., Fig. 4i) and subsequent gravitational instability that is improperly represented in the hydrostatic model used in this study. We also comment on the potential for SI motions to effectively amplify the vertical boundary layer turbulence (and thus Ek) at the front.

Snapshots of overturning circulation qualitatively indicate the onset of SI (Fig. 13a-c). Note that in this solution no discernible ASC develops for $t \lesssim 2 T_i$. Instead, multi-signed overturning cells appear after $\approx 2 T_i$ and align approximately along isopycnals (Fig. 13a). These overturning cells are distinct from the single-signed, larger-scale ASC generated via TTW or T³W (Fig. 2 top).

We diagnose the geostrophic shear production (GSP; Thomas et al. (2013)) to quantitatively identify these motions as SI, where GSP > 0 indicates SI. GSP is defined as:

$$GSP = -\overline{v'w'} \frac{\partial \overline{v_g}}{\partial z}, \qquad (13)$$

where the overbar indicates a horizontal average; the primes denote the horizontal anomaly; and $v_g$ the geostrophic, along-front velocity. We define the horizontal average and anomaly relative to a 3.5 km region surrounding the front, indicated by the dashed grey lines in Fig. 13a-c.

Fig. 13d-e demonstrates that the emergence of the multi-signed overturning circulation cells are associated with GSP > 0 ($t \approx 2.2-3.1 T_i$), indicating that the unstable motions derive energy from the geostrophic vertical shear. The multiple overturning cells associated with the onset of SI result in frontogenesis of multiple fronts (most visible in Fig. 13c and Fig. 4i), but also create negative stratification (not shown) that makes the solution gravitationally unstable. Other solutions
<figure><img src="image_14.png" /><figcaption>FIG. 13. Example of symmetric instability (SI) in a solution with $Ro = 0.5$, $Ek = 10^{-4}$, $Pr = 1$ (Fig. 4i). (a-c): snapshots of the overturning streamfunction ($\Psi$) and isopycnals (time indicated at the top of each panel). (d-e) geostrophic shear production (GSP; Eq. 13) as a function of depth and time (d) and vertically averaged in the upper 71 m (black curve, left y-axis); GSP is computed for the area between the dashed grey vertical lines in (a)-(c), which surround the front. The grey curve in (d) shows the vertical average of the diagnosed vertical eddy viscosity associated with the SI motions ($\nu_{v,SI}$; Eq. 14), corresponding to the right y-axis. The vertical dashed lines in (d-e) indicate the snapshot times in (a-c). Note the approximately along-isopycnal overturning cells in (a-c) – which are structurally distinct from the ASC in Fig. 2 – that occur with $GSP > 0$ (d-e), indicating the onset of SI. Also note the ‘kink’ in the mixed layer isopycnals in (a) ($x \approx 13$ km, $z \approx -40$ m), indicating how SI motions create negative stratification, that leads to gravitational instability that cannot be resolved with the hydrostatic model used in this study.</figcaption></figure>

with low $Ro = 0.25, 0.5$ and low $Ek = 10^{-3}, 10^{-4}$ exhibit variants of this behavior, with onset of SI occurring later for smaller $Ro$. Additionally, solutions bordering this region of the parameter space can exhibit less intense SI motions along with ASC-driven frontogenesis (e.g., Fig. 4e; see Movie S1).

Previous studies (Verma et al. 2019; Wenegrat et al. 2020) which more completely resolve SI turbulence at submesoscale fronts demonstrate a complex interplay between SI, frontogenesis, and frontolysis. Fundamental to this interplay is the amplification of boundary layer turbulence by the SI motions, which can feedback on the front. Following Bachman et al. (2017), we compute a turbulent vertical eddy viscosity associated with SI:

$$ \nu_{v,SI} = \frac{\text{GSP}}{(\partial v_g / \partial z)^2}. \qquad (14) $$
The grey curve in Fig. 13e demonstrates that (the vertically averaged) $v_{v,SI}$ increases to approximately $6 \times 10^{-3} \text{m}^2\text{s}^{-1}$ ($t \approx 2.4 T_i$). This viscosity is an order of magnitude larger than the prescribed $v_{max} = 10^{-4} \text{m}^2\text{s}^{-1}$ (Fig. 1d) for this solution. Taking this peak in $v_{v,SI}$ gives an effective SI Ekman number $Ek_{SI} \approx 10^{-2}$ (assuming the same $h_{sbl}$ as the initial condition). This $Ek_{SI}$ is two orders of magnitude larger than the prescribed $Ek = 10^{-4}$ via $v_v$ in the experimental setup (Sec. 2c). In the paradigm of this idealized study, the amplification of $Ek$ by the effective SI mixing ($v_{v,SI}$) suggests that SI could actually result in frontolysis (i.e., smaller $Ro^2/Ek(\approx 25)$), assuming that $Pr \sim 1$ for SI turbulence (as demonstrated in Bachman et al. (2017)).

The above interpretation is consistent with Wenegrat et al. (2020), who suggest with measurements and LES simulations that shear dispersion driven frontolysis of a front in the Gulf Stream results from SI amplification of the vertical buoyancy mixing. Conversely, Verma et al. (2019) demonstrate – in a 3D LES simulating the spindown of a geostrophic front with initial $Ro = 0.32$ – that the onset of SI results in shear instabilities, with both symmetric and shear instabilities supplying boundary layer turbulence that induces a re-stratifying, larger-scale (frontogenetically favorable) secondary circulation; D’Asaro et al. (2011) also observe SI driving a frontogenetic ASC in measurements. Given the seemingly contradictory role of SI in these studies, the interpretation of an effective SI viscosity (as in Eq. 14) amplifying $Ek$ is potentially too simplistic. A suite of SI-resolving LES simulations would be required to systematically investigate the interactions of SI-induced turbulence, frontogenesis, and frontolysis over the parameter space explored in this study.

Symmetric and gravitational instabilities have been observed at fronts in the real ocean (Thomas et al. 2016; Yu et al. 2019; Peng et al. 2021) with measured mixed layer stratification ($N^2 \sim 10^{-5} - 10^{-7} \text{s}^{-2}$) comparable to our initial condition ($N^2 \sim 10^{-7} \text{s}^{-2}$). While these and other studies (Thomas et al. 2013; Haney et al. 2015; Bachman et al. 2017; Wenegrat et al. 2020; Dong et al. 2021; Chor et al. 2022) evidence typicality of symmetric or gravitational instabilities at submesoscale fronts, their emergence in the present solutions, particularly at low $Ek$, is unconstrained due to the model resolution, hydrostatic assumption, and most importantly, inability of the fixed vertical mixing to respond to these unstable motions. Given the incompleteness of these unstable motions in our simulations, we do not over-interpret solution behavior in the low $Ek$, low $Ro$ portion of the parameter space (Fig. 4i,j,m,n).
## b. Comparisons with past interpretations

The present results reconcile previously competing interpretations (Sec. 1a) of whether and how vertical mixing induces sharpening or weakening of submesoscale fronts, with Fig. 5 (regime collapse on $Ro^2/Ek$) quantitatively summarizing our attempt at a unified paradigm for this problem. These previous interpretations take the separate views that either momentum dynamics or buoyancy dynamics dominate frontal evolution in response to vertical mixing. The former view anticipates vertical mixing induced frontogenesis via TTW at large $Ro$ (McWilliams et al. 2015; Sullivan and McWilliams 2017, 2024) and the latter anticipates vertical mixing induced frontolysis via shear dispersion or vertical diffusion at small $Ro$ and intermediate to large $Ek$ (Crowe and Taylor 2018, 2019). Our simulations demonstrate all of these frontogenetic or frontolytic mechanisms and we map their relative dominance across a broader ($Ek, Ro$) space (Fig. 4) compared to previous individual studies (see Table A1).

The identification of a single parameter (either $Ro^2/Ek$ or $Ro/Ek^{1/2}$; see Sec. 2b) as an approximate governor of vertical mixing induced frontogenesis results from the insight that the competition between cross-front advection of buoyancy and vertical buoyancy diffusion fundamentally controls solution outcome (Fig. 6, 7). Thompson (2000) and Thomas and Ferrari (2008) highlight this same competition, with Thompson (2000) arriving at a different scaling (see their Appendix A) in a semi-geostrophic framework.

The parameters $Ro^2/Ek$ and $Ro/Ek^{1/2}$ equivalently measure these competing buoyancy dynamics. They arise from different scaling choices for the cross-front velocity ($u_{TF08}, u_{B19}$; Sec. 2b), which originate in Barkan et al. (2019) and Thomas and Ferrari (2008), respectively. Application of the scaling in Crowe and Taylor (2018) gives $D_t b \sim Ro/Ek$ (with $Pr = 1$) for their initial transient (more frontogenetic) regime, and interestingly, $D_t b \sim Ro^2/Ek$ for their slow shear dispersion regime. $Ro/Ek$ exhibits some success for mapping maximum $b_x$ (as in Fig. 5a,b) in our solutions, but does not select for the frontal sharpening rates (as in Fig. 5c,d) as well as $Ro^2/Ek$ (and the equivalent $Ro/Ek^{1/2}$; see Fig. 1 in supplemental materials). Fundamentally, $Ro^2/Ek, Ro/Ek^{1/2}$ and $Ro/Ek$ all quantify the advective (frontogenetic) versus diffusive (frontolytic) competition. Importantly, $Ro/Ek^{1/2}$ and $Ro^2/Ek$ emerge from scaling choices that capture relevant regimes of a submesoscale ASC: linear, TTW dynamics that generate or sustain the ASC as a pre-cursor to
849 frontogenesis ($u_{TF08}$; Eq. 3) and nonlinear dominance observed during frontogenesis ($u_{B19}$; Eq. 4; detailed further in Appendix B).
850

851 While past studies refer to vertical mixing induced frontogenesis as 'TTW frontogenesis' (McWilliams et al. 2015), the TTW balance in the present simulations is only valid in frontogenetically inhibited cases with large $Ek$ (Fig. 9c). In our simulations, a linear, transient adjustment creates the frontogenetic secondary circulation (Fig. 9a) that transitions to a nonlinear balance during peak frontogenesis. Of course, this transient adjustment is a consequence of the geostrophic initial condition, however, we note that simulations with an initial condition in TTW balance (not shown) behave analogously.
852
853
854
855
856
857

We primarily focus on early-time solution behavior ($t \lesssim 2$ inertial periods), noting the superinertial nature of submesoscale frontogenesis (Barkan et al. 2019). This early-time focus contrasts with the long-time ($t \gtrsim 10$ inertial periods) focus of Crowe and Taylor (2018, 2019, 2020) (and their theoretical predecessor (Young 1994)). The transience of early-time frontogenesis in the present simulations (Fig. 7 solid curves) is, however, qualitatively consistent with the even more idealized simulations in Crowe and Taylor (2019) (see their Fig. 7); that study initializes a front with zero vertical stratification and employs a free-slip bottom and uniform mixing. We note that Crowe and Taylor (2019) interpret observation of early-time frontogenesis in their simulations to result from nonlinearity, analogous to our results. We observe that momentum advection and convergence dominate our strongest frontogenetic cases (Fig. 9a, 11a), with this frontogenesis occurring for $\lesssim 1$ inertial period. This result corroborates the inviscid, asymptotic model of Barkan et al. (2019), which predicts that the surface convergence of the ASC (as opposed to deformation) dominates the frontogenetic tendency of buoyancy (Fig. 8), divergence (Fig. 11), and vorticity (Fig. 12).
858
859
860
861
862
863
864
865
866
867
868
869
870

When $Pr = 1$, vertical mixing induced frontogenesis does not 'run away' to a singularity, as anticipated in inviscid theory for strain-induced frontogenesis (Hoskins and Bretherton 1972) and some interpretations of TTW frontogenesis (Sec. 6 in McWilliams et al. (2015)). Instead, the temporally fixed vertical buoyancy diffusion acts as the primary inhibitor of frontal sharpening. This frontal weakening by $\kappa_v$ occurs via vertical diffusion at early-time (large $Ek$; Fig. 8c) or via shear dispersion at later time (small $Ek$), where viscously damped inertial oscillations of the cross-frontal flow aid in spreading the previously sharpened front (see Movie S1). While Wenegrat et al. (2020) evidences shear dispersion at submesoscale fronts in the Gulf Stream, the late-time
871
872
873
874
875
876
877
878
879 oscillations (and associated frontal spreading) in our simulations may result artificially from the fixed mixing assumption and 2D posing; the fixed mixing limits the ability of the (weak) vertical mixing to damp inertial oscillations and the 2D posing excludes 3D instabilities (e.g., baroclinic mixed layer, horizontal shear) that may preclude these late-stage (damped) inertial oscillations (discussed further in Sec. 5c).

884 As in Bodner et al. (2019), we attempt to distinguish the roles of $v_v$ and $\kappa_v$ in sharpening or weakening fronts, albeit with a different approach. Bodner et al. (2019) treat the vertical mixing as a first-order correction to inviscid, strain-induced frontogenesis theory (Shakespeare and Taylor 2013), while we prescribe the vertical eddy viscosity and diffusivity in a primitive equation system with no straining. In our posing, the vertical eddy viscosity is necessary for inducing the ASC and initiating frontogenesis. This contrasts with Bodner et al. (2019), who find that $v_v$ weakens the strain-driven frontogenesis. However, we also demonstrate that $v_v$ can actually damp frontogenesis at later time via mixing of the along-front velocity and vorticity (Fig. 12d). We do not observe that $\kappa_v$ enhances frontogenesis, in contrast to Bodner et al. (2019), who observe frontogenetic enhancement by $\kappa_v$ at later-time; although, they note this later-time is beyond the limit of the perturbation approach.

### c. Applicability of interpretations to submesoscale fronts in nature

895 A utility of this study is the prediction of frontal evolution given a measure of frontal strength (Ro) and vertical mixing intensity (Ek), with $Ro^2/Ek$ quantifying the competition between cross-front buoyancy advection and vertical diffusion that approximately governs solution outcome (Fig. 5). If we assume $Pr = v_v/\kappa_v = 1$, this competition could be measured locally at a front as:

$$ \frac{Ro^2}{Ek} \approx \frac{\zeta^2 h_{ml}^2}{\kappa_v f} \quad \text{or} \quad \frac{\delta^2 h_{ml}^2}{\kappa_v f}. \qquad (15) $$

900 While this metric (or its practically equivalent square root, i.e., $Ro/Ek^{1/2}$; see Sec. 2b) could potentially explain whether real submesoscale fronts sharpen or weaken, the exclusion of other intrinsic processes in the idealized posing may limit applicability of our interpretations to the real ocean. These additional processes primarily include the response of vertical mixing to frontal evolution, straining or deformation flows, and 3D instabilities.
The fixed vertical mixing assumption allows us to treat *Ro* and *Ek* as independent parameters in the present idealized framework. In reality, *Ro* and *Ek* are not independent; the boundary layer turbulence (*Ek*) evolves spatially and temporally in response to the frontal evolution (*Ro*) and vice versa. We illustrate the spatial variability of *Ek* at submesoscale fronts and filaments in a realistically configured simulation in Fig. 2 of the supplementary materials, leaving comprehensive investigation of *Ro* and *Ek* dependencies in such simulations for future work. Johnson and Fox-Kemper (2024) demonstrate that the partially resolved boundary layer turbulence in 3D LES of a frontal spin down responds differently to fronts than traditional surface boundary layer turbulence parameterizations (e.g., KPP); this result makes interpretation of *Ek* in realistic simulations, such as the one in Fig. 2 of the supplemental information, more difficult. However, past numerical studies of submesoscale dense filaments – with both partially resolved (Sullivan and McWilliams 2017, 2024) and parameterized (Gula et al. 2014; McWilliams et al. 2015) vertical boundary layer turbulence – consistently demonstrate a horizontal structure in vertical mixing: stronger mixing at the filament center relative to the surrounding, re-stratified regions. In the context of our 2D framework, the spatio-temporal response of the vertical mixing to frontal evolution may alter the (late-time) frontolytic behavior in the present idealizations, which results from a fixed vertical eddy diffusivity and/or inertial oscillations (Sec. 4). Vertical buoyancy mixing that reaches further into the pycnocline, relative to our posing (Fig. 1e), can mix stratified water into the mixed layer and accelerate frontolysis. More generally, spatial structure in the vertical boundary layer turbulence raises questions regarding the utility or most dynamically apt definition of *Ek* at a front or filament

Straining currents can induce or maintain the frontogenetic secondary circulation, separate from vertical momentum mixing via TTW (or transient TTW). This straining can be supplied by mesoscale currents or submesoscale mixed layer eddies (Boccaletti et al. 2007; Zhang et al. 2021). While Bodner et al. (2019) attempt to diagnose the relative roles of straining and vertical mixing during 2D frontogenesis, there remains an open question regarding the role of these processes at different stages in a frontal life-cycle. A less-highlighted, but relevant result of this study is the demonstration that temporally fixed vertical mixing does not induce frontogenesis for initial *Ro* = 0.25 (Fig. 4m-p, excluding the very late-time SI induced frontogenesis in m,n; see Sec. 5a). This result indirectly suggests that either straining currents or vertical mixing response to the front are required to drive a transition from *Ro* ~ 0.1 to *Ro* >> 1. However, this view assumes that
submesoscale fronts ‘start’ with a particular (mesoscale or mixed layer eddy) *Ro* and motivates
clarification on ‘typical’ precursor conditions to frontogenesis.

Along-front, horizontal shear instabilities (Sullivan and McWilliams 2024; Wu et al. 2022; Gula et al. 2014) provide a separate route to frontal erosion or arrest that can preclude the vertical diffusion or shear dispersion frontolytic mechanisms demonstrated in this 2D study. The respective roles of vertical buoyancy diffusion and instabilities (horizontal shear, vertical shear, centrifugal, symmetric, mixed layer baroclinic) in driving frontolysis remain to be systematically quantified, while also noting that some of these instabilities can actually set up the ASC as demonstrated in Verma et al. (2019) and D’Asaro et al. (2011). The expectation is that the competition between the vertical mixing rate ($h_{ml}^2/\kappa_v$) and the growth rate of the instability – both of which compete with (or contribute to; Verma et al. (2019)) frontogenesis by the secondary circulation – determines the dominant frontolytic mechanism.

Recent observations (Swart et al. 2020) and realistic simulations (Sun et al. 2020) demonstrate that strong winds can erode submesoscale fronts, with Swart et al. (2020) suggesting shear dispersion via inertial oscillations as a frontolytic mechanism initiated by winds. While we capture analogous behavior in our simulations, particularly at lower *Ek* (Sec. 4), our results demonstrate that strong vertical buoyancy diffusion also drives frontolysis (at small *Ro*²/*Ek*). This motivates further work to diagnose the relative roles of vertical diffusion and shear dispersion during strong wind events. Absent from these frontolytic paradigms is the consideration of favorably aligned wind-stresses that can drive frontogenesis (Crowe and Taylor 2020) as well as wind-driven across-front buoyancy fluxes that create conditions for symmetric instabilities (Thomas et al. 2013, 2016). Such instabilities will enhance localized mixing ($v_v, \kappa_v$), but may also initiate a transition to 3D, frontolytic instabilities (e.g., horizontal shear instability) as described above. This interplay remains relatively unexplored.

### d. Implications for realistic modeling and parameterization of submesoscale fronts

The idealized demonstration of *Ro*²/*Ek* controlling vertical mixing induced frontogenesis suggests one explanation for the observation of a resolution-dependent transition in frontogenesis (and *Ro*) in realistically configured simulations. Previous studies (Capet et al. 2008; Barkan et al. 2017; Garabato et al. 2022; Srinivasan et al. 2023) collectively indicate that a transition from
‘submesoscale-permitting’ to ‘submesoscale-resolving’ occurs at horizontal resolutions of approximately 500 m; *O* (1 km) simulations do not adequately capture nonlinearity and forward energy fluxes (Garabato et al. 2022; Srinivasan et al. 2023) associated with frontogenesis at higher resolutions. Our results suggest one interpretation for this observed regime transition across resolution: nonlinear frontogenesis induced by vertical mixing can only occur once *Ro*²/*Ek* crosses a threshold. This interpretation relies on an (untested) assumption that (the bulk) *Ek* remains constant across resolution and responds primarily to the atmospheric forcing; that is, the smaller Δ*x* allows for larger *Ro*, given the same *Ek*. However, a population of submesoscale fronts could similarly modify large-scale *Ek* via integrated re-stratification fluxes. In that sense, while our results suggest one paradigm for *Ro*, *Ek* control on observed resolution-dependencies of submesoscale motions in realistic models, further work is required to determine the validity of this hypothesis.

Recently designed parameterizations of submesoscale re-stratification fluxes rely on assumption that the steady-state TTW balance well-predicts either the width of “stable” fronts (Bodner et al. 2023) or the submesoscale secondary circulation and thus re-stratification (Yang et al. 2024). This study demonstrates that non-steady and nonlinear dynamics dominate re-stratification during frontogenesis (Sec. 4). In the present simulations, the TTW balance is only valid in the frontogenetically inhibited cases (Fig. 9c). This result partially supports an interpretation that TTW controls the ‘arrested’ frontal width (Bodner et al. 2023), while noting that vertical buoyancy diffusion primarily drives frontolysis in these large *Ek* cases (Fig. 8c; Fig. 7a black curve). More generally, the present demonstrations of non-steady momentum and buoyancy balances – whether during linear secondary circulation adjustment, nonlinear frontogenesis, or vertical buoyancy diffusion driven frontolysis – motivate consideration of transient dynamics in designing parameterization of submesoscale fluxes, while noting the caveats of the idealization (Sec. 5c).

## 6. Summary and conclusions

This study re-litigates the role of vertical mixing in submesoscale frontogenesis and frontolysis (Sec. 1a) with a suite of idealized simulations (Fig. 4,6; Table 1). These simulations evolve 2D fronts initially in geostrophic balance over a range of initial frontal strengths (*Ro*) and vertical mixing intensities (*Ek*), where the introduction of vertical mixing ($v_v$, $κ_v$) triggers a frontal evolution. Our problem posing (Sec. 2) prescribes an initial surface mixed layer front that is guided by
realism (in *Ro* and stratification); however, in order to explicitly isolate the role of vertical mixing, we artificially hold the vertical eddy viscosity and diffusivity profiles (Fig. 1e) as constant in time.

We observe that vertical mixing can both induce and inhibit frontogenesis (Fig. 4), primarily focusing on the first ≈ 1–2 inertial periods, which corresponds to a typical life-cycle of a submesoscale front. The identification of two practically equivalent controlling parameters ($Ro^2/Ek \sim Ro/Ek^{1/2}$; Sec. 2b) that approximately map regime transitions (Fig. 5) quantitatively signifies our attempt at a common paradigm for this problem (assuming $v_v = \kappa_v$). Both $Ro^2/Ek$ and $Ro/Ek^{1/2}$ measure the competition between cross-front buoyancy advection and vertical diffusion. These ratios reflect particular scaling choices for the cross-front velocity–originating in Barkan et al. (2019) and Thomas and Ferrari (2008), respectively–that capture different aspects of the ageostrophic secondary circulation that drives submesoscale frontogenesis: non-conservative, linear dynamics (i.e., TTW) that can generate or sustain an ASC (Thomas and Ferrari 2008) and conservative, nonlinear dynamics during frontogenesis when $Ro > 1$ (Barkan et al. 2019). We also note the potential applicability of an analogous parameter ($Ro/Ek$) utilizing the scaling of Crowe and Taylor (2018) (see Sec. 5b), despite its theoretical limitation to $Ro < 1$.

The controlling dynamics elucidated in this study blend and update previous interpretations of vertical mixing impacts on submesoscale frontogenesis: turbulent thermal wind (TTW) frontogenesis and shear dispersion or vertical diffusion frontolysis (Sec. Sec. 1a and 5b). For large $Ro^2/Ek$, vertical momentum mixing can induce a transition to nonlinear, convergence dominated (Barkan et al. 2019) frontogenesis via the generation of an ageostrophic secondary circulation (McWilliams et al. 2015; McWilliams 2017; Sullivan and McWilliams 2017, 2024). Conversely, for small $Ro^2/Ek$, vertical buoyancy mixing suppresses frontogenesis via strong vertical diffusion that inhibits frontal sharpening by the secondary circulation (Crowe and Taylor 2018, 2019). This distinction between the generally frontogenetic impact of $v_v$ and solely frontolytic impact of $\kappa_v$ is made explicit with simulations that set $\kappa_v = 0$ (Fig. 6), which all exhibit frontogenesis, remarkably, even at large $Ek$. In all simulations, transient TTW dynamics (Wenegrat and McPhaden 2016; Dauhajre and McWilliams 2018) – as opposed to the usually-invoked diagnostic TTW balance (Garrett and Loder 1981; Gula et al. 2014; McWilliams et al. 2015; McWilliams 2017) – generate the secondary circulation; we observe that TTW balance is valid only in the frontogenetically in-
hibited cases (Fig. 9, Sec. 4) or unrealistically frontogenetic cases with $κ_v = 0$ (Fig. 6), particularly
at large *Ek*.

We expect similar results for submesoscale dense filaments, a common (McWilliams et al. 2015; McWilliams 2017; Sullivan and McWilliams 2017, 2024) and dynamically relevant target for this problem. We note that Sullivan and McWilliams (2024) observe weak-to-no frontogenesis in a 3D LES for an initially weak dense filament that is subject to strong atmospheric cooling; this is qualitatively consistent with the small *Ro*, large *Ek* outcome here.

The present simulations also exhibit late-stage (≥ 1.5 – 2 inertial periods) frontolysis at larger *Ro*²/*Ek* (after early-time frontogenesis) that is qualitatively consistent with shear dispersion (Young and Jones 1991; Crowe and Taylor 2018; Wenegrat et al. 2020; Swart et al. 2020) as well as frontogenesis induced by symmetric instability (Verma et al. 2019) at very small *Ek*, *Ro* (e.g., Fig. 4i). However, we caution interpretation of these regimes due to the fixed vertical mixing assumption that limits solution validity at later-time (Sec. 5a,c). More generally, the assumption of fixed vertical mixing in this study places a fundamental limitation on extrapolating the present interpretations to submesoscale fronts nature, where *Ro* and *Ek* are not independent. Future work can interrogate the applicability of the present idealized framework (e.g., predictive utility of *Ro*²/*Ek*) in more realistic scenarios.
*Acknowledgments.* We thank Leif Thomas for very constructive comments during the review process; in particular, for his comments on scaling choices. We also thank Matthew Crowe and anonymous reviewer for similarly constructive comments and suggestions that improved this manuscript. D.P. Dauhajre thanks Jacob Wenegrat for insightful discussion of the study. D.P. Dauhajre, D. Hypolite and J.C. McWilliams were supported by the National Science Foundation (NSF) grant OCE-2124174 and the Office of Naval Research (ONR) grant N00014-23-1-2812. J. Gula was supported by the French National Agency for Research (ANR) through the project DEEPER (ANR-19-CE01-0002-01), and by ONR grant N00014-23-1-2226.

*Data availability statement.* The simulations in this study are run with the UCLA ROMS code: https://github.com/CESR-lab/ucla-roms. Setup files for the idealized double front configuration are available at: https://github.com/ddauhajre/IdealFrontVmixROMS

# APPENDIX A

## Inventory of past numerical and theoretical studies

Table A1 summarizes the (*Ek*, *Ro*) parameter space of past idealized modeling or theoretical studies that investigate the role of vertical mixing in submesoscale frontogenesis. This table is provided as context for the discussion in Sec. 1 and 5.

TABLE A1. Summary of past numerical and theoretical studies that investigate vertical mixing impacts on frontogenesis. The listed values of *Ro* (of either an initial condition or applicability in theory) and *Ek* are either reported in the studies or estimated here, with blank values indicating that *Ro* or *Ek* are either not reported and/or difficult to estimate (e.g., for Large eddy simulations).

<table><thead><tr><td>Study</td><td>Approach</td><td><i>Ro</i></td><td><i>Ek</i></td><td>Vertical mixing representation</td></tr></thead><tbody><tr><td>Thompson (2000)</td><td>Semi-geostrophic model</td><td>&lt; 0.1</td><td></td><td>constant</td></tr><tr><td>McWilliams et al. (2015)</td><td>2D primitive equation model</td><td>≈ 2</td><td>≈ 0.05 – 0.1</td><td>K-profile parameterization</td></tr><tr><td>McWilliams (2017)</td><td>TTW+Omega equation diagnostics</td><td>≤ 1</td><td>≈ 0.05 – 0.1</td><td>analytical formulation</td></tr><tr><td>Sullivan and McWilliams (2017)</td><td>3D Large eddy simulation</td><td>≥ 1</td><td></td><td>partially resolved</td></tr><tr><td>Crowe and Taylor (2018, 2020)</td><td>Asymptotic theory</td><td>&lt; 1</td><td>O(1), O(<i>Ro</i>)</td><td>constant</td></tr><tr><td>Crowe and Taylor (2019)</td><td>2D Large eddy simulation</td><td>0.1 – 1</td><td>0.01 – 1</td><td>constant</td></tr><tr><td>Bodner et al. (2019)</td><td>Perturbation analysis</td><td>0.4</td><td></td><td>first-order correction to strain theory</td></tr><tr><td>Verma et al. (2019)</td><td>3D Large eddy simulation</td><td>0.32</td><td></td><td>partially resolved</td></tr><tr><td>Sullivan and McWilliams (2024)</td><td>3D Large eddy simulation</td><td>0.3 – 4</td><td></td><td>partially resolved</td></tr></tbody></table>
Application of $u_{TF08}$, $u_{B19}$ scalings

Here, we detail the scaled complete primitive equation systems (Eq. 1) that result from applying both scalings for the cross-front velocity described in Sec. 2b: $u_{TF08}$ (Eq. 3) and $u_{B19}$ (Eq. 4). We follow generally the scaling choices in Barkan et al. (2019), with an interchangeable choice for the cross-front velocity (denoted by $U_*$). This scaling makes an additional assumption of anisotropy of submesoscale fronts ($l/L \ll 1$) in a surface mixed layer of depth $h_{ml}$, where $l$ and $L$ are across- and along-front length scales, respectively. The 2D system (Eq. 1) has no along-front dimension, so $l/L \ll 1$ by construction. The scaling is as follows:

$$x \sim l, \quad z \sim h_{ml} \qquad (B1a)$$

$$v \sim V, \quad u \sim U_*, \quad w \sim U_* \frac{h_{ml}}{l} \qquad (B1b)$$

$$t \sim \frac{l}{U_*}, \quad \phi \sim fVl, \quad b \sim \frac{fVl}{h_{ml}} \qquad (B1c)$$

$$v_v \sim v_0, \quad \kappa_v \sim \kappa_0 \qquad (B1d)$$

where $U_* = u_{T08}$ or $u_{B19}$.

Applying $U_* = u_{TF08}$ gives:

$$\frac{Ro}{Ek^{1/2}} [D_t u] - \frac{1}{Ek^{3/2}} [v] = -\frac{1}{Ek^{3/2}} [\phi_x] + [u_{zz}], \qquad (B2a)$$

$$\frac{Ro}{Ek^{1/2}} [D_t v] + \frac{1}{Ek^{1/2}} [u] = [v_{zz}], \qquad (B2b)$$

$$Pr \frac{Ro}{Ek^{1/2}} [D_t b] = [b_{zz}], \qquad (B2c)$$

$$[\phi_z] = [b], \qquad (B2d)$$

$$[u_x] + [w_z] = 0. \qquad (B2e)$$
Applying $U_* = u_{B19}$ gives:

$$ \frac{Ro^2}{Ek} [D_t u] - \frac{1}{RoEk} [v] = - \frac{1}{RoEk} [\phi_x] + [u_{zz}], \qquad (B3a) $$

1078

$$ \frac{Ro^2}{Ek} [D_t v] + \frac{Ro}{Ek} [u] = [v_{zz}], \qquad (B3b) $$

1079

$$ Pr \frac{Ro^2}{Ek} [D_t b] = [b_{zz}], \qquad (B3c) $$

1080

$$ [\phi_z] = [b], \qquad (B3d) $$

1081

$$ [u_x] + [w_z] = 0. \qquad (B3e) $$

In the above systems, we organize terms so that the emergence of the controlling parameters ($Ro/Ek^{1/2}$ and $Ro^2/Ek$; see Sec. 2b and 3c) is apparent in the momentum and buoyancy equations. However, it is more useful to re-organize the terms in the momentum equations to compare the scaled momentum dynamics resulting from $u_{TF08}$ and $u_{B19}$. We focus on the along-front momentum equations, dividing Eq. B2b by $Ro/Ek^{1/2}$ and Eq. B3b by $Ro^2/Ek$:

1082

1083

1084

1085

1086

$$ [D_t v] + \frac{1}{Ro} [u] = \frac{Ek^{1/2}}{Ro} [v_{zz}], \qquad (B4a) $$

1087

$$ [D_t v] + \frac{1}{Ro} [u] = \frac{Ek}{Ro^2} [v_{zz}], \qquad (B4b) $$

1088

1089

1090

1091

1092

1093

The distinction between the $u_{TF08}$ and $u_{B19}$ scalings arises as expected: the scaled vertical mixing term $[v_{zz}]$ is more dominant, relative to the Coriolis term $(1/Ro[u])$, in Eq. B4a ($u_{TF08}$) than in Eq. B4b ($u_{B19}$). That is, the $u_{TF08}$ scaling identifies the linear dynamics, non-conservative dynamics that generate or sustain an ASC $(1/Ro[u])$, while the $u_{B19}$ scaling identifies the nonlinear (approximately conservative) dynamics that 'run away' with frontogenesis, under an assumption of a pre-existing front and ASC.

1094

# APPENDIX C

1095

## Idealized initial condition and vertical mixing profile
The double front initial condition (Fig. 1) prescribes a 2D buoyancy field ($b(x, z)$) and associated geostrophic (along-front) velocity ($v(x, z)$). Fig. C1 shows the potential vorticity of the initial condition for each Ro. Negative potential vorticity ($q = (v_x + f)b_z - |b_x|^2/f$) in the initial condition – a consequence of quasi-realistic $b_z, b_x$ – leads to the onset of symmetric instability in cases with weak mixing (see Sec. 5a).

The construction of $b(x, z)$ follows McWilliams (2017):

$$b(x, z) = b_0 + N_b^2 (z + H) + \frac{N_0^2}{2} \left[ (1 + B)z - (1 - B) \left( h_{ml}(x) + \lambda^{-1} \log \cosh[\lambda(z + h_{ml}(x))] \right) \right]. \quad (C1)$$

where $N_b^2$ is a minimum background stratification, $N_0^2$ the interior stratification, $\lambda$ a scale of the transition between the surface boundary layer and interior stratification that exhibits a fractional reduction in stratification of B.

The mixed layer depth ($h_{ml}(x)$) sets the double front structure (shape and magnitude of $b_x$):

$$h_{ml}(x) = h_0 + \delta h \left[ \tanh(M_f(x - x_f)) - \tanh(M_f(x + x_f)) \right], \quad (C2)$$

where $h_0$ is the mixed layer depth away from the front; $h_0 + \delta_h$ the mixed layer depth at the front; and $\pm x_f$ the location of the front. We modulate the initial frontal strength (Ro in Fig. 1a-d) via $M_f$ in Eq. C2; we set $M_f = [3.11 \times 10^{-4}, 4.4 \times 10^{-4}, 6.24 \times 10^{-4}, 8.83 \times 10^{-4}]$ m$^{-1}$ for Ro = [0.25, 0.5, 1, 2], respectively.

We set the following (in all simulations) relative to Eq. C1-C2:

$$b_0 = 5 \times 10^{-2} \text{ms}^{-2}, \quad B = 0.025 \quad (C3)$$

$$N_0^2 = 3 \times 10^{-5} \text{s}^{-2}, \quad N_b^2 = 10^{-7} \text{s}^{-2}, \quad \lambda^{-1} = 8 \text{ m} \quad (C4)$$

$$h_0 = 60 \text{ m}, \quad \delta_h = 15 \text{ m}, \quad x_f = 12.8 \text{ km} \quad (C5)$$

The prescribed vertical mixing profile ($v_v = \kappa_v$) is cubic and non-zero only above a threshold mixed layer depth ($h_{ml}^* = 70$ m):

$$v_v(z) = v_{\max} \frac{z'(1-z')^2}{0.14805} \quad (C6)$$
<figure><img src="image_15.png" /><figcaption>FIG. C1. Potential vorticity $q = (v_x + f)b_z - |b_x|^2/f$ for the four Ro double front initial conditions, as in Fig. 1.</figcaption></figure>

where $z' = \frac{\eta-z}{h_{ml}^*}$ and $v_{max}$ sets the magnitude (and thus $Ek$). For $f = 10^{-4}$ and $h_{ml} = 70$ m, we set $v_{max} = [10^{-4}, 10^{-3}, 10^{-2}, 10^{-1}] \text{ m}^2\text{s}^{-1}$ to give $Ek = 1.14 \times [10^{-4}, 10^{-3}, 10^{-2}, 10^{-1}]$, where the 1.14 factor comes from the vertical average of $\frac{v_v(z)}{v_{max}}(= 0.56)$ in the upper 70 m.

## References

[1117] Bachman, S., B. Fox-Kemper, J. Taylor, and L. Thomas, 2017: Parameterization of frontal symmetric instabilities. i: Theory for resolved fronts. *Ocean Modelling*, **109**, 72–95, https://doi.org/10.1016/j.ocemod.2016.12.003.

[1120] Bachman, S. D., and J. R. Taylor, 2016: Numerical simulations of the equilibrium between eddy-induced restratification and vertical mixing. *Journal of Physical Oceanography*, **46** (3), https://doi.org/10.1175/JPO-D-15-0110.1.

[1123] Barkan, R., J. C. McWilliams, A. F. Shchepetkin, M. J. Molemaker, L. Renault, A. Bracco, and J. Choi, 2017: Submesoscale dynamics in the northern gulf of mexico. part i: Regional and seasonal characterization and the role of river outflow. *Journal of Physical Oceanography*, **47** (9), 2325–2346, https://doi.org/10.1175/JPO-D-17-0035.1.

[1127] Barkan, R., M. Molemaker, K. Srinivasan, J. McWilliams, and E. D'Asaro, 2019: The role of horizontal divergence in submesoscale frontogenesis. *J. Phys. Oceanogr.*, **49**, 1593–1618, https://doi.org/10.1175/JPO-D-18-0162.1.

[1130] Boccaletti, G., R. Ferrari, and B. Fox-Kemper, 2007: Mixed layer instabilities and restratification. *J. Phys. Oceanogr.*, **37**, 2228–2250, https://doi.org/10.1175/JPO3101.1.
1132 Bodner, A. S., B. Fox-Kemper, L. Johnson, L. P. V. Roekel, J. C. McWilliams, P. P. Sullivan, P. S. Hall, and J. Dong, 2023: Modifying the mixed layer eddy parameterization to include frontogenesis arrest by boundary layer turbulence. *Journal of Physical Oceanography*, **53 (1)**, 323–339, https://doi.org/10.1175/JPO-D-21-0297.1.

1136 Bodner, A. S., B. Fox-Kemper, L. P. VanRoekel, J. C. McWilliams, and P. P. Sullivan, 2019: A perturbation approach to understanding the effects of turbulence on frontogenesis. *Journal of Fluid Mechanics*, **883**, A25, https://doi.org/10.1017/jfm.2019.804.

1139 Buckingham, C. E., N. S. Lucas, S. E. Belcher, T. P. Rippeth, A. L. M. Grant, J. Le Sommer, A. O. Ajayi, and A. C. Naveira Garabato, 2019: The contribution of surface and submesoscale processes to turbulence in the open ocean surface boundary layer. *Journal of Advances in Modeling Earth Systems*, **11 (12)**, 4066–4094, https://doi.org/10.1029/2019MS001801.

1143 Capet, X., J. McWilliams, M. Molemaker, and A. Shchepetkin, 2008: Mesoscale to Submesoscale Transition in the California Current System. Part I: Flow Structure, Eddy Flux, and Observational Tests. *J. Phys. Oceanogr.*, **38**, 29–43, https://doi.org/10.1175/2007JPO3671.1.

1146 Carpenter, J. R., A. Rodrigues, L. K. P. Schultze, L. M. Merckelbach, N. Suzuki, B. Baschek, and L. Umlauf, 2020: Shear instability and turbulence within a submesoscale front following a storm. *Geophysical Research Letters*, **47 (23)**, e2020GL090365, https://doi.org/10.1029/2020GL090365.

1150 Chor, T., J. O. Wenegrat, and J. Taylor, 2022: Insights into the mixing efficiency of submesoscale centrifugal-symmetric instabilities. *Journal of Physical Oceanography*, **52 (10)**, 2273–2287, https://doi.org/10.1175/JPO-D-21-0259.1.

1153 Contreras, M., L. Renault, and P. Marchesiello, 2023: Understanding energy pathways in the gulf stream. *Journal of Physical Oceanography*, **53 (3)**, 719–736, https://doi.org/10.1175/JPO-D-22-0146.1.

1156 Crowe, M. N., and J. R. Taylor, 2018: The evolution of a front in turbulent thermal wind balance. part 1. theory. *Journal of Fluid Mechanics*, **850**, 179–211, https://doi.org/10.1017/jfm.2018.448.
1158

Crowe, M. N., and J. R. Taylor, 2019: The evolution of a front in turbulent thermal wind balance.

1159

part 2. numerical simulations. *Journal of Fluid Mechanics*, **880**, 326–352, https://doi.org/10.

1160

1017/jfm.2019.688.

1161

Crowe, M. N., and J. R. Taylor, 2020: The effects of surface wind stress and buoyancy flux

1162

on the evolution of a front in a turbulent thermal wind balance. *Fluids*, **5 (2)**, https://doi.org/

1163

10.3390/fluids5020087.

1164

Damien, P., D. Bianchi, J. C. McWilliams, F. Kessouri, C. Deutsch, R. Chen, and L. Renault,

1165

2023: Enhanced biogeochemical cycling along the u.s. west coast shelf. *Global Biogeochemical*

1166

Cycles, **37 (1)**, e2022GB007572, https://doi.org/https://doi.org/10.1029/2022GB007572.

1167

D'Asaro, E., C. Lee, L. Rainville, R. Harcourt, and L. Thomas, 2011: Enhanced turbulence

1168

and energy dissipation at ocean fronts. *Science*, **332**, 318–322, https://doi.org/10.1126/science.

1169

1201515.

1170

Dauhajre, D. P., and J. McWilliams, 2018: Diurnal evolution of submesoscale front and filament

1171

circulations. *J. Phys. Oceanogr.*, **48**, 2343–2361, https://doi.org/10.1175/JPO-D-18-0143.1.

1172

Dauhajre, D. P., J. McWilliams, and Y. Uchiyama, 2017: Submesoscale Coherent Structures on the

1173

Continental Shelf. *J. Phys. Oceanogr.*, **47**, 2949–2976, https://doi.org/10.1175/JPO-D-16-0270.

1174

1.

1175

Dong, J., B. Fox-Kemper, H. Zhang, and C. Dong, 2021: The scale and activity of symmetric

1176

instability estimated from a global submesoscale-permitting ocean model. *Journal of Physical*

1177

*Oceanography*, **51 (5)**, 1655–1670, https://doi.org/10.1175/JPO-D-20-0159.1.

1178

D'Asaro, E. A., and Coauthors, 2018: Ocean convergence and the dispersion of flotsam.

1179

*Proceedings of the National Academy of Sciences*, **115 (6)**, 1162–1167, https://doi.org/

1180

10.1073/pnas.1718453115.

1181

Fahlbusch, J. A., D. E. Cade, E. L. Hazen, M. L. Elliott, B. T. Saenz, J. A. Goldbogen, and

1182

J. Jahncke, 2024: Submesoscale coupling of krill and whales revealed by aggregative lagrangian

1183

coherent structures. *Proc. R. Soci. B.*, **291 (20232461)**, https://doi.org/10.1098/rspb.2023.2461.
1184 Farrar, J. T., and Coauthors, 2020: S-mode: The sub-mesoscale ocean dynamics experiment. *IGARSS 2020 - 2020 IEEE International Geoscience and Remote Sensing Symposium*, 3533–3536, https://doi.org/10.1109/IGARSS39084.2020.9323112.

1187 Fox-Kemper, B., R. Ferrari, and R. Hallberg, 2008: Parameterization of mixed layer eddies. part i: Theory and diagnosis. *Journal of Physical Oceanography*, **38** (6), 1145 – 1165, https://doi.org/10.1175/2007JPO3792.1.

1190 Fox-Kemper, B., and Coauthors, 2011: Parameterization of mixed layer eddies. iii: Implementation and impact in global ocean climate simulations. *Ocean Modelling*, **39** (1), 61–78, https://doi.org/10.1016/j.ocemod.2010.09.002.

1193 Freilich, M. A., G. Flierl, and A. Mahadevan, 2022: Diversity of growth rates maximizes phytoplankton productivity in an eddying ocean. *Geophysical Research Letters*, **49** (3), e2021GL096180, https://doi.org/10.1029/2021GL096180.

1196 Garabato, A. C. N., X. Yu, J. Callies, R. Barkan, K. L. Polzin, E. E. Frajka-Williams, C. E. Buckingham, and S. M. Griffies, 2022: Kinetic energy transfers between mesoscale and submesoscale motions in the open ocean’s upper layers. *Journal of Physical Oceanography*, **52** (1), 75 – 97, https://doi.org/10.1175/JPO-D-21-0099.1.

1200 Garrett, C., and J. Loder, 1981: Dynamical aspects of shallow sea fronts. *Phil. Trans. R. Soc. Lond.*, **302**, 563–581, https://doi.org/10.1098/rsta.1981.0183.

1202 Gula, J., M. Molemaker, and J. McWilliams, 2014: Submesoscale cold filaments in the Gulf Stream. *J. Phys. Oceanogr.*, **44**, 2617–2643, https://doi.org/10.1175/JPO-D-14-0029.1.

1204 Gula, J., J. Taylor, A. Shcherbina, and A. Mahadevan, 2022: Chapter 8 - Submesoscale processes and mixing. *Ocean Mixing*, M. Meredith, and A. Naveira Garabato, Eds., Elsevier, 181–214, https://doi.org/10.1016/B978-0-12-821512-8.00015-3.

1207 Haney, S., B. Fox-Kemper, K. Julien, and A. Webb, 2015: Symmetric and geostrophic instabilities in the wave-forced ocean mixed layer. *Journal of Physical Oceanography*, **45** (12), 3033 – 3056, https://doi.org/10.1175/JPO-D-15-0044.1.
Heaps, N. S., 1972: Estimation of Density Currents in the Liverpool Bay Area of the Irish Sea. *Geophysical Journal International*, **30 (4)**, 415–432, https://doi.org/10.1111/j.1365-246X.1972.tb05825.x.

Hoskins, B. J., 1974: The role of potential vorticity in symmetric stability and instability. *Quarterly Journal of the Royal Meteorological Society*, **100 (425)**, 480–482, https://doi.org/10.1002/qj.49710042520.

Hoskins, B. J., 1982: The mathematical theory of frontogenesis. *Ann. Rev. Fluid Mech.*, **14**, 131–151, https://doi.org/10.1146/annurev.fl.14.010182.001023.

Hoskins, B. J., and F. P. Bretherton, 1972: Atmospheric frontogenesis models: Mathematical formulation and solution. *Journal of Atmospheric Sciences*, **29 (1)**, 11–37, https://doi.org/10.1175/1520-0469(1972)029<0011:AFMMFA>2.0.CO;2.

Johnson, L., and B. Fox-Kemper, 2024: Modification of boundary layer turbulence by submesoscale flows. *Flow*, **4**, E20, https://doi.org/10.1017/flo.2024.17.

Johnson, L., C. M. Lee, E. A. D'Asaro, L. Thomas, and A. Shcherbina, 2020a: Restratification at a california current upwelling front. part i: Observations. *Journal of Physical Oceanography*, **50 (5)**, 1455–1472, https://doi.org/10.1175/JPO-D-19-0203.1.

Johnson, L., C. M. Lee, E. A. D'Asaro, J. O. Wenegrat, and L. N. Thomas, 2020b: Restratification at a california current upwelling front. part ii: Dynamics. *Journal of Physical Oceanography*, **50 (5)**, 1473–1487, https://doi.org/10.1175/JPO-D-19-0204.1.

Johnston, T. M. S., D. L. Rudnick, and E. Pallàs-Sanz, 2011: Elevated mixing at a front. *Journal of Geophysical Research: Oceans*, **116 (C11)**, https://doi.org/10.1029/2011JC007192.

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

Lemarie, F., J. Kurian, A. Shchepetkin, M. Molemaker, F. Colas, and J. McWilliams, 2012: Are there inescapable issues prohibiting the use of terrain-following coordinates in climate models? *Ocean Modelling*, **42**, 57–79, https://doi.org/10.1016/j.ocemod.2011.11.007.
1238

Lentz, S. J., 2022: Turbulent thermal-wind-driven coastal upwelling: Current observations and dynamics. *Journal of Physical Oceanography*, **52 (12)**, 2909 – 2921, https://doi.org/10.1175/JPO-D-22-0063.1.

1239

1240

1241

Levy, M., R. Ferrari, P. Franks, A. Martin, and P. Riviere, 2012: Bringing physics to life at the submesoscale. *Geophys. Res. Lett.*, **39**, 1–13, https://doi.org/10.1029/2012GL052756.

1242

1243

Li, D., 2019: Turbulent prandtl number in the atmospheric boundary layer - where are we now? *Atmospheric Research*, **216**, 86–105, https://doi.org/10.1016/j.atmosres.2018.09.015.

1244

1245

1246

Lévy, M., P. Franks, and K. Smith, 2018: The role of submesoscale currents in structuring marine ecosystems. *Nature Communications*, **9 (4758)**, https://doi.org/10.1038/s41467-018-07059-3.

1247

1248

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

1249

1250

McWilliams, J., 2017: Submesoscale surface fronts and filaments: secondary circulation, buoyancy flux, and frontogenesis. *J. Fluid Mech.*, **823**, 391 – 432, https://doi.org/10.1017/jfm.2017.294.

1251

1252

McWilliams, J., 2021: Oceanic frontogenesis. *Ann. Rev. Mar. Sci., in press.*, **13**, 227–53, https://doi.org/10.1146/annurev-marine-032320-120725.

1253

1254

McWilliams, J., J. Gula, J. Molemaker, L. Renault, and A. Shchepetkin, 2015: Filament frontogenesis by boundary layer turbulence. *J. Phys. Oceanogr.*, **45**, 1988–2005, https://doi.org/10.1175/JPO-D-14-0211.1.

1255

1256

1257

Molemaker, M. J., J. C. McWilliams, and X. Capet, 2010: Balanced and unbalanced routes to dissipation in an equilibrated eady flow. *Journal of Fluid Mechanics*, **654**, 35–63, https://doi.org/10.1017/S0022112009993272.

1258

1259

1260

Nagai, T., A. Tandon, and D. Rudnick, 2006: Two-dimensional ageostrophic secondary circulations at ocean fronts due to vertical mixing and large-scale deformation. *J. Geophys. Res.*, **111**, https://doi.org/10.1029/2005JC002964.

1261

1262
1263

Peng, J.-P., J. Dräger-Dietel, R. P. North, and L. Umlauf, 2021: Diurnal variability of frontal dynamics, instability, and turbulence in a submesoscale upwelling filament. *Journal of Physical Oceanography*, **51 (9)**, 2825–2843, https://doi.org/https://doi.org/10.1175/JPO-D-21-0033.1.

1264

1265

1266

Qu, L., and Coauthors, 2022: Rapid vertical exchange at fronts in the northern gulf of mexico. *Nature Communications*, **13 (5624)**, 1455–1472, https://doi.org/10.1038/s41467-022-33251-7.

1267

1268

Shakespeare, C., and J. Taylor, 2013: A generalized mathematical model of geostrophic adjustment and frontogenesis: uniform potential vorticity. *J. Fluid Mech.*, **736**, 366 – 413, https://doi.org/10.1017/jfm.2013.526.

1269

1270

1271

Shchepetkin, A., and J. McWilliams, 2005: The Regional Oceanic Modeling System (ROMS): A split-explicit, free-surface, topography-following-coordinate ocean model. *Ocean Modell.*, **9**, 347–404, https://doi.org/10.1016/j.ocemod.2004.08.002.

1272

1273

1274

Siegelman, L., P. Klein, P. Rivière, A. F. Thompson, H. S. Torres, M. Flexas, and D. Menemenlis, 2020: Enhanced upward heat transport at deep submesoscale ocean fronts. *Nature Geoscience*, **13**, 50 – 55, https://doi.org/doi.org/10.1038/s41561-019-0489-1.

1275

1276

1277

Srinivasan, K., R. Barkan, and J. C. McWilliams, 2023: A forward energy flux at submesoscales driven by frontogenesis. *Journal of Physical Oceanography*, **53 (1)**, 287 – 305, https://doi.org/10.1175/JPO-D-22-0001.1.

1278

1279

1280

Su, Z., J. Wang, P. Klein, A. F. Thompson, and D. Menemenlis, 2018: Ocean submesoscales as a key component of the global heat budget. *Nature Communication*, **9 (775)**, https://doi.org/10.1038/s41467-018-02983-w.

1281

1282

1283

Sullivan, P., and J. McWilliams, 2017: Frontogenesis and frontal arrest for a dense filament in the oceanic surface boundary layer. *J. Fluid Mech.*, **837**, 341 – 380, https://doi.org/10.1017/jfm.2017.833.

1284

1285

1286

Sullivan, P. P., and J. C. McWilliams, 2024: Oceanic frontal turbulence. *Journal of Physical Oceanography*, **54 (2)**, 333 – 358, https://doi.org/10.1175/JPO-D-23-0033.1.

1287

1288

Sun, D., and Coauthors, 2020: Diurnal cycling of submesoscale dynamics: Lagrangian implications in drifter observations and model simulations of the northern gulf of mexico. *Journal of Physical Oceanography*, **50 (6)**, 1605 – 1623, https://doi.org/10.1175/JPO-D-19-0241.1.

1289

1290
1291 Swart, S., M. D. du Plessis, A. F. Thompson, L. C. Biddle, I. Giddy, T. Linders, M. Mohrmann, and S.-A. Nicholson, 2020: Submesoscale fronts in the antarctic marginal ice zone and their response to wind forcing. *Geophysical Research Letters*, **47** (6), e2019GL086649, https://doi.org/10.1029/2019GL086649.

Taylor, J. R., 2016: Turbulent mixing, restratification, and phytoplankton growth at a submesoscale eddy. *Geophysical Research Letters*, **43** (11), 5784–5792, https://doi.org/10.1002/2016GL069106.

Taylor, J. R., and A. F. Thompson, 2023: Submesoscale dynamics in the upper ocean. *Annual Review of Fluid Mechanics*, **55** (1), 103–127, https://doi.org/10.1146/annurev-fluid-031422-095147.

Thomas, L., and R. Ferrari, 2008: Friction, frontogenesis, and the stratification of the surface mixed layer. *Journal of Physical Oceanography*, **38** (11), 2501–2518, https://doi.org/10.1175/2008JPO3797.1.

Thomas, L., A. Tandon, and A. Mahadevan, 2008: Submesoscale Processes and Dynamics. *Ocean Modeling in an Eddying Regime, Geophysical Monograph Series*, **177**, 17–38, https://doi.org/10.1029/177GM04.

Thomas, L. N., J. R. Taylor, E. A. D'Asaro, C. M. Lee, J. M. Klymak, and A. Shcherbina, 2016: Symmetric instability, inertial oscillations, and turbulence at the gulf stream front. *Journal of Physical Oceanography*, **46** (1), 197–217, https://doi.org/10.1175/JPO-D-15-0008.1.

Thomas, L. N., J. R. Taylor, R. Ferrari, and T. M. Joyce, 2013: Symmetric instability in the gulf stream. *Deep Sea Research Part II: Topical Studies in Oceanography*, **91**, 96–110, https://doi.org/10.1016/j.dsr2.2013.02.025.

Thompson, L., 2000: Ekman layers and two-dimensional frontogenesis in the upper ocean. *Journal of Geophysical Research: Oceans*, **105** (C3), 6437–6451, https://doi.org/10.1029/1999JC900336.

Verma, V., H. T. Pham, and S. Sarkar, 2019: The submesoscale, the finescale and their interaction at a mixed layer front. *Ocean Modelling*, **140**, 101–100, https://doi.org/10.1016/j.ocemod.2019.05.004.
1319 Wang, T., R. Barkan, J. McWilliams, and M. J. Molemaker, 2021: Structure of submesoscale fronts of the mississippi river plume. *J. Phys. Oceanogr.*, **51**, 1113–1131, https://doi.org/10.1175/JPO-D-20-0191.1.

1322 Wenegrat, J., and M. McPhaden, 2016: A simple analytical model of the diurnal Ekman layer. *J. Phys. Oceanogr.*, **46**, 2877–2894, https://doi.org/10.1175/JPO-D-16-0031.1n.

1324 Wenegrat, J. O., L. N. Thomas, M. A. Sundermeyer, J. R. Taylor, E. A. D'Asaro, J. M. Klymak, R. K. Shearman, and C. M. Lee, 2020: Enhanced mixing across the gyre boundary at the gulf stream front. *Proceedings of the National Academy of Sciences*, **117 (30)**, 17607–17614, https://doi.org/10.1073/pnas.2005558117.

1328 Wu, K., H. Cao, and G. Liao, 2022: Submesoscale frontal waves and instabilities driven by sheared flows. *Deep Sea Research Part II: Topical Studies in Oceanography*, **202**, 105145, https://doi.org/10.1016/j.dsr2.2022.105145.

1331 Yang, P., Z. Jing, H. Yang, and L. Wu, 2024: A scale-aware parameterization of restratification effect of turbulent thermal wind balance. *Journal of Physical Oceanography*, https://doi.org/10.1175/JPO-D-23-0169.1.

1334 Young, W. R., 1994: The subinertial mixed layer approximation. *Journal of Physical Oceanography*, **24 (8)**, 1812–1826, https://doi.org/10.1175/1520-0485(1994)024(1812:TSMLA)2.0.CO; 2.

1337 Young, W. R., and S. Jones, 1991: Shear dispersion. *Physics of Fluids A: Fluid Dynamics*, **3 (5)**, 1087–1101, https://doi.org/10.1063/1.858090.

1339 Yu, X., R. Barkan, and A. N. Garabato, 2024: Intensification of submesoscale frontogenesis and forward energy cascade driven by upper-ocean convergent flows. *Nat. Commun.*, **15 (9214)**, https://doi.org/10.1038/s41467-024-53551-4.

1342 Yu, X., A. C. Naveira Garabato, A. P. Martin, D. Gwyn Evans, and Z. Su, 2019: Wind-forced symmetric instability at a transient mid-ocean front. *Geophysical Research Letters*, **46 (20)**, 11281–11291, https://doi.org/10.1029/2019GL084309.
Zhang, J., Z. Zhang, and B. Qiu, 2023: Parameterizing submesoscale vertical buoyancy flux by simultaneously considering baroclinic instability and strain-induced frontogenesis. *Geo-physical Research Letters*, **50** (8), e2022GL102292, https://doi.org/10.1029/2022GL102292.

Zhang, Z., X. Zhang, B. Qiu, W. Zhao, C. Zhou, X. Huang, and J. Tian, 2021: Submesoscale currents in the subtropical upper ocean observed by long-term high-resolution mooring arrays. *Journal of Physical Oceanography*, **51** (1), 187–206, https://doi.org/10.1175/JPO-D-20-0100.1.