Modification of boundary layer turbulence by submesoscale flows

Abstract

The stirring and mixing of heat and momentum in the ocean surface boundary layer (OSBL) are dominated by 1 to 10 km fluid flows – too small to be resolved in global and regional ocean models. Instead, these processes are parametrized. Two main parametrizations include vertical mixing by surface-forced metre-scale turbulence and overturning by kilometre-scale submesoscale frontal flows and instabilities. In present models, these distinct parametrizations are implemented in tandem, yet ignore meaningful interactions between these two scales that may influence net turbulent fluxes. Using a large-eddy simulation of frontal spin down resolving processes at both scales, this work diagnoses submesoscale and surface-forced turbulence impacts that are the foundation of OSBL parametrizations, following a traditional understanding of these flows. It is shown that frontal circulations act to suppress the vertical buoyancy flux by surface forced turbulence, and that this suppression is not represented by traditional boundary layer turbulence theory. A main result of this work is that current OSBL parametrizations excessively mix buoyancy and overestimate turbulence dissipation rates in the presence of lateral flows. These interactions have a direct influence on the upper ocean potential vorticity and energy budgets with implications for global upper ocean budgets and circulation.

Impact Statement

The ocean surface layer contains fluid flows that play an essential role in the communication between the atmosphere and ocean. Two small-scale processes that are routinely approximated in global models include turbulent mixing due to atmospheric forcing and turbulent circulations from small-scale currents. Interactions between surface-forced and current-forced turbulence indicate that our modelling approaches need updating.

1. Introduction

Turbulence in the ocean surface boundary layer (OSBL) involves a multitude of dynamics that drive fluid motions on a range of scales. This turbulence defines the boundary layer and therefore the upper ocean's role in transferring heat and momentum between the atmosphere and the ocean. Surface fluxes can drive fine-scale turbulence ($O(1\,{\rm cm})$–$O(100\,{\rm m})$) in the OSBL as energy injection from winds, waves and cooling that compete with stratification by solar warming and freshwater fluxes (Belcher et al. 2012; Fox-Kemper, Johnson & Qiao 2021). In the presence of strong lateral density variability, submesoscale fronts ($O(50\,{\rm m})$–$O(2\,{\rm km})$; Bodner et al. 2023) and baroclinic mixed-layer instabilities ($O(500\,{\rm m})$–$O(50\,{\rm km})$; Dong et al. 2020) can restratify the upper ocean as dynamics release available potential energy to convert horizontal gradients into vertical ones (e.g. Fox-Kemper, Ferrari & Hallberg 2008; Johnson et al. 2020). Additionally, a range of instabilities associated with $O(1)$ Rossby number submesoscale flows are known to drive energy in shear flows towards dissipative scales (Taylor & Ferrari 2009; Thomas et al. 2013) producing turbulence comparable in length scale ($O(20\,{\rm m})$–$O(500\,{\rm m})$; Dong et al. 2021) but distinct from surface-forced turbulence. The dynamic range of submesoscale and turbulent motions cannot be captured by global, regional and even many submesoscale-permitting process models and therefore the unresolved motions are parametrized to a varying extent. The ability to parametrize transport by unresolved physics rests on averaging over the field of eddies, whether it be surface-forced boundary layer turbulent eddies (e.g. Troen & Mahrt 1986; Large, McWilliams & Doney 1994; Umlauf & Burchard 2003; Reichl & Hallberg 2018), or larger, flatter eddies associated with baroclinic instabilities of mesoscale or submesoscale lateral density gradients (Gent & Mcwilliams 1990; Fox-Kemper & Ferrari 2008; Bachman, Fox-Kemper & Bryan 2020; Bodner & Fox-Kemper 2020). These parametrizations have an important impact on how the upper ocean is represented in global models (Fox-Kemper et al. 2011; Li et al. 2016; Fox-Kemper et al. 2019; Dong et al. 2020).

The theory of each parametrization targets a specific dynamical regime and isolates the time and length scales over which to average the processes to arrive at meaningful transport estimates. It is routine for global and regional circulation models to include separate parametrizations for boundary layer eddies and submesoscale eddies. Yet, implementing these individual processes in tandem typically ignores cross-scale interactions such as suppression of one scale of eddy by the eddies of the other scale. Note that this eddy–eddy effect is different from the effect that a parametrization of one scale may have on the eddies of the other scale. For example, a parametrization of submesoscale restratification (Fox-Kemper et al. 2008) will reduce boundary layer depth and thus also reduce boundary layer eddy mixing because mixing strength scales with boundary layer depth (Monin & Obukhov 1954). In the other direction, more surface forcing of a boundary layer parametrization will deepen the mixed layer and make submesoscale instabilities appear at larger scales. Here we focus on the eddy–eddy interaction in a multiscale simulation, not the parametrization–eddy interaction that could be carried out with a single scale resolved and the other parametrized. This work addresses how the foundational approximations for surface forced boundary layer turbulence and submesoscale flows are modified when the two co-exist and interact in the OSBL.

Understanding nonlinear interactions between different eddying flows is a challenge due to the computational resources required for a large domain and small-grid-scale simulations. A computationally approachable way to explore the role of boundary layer turbulence on submesoscale motions is through simulations that resolve submesoscale instabilities, yet parametrize boundary layer mixing (e.g. Mahadevan, Tandon & Ferrari 2010; Wyngaard 2010; Callies & Ferrari 2018; Dauhajre & McWilliams 2018; Wenegrat et al. 2020), and has revealed how the evolution of submesoscale fronts is intimately linked to the turbulent environment they exist in. Yet these works are limited in their ability to understand how submesoscale flows impact fine-scale turbulence not resolved by such simulations. To do this, studies employ large-eddy simulation (LES), with small enough grid scales $(O(1\,{\rm m}))$ to resolve boundary layer eddies, yet large enough to simulate baroclinic waves and frontal instabilities $(O(10\,{\rm km}))$. These simulations have been used to investigate the fine-scale motions associated with submesoscale flows (Taylor & Ferrari 2009; Skyllingstad & Samelson 2012; Verma, Pham & Sarkar 2019), or how submesoscale fronts interact with wind-driven (Hamlington et al. 2014; Skyllingstad, Duncombe & Samelson 2017; Whitt, Lévy & Taylor 2019) or convective (Taylor 2016; Skyllingstad et al. 2017; Taylor, Smith & Vreugdenhil 2020; Verma, Pham & Sarkar 2022) turbulence. Many of these studies identify the role of turbulence in enhancing or driving submesoscale circulation and the transfer of energy across spatial scales. Under a different lens, this work will revisit simulations presented in Hamlington et al. (2014) to understand how the interactions between fine-scale and submesoscale flows modify current idealized frameworks and state-of-the-art parametrizations that treat them separately. It will be shown that OSBL parametrizations can overestimate turbulence in the presence of submesoscale flows.

Section 2 overviews the LES and data used in the paper. Section 3 presents the multi-level Reynolds decomposition to separate different scales of motion. Section 4 reviews the buoyancy budget and evaluates how well different scales of transport relate to current scalings. These results are expanded to explore the impact of approximating turbulent flows using current state-of-the-art turbulence parametrizations on the upper ocean potential vorticity (PV) budget in § 5, and the impact on the dissipation rate of kinetic energy in § 6. Section 7 discusses these results in the context of current work along with implications for steps forward.

2. Simulations and data summary

This work utilizes an LES of a frontal spin down under uniform wind stress and without waves (or Stokes drift) as described in Hamlington et al. (2014). Other studies have revisited passive tracers mimicking biological tracer transport (Smith, Hamlington & Fox-Kemper 2016), frontogenesis (Suzuki & Fox-Kemper 2016) and PV spectra (Bodner & Fox-Kemper 2020) in these runs; this analysis focuses on the standard diagnostic practices distinguishing turbulent ‘mixing’ from submesoscale ‘overturning’ underpinning parametrizations affecting buoyancy, energy and PV.

The horizontal doubly periodic domain is $20\,{\rm km}\times 20\,{\rm km}$ with 4.9 m resolution, and the vertical extent in 160 m with 1.25 m resolution. The simulation contains a warm filament with a uniform surface wind stress of $\tau = 0.025\,{\rm N}\,{\rm m}^{-2}$ aligned at a $30^{\circ }$ angle from the geostrophic flow. The analysis includes data from days 10 and 11 of the simulation. It is convenient to divide the simulation spatially into three regions: (1) NOFRONT, (2) STABLE and (3) UNSTABLE (figure 1). The NOFRONT region represents the neutral boundary layer, where turbulent motions are driven by surface wind stress and deliver modest entrainment. In the STABLE region, the constant and uniform, partly upfront, wind stress induces an Ekman transport to the right of the wind direction. This Ekman flow transports warm water over the cold side of the front (upfront wind/stable), thereby stratifying the boundary layer. In the UNSTABLE region, Ekman flow transports cold water over the warm side of the front, causing boundary layer convection from the Ekman buoyancy flux (i.e. downfront winds; Thomas & Lee 2005). This front develops unstable baroclinic instabilities that drive restratifying circulations within the boundary layer. The simulations were designed such that eddy restratification in this UNSTABLE region ($Q \sim 25\,{\rm W}\,{\rm m}^{-2}$) is larger than Ekman-driven convection ($Q \sim -13\,{\rm W}\,{\rm m}^{-2}$) and similar in magnitude to the vertical heat flux in the NOFRONT region from shear-driven mixing and entrainment ($Q\sim -10\,{\rm W}\,{\rm m}^{-2}$). Due to the finite-amplitude baroclinic waves, the Ekman-induced boundary layer convection impacts 20 % of the UNSTABLE region. Realistic forcing conditions might strengthen any one mechanism to dominate over the others (Haney et al. 2012). Note that the UNSTABLE region may also possess forced symmetric instabilities (Bachman et al. 2017). For the time window of data analysed, the STABLE and UNSTABLE regions are more stratified than the NOFRONT region, because of Ekman and mixed-layer instability restratification, respectively. On average, mixed-layer depths are about 50, 45 and 40 m for the NOFRONT, UNSTABLE and STABLE regions, respectively. The fully developed eddy field and restratifying front at day 10 can be seen in figure 1.

Figure 1. Surface temperature at day 10 of the H14 runs. (a) The entire model domain is separated into regimes: homogeneous turbulence region (NO FRONT, orange), front with stabilizing winds (STABLE, green), front with unstable winds and baroclinic instability (UNSTABLE,blue). (b) A $2\,{km}\times 2\,{km}$ close-up of the UNSTABLE region. The dashed area is a $400\,{m}\times 400\,{m}$ box, where 400 m length scale is the transition between quasi-2-D and 3-D flows (see H14). The colourbar is the same for both panels.

Fluid motions in each of these domains are influenced by the dynamics on all scales, with velocity spectral slopes consistent with two-dimensional (2-D) circulation involving frontal velocity jumps at larger scales ($k^{-2}$ slope), transitioning to three-dimensional (3-D) turbulence at ${\sim }400\,{\rm m}$ (Hamlington et al. 2014). Though velocity spectra support 3-D turbulence below 400 m, a detailed look at surface temperature (figure 1) reveals rich horizontal structure at the $O({\rm m})$ scale, thus challenging the assumption of horizontal homogeneity that often sets the prerequisite for OSBL turbulence physics (Monin & Obukhov 1954). The goal is to evaluate how these fine-scale flows compare with buoyancy fluxes approximated by current OSBL parametrizations. To do this, three different OSBL parametrizations are implemented using the General Ocean Turbulence Modelling (GOTM; Burchard, Bolding & Villarreal 1999) framework (hereafter referred to as one-dimensional (1-D)) as summarized in table 1.

Table 1. List of parametrizations used in this study.

The ingredients of extant submesoscale parametrizations (Fox-Kemper et al. 2008; Bachman et al. 2017) are diagnosed by scale separation.

3. Separation of scales

Turbulent flows are often expressed in terms of their statistical mean through the Reynolds-averaged equations, whereby instantaneous values are represented as the sum of an average value (in space or time) and fluctuations from the mean. The interpretation of the turbulence then is subject to the resolution of the fluctuations and definition of the mean fields. Here, a multi-level Reynolds decomposition is utilized to separate submesoscale turbulence from fine-scale turbulence. The first decomposition for any variable $c$ is ${c} = \tilde {c} + c'$, where $\tilde {c}$ represents a square boxcar average over a to-be-determined length scale to be permitted as submesoscale, and $c'$ are the finer-than-submesoscale fluctuations from that mean. The submesoscale can also be averaged further along the front indicated by an overbar, $\bar {c}$. Thus, $\tilde {c} = {c^a} + c^{s}$, where ${c^a}$ is the domain average of $\tilde {c}$ in the along-front direction (hereafter referred to as along-front average). Here $c^{s}$ are the submesoscale fluctuations from the along-front mean. Combined, the decomposition becomes

(3.1a–c)\begin{equation} c = \tilde{c} + c' = c^a + c^s + c',\quad \tilde{c} = c^a + c^s,\quad \overline{\tilde{c}} = c^a. \end{equation}

Note that $c'$ is not the fluctuation term in a traditional along-front Reynolds decomposition, $c - \bar {c} \ne c'$; rather, it is the fluctuation from a submesoscale-permitting resolution as in (3.1a–c) – similar to the role that a boundary layer mixing parametrization would play in a submesoscale-permitting model (e.g. Gula, Molemaker & McWilliams 2014; Su et al. 2018). On the other hand, the submesoscale-permitting mean, $\tilde {c}$, may be analogous to a coarse resolution regional or process model that resolves submesoscale turbulence but not fine-scale motions, which is parametrized. Making this scale separation specific then rests on the definition of the minimum submesoscale length permitted in $\tilde {c}$. A sensible definition is the transition scale from 2-D to 3-D turbulence scaling. Our choice is 396.9 m to accommodate grid resolution and is based on the 400 m flow transition for this simulation described in Hamlington et al. (2014) (horizontal kinetic energy spectral slope) and Bodner & Fox-Kemper (2020) (PV spectral slope).

The along-front and submesoscale means are both linear operators, and they are furthermore idempotent Reynolds averages, which distinguishes them from more generic ‘linear filters’ more commonly used in LES (e.g. Fox-Kemper & Menemenlis 2008; Storer et al. 2022). Differentiation and integration are also linear operators, but they do not commute with averaging generally. A further discussion of the specific operators used, their limitations and estimates of error resulting from these choices is in Appendix A.

The multi-level Reynolds decomposition for a Boussinesq fluid is thus

(3.2)\begin{equation} \frac{\partial }{\partial t}(b^a + b^s + b') = {-}\frac{\partial}{\partial x_j}[(u_j^a + u_j^s + u'_j)(b^a + b^s + b')] -\mathcal{V}_{ij}^{b}, \end{equation}

(3.3)\begin{align} \frac{\partial }{\partial t}(u_i^a + u_i^s + u'_i) & = {-}\frac{\partial}{\partial x_j}[(u_j^a + u_j^s + u'_j) (u_i^a + u_i^s + u'_i)] -\frac{1}{\rho_o}\frac{\partial}{\partial x_i}(p^a + p^s + p') \nonumber\\ &\quad - 2\epsilon_{ikj}\varOmega_j(u_k^a + u_k^s + u'_k) - \mathcal{V}_{ij}^{m}. \end{align}

Here, $b$ is buoyancy, $u$ is velocity, $p$ is pressure, $\varOmega$ is planetary rotation rate and $\mathcal {V}_{ij}^{b}$ and $\mathcal {V}_{ij}^{m}$ are the divergence of the diffusive buoyancy and viscous momentum fluxes, respectively. Operating on (3.2) and (3.3), the submesoscale Reynolds average, $\widetilde {()}$, is taken first (noting that all terms linear in primes vanish), followed by the along-front average, $\overline {()}$ (noting that all terms linear in submesoscale variables vanish), yielding the multi-level Reynolds-averaged equations for along-front mean of buoyancy and momentum:

(3.4)\begin{equation} \frac{\partial b^{a}}{\partial t} = {-}\frac{\partial}{\partial x_j}[u^{a}_j b^{a}] - \frac{\partial}{\partial x_j} \overline{\widetilde{b^{s} u^{s}_j}} - \frac{\partial}{\partial x_j} \overline{\widetilde{b' u'_j}} - \overline{\widetilde{\mathcal{V}_{ij}^{b}}} \end{equation}

and

(3.5)\begin{equation} \frac{\partial u^{a}_i}{\partial t} = {-}\frac{\partial }{\partial x_j}[u^{a}_j u^{a}_i] -\frac{1}{\rho_o}\frac{\partial p^{a}}{\partial x_i} - 2\epsilon_{ikj}\varOmega_ju^{a}_k + \frac{\partial}{\partial x_j}\overline{\widetilde{u^{s}_i u^{s}_j}} - \frac{\partial}{\partial x_j} \overline{\widetilde{u'_i u'_j}} - \overline{\widetilde{\mathcal{V}_{ij}^{m}}}. \end{equation}

The covariance terms capture the transport of properties at the different scales, separating clearly the transport by submesoscale turbulence from transport due to fine-scale motions.

This Reynolds-averaging approach is more similar to the scale-separation assumption typical of ocean climate models (i.e. scale-separated Reynolds averaging) rather than cascade-filtering approaches more typical of LES studies (see e.g. Fox-Kemper & Menemenlis 2008; Aluie, Hecht & Vallis 2018). The reason for the choice of Reynolds averaging is to explore whether the approach taken in ocean climate modelling is deficient when considering the turbulence of both the submesoscale (e.g. Fox-Kemper et al. 2008) and the fine scale (e.g. Large et al. 1994), potentially including resolving the submesoscale but not the fine scale (e.g. Gula et al. 2014; Su et al. 2018). Note that these scales are neighbouring, but with spectrally distinct slopes in this simulation and recognizable as quasi-2-D and 3-D cascades (Hamlington et al. 2014; Bodner & Fox-Kemper 2020). Equation (3.4) is similar to the triple decomposition used with observations to evaluate along isopycnal tracer dispersion in the North Atlantic (Ferrari & Polzin 2005; Smith & Ferrari 2009) and Southern Ocean (Garabato et al. 2016) stratified ocean interior, but here the decomposition is used in the context of the OSBL. This Reynolds decomposition will be used to evaluate the role of multi-scale transport in the upper ocean buoyancy budget in a form that lends itself to being compared with current subgridscale parametrizations. The Reynolds fluxes that arise from the nonlinearity of the advective term transport and stir tracers. These processes are distinct from mixing that acts to homogenize tracers. Yet delineating between the two becomes convoluted when considering how they are represented by parametrizations, each of which treats these terms differently. In this paper, the term mixing is consistent with the OSBL parametrization convention, where the stirring by boundary layer eddies is conceptualized as a mixing term.

Submesoscale motions manifest in the buoyancy and energy budgets and are linked to the fluid's PV. In terms of the buoyancy budget and energetics (Fox-Kemper et al. 2008, 2011), submesoscales impact stratification and enhance the transfer of energy towards dissipative scales (Capet et al. 2008; Taylor & Ferrari 2010; Thomas & Taylor 2010; McWilliams 2016). Additionally, PV is a dynamically relevant tracer to understand and identify submesoscale processes (Bodner & Fox-Kemper 2020; Johnson et al. 2020), especially symmetric instability (Hoskins 1974; Thomas et al. 2013; Haney et al. 2015; Bachman et al. 2017). The following sections evaluate the roles of submesoscales and turbulent scales in the tendencies of upper ocean buoyancy and PV as well as in the dissipation of kinetic energy.

4. The buoyancy budget

4.1 The transformed Eulerian mean with mixing

When averaged in the along-front direction and neglecting viscous dissipation, the Reynolds-averaged equation for buoyancy (3.4) can be rewritten as

(4.1)\begin{equation} \frac{\partial b^{a}}{\partial t} = {-} \bar{v}M^2 - \bar{w}N^2 - \frac{\partial}{\partial y} \overline{b^{s} v^{s}} - \frac{\partial}{\partial z} \overline{w^{s}b^{s}} - \frac{\partial}{\partial y} \overline{v'b'} - \frac{\partial}{\partial z} \overline{w'b'}, \end{equation}

where $M^2 = \partial b^{a}/\partial y$ and $N^2 = \partial b^{a}/\partial z$ are horizontal and vertical gradients of the along-front averaged buoyancy, respectively. The mean circulation (first two terms of (4.1)) is driven mostly by the injection of horizontal vorticity by the winds, resulting in an Ekman transport to the right of the wind stress. This Ekman overturning is represented by the Eulerian mean streamfunction, $\varPsi ^{a}$, where $v^a=\partial \varPsi ^{a}/\partial z$ and $w^a=-\partial \varPsi ^{a}/\partial y$. In the stable region, Ekman overturning restratifies the front, as differential advection moves warmer water (less dense) over the cold (more dense) side of the front. In the UNSTABLE region, Ekman transport advects dense (cold) water over light (warm), resulting in convection. This destabilizing Ekman overturning is opposed by eddy-driven overturning that acts to restratify the front (third and fourth terms in (4.1)). The eddy overturning is diagnosed by the Held & Schneider (1999) streamfunction:

(4.2)\begin{equation} \varPsi^{s} = {-}\frac{\overline{w^sb^s}}{M^2}. \end{equation}

Recognizing the eddy-induced transport slope ($S$) and the isopycnal slope ($I$) as

(4.3a,b)\begin{equation} S\equiv \frac{\overline{w^sb^s}}{\overline{v^sb^s}},\quad I\equiv \frac{-M^2}{N^2}, \end{equation}

an eddy forcing term (Held & Schneider 1999) can then be defined as

(4.4)\begin{equation} T_{ * } \equiv \overline{v^sb^s}\left(1-\frac{S}{I}\right). \end{equation}

If more tracers were present (per Smith et al. 2016) and used in the analysis of the submesoscale transport, a method like that used in Bachman, Fox-Kemper & Bryan (2015) and Bachman et al. (2020) might reveal that the eddy-induced transport has more structure than just a streamfunction.

Finally, under the assumption that fine-scale mixing is dominated by 3-D turbulence that is anisotropic according to the surface orientation, the third term on the right-hand side of (3.4) can be written using k-theory as an eddy diffusivity, such that $\overline {w'b'} = -\kappa _v N^2$ and $\overline {v'b'} = -\kappa _h M^2$. This representation assumes that fine-scale transport occurs down-gradient, an assumption that is known to break down in the presence of non-local transport (e.g. Large et al. 1994).

The buoyancy budget can now be written as

(4.5) \begin{equation} \frac{\partial b}{\partial t} = \underbrace{\boldsymbol{\nabla} \boldsymbol{\cdot}(\varPsi^a \times \boldsymbol{\nabla} b^a)}_{{\textit{Ekman}}} + \underbrace{\boldsymbol{\nabla}\boldsymbol{\cdot}(\varPsi^s \times \boldsymbol{\nabla} b^a)}_{{\textit{Eddies}}} + \underbrace{\frac{\partial}{\partial z}(\kappa_v N^2)}_{{\textit{Turbulence}}} + \underbrace{\frac{\partial}{\partial y}(T_{ * } + \kappa_h M^2)}_{{\textit{Residual}}}. \end{equation}

We note that the form of these diagnosed quantities is not guaranteed to be the most meaningful but is chosen to match the parametrization forms presently in use. A non-local convective transport (Large et al. 1994) or along-isopycnal submesoscale diffusivity (Redi 1982; Bachman et al. 2015) might also be added as common parametrizations, but they are less commonly used and cannot be reliably diagnosed from the particular LES set-up.

Combining Ekman and eddy streamfunctions, $\varPsi ^{tem} = \varPsi ^{a} + \varPsi ^{s}$ (e.g. Held & Schneider 1999; Marshall & Radko 2003), transforms (4.5) into the classic transform Eulerian mean equations modified to include fine-scale mixing. The inclusion of the turbulence term is particularly important in the near-surface ocean where strong surface forcing can cause fine-scale advective fluxes that compete with submesoscale and mean transports. This work will focus on the contribution of the first three terms, Ekman, eddies and turbulence, to stratification. Overall, along-front frontogenesis is not expected to be strong in this configuration due to a lack of large-scale strain, but it is included in the Ekman term if present. Submesoscale frontogenesis (Suzuki & Fox-Kemper 2016) involves fronts that deviate from the along-channel direction and thus are likely to be strongest in the eddies term.

4.2 Scaling buoyancy fluxes

Focusing on stratification, it is common to compare the relative strength of the vertical buoyancy flux associated with the first three terms on the right-hand side of (4.5). These scalings are

(4.6) \begin{equation} \underbrace{w^{a}b^{a} \sim \frac{\tau\times\hat{k}}{\rho}M^2}_{{\textit{Ekman buoyancy flux}}}\quad \underbrace{\overline{w^{s}b^{s}} \sim C_e\frac{M^4H^2}{f}}_{{\textit{Mixed-layer eddies}}} \quad \underbrace{\overline{w'b'} \sim{-}\kappa N^2}_{{\textit{k-theory turbulence}}}. \end{equation}

The first term (Ekman buoyancy flux (EBF); Thomas & Lee 2005) represents differential advection of horizontal buoyancy gradients by Ekman transport (figure 2a,d), causing a destabilizing buoyancy flux in the UNSTABLE region and a restratifying buoyancy flux in the STABLE region. The second term is specific to the UNSTABLE region, where baroclinic instability drives a restratifying submesoscale eddy overturning (mixed-layer eddies (MLE); Fox-Kemper et al. 2008), as in figure 2(a–c). Previous work has evaluated the trade-off between destratifying Ekman transport and restratifying submesoscale eddies using the EBF and MLE scalings (e.g. Mahadevan et al. 2010; Taylor & Ferrari 2011). Finally, the last term is relevant to all regions, as wind stress injects energy into fine-scale motions, approximated in k-theory as an isotropic, down-gradient eddy diffusivity, $\kappa$ (e.g. Large et al. 1994; Umlauf & Burchard 2003; Reichl & Hallberg 2018).

Figure 2. Overturning streamfunctions. (a) Schematic representing the different components of the buoyancy budget. Blue arrows, Ekman overturning; grey dashed line, submesoscale eddy overturning; light-grey curls, wind-driven mixing. (b) The sum of the Ekman streamfunction ($\varPsi ^a$ calculated from along-front velocities) and the Held and Schnieder streamfunction ($\varPsi ^s$ in (4.2)). Note that the positive circulation in the UNSTABLE region is consistent with dominating eddy restratification. (c) Held and Schnieder streamfunction ($\varPsi ^s$), only. (d) Ekman overturning ($\varPsi ^a$), only.

Values for $\varPsi$ (equation (4.2)) and an effective mixing coefficient, $\kappa _{eff}=\overline {w'b'}/N^2$, can be calculated from the numerical simulation directly and compared with current scalings and parametrizations that also estimate the vertical structure of these effects (see figure 3). Scalings are estimated using the mean fields on days 10 and 11 (i.e. well-developed circulations), not the initial field, which had ${\sim }5$ times stronger horizontal buoyancy gradient initially. The focus here is to assess how well parametrizations for $\kappa$, derived from theory in a homogeneous boundary layer, represent fine-scale vertical transport of buoyancy, $\overline {w'b'}$. Yet it is useful to also evaluate how the other scalings agree with the mean wind-driven and submesoscale overturning.

Figure 3. Buoyancy fluxes for different regions across the front, and colours across panels follow the legend in (b). (a) Submesoscale buoyancy flux $\overline {w^sb^s}$ as in (4.2). The only notable flux is in the UNSTABLE region and agrees with $\varPsi _{MLE}$ of Fox-Kemper et al. (2008). (b) Turbulent buoyancy flux across regions. Negative $\overline {w'b'}$ for the NOFRONT and STABLE regions is consistent with up-gradient transport and positive diffusivity in figure 4(a). The UNSTABLE region has a net positive $\overline {w'b'}$ at depths less than 15 m in the presence of stable stratification, consistent with a counter-gradient flux or frontal overturning. Averaging only in regions where lateral stratification is particularly strong and UNSTABLE (purple) or where PV is negative (magenta) increases the degree of counter-gradient flux. (c) Mixed layer averaged buoyancy fluxes for mean, submesoscale and turbulent fields across domains. Dashed lines mark the classic scalings (4.6) for EBF and MLE.

As anticipated, the submesoscale vertical buoyancy flux is active in the presence of baroclinic waves (UNSTABLE) and negligible in the other regions (figures 3a,c and 2c). In the UNSTABLE region, MLE restratification scales as expected, with a vertical profile of $\overline {w^sb^s}$ similar to that predicted ($\varPsi _{MLE}$) by Fox-Kemper et al. (2008) and a vertically averaged $\overline {w^sb^s}$ ($4.4\times 10^{-10}\,{\rm m}^2\,{\rm s}^{-3}$) within 6 % of the $\varPsi _{MLE}$ prediction.

The interpretation of mean $\overline {w^ab^a}$ through EBF is more nuanced. The weakening of the fronts during spin down implies that EBF scales are stronger than MLE for the data analysed. This stronger EBF is not found in $\overline {w^ab^a}$ (figure 3c). In UNSTABLE, $\overline {w^ab^a}$ ($4.4\times 10^{-10}\,{\rm m}^2\,{\rm s}^{-3}$) is similar and opposite to MLE, while in STABLE $\overline {w^ab^a}$ is smaller yet ($2.4\times 10^{-10}\,{\rm m}^2\,{\rm s}^{-3}$), balancing $\overline {w'b'}$. The scaling for EBF is rooted in horizontal advection of buoyancy by Ekman transport and so can also be compared with vertically averaged over the Ekman depth, which yields $6.5\times 10^{-10}\,{\rm m}^2\,{\rm s}^{-3}$ in UNSTABLE and $1.6\times 10^{-9}\,{\rm m}^2\,{\rm s}^{-3}$ in STABLE. Some of the differences may be the result of geostrophic stress excluded from the classical EBF scaling (Wenegrat & McPhaden 2016). Overall, the asymmetry in horizontal and vertical transport highlights nonlinear interactions that complicate the interpretation of EBF in well-developed flows. Ultimately, the interplay between horizontal vorticity injection by winds and restratification by the resulting overturning alters the Ekman overturning streamfunction.

Though surface forcing is uniform throughout, fine-scale turbulent fluxes $\overline {w'b'}$ vary across regimes as mean Ekman overturning and submesoscale flows interact with transport by wind-driven boundary layer eddies (figure 3b). This multi-scale interaction has implications for effective $\kappa _{eff}$ in each region (figure 4a). The smallest $\overline {w'b'}$ magnitude is in NOFRONT, absent of lateral buoyancy advection. Using the Monin & Obukhov (1954) similarity theory to arrive at non-dimensional $\kappa _{eff}$ in NOFRONT scales reasonably well with traditional 1-D boundary layer turbulence parametrizations, being within 20 % of ePBL and GLS, but with KPP exceeding diagnosed mixing by 60 %.

Figure 4. Effective diffusivity from different regions. (a) Turbulent diffusivity, $\kappa _{eff}$, is estimated as $\overline {w'b'}/N^2$ and is compared against $\kappa$ produced by three OSBL parametrizations, k-eps, ePBL and KPP. Note that in UNSTABLE, strong $M^2$ and negative PV regions,these parametrizations produce the wrong sign of diffusivity and in the STABLE region they overestimate the diffusivity. (b) Mixed-layer averaged $\kappa _{eff}$ across domains, normalized by $H$ and $u_*$.

The agreement between 1-D and $\kappa _{eff}$ falls off in the frontal regions. In STABLE, turbulent fluxes are negative and larger in magnitude than in the NOFRONT region as turbulent eddies drive down warm water advected near the surface by Ekman shear. The larger magnitude $\overline {w'b'}$ than NOFRONT is expected as enhanced stratification implies there is more buoyancy to flux. Yet the Monin & Obukhov (1954) non-dimensional mixing $\kappa _{eff}$ in this region is smaller than NOFRONT and 1-D, so turbulence is suppressed in the retratifying STABLE region compared to boundary layer scalings. This may seem obvious as lateral restratification is expected to suppress turbulence, but this suppression is stronger than Monin & Obukhov (1954) scaling predicts. These results are not trivial, particularly for similarity-based parametrizations where surface fluxes are intrinsic to the mixing amplitude and do not include information about lateral fluxes and lateral gradients – note that Monin & Obukhov (1954) theory assumes horizontal homogeneity which prohibits both. This implies that suppressed turbulence by lateral restratification is not well represented in 1D models and that current turbulence parametrizations are systematically over-mixing surface layers in the presence of lateral gradients and flows.

Turbulent buoyancy flux with NOFRONT and STABLE are both negative (positive $\kappa _{eff}$), consistent with the classic k-theory assumption that turbulence can be approximated as an eddy diffusivity mixing a flux down-gradient. This assumption breaks down in the UNSTABLE region, as $\overline {w'b'}$ becomes positive in the presence of positive mean $N^2$. The positive $\overline {w'b'}$ is not uniform across the region, but $\widetilde {w'b'}$ varies spatially ranging from $-2\times 10^{-10}$ to $1\times 10^{-9}$. Isolating regions of the strongest buoyancy gradients ($|\boldsymbol {\nabla }_H b|$ stronger than the initial front $M^2_0=2.1\times 10^{-8}\,{\rm s}^{-2}$) and regions of negative PV confirm that fine-scale circulations within the sharpest portions of the front are transporting buoyancy upward and against mean stratification, thereby resulting in negative $\kappa _{eff}$. A negative diagnosed $\kappa _{eff}$ may be associated with (non-local) convective instabilities or a result of frontal restratifying overturning circulations (as in the lateral-instability-driven circulations of Sullivan & McWilliams (2019)) which might be represented by the Held and Schnieder streamfunction, but here this effect is occurring on scales finer than the submesoscale. Therefore, the near-zero $\overline {w'b'}$ in the UNSTABLE region does not imply that fine-scale fluxes are not important; rather, fine-scale fluxes are transporting significant buoyancy up-gradient to counterbalance the transport down-gradient by surface-forced boundary layer eddies. This fine-scale circulation is not captured by coarser grain models that implement submesoscale parametrizations alongside OSBL parametrizations and again suggests that these models are over-mixing buoyancy in the presence of lateral gradients and flows.

5. The PV fluxes

The full Ertel PV, $q$, is defined as

(5.1)\begin{equation} q = (f\delta_{i3} + \omega_i)\frac{\partial b}{\partial x_i}, \end{equation}

where $\omega _i = \varepsilon _{ijk}({\partial u_k}/{\partial x_j})$ and subscript index $3$ is taken to be the rotation axis direction (vertical).

The nonlinearity that arises from the correlation between velocity gradients and buoyancy gradients implies that $\overline {(f\delta _{i3} + \omega '_i) \partial b'/ \partial x_i} \neq 0$, and that small-scale gradients impact mean PV. In a turbulent regime, PV at the small scale has a different behaviour from those expected from geophysical fluids (Bodner & Fox-Kemper 2020) as small-scale correlations create noisy fluctuations that dominate mean PV even on larger scales. Bodner & Fox-Kemper (2020) suggest that pre-filtered buoyancy and momentum be used to define PV in LES and that turbulent scale momentum and buoyancy fluxes be contained as a turbulent flux divergence. Therefore, PV relevant to geophysical flows is defined by submesoscale or larger fields, but the tendency of PV may be altered by fine-scale processes. The following section adopts the approach of Bodner & Fox-Kemper (2020) for the Reynolds-averaged equations, thus adapting it to be comparable with boundary layer parametrizations.

When defining the submesoscale PV ($\tilde {q}$), the Reynolds average includes the submesoscale-permitting fields only, $\tilde {c}$, not yet averaged in the along-front direction. A full derivation of the PV in the multi-Reynolds decomposition, the mean-eddy PV (MEPV), can be found in Appendix B. The submesoscale-permitting PV is

(5.2) \begin{equation} \tilde{q} = (f\delta_{k3} + \tilde{\omega}_i)\frac{\partial \tilde{b}}{\partial x_i}. \end{equation}

In the full MEPV decomposition (Appendix B), $\overline {\omega ^s_i\partial b^s/\partial x_i}$ contributes 1 %–5 % of $\tilde {q}$, yet it will be shown that submesoscale fluctuations play a leading-order role in $\tilde {q}$ tendencies.

Following Bodner & Fox-Kemper (2020), derivation of the submesoscale-averaged PV tendency equation begins with combining the turbulent transport terms into frictional and diffusive fluxes by defining

(5.3)\begin{gather} \mathcal{F}^{ + }_i = {-}\frac{\partial}{\partial x_j}(\widetilde{u_i'u_j'}) - \mathcal{V}_{ij}^{m}; \end{gather}

(5.4) \begin{gather}\mathcal{D}^{ + } = {-}\frac{\partial}{\partial x_j}(\widetilde{b'u_j'}) - \mathcal{V}_{ij}^{b}. \end{gather}

Note that (5.3) and (5.4) assume that Reynolds-averaging to the submesoscale-permitting scale ($\tilde {\cdot }$) commutes with differentiation. The uncertainty in (5.3) and (5.4) implied by the fine-scale variations in the choice of boundary locations for the $\tilde {\cdot }$ average can be estimated from Leibniz's theorem and gives an error estimate for the horizontal derivatives that are an order of magnitude larger than the signal (see Appendix A). Nonetheless, the Reynolds-averaged expression is most analogous to the parametrized form for turbulence solved by submesoscale-permitting simulations. Evaluating PV in this framework, in light of these uncertainties, allows a comparison between how PV is modelled in larger grid-scale ocean simulations to the $\tilde {q}$ in the LES.

The $\tilde {q}$ tendency equation can be found by multiplying the $\tilde {u}$ evolution equation by $\partial \tilde {b}/\partial x^i$, and multiplying the $\tilde {b}$ evolution equation by $(f\delta _{k3} + \widetilde {\omega }_i)$ (Appendix B). Combining the two and rearranging gives the PV evolution in flux form:

(5.5) \begin{equation} \frac{\partial}{\partial t}\tilde{q} = {-}\frac{\partial}{\partial x_i} \left[\underbrace{\tilde{u}_i \widetilde{q}}_{{\text{ADV}}} + \underbrace{\epsilon_{ikj}F_j^{ + }\frac{\partial \tilde{b}}{\partial x_k}}_{{\text{FRIC}}} - \underbrace{(f\delta_{k3} + \widetilde{\omega})\mathcal{D}^{ + }}_{{\text{DIA}}}\right]. \end{equation}

The advective flux (ADV) includes correlations between submesoscale currents and submesoscale gradients that define $\tilde {q}$. Similarly, turbulent scale motions interact with submesoscale buoyancy gradients through the friction (FRIC) and diabatic (DIA) flux terms.

The divergences of the ADV, FRIC and DIA terms are estimated from the upper 30 m in the LES regions, to avoid noise near the mixed-layer base. The magnitude of these terms for each region is shown in figure 5.

Figure 5. The MEPV tendency equation terms separated by region as described in (5.5). Green is $\partial /\partial x_i$ (ADV), blue is $\partial /\partial x_i$ (FRIC), purple is $\partial /\partial x_i$ (DIA)and dark grey is the sum of all three terms. Light blue and light purple represent the parametrized equivalent of FRIC and DIA, where the LES turbulent fluxes have been replaced with parametrized fluxes estimated from 1-D models (see table 1). Light grey is the sum of the advective term (green) and the parametrized terms (light blue and light purple).

In the STABLE region, ADV increases PV, consistent with Ekman overturning which drives restratification, while DIA competes with ADV to decrease PV. The process underlying DIA's reduction of PV is dominated by submesoscale correlations between the vertical buoyancy flux divergence and vertical vorticity as wind-driven turbulent motions homogenize buoyancy in the presence of geostrophic frontal flow.

In the UNSTABLE region, ADV, FRIC and DIA increase PV, consistent with the dominance of MLE restratification over destratifying Ekman overturning in the buoyancy budget discussed in § 4. Small-scale DIA transport increases PV as positive vertical buoyancy flux interacts with vertical vorticity gradients near sharp fronts. The tendency for fine-scale turbulence to increase PV is reminiscent of the injection of PV into the mixed layer (ADV) by baroclinic instability as described by Boccaletti, Ferrari & Fox-Kemper (2007) as well as due to frictional geostrophic stress (FRIC) at the surface (Wenegrat & McPhaden 2016) and secondary instabilities such as symmetric instability (Thomas et al. 2013).

In § 4, the diffusivities ($\kappa _{eff}$) in the unstable and stable regions are markedly different from that predicted by current state-of-the-art OSBL parametrizations. It is therefore of interest to explore how approximating turbulent fluxes with an OSBL parametrization might impact PV tendency (figure 5, light shading for parametrizations). Parametrized turbulence from GOTM can easily be used to replace the turbulent flux divergences in (5.3) and (5.4) with the following (assuming isotropic mixing on fine scales):

(5.6)\begin{gather} \mathcal{F}^{ + }_j \rightarrow \mathcal{F}^{param} = \frac{\partial}{\partial x_j}\left(\nu_H \frac{\partial \tilde{u}^a_i}{\partial x_j}\right), \end{gather}

(5.7)\begin{gather}\mathcal{D}^{ + } \rightarrow \mathcal{D}^{param} = \frac{\partial}{\partial x_j}\left(\kappa_H \frac{\partial \tilde{b}^a}{\partial x_j}\right), \end{gather}

where $\nu _H$ and $\kappa _H$ are from the GOTM 1-D models’ mixed-layer eddy viscosity and diffusivity, respectively. The mixing coefficients $\nu _H$ and $\kappa _H$ are specific to each submesoscale average grid and estimated using the average of the non-dimensionalized 1-D $\kappa$, i.e. taking the average of $\kappa _{\theta } {u_*}^{-1} H^{-1}$ over the ensemble of 1-D models in figure 4. The dimensional diffusivity and viscosity are then restored (per Monin & Obukhov 1954), by multiplying the average by the local $u_* H$ for each submesoscale bin. This PARAM estimate is analogous to what would occur in a submesoscale-permitting model that uses parametrized turbulence in the form $\kappa _c = \mu v_t l$, where $\mu$ is a non-dimensional coefficient, $v_t$ is the turbulent velocity scale (proportional to $u_*$) and $l$ is a typical turbulence length scale proportional to the depth of the boundary layer (i.e. $H$) (Large et al. 1994; Tennekes & Lumley 2018). Further assumptions are isotropy and that the Prandtl number is assumed equal to one, $Pr=\nu _H/\kappa _H=1$, so $\nu _V=\nu _H=\kappa _V=\kappa _H$. With these assumptions and (5.3)–(5.4) and (5.6)–(5.7), the PV fluxes of the parametrizations are found. The vertical gradients dominate the flux divergence, accounting for 99 % of the fluxes.

In the STABLE region, $\mathcal {F}^{param}$ and $\mathcal {D}^{param}$ are of opposite sign, but $\mathcal {D}^{param}$ and their sum are much larger in magnitude than $\mathcal {F}^{+}$ and $\mathcal {D}^{+}$ (figure 5). The dominant signal is the overestimation of $\mathcal {D}^{param}$, consistent with the results in § 4, where the resolved turbulent flux of buoyancy in the STABLE region was suppressed compared with parametrized fluxes and Monin–Obukhov scaling. However, in the dimensionless diffusivities, this suppression was only a factor of two (figure 4a), whereas $\mathcal {D}^{param}$ is more than four times $\mathcal {D}^+$ indicating that the averaged gradients and covariation over a submesoscale grid scale lead to further errors.

In the UNSTABLE region, $\mathcal {F}^{param}$ and $\mathcal {D}^{param}$ are of the same sign, but both are larger in magnitude and of opposite sign to their turbulence diagnosis partners $\mathcal {F}^{+}$ and $\mathcal {D}^{+}$ (figure 5). Again, $\mathcal {D}^{param}$ dominates the parametrized PV tendency terms, as parametrized fluxes mix buoyancy down-gradient. The positive $\mathcal {D}^{+}$ from fine-scale PV injection into the mixed layer is not captured by the parametrized mixing $\mathcal {D}^{param}$, thereby inaccurately representing PV changes in the presence of unstable flows. Therefore, in both STABLE and UNSTABLE regions, but for different reasons, parametrized mixing tends to drive PV towards negative values – in the opposite direction or faster than implied by the resolved LES fine-scale turbulent fluxes. Over the whole domain, this implies that the parametrizations predict a completely inaccurate PV budget, which should be a caution for using even the sign of PV fluxes as inferred from parametrizations.

6. Dissipation of kinetic energy

Instabilities of submesoscale frontal flows are known to transfer energy to smaller scales, thereby enhancing the dissipation of kinetic energy. Terms in the turbulent or eddy kinetic energy equations, often used to gain insight into transfer pathways for energy towards smaller scales, are difficult to estimate locally (Cao, Fox-Kemper & Jing 2021) and therefore are not calculated here. Yet the dissipation rate of kinetic energy, $\epsilon$, can be averaged over domains and is used here to understand how energy is dissipated in the presence of submesoscale fronts. Following Sullivan, McWilliams & Moeng (1994) and Bodner & Fox-Kemper (2020), the grid-scale dissipation rate, $\epsilon$, is calculated from the subgrid turbulent kinetic energy ($k$) and modelled turbulent length scale that depends on $k$ and grid-scale stratification. Profiles of $\bar {\epsilon }$ (figure 6a) from each region are compared with that estimated from GLS (Umlauf & Burchard 2003). Near-surface deviations between LES and GLS are expected partly due to discretized implementation of the boundary conditions and partly due to parametrization inaccuracy.

Figure 6. (a) Profiles of turbulent energy dissipation averaged over different regions of the domain, $\bar {\epsilon }$, including NOFRONT (orange), STABLE (green), UNSTABLE (blue), UNSTABLE region with $M^2>M_o^2$ (purple), UNSTABLE region with ${PV}<0$ (magenta) and UNSTABLE region with $H/h<0.6$ (red dashed). General length scale $\epsilon$ is included for reference (grey dashed line). (b) Depth and region averages of $\bar {\epsilon }$ with colours corresponding to lines in (a). The grey dashed line in the average $\epsilon$ across the domain.

Below this surface deviation, $\bar {\epsilon }$ in GLS is similar to $\bar {\epsilon }$ in NOFRONT, in line with the agreement between NOFRONT $\kappa _\textrm {eff}$ and $\kappa$ from 1-D models in § 4. In UNSTABLE, vertically integrated LES $\bar {\epsilon }$ is 26 % larger than in NOFRONT (figure 6b), reinforcing that submesoscale instabilities enhance dissipation (Thomas & Taylor 2010; D'Asaro et al. 2014). Conversely, dissipation rate is suppressed by lateral advection and its restratification in STABLE (figure 6b; Taylor 2016; Taylor et al. 2020), with a vertically integrated $\bar {\epsilon }$ that is 20 % less than that in NOFRONT. This compensation between enhanced and suppressed dissipation results in the entire domain having just 1 % larger $\bar {\epsilon }$ than the NOFRONT boundary layer turbulence generated by wind input into a homogeneous ocean. The overall dissipation is thus surprisingly similar to that if the submesoscale effects were neglected.

7. Discussion

To further understand the impact of different flow regimes, values of $\widetilde {w'b'}$ and $\tilde {\epsilon }$ are plotted against the submesoscale PV, $\tilde {q}$ from (5.2), and the submesoscale advection (ADV) of buoyancy by horizontal shear, $-\partial \tilde {u}_i/\partial z$ $\partial \tilde {b}/\partial x_i$ (figure 7). In regions of negative PV, flow may be subject to instability, including gravitational and symmetric instabilities. Differential advection of buoyancy gradients impacts stratification, with positive ADV increasing stratification and negative ADV reducing stratification. The flow in NOFRONT occupies a narrow region of this parameter space, with homogeneous values of $\epsilon$ and $\widetilde {w'b'}$. It was shown that the transport of buoyancy (§ 4) and dissipation of energy (§ 6) in NOFRONT agree with estimates from classic boundary layer turbulence theory and are represented by turbulence parametrizations. The rest of this section will focus on STABLE and UNSTABLE, and how these regions differ from NOFRONT and classic boundary layer turbulence theory.

Figure 7. The correlation between PV (5.1) and the horizontal advection (ADV) of buoyancy by the horizontal shear, $-{\rm d}\tilde {u}_i/{\rm d} z \,{\rm d}\tilde {b}/{\rm d} x_i$ (i.e. a horizontal advective term in the stratification tendency equation). The lower-left quadrant is destratifying flow with negative PV, and the upper-right quadrant is stratifying flow with positive PV. Dots are coloured by magnitude of $\widetilde {w'b'}$ and $\tilde {\epsilon }$: (a) NOFRONT $\tilde {\epsilon }$, (b) STABLE $\tilde {\epsilon }$, (c) UNSTABLE $\tilde {\epsilon }$, and parameter space of STABLE and NOFRONT as in (a,b); (d) NOFRONT $\widetilde {w'b'}$, (e) STABLE $\widetilde {w'b'}$, (f) UNSTABLE $\widetilde {w'b'}$, and parameter space of STABLE and NOFRONT as in (d,e).

7.1 Stable region

In STABLE, ADV tends to restratify and tends towards positive PV within this region, with a narrow range of $\tilde {\epsilon }$ and positive $\widetilde {w'b'}$. The strongest $\widetilde {w'b'}$ occurs in positive PV flows where strong stratification from lateral advection is mixed by wind-driven turbulence. Strong PV and enhanced $\widetilde {w'b'}$ are coincident with attenuated $\tilde {\epsilon }$ as dissipation of kinetic energy is suppressed in regions of lateral advection and its restratification. Both FRIC and DIA tend to reduce the positive PV in this region (as predicted but overestimated by parametrizations; figure 5).

7.2 Unstable region

The flow in UNSTABLE exhibits larger fluctuations across phase space and in value and sign compared with the other regions, highlighting the heterogeneity of unstable flows and the wide range of processes that contribute to the along-front and domain average under downfront winds. Positive ADV tendency co-occurs with the strongest upward flux of buoyancy and smaller values of $\epsilon$: intensive, intermittent restratification with weak turbulence dissipation. Negative ADV tendency and, particularly, though not exclusively, negative PV are associated with the strongest downward/destratifying flux of buoyancy and the largest rates of energy dissipation: intermittent convective mixing, entrainment and overturning within the submesoscale features. In summary, restratifying, positive ADV wins out in the resolved fluxes, counter to that predicted using parametrized fluxes (figure 5), and points to the importance of small scales for vertical heat transports (Cao, Jing & Fox-Kemper 2024).

A positive buoyancy flux in the presence of stable stratification implies that small-scale motions act to mix cold water down (i.e. convection) or include 3-D circulations that increase stratification. The resulting negative effective $\kappa$ (figure 4) indicates that the fine-scale flow is not isotropic turbulence which would flux down the gradient. Similarly, baroclinic instabilities wrap up to increase local and overall stratification. In the same regions, symmetric instability (SI) develops slantwise convection that extracts PV from the thermocline and increases mixed-layer stratification and PV. Although the along-isopycnal structures commonly associated with SI (e.g. Haney et al. 2015) are not discernable in this complex flow, the PV and buoyancy fluxes resemble this phenomenon. However, regions of the strongest stratification have lower energy dissipation rates, which is not consistent with the present understanding of ongoing SI (Taylor & Ferrari 2009).

Conditions for SI arise as downfront winds drive down PV at ocean fronts, such that $qf<0$. Though SI is known to enhance dissipation rates (Taylor & Ferrari 2010; Thomas & Taylor 2010), it is intermittent and highly localized. The Stone (1966) scaling for the largest unstable SI provides a scale of 50–100 m under the conditions in this simulation, thus over 10 collocation points per instability are available. Though smaller symmetrically unstable modes are predicted to grow, transport is typically dominated by larger modes due to their size, with less dissipation rate, leading to higher energy content. Taylor & Ferrari (2009) show that a pathway from SI to smaller, secondary instability turbulence is important for energy dissipation. The $4.9\,{\rm m}\times 1.25\,{\rm m}$ resolution here is unlikely to capture these secondary instabilities, but their dissipative effect is parametrized by the LES closure (Sullivan et al. 1994). The instability can be active when the ratio between the convective depth (due to unstable Ekman transport), $\tilde {h}$, and the mixed-layer depth, $\tilde {H}$, is small, or when $h/H<<1$. This ratio is found to obey (Taylor & Ferrari 2010)

(7.1)\begin{equation} \left(\frac{\tilde{h}}{\tilde{H}}\right)^4 - c^3 \left(1-\frac{\tilde{h}}{\tilde{H}}\right)^3\left[\frac{u_ * ^3}{|{\rm \Delta} \tilde{u}_g|^2}\cos\theta\right]^2. \end{equation}

The SI parametrization put forth by Bachman et al. (2017) requires a threshold of $h/H<0.9$ for SI to be considered important, with stronger SI occurring as $h/H\rightarrow 0$. In the UNSTABLE region, $\tilde {h}/\tilde {H}<0.9$ occupies 26 % of the domain, whereas $\tilde {h}/\tilde {H}<0.6$ occupies 1 % of the UNSTABLE domain. Isolating regions $\tilde {h}/\tilde {H}<0.9$ does not significantly change $\epsilon$ compared with NOFRONT, but the 1 % of the region with $\tilde {h}/\tilde {H}<0.6$ has 86 % higher dissipation than NOFRONT. Overall, there is no correlation between $\epsilon$ and ADV or PV.

8. Conclusions and recommendations

While Monin & Obukhov (1954) and 1-D parametrizations do a reasonable job of approximating the turbulent mixing in terms of diffusion, energy dissipated and PV sources in between fronts, there are major discrepancies with these theories in both an upfront-wind/Ekman-stabilized front and in a downfront-wind/Ekman-destabilized front, even when separately accounting for the effects of submesoscale mixed-layer eddy restratification. A main message from this work is that surface-forced boundary layer turbulence is modified – sometimes even reversed in net effect – by the presence of submesoscale flows, with implications for how regional and global circulation models are able to simulate the distribution of heat and momentum in the near-surface ocean. In the STABLE region, buoyancy fluxes and dissipation rates are suppressed compared with NOFRONT and OSBL parametrizations, while in the UNSTABLE region, buoyancy fluxes are suppressed, but dissipation rates are enhanced.

In the UNSTABLE region, a positive mean buoyancy flux and negative $\kappa _{eff}$ suggest that fine-scale circulations transport dense fluids up-gradient. From a broad perspective, this mixing resembles the non-local fluxes embedded in KPP used to simulate convection during unstable conditions. In KPP, the non-local term is meant to capture the transport of cold water from the surface to depth while ignoring smaller gradients along the way; here, the negative $\kappa _{eff}$ represents fine-scale overturning due to a combination of surface cooling, convection from Ekman transport and submesoscale overturning circulations. While transport equivalent to a net overturning on the submesoscale is expected at submesoscales (Fox-Kemper et al. 2008), this diagnosis finds non-negligible fine-scale up-gradient transport that competes with and often overtakes boundary layer mixing. An examination of PV generated by turbulence in this UNSTABLE region suggests that present parametrizations fail in sign (destruction rather than creation) and magnitude (too large of an effect). This region has the most variations in buoyancy and energy dissipation – indicating a heterogenous collection of restratifying flow in high-PV, counter-gradient, weakly energetic fine-scale subregions and destratifing flow in low-PV, strongly energetic plumes. Overall, restratification and positive PV injection dominate in this region in this simulation, which is opposite to that predicted by parametrizations.

In the STABLE region, turbulence is suppressed, more so than the Monin & Obukhov (1954) expectations due to the increased stratification from lateral advection and the accompanying shallower boundary layer. Consistently, in this region, the PV tendency and dissipation rates predicted by 1-D scalings are overestimated compared with those diagnosed in the simulation.

The Reynolds-averaged approach for estimating PV fluxes adopted here lends itself to be compared with what is assumed in numerical schemes for simulations with grid scales too large to resolve fine-scale motions and instead use boundary layer parametrizations. This can be seen by Leibniz's theorem (see Appendix A) where the Reynolds-averaged approach results in a loss of information compared with the filtered approach presented by Bodner & Fox-Kemper (2020). The propagation of uncertainty from Leibniz's theorem suggests that errors in the frictional and diabatic PV flux terms are orders of magnitude larger than the fluxes themselves. This questions our ability to understand the role of parametrized down-gradient turbulent fluxes in PV tendencies in the presence of submesoscale flows.

Submesoscale flows alter the kinetic energy dissipation rate, with $\epsilon$ suppressed in regions of lateral restratification – by large-scale fronts or intermittent submesoscales – and enhanced in the presence of negative PV and destratifying circulations. While much work has been done to evaluate the enhanced dissipation and mixing due to submesoscales, in this simulation, enhanced dissipation at the downwind front is balanced by turbulence suppression at the upwind front such that the average dissipation is indistinguishable on average over the whole domain from that predicted by boundary layer scalings. More work is needed to understand the interaction of submesoscales and turbulence across parameter space (Taylor & Thompson 2023).

Finally, this simulation is in a modest part of parameter space for the most energetic instabilities, with $M^2\sim 10^{-8}\,{\rm s}^{-2}$, which is more than an order of magnitude less than observed during frontal process studies (e.g. $M^2\sim 10^{-6}\,{\rm s}^{-2}$; D'Asaro et al. 2014; Johnson et al. 2020). Other recent simulations (e.g. Sullivan & McWilliams 2019) have much stronger fronts but domains that are too small to resolve mixed-layer eddies and the intricate processes as the front destabilizes. More work is needed to extend these results to more energetic frontal cases with varying surface fluxes.