The Effects of Surface Wind Stress and Buoyancy Flux on the Evolution of a Front in a Turbulent Thermal Wind Balance

Abstract

Here we consider the effects of surface buoyancy flux and wind stress on a front in turbulent thermal wind (TTW) balance using the framework of Crowe and Taylor (2018). The changes in the velocity and density profiles induced by the wind stress and buoyancy flux interact with the TTW and can qualitatively change the evolution of the front. In the absence of surface-forcing, Crowe and Taylor (2018) found that shear dispersion associated with the TTW circulation causes the frontal width to increase. In many cases, the flow induced by the surface-forcing enhances the spreading rate. However, if the wind stress drives a cross-front flow which opposes the frontal buoyancy gradient or the buoyancy flux drives an unstable stratification, it is possible to obtain an up-gradient cross-front buoyancy flux, which can act to sharpen the front. In certain conditions, an equilibrium state develops where the tendency for the TTW circulation to spread the front is balanced by the frontogenetic tendency of the surface forces. We use numerical solutions to a nonlinear diffusion equation in order to test these predictions. Finally, we describe the connection between surface-forcing and vertical mixing and discuss typical parameters for mid-ocean fronts.

1. Introduction

Fronts, or regions of large horizontal density gradient are important features of the upper ocean and atmosphere where they play an important role in the exchange of heat, carbon, and other important tracers. Fronts are highly anisotropic features where the cross-front width (in the direction of the density gradient) is small compared to the along-front length. Fronts can be characterised by their cross-front width which can range from tens of kilometers for large fronts including the Gulf Stream, Kuroshio, and fronts associated with the Antarctic Circumpolar Current (ACC), to scales of meters for fronts associated with freshwater plumes and gravity currents. The focus here will be on relatively large scale fronts where the Rossby number is small and the frontal width is large compared to the mixed layer depth. A typical front is shown in Figure 1.

The dominant balance in many large-scale frontal systems is between the geostrophic shear and horizontal density gradients, resulting in a balanced state with an along-front ‘thermal wind’ flow [1]. Elevated levels of small-scale turbulence in the ocean surface mixed layer act to disrupt this balance by driving a circulation around the front which can in turn lead to an intensification of the front (frontogenesis) or spreading of the front (frontolysis) [2,3]. Including the effects of vertical mixing in the horizontal momentum equations leads to a state of ‘turbulent thermal wind’ (TTW) or ‘generalised Ekman’ balance, which has been seen in both numerical models of ocean fronts [3,4,5,6] and observational data [7].

Frontogenesis is the process by which fronts form and sharpen. Externally imposed strain flows are often used to model the squeezing of fronts by large-scale currents and analytic work by Hoskins and Bretherton [8] found that for an inviscid front, this external strain flow can lead to a finite-time discontinuity in the surface density. Physically, we expect this collapse to be prevented by the effects of enhanced turbulence and mixing in the frontal region [9].

Surface wind and buoyancy forcing can have a strong influence on the evolution of ocean fronts. For example, Thomas and Lee [10] showed that when the surface wind stress is aligned with the thermal wind (in the ‘down-front’ direction), the resulting Ekman flow destabilises the water column resulting in convection, which, combined with Ekman pumping, drives a frontogenetic secondary circulation. This secondary circulation can spontaneously form very sharp fronts and also results in large vertical velocities. Thomas and Ferrari [11] examined the relative influence of wind-driven circulation and horizontal strain flow on the secondary circulation and stratification at ocean fronts and found that wind-driven circulation can dominate strain-driven frontogenesis for typical ocean conditions.

Surface buoyancy fluxes and wind stresses also act to destabilise vertical density profiles leading to small scale instabilities and generating mixed layer turbulence [12] however, observations by D’Asaro et al. [13] and numerical work by Capet et al. [14] suggests that boundary layer turbulence can be significantly enhanced by the presence of fronts and is further increased by the intensification of surface buoyancy gradients during frontogensis [9]. Surface fluxes and wind stresses have been used to drive turbulent mixing in many large-eddy simulations (LES) of ocean fronts [9,15,16] where the properties of sub-grid scale turbulence are parametrised in terms of an eddy viscosity.

The secondary circulation driven by small-scale turbulence can have frontogenetic effects as shown by McWilliams [17] who considered the frontogenetic tendency of an idealised front in the presence of strong turbulent mixing. The secondary circulation resulting from TTW balance drives a cross-front differential advection which leads to a sharpening of the horizontal surface buoyancy gradient on the dense side of the front. However, in Crowe and Taylor [1] we showed that the TTW circulation can also lead to frontal spreading through shear dispersion in the case of a small Rossby number. The cross-front shear associated with the TTW circulation acts to project the vertical mixing horizontally, leading to a nonlinear diffusion equation for the depth-averaged density field. This equation can be solved analytically in terms of a spreading similarity solution.

The central aim of this paper is to study the competition between frontogenesis and frontolysis associated with vertical mixing and surface-forcing in the form of an applied surface wind stress and surface buoyancy flux. This leads to the following question: for a prescribed level of mixing and surface forcing, does the front evolve towards an equilibrium state, and if so, what controls the associated frontal width and horizontal buoyancy gradient? An answer to this question could provide a basis for improving parametrisations of submesoscale processes where the horizontal buoyancy gradient is an important input parameter that is often under-resolved in ocean models (e.g., [18,19]).

To address these questions, we extend the model of Crowe and Taylor [1], where turbulence is represented by vertical mixing, to include a surface wind stress and buoyancy flux imposed as boundary conditions. This system is solved analytically in the limit of small Rossby number in Section 3 in terms of a background buoyancy field. The evolution of the background buoyancy field is examined in Section 4 for different values of wind stress and heat flux. We find that positive heat fluxes lead to a spreading of the front while negative heat fluxes may lead to sharpening if they are sufficiently strong to overcome the spreading effects of horizontal mixing. Similarly, wind stresses can cause either spreading or sharpening according to the strength and direction of the wind. These different cases are examined in Section 5 using numerical simulations of the buoyancy evolution equation.

In Section 6, we use the mixing length arguments of Taylor and Ferrari [16] and Enriquez and Taylor [20] to link the vertical diffusivity to the strength of the surface-forcing and hence describe the frontal evolution in terms of the external forcing and buoyancy gradient. Finally, in Section 7, we use typical parameters to determine if equilibrium fronts are possible physically and on what scales they might occur.

2. Governing Equations

Here we consider a three-dimensional front in a geometry bounded from above and below by flat, rigid horizontal surfaces separated by a distance H. We assume that the fluid is horizontally infinite and rotating about the z-axis with Coriolis parameter f. We define $\Delta b$ to be a typical buoyancy difference across the front and L to be a typical horizontal lengthscale with the assumption that the depth of the fluid is thin and hence $H\ll L$.

We invoke the Boussinesq approximation and use a linear equation of state to represent density by a single scalar equation for buoyancy. The effects of turbulence are modelled by a constant viscosity, $\nu$, and diffusivity, $\kappa$. Note that we could follow Crowe and Taylor [1] and take $\nu$ and $\kappa$ to be depth-dependent and still solve the resulting system. A depth-dependent $\nu$ and $\kappa$ would not qualitatively change our method or solution and would lead only to a modification of the coefficients in our solution. Since the coefficients would have to be calculated numerically for general $\nu =\nu \left(z\right)$ and $\kappa =\kappa \left(z\right)$, here we will assume depth-independent values for ease of calculation and simplify the presentation [21,22].

The effects of wind and surface heat fluxes are represented by a two-dimensional wind stress vector, ${\tau }^{*}$, and a surface heat flux $Q^{*}$. Note that we could take ${\tau }^{*}$ and $Q^{*}$ to depend on $(x,y)$ and still obtain solutions, although here we consider constant values for simplicity.

We now non-dimensionalise the governing equations using the horizontal velocity scale, $U=\Delta bH/\left(fL\right)$, vertical velocity scale $W=UH/L=\Delta b{H}^2/\left(f{L}^2\right)$, pressure scale $P=fUL=\Delta bH$, and timescale $T=L/U=f{L}^2/(H\Delta b)$. We define the Rossby number, $\mathrm{Ro}=U/\left(fL\right)$, using the geostrophic shear, $U/H=\partial b/\partial x/f=\Delta b/\left(fL\right)$ and write the aspect ratio, $H/L$, as $ϵ$. The Ekman number and Prantl number are $\mathrm{E}=\nu /\left(f{H}^2\right)$ and $\mathrm{Pr}=\nu /\kappa$ respectively. The dimensionless wind stress and buoyancy flux are denoted $\tau$ and Q respectively and are given in terms of dimensional parameters in Section 3. This gives the following non-dimensional equations [23]:

$$\begin{array}{cc}\hfill \mathrm{Ro}\frac{Du}{Dt}-v& =-\frac{\partial p}{\partial x}+\mathrm{E}{▽}_{ϵ}^{2}u,\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill \mathrm{Ro}\frac{Dv}{Dt}+u& =-\frac{\partial p}{\partial y}+\mathrm{E}{▽}_{ϵ}^{2}v,\hfill \end{array}$$

(1b)

$$\begin{array}{cc}\hfill \mathrm{Ro}{ϵ}^2\frac{Dw}{Dt}& =-\frac{\partial p}{\partial z}+b+{ϵ}^2\mathrm{E}{▽}_{ϵ}^{2}w,\hfill \end{array}$$

(1c)

$$\begin{array}{cc}\hfill \mathrm{Ro}\frac{Db}{Dt}& =\frac{\mathrm{E}}{\mathrm{Pr}}{▽}_{ϵ}^{2}b,\hfill \end{array}$$

(1d)

$$\begin{array}{cc}\hfill \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}& =0,\hfill \end{array}$$

(1e)

with top and bottom boundary conditions

$$\begin{array}{cc}\hfill \mathrm{E}\frac{\partial {\mathbf{u}}_H}{\partial z}& =\tau ,\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill \frac{\mathrm{E}}{\mathrm{Pr}}\frac{\partial b}{\partial z}& =Q,\hfill \end{array}$$

(2b)

at $z=1/2$ and $z=-1/2$. The parameters are given in

Table 1. Definitions of the dimensionless parameters and their values for buoyancy difference, $\Delta b$, Coriolis parameter, f, horizontal lengthscale, L, vertical lengthscale, H and dimensional viscosity and diffusivity, $\nu$ and $\kappa$.

Parameter	Rossby No.	Ekman No.	Prandtl No.	Aspect Ratio
Symbol	$\mathrm{Ro}$	$\mathrm{E}$	$\mathrm{Pr}$	$ϵ$
Definition	$ϵ\Delta b/f^{2}L$	$\nu /f{H}^2$	$\nu /\kappa$	$H/L$

and

$$\frac{D}{Dt}=\frac{\partial }{\partial t}+u\frac{\partial }{\partial x}+v\frac{\partial }{\partial y}+w\frac{\partial }{\partial z},$$

(3)

denotes the advective or total derivative and

$${▽}_{ϵ}^2=\frac{{\partial }^2}{\partial z^2}+{ϵ}^2\left[\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}\right].$$

(4)

denotes the rescaled Laplacian operator. The horizontal velocity, ${\mathbf{u}}_H$, is given by $(u,v,0)$ and the stress is a 2D vector given by $\tau =({\tau }_x,{\tau }_y,0)$. Our nondimensional setup is shown in Figure 2 for a typical isolated front.

Note that we have applied opposing wind stress and buoyancy flux to each surface due to the change in the direction of the surface normal vector. For stress-driven flow in a shallow sea, we might anticipate that the bottom stress will oppose the surface stress at a steady state. However, rather than representing any physical processes, our approach is primarily for mathematical convenience as it makes the problem vertically symmetric and eliminates any net translation of the front. We expect that having forcing at the top and bottom likely enhances the effect of the forcing by about a factor of two for a given forcing strength, relative to what would happen if the forcing was only applied to one boundary. A similar choice was made in the simulations of convection reported in Callies and Ferrari [24].

3. Analytic Solution

We consider the limit of small Rossby number and aspect ratio, $ϵ=O\left(\mathrm{Ro}\right)$, and assume that the Ekman number and Prandtl number are order 1. We now expand the fields, $\mathbf{u}$ and b, in powers of $\mathrm{Ro}$ as

$$(\mathbf{u},b)=({\mathbf{u}}_0,b_0)+\mathrm{Ro}({\mathbf{u}}_1,b_1)+\dots ,$$

(5)

and expand Equations (1) and (2) in powers of $\mathrm{Ro}$. We also use a multiple timescale approach to expand the time derivatives as

$$\frac{\partial }{\partial t}=\frac{1}{\mathrm{Ro}}\frac{\partial }{\partial \mathcal{T}}+\frac{\partial }{\partial t}+\mathrm{Ro}\frac{\partial }{\partial T},$$

(6)

for fast timescale, $\mathcal{T}=t/\mathrm{Ro}$, intermediate timescale, t, and slow timescale, $T=\mathrm{Ro}t$. We can now consider the resulting equations at each order in $\mathrm{Ro}$. For simplicity we consider the case of constant $\nu$ and $\kappa$. Solutions can be found for $\nu$ and $\kappa$ as general functions of z (see [1]). Depth-averaged quantities are written as $\overline{f}$ and departures from the depth-average are defined as $f^{\prime }=f-\overline{f}$. We consider a vertical domain $z\in [-0.5,0.5]$ so the depth average is defined as

$$\overline{f}={\int }_{-1/2}^{1/2}fdz.$$

(7)

The dimensionless heat flux is given by

$$Q=\frac{\alpha g{Q}^{*}}{c_p{\rho }_{0}fHL{M}^2},$$

(8)

and the dimensionless wind stress magnitude by

$$\tau =\frac{{\tau }^{*}}{{\rho }_0{H}^2{M}^2},$$

(9)

where $Q^{*}$ and ${\tau }^{*}$ are the dimensional heat flux and wind stress magnitude respectively. Here $\alpha$ is the heat capacity, ${\rho }_0$ is a reference density, g is the gravitational acceleration, and $c_p$ is the specific heat capacity. Using typical parameters of $H=100\mathrm{m}$, $L=10\mathrm{km}$, $M^2=\Delta b/L={10}^{-8}{\mathrm{s}}^{-2}$, ${\tau }^{*}=0.1{\mathrm{Nm}}^{-2}$, $Q^{*}=100{\mathrm{Wm}}^{-2}$, and $f={10}^{-4}{\mathrm{s}}^{-1}$ we find that $Q\sim 0.01$ and $\tau \sim 1$, and motivated by this, we let $Q=\mathrm{Ro}{Q}_1$ and $\tau ={\tau }_0$.

3.1. O(1) Balance

The leading order buoyancy equation is

$$\frac{\partial b_0}{\partial \mathcal{T}}=\frac{\mathrm{E}}{\mathrm{Pr}}\frac{{\partial }^2{b}_0}{\partial z^2},$$

(10)

with boundary conditions

$$\frac{\mathrm{E}}{\mathrm{Pr}}\frac{\partial b_0}{\partial z}=0,$$

(11)

at $z=\pm 1/2$. Assuming that any transients acting on the fast timescale have decayed, the solution is

$$b_0=b_0(x,y,t,T),$$

(12)

which is vertically homogeneous. The leading order vertical momentum equation is hydrostatic balance

$$\frac{\partial p_0}{\partial z}=b_0,$$

(13)

and hence

$$p_0={\overline{p}}_0(x,y,t,T)+b_0(x,y,t,T)z.$$

(14)

The leading order horizontal momentum equation is

$$\frac{\partial {\mathbf{u}}_{H0}}{\partial \mathcal{T}}+\mathbf{k}\times {\mathbf{u}}_{H0}=-{▽}_H{p}_0+\mathrm{E}\frac{{\partial }^2{\mathbf{u}}_{H0}}{\partial z^2},$$

(15)

which, upon ignoring fast time transients and substituting for $p_0$, gives

$$\mathbf{k}\times {\mathbf{u}}_{H0}=-{▽}_H{\overline{p}}_0-z{▽}_H{b}_0+\mathrm{E}\frac{{\partial }^2{\mathbf{u}}_{H0}}{\partial z^2},$$

(16)

with boundary conditions

$$\mathrm{E}\frac{\partial {\mathbf{u}}_{H0}}{\partial z}={\tau }_0,$$

(17)

at $z=\pm 1/2$. This system can be solved (see Appendix A) to obtain the solution

$$\begin{array}{c}\hfill {\mathbf{u}}_{H0}=-{▽}_H\times \left({\overline{p}}_0\mathbf{k}\right)+\frac{1}{\sqrt{\mathrm{E}}}\left[\left(K_0+\frac{z}{\sqrt{\mathrm{E}}}\right){\tau }_0-K_0^{″}\mathbf{k}\times {\tau }_0\right]\\ \hfill -\sqrt{\mathrm{E}}\left[K_0^{″}{▽}_H{b}_0+K_0\mathbf{k}\times {▽}_H{b}_0\right],\end{array}$$

(18)

where $K_0$ describes the modification of the thermal wind velocity by vertical mixing and is defined in Equation (A7). Note that ${\overline{p}}_0$ acts as a streamfunction for the depth-independent component of velocity, as in geostrophic balance. The first term in square brackets describes the stress-driven components of the velocity with components both parallel and perpendicular to the surface stress. The second term in square brackets describes the buoyancy-driven flow and also has components parallel and perpendicular to the direction of the buoyancy gradient. In the case of no wind stress, ${\tau }_0=\mathbf{0}$, and no depth-independent flow, ${\overline{p}}_0$, the cross-front and along front velocities are given respectively by

$$u=-\sqrt{\mathrm{E}}{K}_0^{″}{▽}_H{b}_0\mathrm{and}v=-\sqrt{\mathrm{E}}{K}_0{▽}_H{b}_0,$$

(19)

where we note that $-\sqrt{\mathrm{E}}{K}_0$ reduces to the linear ‘thermal wind’ profile, z, in the case of $\mathrm{E}\to 0$ [1]. By mass conservation we can calculate the vertical velocity as

$$w_0=K_0^{‴}{▽}_H·{\tau }_0-K_0^{\prime }\mathbf{k}·[{▽}_H\times {\tau }_0]+\mathrm{E}{K}_0^{\prime }{▽}_H^2{b}_0,$$

(20)

which reduces to

$$w_0=\mathrm{E}{K}_0^{\prime }{▽}_H^2{b}_0,$$

(21)

in the case where ${\tau }_0$ is homogeneous.

3.2. $O\left(\mathrm{Ro}\right)$ Balance

The $O\left(\mathrm{Ro}\right)$ buoyancy equation is

$$\frac{\partial b_0}{\partial t}+\frac{\partial b_1}{\partial \mathcal{T}}+{\mathbf{u}}_{H0}·{▽}_H{b}_0=\frac{\mathrm{E}}{\mathrm{Pr}}\frac{{\partial }^2{b}_1}{\partial z^2},$$

(22)

with boundary conditions

$$\frac{\mathrm{E}}{\mathrm{Pr}}\frac{\partial b_1}{\partial z}=Q_1,$$

(23)

at $z=\pm 1/2$. We begin by assuming that the fast time transients have decayed and depth average Equation (22) to get

$$\frac{\partial b_0}{\partial t}+{\overline{\mathbf{u}}}_{H0}·{▽}_H{b}_0=0.$$

(24)

Subtracting Equation (24) from Equation (22) gives an equation for $b_1^{\prime }$,

$${\mathbf{u}}_{H0}^{\prime }·{▽}_H{b}_0=\frac{\mathrm{E}}{\mathrm{Pr}}\frac{{\partial }^2{b}_1^{\prime }}{\partial z^2}.$$

(25)

The depth-dependent horizontal velocity is

$${\mathbf{u}}_{H0}^{\prime }=\frac{1}{\sqrt{\mathrm{E}}}\left[\left(K_0+\frac{z}{\sqrt{\mathrm{E}}}\right){\tau }_0-K_0^{″}\mathbf{k}\times {\tau }_0\right]-\sqrt{\mathrm{E}}\left[K_0^{″}{▽}_H{b}_0+K_0\mathbf{k}\times {▽}_H{b}_0\right],$$

(26)

and hence Equation (25) becomes

$$\begin{array}{c}\hfill \frac{\mathrm{E}}{\mathrm{Pr}}\frac{{\partial }^2{b}_1^{\prime }}{\partial z^2}=\frac{1}{\sqrt{\mathrm{E}}}\left[\left(K_0+\frac{z}{\sqrt{\mathrm{E}}}\right){\tau }_0-K_0^{″}\mathbf{k}\times {\tau }_0\right]·{▽}_H{b}_0-\sqrt{\mathrm{E}}{K}_0^{″}{\left|{▽}_H{b}_0\right|}^2,\end{array}$$

(27)

which can be integrated to give

$$\begin{array}{c}\hfill b_1^{\prime }=\frac{\mathrm{Pr}{Q}_1}{\mathrm{E}}z-\frac{\mathrm{Pr}}{\sqrt{\mathrm{E}}}\left[K_0^{″}{\tau }_0+K_0\mathbf{k}\times {\tau }_0\right]·{▽}_H{b}_0-\mathrm{Pr}\sqrt{\mathrm{E}}{K}_0{\left|{▽}_H{b}_0\right|}^2.\end{array}$$

(28)

Recall that $K_0$ is a function of $z/\sqrt{E}$ as described in (A7). The first term in $b_1^{\prime }$ describes the stratification induced by the surface buoyancy flux, the second term represents changes in buoyancy induced by wind forcing, and the third term represents the stratification induced by the TTW circulation [1].

4. The Evolution of the Depth-Averaged Buoyancy

In order to determine the evolution of the buoyancy valid to $O\left(\mathrm{Ro}\right)$ we solve for the leading order component of the buoyancy field, $b_0$, and can then use this in Equation (28) to find the $O\left(\mathrm{Ro}\right)$ contribution. Following [1] we derive an equation in terms of $b_0$ only by depth averaging the buoyancy evolution equation and substituting for ${\mathbf{u}}_{H0}$ and $b_1^{\prime }$. Depth averaging the full buoyancy equation (see Crowe and Taylor [1]) gives

$$\mathrm{Ro}\left[\frac{\partial \overline{b}}{\partial t}+{\overline{\mathbf{u}}}_H·{▽}_H\overline{b}+{▽}_H·\overline{{\mathbf{u}}_H^{\prime }{b}^{\prime }}\right]=\frac{{ϵ}^2\mathrm{E}}{\mathrm{Pr}}{▽}_H^2\overline{b},$$

(29)

which to leading order, $O\left(\mathrm{Ro}\right)$, reduces to Equation (24) and can be written using the depth-averaged velocity as

$$\frac{\partial b_0}{\partial t}+J\left({\overline{p}}_0,b_0\right)=0,$$

(30)

where the depth-averaged pressure acts as a stream function for the depth-averaged velocity. In order to determine the full evolution another equation is needed to determine ${\overline{p}}_0$. In Crowe and Taylor [25] we used a depth-averaged vorticity equation to close the system and obtain two coupled equations in ${\overline{p}}_0$ and $b_0$. Here we instead assume that we can neglect the Jacobian term, $J({\overline{p}}_0,b_0)$. One justification for this would be to assume that the characteristic along-front (y) length scales are large. In this case $b_0$ does not depend on the intermediate timescale, t and so we must look to the $O\left({\mathrm{Ro}}^2\right)$ evolution equation.

The $O\left({\mathrm{Ro}}^2\right)$ depth-averaged buoyancy equation is

$$\frac{\partial {\overline{b}}_1}{\partial t}+\frac{\partial b_0}{\partial T}+J({\psi }_1,b_0)+J({\overline{p}}_0,{\overline{b}}_1)+{▽}_H·\overline{\left[{\mathbf{u}}_{H0}^{\prime }{b}_1^{\prime }+{\mathbf{u}}_{H1}^{\prime }{b}_0^{\prime }\right]}=\frac{{ϵ}^2\mathrm{E}}{{\mathrm{Ro}}^2\mathrm{Pr}}{▽}_H^2{b}_0,$$

(31)

where ${\psi }_1$ is a streamfunction for the depth-averaged flow ${\overline{\mathbf{u}}}_{H1}$. We note that $b_0^{\prime }=0$ and by symmetry we can assume that ${\overline{b}}_1$ is zero to obtain an equation for the slow evolution of $b_0$

$$\frac{\partial b_0}{\partial T}+J({\psi }_1,b_0)+{▽}_H·\overline{{\mathbf{u}}_{H0}^{\prime }{b}_1^{\prime }}=\frac{{ϵ}^2\mathrm{E}}{{\mathrm{Ro}}^2\mathrm{Pr}}{▽}_H^2{b}_0.$$

(32)

Again, we assume that we can neglect the Jacobian term to obtain

$$\frac{\partial b_0}{\partial T}+{▽}_H·\overline{{\mathbf{u}}_{H0}^{\prime }{b}_1^{\prime }}=\frac{{ϵ}^2\mathrm{E}}{{\mathrm{Ro}}^2\mathrm{Pr}}{▽}_H^2{b}_0.$$

(33)

Note that while we have neglected all Jacobian terms, we retain the y derivatives in the flux term and the diffusion term of Equation (33). This is to give as much generality as possible since there is a much wider class of problems for which we can neglect the Jacobian derivatives than for which we can neglect the along-front second order derivatives. For example, in the case of the radially symmetric front considered by Shakespeare [26], the Jacobian terms would vanish. However, the along-front derivatives in the flux term and diffusion term would describe the frontal curvature which can have significant effects on the structure of the front [26]. Similarly, in the case of a fully 3D system, baroclinic instability can generate modes with small along-front scales. For these small y scales, the high number of derivatives in the flux term of Equation (33) can cause it to dominate the Jacobian term and have significant effects on the formation of the instability [25]. We will later neglect all y dependence by assuming that any dependence on the along-front scale y is weak due to along-front scales being much larger than cross-front scales. While this is a common approximation when studying frontal problems, it does prevent the growth of baroclinic instability [25,27]. The limitations of this approximation will be discussed in Section 8.

Equation (33) does not involve the $O\left(\mathrm{Ro}\right)$ contribution to the horizontal velocity, ${\mathbf{u}}_{H1}$, and hence (33) is a closed system involving $b_0$ and known quantities. The term ${\mathbf{u}}_{H0}^{\prime }{b}_1^{\prime }$ can be written as

$$\begin{array}{c}\hfill {\mathbf{u}}_{H0}^{\prime }{b}_1^{\prime }=\left(\frac{1}{\sqrt{\mathrm{E}}}\left[\left(K_0+\frac{z}{\sqrt{\mathrm{E}}}\right){\tau }_0-K_0^{″}\mathbf{k}\times {\tau }_0\right]-\sqrt{\mathrm{E}}\left[K_0^{″}{▽}_H{b}_0+K_0\mathbf{k}\times {▽}_H{b}_0\right]\right)\\ \hfill \times \left(\frac{\mathrm{Pr}{Q}_1}{\mathrm{E}}z-\frac{\mathrm{Pr}}{\sqrt{\mathrm{E}}}\left[K_0^{″}{\tau }_0+K_0\mathbf{k}\times {\tau }_0\right]·{▽}_H{b}_0-\mathrm{Pr}\sqrt{\mathrm{E}}{K}_0{\left|{▽}_H{b}_0\right|}^2\right),\end{array}$$

(34)

using the results for ${\mathbf{u}}_{H0}^{\prime }$ and $b_1^{\prime }$. Depth-averaging this result gives

$$\begin{array}{c}\hfill \overline{{\mathbf{u}}_{H0}^{\prime }{b}_1^{\prime }}=-\mathrm{Pr}[Q_1{\mathbf{C}}_{10}·{\tau }_0+Q_1{\mathbf{C}}_{01}·{▽}_H{b}_0+{\mathbf{C}}_{21}·\left[{\tau }_0{\tau }_0{▽}_H{b}_0\right]\\ \hfill +{\mathbf{C}}_{12}·\left[{\tau }_0{▽}_H{b}_0{▽}_H{b}_0\right]+{\mathbf{C}}_{03}·{▽}_H{b}_0{\left|{▽}_H{b}_0\right|}^2],\end{array}$$

(35)

where the tensors ${\mathbf{C}}_{10}$, ${\mathbf{C}}_{01}$, ${\mathbf{C}}_{21}$, ${\mathbf{C}}_{12}$ and ${\mathbf{C}}_{03}$ are defined in Appendix B and are functions only of $\mathrm{E}$. Note that ${\mathbf{C}}_{21}$ and ${\mathbf{C}}_{12}$ are fourth-rank tensors so here we use “·” to denote the contraction ${\left[\mathbf{A}·\left[\mathbf{bcd}\right]\right]}_i=A_{ijkl}{b}_j{c}_k{d}_l$. Substituting this result into Equation (33) gives the evolution equation for the depth-averaged buoyancy,

$$\begin{array}{c}\hfill \frac{\partial b_0}{\partial T}=\mathrm{Pr}{▽}_H·\left(Q_1{\mathbf{C}}_{10}·{\tau }_0+Q_1{\mathbf{C}}_{01}·{▽}_H{b}_0+{\mathbf{C}}_{21}·\left[{\tau }_0{\tau }_0{▽}_H{b}_0\right]\right.\\ \hfill \left.+{\mathbf{C}}_{12}·\left[{\tau }_0{▽}_H{b}_0{▽}_H{b}_0\right]+{\mathbf{C}}_{03}·{▽}_H{b}_0{\left|{▽}_H{b}_0\right|}^2\right)+\frac{{ϵ}^2\mathrm{E}}{{\mathrm{Ro}}^2\mathrm{Pr}}{▽}_H^2{b}_0.\end{array}$$

(36)

Note that Equation (36) is valid for ${\tau }_0={\tau }_0(x,y)$ and $Q_1=Q_1(x,y)$ however, in the case of constant ${\tau }_0$ and $Q_1$ the term $Q_1{\mathbf{C}}_{10}·{\tau }_0$ is constant and hence can be neglected.

Physically, the leading-order velocity perturbation $u^{\prime }$ consists of a TTW component, $u_{TTW}\propto {▽}_H{b}_0$, plus an Ekman velocity driven by the wind stress, $u_{WS}\propto {\tau }_0$. Since the leading order velocity controls the vertical buoyancy gradient, $N^2$, through an advection–diffusion balance, $N^2$ is determined by the advection of the background horizontal buoyancy gradient and the surface buoyancy flux. Therefore $N^2$ consists of three terms, $N_{TTW}^2\propto {\left({▽}_H{b}_0\right)}^2$ corresponding to the stratification maintained by TTW flow, $N_{WS}^2\propto {\tau }_0{▽}_H{b}_0$ corresponding to the stratification maintained by the wind-driven Ekman flow and $N_Q^2\propto Q_1$ describing the buoyancy flux induced stratification.

From Equation (33) we can see that the correlation between the stratified component of the buoyancy and the leading order velocity drives the evolution of the leading order buoyancy field via shear dispersion [1]. This correlation term describes the horizontal flux of depth-averaged buoyancy and consists of six terms due to the advection of each of the stratification terms $N_{TTW}^2$, $N_{WS}^2$ and $N_Q^2$ by the velocity terms $u_{TTW}$ and $u_{WS}$. We therefore have terms proportional to ${\left({▽}_H{b}_0\right)}^3$, ${\left({▽}_H{b}_0\right)}^2{\tau }_0$, ${▽}_H{b}_0{\tau }_0^2$, ${▽}_H{b}_0{Q}_1$ and $Q_1{\tau }_0$ as shown in Equation (36).

4.1. The y Independent Case

Our results can be simplified further if we neglect the remaining y dependence from this point onward by assuming that any dependence on the along-front distance y is small. In this case, Equation (36) simplifies to

$$\begin{array}{c}\hfill \frac{\partial b_0}{\partial T}=\mathrm{Pr}\frac{\partial }{\partial x}\left[C_0(Q_1,{\tau }_0)+C_1(Q_1,{\tau }_0)\frac{\partial b_0}{\partial x}+C_2\left({\tau }_0\right){\frac{\partial b_0}{\partial x}}^2+C_3{\frac{\partial b_0}{\partial x}}^3\right]+D_H\frac{{\partial }^2{b}_0}{\partial x^2},\end{array}$$

(37)

where

$$D_H=\frac{{ϵ}^2\mathrm{E}}{{\mathrm{Ro}}^2\mathrm{Pr}},$$

(38)

and the $C_i$ are given in Appendix B. $C_0$ is linear in ${\tau }_0$ and $Q_1$, $C_1$ consists of a part linear in $Q_1$ plus a part quadratic in ${\tau }_0$ and $C_2$ is linear in ${\tau }_0$. Note that $C_3$ does not depend on the wind stress or buoyancy flux and hence is the same as in the unforced case [1]. Spatial variations in the wind stress and buoyancy flux can cause the buoyancy field to evolve from an initially constant buoyancy through the $C_0$ term which we anticipate could drive frontogenesis. Note that when any one coefficient is dominant we can find an approximate solution to Equation (37) using the similarity solutions in Appendix C.

For simplicity, we now reduce to the case where the wind stress, ${\tau }_0$, and buoyancy flux, $Q_1$, are spatially uniform and hence constant. The limitations of this assumption are discussed in Section 8. Equation (37) becomes

$$\frac{\partial b_0}{\partial T}=\mathrm{Pr}\frac{\partial }{\partial x}\left[C_1(Q_1,{\tau }_0)\frac{\partial b_0}{\partial x}+C_2\left({\tau }_0\right){\frac{\partial b_0}{\partial x}}^2+C_3{\frac{\partial b_0}{\partial x}}^3\right]+D_H\frac{{\partial }^2{b}_0}{\partial x^2},$$

(39)

and we write

$${\tau }_0={\tau }_0(cos\theta ,sin\theta ,0),$$

(40)

where ${\tau }_0>0$ describes the magnitude of the wind stress and $\theta$ describes the angle from the positive cross-front direction in the right-handed sense. The coefficients, $C_1$, $C_2$ and $C_3$ can now be written as

$$C_1=Q_1\left[\overline{K_0^{\prime 2}}+\overline{K_0^{‴2}}\right]+\frac{1}{\mathrm{E}}\left[\overline{K_0^{‴2}}{cos}^2\theta +2\overline{K_0^{″2}}cos\theta sin\theta +\overline{K_0^{\prime 2}}{sin}^2\theta \right]{\tau }_0^2,$$

(41)

$$C_2=-2\left[\overline{K_0^{″2}}cos\theta +\overline{K_0^{\prime 2}}sin\theta \right]{\tau }_0,$$

(42)

and

$$C_3=\mathrm{E}\overline{K_0^{\prime 2}},$$

(43)

using the results from Appendix B. We now group the buoyancy flux term from $C_1$ with $D_H$ to write Equation (39) as

$$\frac{\partial b_0}{\partial T}=\mathrm{Pr}\frac{\partial }{\partial x}\left[C_1^{\prime }(E,\theta ){\tau }_0^2\frac{\partial b_0}{\partial x}+C_2^{\prime }(E,\theta ){\tau }_0{\frac{\partial b_0}{\partial x}}^2+C_3{\frac{\partial b_0}{\partial x}}^3\right]+D\frac{{\partial }^2{b}_0}{\partial x^2},$$

(44)

where

$$C_1^{\prime }(E,\theta )=\frac{1}{\mathrm{E}}\left[\overline{K_0^{‴2}}{cos}^2\theta +2\overline{K_0^{″2}}cos\theta sin\theta +\overline{K_0^{\prime 2}}{sin}^2\theta \right],$$

(45)

$$C_2^{\prime }(E,\theta )=-2\left[\overline{K_0^{″2}}cos\theta +\overline{K_0^{\prime 2}}sin\theta \right],$$

(46)

and

$$D=D_H+D_Q=\frac{{ϵ}^2\mathrm{E}}{{\mathrm{Ro}}^2\mathrm{Pr}}+\mathrm{Pr}{Q}_1\left[\overline{K_0^{\prime 2}}+\overline{K_0^{‴2}}\right].$$

(47)

Note that $C_1^{\prime }$ and $C_3$ are positive though $C_2$ and D may be negative.

4.2. Frontal Evolution

We now study the evolution of a simple y independent front and find that wind stress and buoyancy flux may drive frontal spreading or sharpening depending on the parameter values. W begin by writing Equation (44) as

$$\frac{\partial b_0}{\partial T}={\kappa }_0\frac{{\partial }^2{b}_0}{\partial x^2},$$

(48)

where

$${\kappa }_0=\mathrm{Pr}\left[C_1^{\prime }(E,\theta ){\tau }_0^2+2C_2^{\prime }(E,\theta ){\tau }_0\frac{\partial b_0}{\partial x}+3C_3\left(E\right){\frac{\partial b_0}{\partial x}}^2\right]+D.$$

(49)

If ${\kappa }_0$ is constant, then Equation (48) reduces to the linear diffusion equation, although since ${\kappa }_0$ depends on the horizontal buoyancy gradient it will generally not be constant and Equation (48) will be nonlinear.

It is difficult to obtain exact solutions to Equation (48) due to the nonlinearity in $b_0$ (implicit through ${\kappa }_0$). In Section 5, we will discuss numerical solutions to this equation. However, we can use insights from a linear diffusion equation with constant diffusivity to qualitatively predict the evolution of the front. In the linear diffusion equation, when ${\kappa }_0$ is a positive constant, the solutions will exhibit frontolysis through diffusive spreading. On the other hand, when ${\kappa }_0<0$, the linear diffusion equation is ill-posed and solutions will exhibit sensitive dependence on the initial conditions and frontogenesis will occur through the development of step discontinuities in $b_0\left(x\right)$.

Using the analogy of the linear diffusion equation we anticipate that frontolysis will occur when ${\kappa }_0>0$ for all x. Note that the nonlinearity introduced through ${\kappa }_0$ and its dependence on x mean that the behaviour is not guaranteed to hold. For example, Taylor and Zhou [28] discussed a mechanism whereby steps can develop in a one-dimensional diffusion equation when the diffusivity depends on the buoyancy gradient. Nevertheless, our numerical solutions to Equation (48) always exhibit spreading when ${\kappa }_0>0$.

It is therefore of interest to determine the sign of ${\kappa }_0$ under various conditions, and this can be done by comparing the sizes of the terms in Equation (48). If we assume that the front is oriented such that $\partial b_0/\partial x\ge 0$ then ${\kappa }_0$ can only be negative if $C_2^{\prime }{\tau }_0$ or D are negative and sufficiently large compared to the other coefficients. Figure 3 shows the coefficients $C_1^{\prime }$ and $C_2^{\prime }$ as functions of $\theta$ and $\mathrm{E}$. Recall that $C_3$ is a positive function that is linearly proportional to the Ekman number. The coefficient $C_1^{\prime }$ is positive for all values of $\theta$ and E while the largest amplitude negative values of $C_2^{\prime }$ occur for small E and $0<\theta <{180}^{\circ }$. Note that $C_1^{\prime }\sim {\mathrm{E}}^{-1}$, $C_2^{\prime }\sim 1$ and $C_3\sim \mathrm{E}$ for small $\mathrm{E}$, hence for small Ekman number and ${\tau }_0=O\left(1\right)$ we expect the front to spread through diffusion since the $C_1^{\prime }{\tau }_0^2$ term will be dominant in that limit. For intermediate values of $\mathrm{E}$ it may be possible for a negative $C_2^{\prime }{\tau }_0$ or D to temporarily sharpen the frontal gradients. However, since the front sharpens, the cubic term in Equation (44) will become large and, since $C_3>0$, the gradient will not increase indefinitely.

From Equation (49), we see that ${\kappa }_0=0$ when

$$\frac{\partial b_0}{\partial x}=B_c^{\pm },$$

(50)

where

$$B_c^{\pm }=\frac{-C_2^{\prime }{\tau }_0\pm \sqrt{C_2^{\prime 2}{\tau }_0^2-3C_3(C_1^{\prime }{\tau }_0^2+D/\mathrm{Pr})}}{3C_3},$$

(51)

assuming that these roots are real. Negative values of ${\kappa }_0$ occur when $B_c^{-}<\partial b_0/\partial x<B_c^{+}$ and ${\kappa }_0$ is positive when $\partial b_0/\partial x<B_c^{-}$ and $\partial b_0/\partial x>B_c+$.

In many cases, we expect fronts to evolve towards an equilibrium state with a buoyancy gradient $\partial b/\partial x=B_c^{+}$. This follows if we first assume that ${\kappa }_0>0$ leads to the spreading of the front so that $\partial b_0/\partial x$ decreases in time. If the buoyancy gradient $\partial b/\partial x>B_c^{+}$, then we anticipate that the buoyancy gradient will decrease until $\partial b/\partial x=B_c^{+}$ where ${\kappa }_0=0$ and an equilibrium state is achieved. On the other hand, if ${\kappa }_0<0$, existing fronts can intensify through reverse diffusion and $\partial b_0/\partial x$ will increase until $\partial b/\partial x=B_c^{+}$. In both cases, $\partial b/\partial x=B_c^{+}$ acts as an attracting equilibrium state. Figure 4, which is discussed below, illustrates this argument in the case of no surface heat flux.

To make this idea more precise, consider an initial non-dimensional buoyancy gradient, ${\left.\partial b/\partial x\right|}_{x=0,t=0}=1$. If there are no positive real roots, then ${\kappa }_0$ is initially positive and remains positive throughout the evolution, and we anticipate that the front will spread. If there are positive real roots then we have three possible cases. Firstly if $B_c^{-}>1$ then ${\kappa }_0$ is initially positive and the front spreads indefinitely with ${\kappa }_0$ increasing throughout the evolution. Secondly if $B_c^{-}<1<B_c^{+}$ then ${\kappa }_0$ is initially negative and the front sharpens towards $\partial b_0/\partial x=B_c^{+}$; as this critical gradient is approached ${\kappa }_0\to 0$ so the front will sharpen towards an equilibrium gradient. Finally if $B_c^{+}<1$ then initially ${\kappa }_0>0$ so the front will spread towards $\partial b_0/\partial x=B_c^{+}$ and again approach an equilibrium width.

Consider a situation where $C_2^{\prime }<0$ and the second term in Equation (48) causes ${\kappa }_0<0$. Physically, this can be interpreted as a wind-driven Ekman flow that opposes the cross-front TTW flow in such a way that the overall buoyancy flux is reversed, resulting in an up-gradient flux. Similarly, when $D<0$ (the diffusivity associated with the last term in Equation (44)) and the surface cooling is sufficiently strong to create regions of unstable stratification also leads to an up-gradient flux of buoyancy.

Note that it is possible for density inversions to occur in cases where the wind stress opposes the cross-front TTW flow or when the buoyancy flux is negative. Theoretical studies by Thomas and Rhines [29] and Benthuysen and Thomas [30] using diffusive parametrisations have also noted these inversions. Below we will consider two limits where the solutions can be simplified, specifically the limits of a weak surface heat flux and a weak wind stress.

4.3. Weak Surface Heat Flux

We begin by considering the case of a weak surface heat flux (and small D in Equation (44)) where the effects of buoyancy flux and horizontal diffusion can be neglected. Figure 5 shows ${\kappa }_0$ as a function of $\theta$ and $log\mathrm{E}$ (the base 10 logarithm is used for the figures throughout) for $\partial b_0/\partial x=1$ and $D=0$. Panels (a)–(d) correspond to four different values of the wind stress, ${\tau }_0$ and the white contour marks ${\kappa }_0=0$. The size and shape of the region where ${\kappa }_0<0$ is strongly dependent on ${\tau }_0$. For a particular Ekman number, ${\kappa }_0<0$ typically occurs for a limited range of wind stress orientations. This is a consequence of the coupling between the along-front and cross-front velocities through the vertical mixing and rotation terms. When the wind stress is weak (and ${\tau }_0$ is small), the largest negative value of ${\kappa }_0$ occurs for small Ekman number and when $\theta \approx {90}^{\circ }$, i.e., when the wind stress is aligned with the thermal wind in a down-front orientation. In this orientation, the cross-front Ekman transport acts to advect water from the dense side of the front over the top of the water from the light side of the front [10].

Conversely, when the wind stress is strong, the largest negative value of ${\kappa }_0$ occurs for large Ekman numbers. In this limit when vertical mixing is strong, the flow is more closely aligned with the wind stress and the strongest cross-front flow occurs when the wind is closer to the cross-front direction. As a result, the largest negative value of ${\kappa }_0$ occurs for smaller values of $\theta$, when the wind stress is more closely aligned with the buoyancy gradient.

For $D=0$, the diffusivity, ${\kappa }_0$, can be written

$${\kappa }_0=Pr{\tau }_0^2\left[C_1^{\prime }+2C_2^{\prime }\left[\frac{1}{{\tau }_0}\frac{\partial b_0}{\partial x}\right]+3C_3{\left[\frac{1}{{\tau }_0}\frac{\partial b_0}{\partial x}\right]}^2\right].$$

(52)

Note that the sign of ${\kappa }_0$ is a function only of three parameters; $\mathrm{E}$, $\theta$ and $[\partial b_0/\partial x]/{\tau }_0$. Figure 5 shows ${\kappa }_0/\left[\mathrm{Pr}{\tau }_0^2\right]$ as a function of $[\partial b_0/\partial x]/{\tau }_0$ for $D=0$ and $\theta$ = 45°, $135$°. These values of $\theta$ were chosen due to having regions of both positive and negative ${\kappa }_0$. For the purpose of discussion, points A–D correspond to particular initial conditions. For cases starting at points A and B, ${\kappa }_0>0$ initially, and the front will spread such that $d{b}_0/dx<0$ as indicated by the arrows. If E is constant, we anticipate that the fronts in these cases will spread indefinitely since ${\kappa }_0$ never reaches zero. In case C ${\kappa }_0>0$ initially. However, as $d{b}_0/dx$ decreases for a constant E, ${\kappa }_0$ will reach zero at an equilibrium buoyancy gradient. Case D has ${\kappa }_0<0$ initially and we anticipate buoyancy gradients to sharpen until ${\kappa }_0=0$.

Following the thought experiment described above, we might anticipate that fronts will evolve until they reach a state of equilibrium where ${\kappa }_0=0$. The equilibrium buoyancy gradient can be found by setting ${\kappa }_0=0$ in Equation (52) and solving for $\partial b_0/\partial x$. The equilibrium value of $\partial b_0/\partial x$, denoted $B_c^{+}$ can then be written

$$B_c^{+}=\left[\frac{-C_2^{\prime }+\sqrt{C_2^{\prime 2}-3C_1^{\prime }{C}_3}}{3C_3}\right]{\tau }_0.$$

(53)

Note that the equilibrium gradient is proportional to the wind stress, so stronger wind stress will sustain a stronger equilibrium gradient. Figure 6 shows $log[B_c^{+}/{\tau }_0]$ as a function of $\theta$ and $log\mathrm{E}$. The white region corresponds to the range of parameters where an equilibrium buoyancy gradient does not exist. Note that the equilibrium buoyancy gradient decreases with increasing Ekman number. This suggests that vertical mixing counteracts wind-driven frontogenesis and that this balance holds for weaker buoyancy gradients when mixing is strong.

4.4. Weak Wind Stress

For small ${\tau }_0$, the evolution of the leading order buoyancy, Equation (44), reduces to

$$\frac{\partial b_0}{\partial T}=\mathrm{Pr}{C}_3\frac{\partial }{\partial x}\left[{\frac{\partial b_0}{\partial x}}^3\right]+D\frac{{\partial }^2{b}_0}{\partial x^2}.$$

(54)

This is the nonlinear Erdogan–Chatwin equation [31,32] discussed in [1] with the horizontal diffusivity, $D_H$, replaced by the combined diffusivity, $D=D_H+D_Q$. The surface buoyancy flux modifies the coefficient of the linear diffusion term and hence increases the spreading rate in the case of surface heating and reduces the spreading rate for surface cooling.

When D is positive ${\kappa }_0$ is strictly positive and is given by

$${\kappa }_0=D+3\mathrm{Pr}{C}_3{\frac{\partial b_0}{\partial x}}^2.$$

(55)

In this case, Equation (54) can be solved using similarity solutions and describes a spreading front [1,32]. Furthermore, if D is large compared with $\mathrm{Pr}{C}_3$, Equation (54) reduces to a linear diffusion equation with the similarity solution

$$b_0=erf\left[\frac{x}{\sqrt{4D(T+T_0)}}\right].$$

(56)

Over long times, the linear terms in Equation (54) dominate the nonlinear term and the front spreads as $T^{1/2}$ and approaches this similarity solution.

If D is negative, corresponding to surface cooling that is sufficiently strong to overcome the spreading effects of horizontal diffusion, we may write

$${\kappa }_0=-\left|D\right|+3\mathrm{Pr}{C}_3{\frac{\partial b_0}{\partial x}}^2,$$

(57)

so ${\kappa }_0/\left|D\right|$ may be written as a function of $[\partial b_0/\partial x]/\sqrt{\left|D\right|}$ as

$$\frac{{\kappa }_0}{\left|D\right|}=-1+3\mathrm{Pr}{C}_3{\left(\frac{1}{\sqrt{\left|D\right|}}\frac{\partial b_0}{\partial x}\right)}^2,$$

(58)

where we recall that $C_3=C_3\left(\mathrm{E}\right).$ Therefore, there will always exist a region where ${\kappa }_0<0$ and we have equilibrium buoyancy gradient (where ${\kappa }_0=0$) of

$$B_c^{+}=\sqrt{\frac{\left|D\right|}{3C_3\mathrm{Pr}}}.$$

(59)

In this case we have $B_c^{-}=-B_c^{+}<0$ so, as discussed above (in the text below Equation (51)), we do not anticipate indefinite frontolysis. Instead, we expect that the front will either spread ($B_c^{+}<1$) or sharpen ($B_c^{+}>1$) until $\partial b_0/\partial x$ reaches the equilibrium value, $B_c^{+}$. Figure 7a shows ${\kappa }_0$ as a function of negative D and $\mathrm{E}$ for $\mathrm{Pr}=1$ and $\partial b_0/\partial x=1$. The region where ${\kappa }_0<0$ occurs to the left of the white line. Figure 7b shows ${\kappa }_0/\left|D\right|$ as a function of $[\partial b_0/\partial x]/\sqrt{\left|D\right|}$ and $\mathrm{E}$ for $\mathrm{Pr}=1$ and $D<0$. Again the white contour denotes ${\kappa }_0=0$ and hence corresponds to the equilibrium gradient $B_c^{+}$. Similarly to Figure 5, point E corresponds to a front spreading towards equilibrium whereas point F corresponds to a front sharpening towards equilibrium.

5. Numerical Simulations

In this section, we describe numerical solutions of the evolution equation for the leading order buoyancy, Equation (44). The objective of these simulations is to test our predictions about frontal spreading or sharpening based on the sign of ${\kappa }_0$ and to examine the behavior in more detail. We use a pseudo-spectral method to evaluate the spatial derivatives and a fourth-order Runge-Kutta scheme for the time stepping. We assume that the leading order buoyancy, $b_0$, is a function of x and t and we use a spatial grid with $N_x=256$ grid points.

We use a numerical domain of nondimensional width $2l$ with $x\in [-l,l]$ and an initial buoyancy profile of

$$b_0(x,0)=erf\left(\frac{\sqrt{\pi }x}{2}\right),$$

(60)

so

$$\frac{\partial b_0}{\partial x}(0,0)=1.$$

(61)

Note that the horizontal lengthscale L used to nondimensionalise $(x,y)$ has been arbitrary up to this point; setting $\partial b_0/\partial x=1$ in the center of the front is equivalent to choosing L to be the frontal width $L=M^2/\Delta b$ so is not an assumption which imposes any restrictions on our method. In order to make the domain periodic and use a pseudo-spectral method we subtract a linear profile from $b_0$ and solve for ${\tilde{b}}_0$ where

$${\tilde{b}}_0(x,T)=b_0(x,T)-\frac{x}{l}.$$

(62)

We now present three different cases of sharpening and spreading fronts that can arise from different choices of parameters. Figure 4 is used to help select the parameters for the different cases.

For each simulation, we plot $b_0$ as a function of x and T showing any spreading, $b_0$ as functions of x for several values of T showing the shape of the cross-front buoyancy profile, ${\kappa }_0$ as a function of x for several T values to show the cross-front variation in spreading rate, and ${\kappa }_0$ as a function of $\partial b_0/\partial x$ to see how the spreading rate changes as the front evolves and if there exists an equilibrium buoyancy gradient, $B_c^{+}$.

5.1. A Spreading Front—Indefinite Spreading

We begin by considering a case where there is no real positive equilibrium gradient, $B_c^{+}$. This case has a surface wind stress partially opposing the thermal wind and no surface heat flux. Specifically, the parameters are $\mathrm{E}=0.1$, ${\tau }_0=0.1$, $\theta =225$°, $\mathrm{Pr}=1$, $Q_1=0$ and $D_H=2.5\times {10}^{-4}$. The corresponding coefficients are $C_1^{\prime }{\tau }_0^2=0.041$, $C_2^{\prime }{\tau }_0=0.059$, $C_3=0.021$ and $D=2.5\times {10}^{-4}$. The numerical solution is shown in Figure 8. The horizontal diffusivity can be much different inside the front than it is towards the edges due to the dependence of ${\kappa }_0$ on $\partial b_0/\partial x$. In this case, ${\kappa }_0$ is larger inside the front than at the edges, leading to an approximately linear cross-front buoyancy profile near $x=0$. We note that it is possible to have cases with $C_2^{\prime }<0$ where ${\kappa }_0$ is smaller inside the front than at the edges.

Qualitatively similar behavior also occurs in the absence of wind stress and a small or zero heat flux. When both the wind stress and the surface heat flux are zero, the solution limits to that described in [1]. In these cases the cubic nonlinear term in Equation (44) is dominant inside the front and the solution resembles the similarity solution, $F_3$ (see Appendix C).

5.2. A Spreading Front—Equilibrium Width

We now consider the case where there exists a real, positive equilibrium gradient with $B_c<1$. Specifically, consider a wind oriented ${25}^{\circ }$ away from the cross-front direction in the absence of a surface heat flux with parameters $\mathrm{E}=0.1$, ${\tau }_0=0.1$, $\theta =115$°, $\mathrm{Pr}=1$, $Q_1=0$ and $D_H=2.5\times {10}^{-4}$. The corresponding coefficients are $C_1^{\prime }{\tau }_0^2=0.0050$, $C_2^{\prime }{\tau }_0=-0.020$, $C_3=0.021$ and $D=2.5\times {10}^{-4}$ and the numerical solution is shown in Figure 9. Near $x=0$ the front spreads towards a constant gradient. Initially, there are regions in which ${\kappa }_0<0$, where the gradient sharpens to match that of the interior region, forming cusps between the spreading edges and the approximately linear interior profile. Once the cross-front buoyancy profile reaches the equilibrium gradient, $B_c^{+}$, we expect it to remain constant in time near $x=0$.

5.3. Sharpening Front

Finally, we consider the case of a real, positive equilibrium buoyancy gradient with $B_c^{+}>1$ with parameters $\mathrm{E}=0.1$, ${\tau }_0=0.5$, $\theta =124$°, $\mathrm{Pr}=1$, $Q_1=0$ and $D_H=2.5\times {10}^{-4}$. The corresponding coefficients are $C_1^{\prime }{\tau }_0^2=0.041$, $C_2^{\prime }{\tau }_0=-0.058$, $C_3=0.021$ and $D=2.5\times {10}^{-4}$ and the numerical solution is shown in Figure 10. At the centre of the front, ${\kappa }_0<0$ initially, and the cross-front buoyancy profile sharpens into a series of steps. Towards the edges of the front ${\kappa }_0>0$ and the front spreads.

We note that due to the ill-posed nature of backwards diffusion, numerical solutions of Equation (44) with ${\kappa }_0<0$ in some regions will not be accurate. In particular, spikes in the buoyancy gradient start to develop, and increasing the grid resolution will only result in more small-scale oscillations as energy is transferred into the highest frequency modes. However, we expect that the qualitative behaviour will be similar to that observed here with a sharpened frontal region and spreading edges.

6. Wind Stress and Buoyancy Flux Driven Turbulence

So far, we have assumed that the Ekman number describing the strength of the turbulence is independent of the surface wind stress and buoyancy flux. In reality, we expect the strength of boundary layer turbulence to be governed by surface-forcing so the Ekman number will depend on $\tau$ and Q [16,20]. We now consider the cases where the turbulence is wind and heat flux-driven, and relate the Ekman number and spreading coefficients to the strength of the forcing.

6.1. Wind Stress

Assuming the the wind stress is dominant, the Ekman number can be related to the wind stress using the scaling described by Enriquez and Taylor [20]

$$\mathrm{E}\simeq \frac{C_{\tau }}{fH}\sqrt{\frac{{\tau }^{*}}{{\rho }_0}},$$

(63)

with an empirical scaling constant $C_{\tau }\approx 0.02$, mixed layer depth H, Coriolis parameter f, dimensional stress ${\tau }^{*}$ and density ${\rho }_0$. We note that Equation (63) was derived for a different flow configuration to the one used here, in particular forcing was only applied at the top surface. We expect the scaling $\nu \sim H\sqrt{{\tau }^{*}/{\rho }_0}$ (and hence $\mathrm{E}\sim \sqrt{{\tau }^{*}/{\rho }_0}/\left(fH\right)$) to still hold for our configuration on dimensional grounds though the coefficient will likely be different. However, as Equation (63) will be purely used for illustrative purposes and any estimates are qualitative we retain the coefficient from Enriquez and Taylor [20].

A typical maximum wind stress of ${\tau }^{*}=0.1{\mathrm{Nm}}^{-2}$ corresponds to $\mathrm{E}\simeq 0.1$. From Equation (9) we note that the non-dimensional wind stress, $\tau$, depends on the horizontal buoyancy gradient, $M^2$, so using Equation (63) we have

$$\mathrm{E}\simeq \frac{C_{\tau }M}{f}\sqrt{{\tau }_0},$$

(64)

and we note that because of the choice of nondimensional variables, $C_{\tau }M/f$ and ${\tau }_0$ depend on $M^2$, while $\mathrm{E}$ does not. In the absence of a buoyancy flux the coefficients $C_1$, $C_2$ and $C_3$ scale as

$$C_1\sim M^{-1}{\tau }_0^{3/2},C_2\sim {\tau }_0,C_3\sim M{\tau }_0^{1/2},$$

(65)

for small $\mathrm{E}$; these results can be derived using the Ekman number dependence given in Appendix B. For large ${\tau }_0$, $C_1$ will be dominant and the front will spread indefinitely, similarly for very small ${\tau }_0$ the coefficient $C_3$ will be dominant and the front will also spread. Therefore, for frontal sharpening or spreading towards an equilibrium width, we need an intermediate value of the wind stress.

6.2. Buoyancy Flux

Similarly, we can relate the Ekman number to convection driven by a surface heat flux, $Q^{*}$, in the case of surface cooling. As described in Taylor and Ferrari [16], the turbulent Ekman number scales according to

$$\mathrm{E}\simeq \frac{C_Q}{f}{\left(\frac{|Q^{*}|\alpha g}{c_p{\rho }_0{H}^2}\right)}^{1/3},$$

(66)

where H is the depth of the convective mixed layer, ${\rho }_0$ is the water density, $c_p$ is the heat capacity, $\alpha$ is the thermal expansion coefficient, $C_Q=0.1$ is an empirical scaling constant and g is gravitational acceleration. Using a typical heat flux of $|Q^{*}|=100{\mathrm{Wm}}^{-2}$ we obtain an Ekman number of $E\sim 0.1$. Using Equation (8) and $Q=\mathrm{Ro}{Q}_1$ the Ekman number can be written in terms of $Q_1$ as

$$\mathrm{E}\simeq C_Q{\left(\frac{M^4\left|Q_1\right|}{f^4}\right)}^{1/3},$$

(67)

where again we note that $M^4/f^4$ and $Q_1$ both depend on $M^2$ while $\mathrm{E}$ does not. Since a positive heat flux does not generate convection, we expect that in the case of positive $Q^{*}$ the Ekman number will be primarily driven by the wind stress.

We now return to the case of arbitrary heat flux and consider the magnitude of the buoyancy flux induced spreading/sharpening compared to the spreading/sharpening effects of horizontal diffusion and wind stress. The diffusivity associated with the surface buoyancy flux for both positive and negative heat flux, $D_Q$ in Equation (44), scales as

$$D_Q=\mathrm{Pr}{Q}_1\left[\overline{K_0^{\prime 2}}+\overline{K_0^{″\prime 2}}\right]\sim Q_1\sim \frac{\alpha gf{Q}^{*}}{c_p{\rho }_0{H}^2{M}^4},$$

(68)

while the explicit horizontal diffusivity added to the buoyancy equation scales as

$$D_H=\frac{{ϵ}^2\mathrm{E}}{{\mathrm{Ro}}^2\mathrm{Pr}}\sim \frac{f^4}{M^4}\mathrm{E},$$

(69)

where $\mathrm{E}$ is driven by either the wind stress or heat flux, $\mathrm{Ro}=\left(H{M}^2\right)/\left(L{f}^2\right)$ and $\mathrm{Pr}=1$. Therefore, the ratio $D_Q/D_H$ can be written

$$\frac{D_Q}{D_H}\sim \frac{\alpha g{Q}^{*}}{c_p{\rho }_0{H}^2{f}^3\mathrm{E}}.$$

(70)

Typical parameters, e.g., $\mathrm{E}\sim 0.1$, $f={10}^{-4}{\mathrm{s}}^{-1}$, $c_p=4.18\times {10}^3{\mathrm{JKg}}^{-1}{\mathrm{K}}^{-1}$, $\alpha =1.67\times {10}^{-4}{\mathrm{K}}^{-1}$, $Q^{*}=\pm 100{\mathrm{Wm}}^{-2}$ and $H=100\mathrm{m}$ give

$$\left|\frac{D_Q}{D_H}\right|\sim 100.$$

(71)

This implies that in the absence of significant wind stress, the influence of the surface buoyancy flux on the evolution of the front is expected to be larger than that of horizontal mixing.

The constant term in ${\kappa }_0$ given by $Pr{C}_1^{\prime }{\tau }_0^2$ (see Equation (49)) describes horizontal diffusion due to the wind stress. Therefore the relative strength of the wind and surface heat flux-driven diffusion is given by the ratio

$$\frac{Pr{C}_1^{\prime }{\tau }_0^2}{D_Q}\sim \frac{Pr{C}_1^{\prime }{c}_p{\tau }_0^{*2}}{{\rho }_0{H}^2\alpha gf{Q}^{*}}\sim \pm C_1^{\prime },$$

(72)

using ${\tau }^{*}=0.1{\mathrm{Nm}}^{-2}$. In general $C_1^{\prime }>1$ for small $\mathrm{E}$, so we anticipate that the wind stress will play a larger role than $D_Q$ in the evolution of the front. Therefore, we predict the surface cooling will primarily act to modulate mixed layer turbulence and hence the Ekman number through Equation (66). Additionally surface heating or cooling may have a spreading or sharpening effect, respectively in the case of very weak winds.

7. Equilibrium Width Due to Stress-Ekman Balance

In this section, we discuss the scaling for the dimensional equilibrium horizontal buoyancy gradient and the implications for typical frontal widths. As discussed in the introduction, these results have implications for parameterisations of submesoscale processes which typically require the horizontal buoyancy gradient as a parameter (e.g., [19,33]).

The expression for the equilibrium horizontal buoyancy gradient given in Equation (53) can be written in terms of dimensional variables,

$$\frac{\partial b^{*}}{\partial x^{*}}=\mathcal{C}(\mathrm{E},\theta ,\mathcal{D})\frac{{\tau }^{*}}{{\rho }_0{H}^2},$$

(73)

where

$$\mathcal{C}(\mathrm{E},\theta ,\mathcal{D})=\frac{-C_2^{\prime }(\mathrm{E},\theta )+\sqrt{{\left(C_2(\mathrm{E},\theta )\right)}^2-3\left(C_1^{\prime }(\mathrm{E},\theta )+\mathcal{D}\right)C_3\left(\mathrm{E}\right)}}{3C_3\left(\mathrm{E}\right)},$$

(74)

and

$$\mathcal{D}=\frac{D_H+D_Q}{Pr{\tau }_0^2}=\left(\frac{{\rho }_0^2{H}^4}{{\tau }^{*2}}\right)\left(\frac{f^4\mathrm{E}}{P{r}^2}+\frac{\alpha gf{Q}^{*}}{c_p{\rho }_0{H}^2}\left[\overline{K_0^{\prime 2}}+\overline{K_0^{‴2}}\right]\right).$$

(75)

Here $\mathcal{D}/C_1^{\prime }$ describes the strength of the combined horizontal diffusivity due to horizontal mixing and the surface heat flux compared with the strength of the diffusivity associated with the wind stress. We expect this ratio to be small as explained in Section 6.2.

For illustration, consider a simple case where $\mathcal{D}=0$ and $\theta$ = 90°, i.e., no explicit horizontal diffusion in the buoyancy equation, no surface heat flux, and a wind stress aligned with the horizontal buoyancy gradient (cross-front wind stress). In this case the coefficients in Equation (48) are

$$(C_1^{\prime },C_2^{\prime },C_3)=(1/\mathrm{E},-2,\mathrm{E})\overline{K_0^{\prime 2}},$$

(76)

and hence the function in Equation (73) reduces to

$$\mathcal{C}=\frac{1}{\mathrm{E}},$$

(77)

and therefore the dimensional equilibrium buoyancy gradient is

$$\frac{\partial b^{*}}{\partial x^{*}}=\frac{{\tau }^{*}}{{\rho }_0{H}^2\mathrm{E}}=\frac{f{\tau }^{*}}{{\rho }_0\nu },$$

(78)

since $\mathrm{E}=\nu /\left(f{H}^2\right)$ for turbulent diffusivity $\nu$. The value for $\nu$ (and hence $\mathrm{E}$) may be taken from observations or estimated via Equation (64) or Equation (66) for the cases of surface wind stress driven turbulence or surface cooling driven turbulence respectively. Typical values, e.g., $H=100\mathrm{m}$, $\mathrm{E}=0.1$ and ${\tau }^{*}=0.1{\mathrm{Nm}}^{-2}$ give an equilibrium buoyancy gradient of

$$\partial b^{*}/\partial x^{*}={10}^{-7}{\mathrm{s}}^{-2}.$$

(79)

This value is well-within the typical range for submesoscale fronts in the open ocean (e.g., [34]). For these parameters that are typical of submesoscale fronts, using our definition, the Rossby number in this case is

$$Ro=\frac{H}{f^{2}L}\frac{\partial b^{*}}{\partial x^{*}}\approx 0.1,$$

(80)

using a horizontal lengthscale of $L=10$ km and $H/L=0.01$. Therefore, our analysis for the equilibrium frontal width can still apply to submesoscale fronts in a parameter range where $Ro$ is small enough for our asymptotic analysis to be valid. However, for weaker mixing or stronger wind stress, the sharpening effects of the wind stress dominate over the buoyancy-driven spreading and the front sharpens towards an equilibrium state with a larger Rossby number where our model equations might break down.

Jindasa et al. [35] found that frontogenesis can sharpen a submesoscale front to scales of only a few meters, in this case, large Rossby number effects dominate rotation and the front behaves like a gravity current. While our model is not valid for large Rossby numbers, it is possible that the mechanism of frontal sharpening due to the coupling of turbulent mixing and wind stress described here can act to sharpen fronts to a point where finite Rossby number effects take over. Therefore our analysis may provide a framework for examining when extreme sharpening can happen.

Note that an equilibrium frontal width will only exist if the Ekman flow opposes the cross-front TTW flow. If the TTW flow is initially stronger (i.e., the initial buoyancy gradient is stronger than the equilibrium buoyancy gradient), then the front will spread towards an equilibrium, whereas if the wind stress-driven flow is stronger, the front will sharpen towards an equilibrium. In the latter case, sharpening may lead to small regions of sharp buoyancy gradient resulting in a step-like series of small fronts, rather than a narrowing of the whole frontal region.

8. Conclusions and Discussion

We have considered an analytical model of a mixed layer front forced by surface fluxes of buoyancy and momentum. Using an asymptotic expansion in the Rossby number we determined the leading-order velocity and buoyancy perturbations. By calculating the depth-averaged buoyancy fluxes, we then showed that the leading order buoyancy satisfies a nonlinear diffusion equation where the diffusivity depends on the surface-forcing and the strength of the mixed layer turbulence.

Depending on the direction of alignment between the front and the wind stress we found that it is possible for the front to sharpen or to spread. Sharpening occurs when the wind stress drives a cross-front flux that opposes and overcomes the TTW-driven cross-front flux and leads to a horizontal convergence in the buoyancy field. Spreading occurs when the stress-driven flux either reinforces the TTW flux or opposes but does not exceed it and the front then spreads via shear dispersion [1,21,36].

We find that surface heating can also act to spread the front due to the upward advection of less buoyancy fluid by the vertical velocity. Surface cooling may lead to frontal sharpening if it is sufficiently strong to overcome the spreading effect of horizontal mixing. Additionally, surface cooling may act to enhance the mixed layer turbulence which will lead to a faster horizontal spreading via shear dispersion due to an increased Ekman number [1].

Note that since we used a diffusive parametrisation for mixed layer turbulence and did not consider the effects of gravitational instability or convective overturning, our model can develop a gravitationally unstable density profile in the cases of surface cooling or strong wind stress opposing the cross-front TTW velocity. These inversions have been noted in previous theoretical studies [29,30], although modelling studies have found that symmetric instability is effective at re-stratifying the mixed layer beneath a ‘convective layer’ [4]. We have not considered the influence of symmetric instability (or other submesoscale instabilities) here.

For intensifying fronts, we find that the increase in the horizontal buoyancy gradient increases the frontolytic tendency of the cross-front TTW velocity. Eventually, this balances the frontogenetic tendency associated with surface-forcing, and the front approaches a balanced state with an equilibrium buoyancy gradient. Similarly, fronts may spread towards an equilibrium width if the frontolytic effect of the cross-front TTW velocity is stronger than the frontogenetic effect of the surface forcing. Due to the nature of backward diffusion, any sharpening will likely form small steps in the buoyancy profile rather than a significant narrowing of the frontal region.

The equilibrium state that we find from the asymptotic theory and numerical solutions when frontogenesis associated with the surface-forcing is balanced by spreading associated with the TTW circulation, may provide a means of maintaining a particular frontal width. Additionally, when surface-forcing induces frontolysis, this effect could act to counteract frontogenesis associated with other processes such as an external strain flow or the TTW secondary circulation [8,17].

Using typical parameters for ocean fronts we have estimated values for the equilibrium buoyancy gradient. Our prediction is of a similar magnitude to the gradients observed in strong submesoscale fronts and much stronger than the buoyancy gradients predicted by global ocean models due to these models not fully resolving submesoscale features. Therefore, our predictions for the scaling of the equilibrium buoyancy gradient may be useful for parametrising the frontal strength in models that do not resolve these scales.

Throughout our analysis, we have made several simplifying assumptions. Some of these can be easily relaxed, while others impose limitations on our method. Firstly, we have assumed that the turbulent viscosity and diffusivity are constant in both space and time. It is possible to use depth-dependent profiles for $\nu$ and $\kappa$ using the approach given in Crowe and Taylor [1]. This would not qualitatively influence solution beyond modifying the coefficients, $C_i$. Physically, we might also expect cross-front and time variations in $\nu$ and $\kappa$ due to the large-scale properties of the front [4,37]. However, the extent of these variations or their influence on frontal dynamics remains unclear. Time-varying viscosity and diffusivity may lead to feedback between the frontal evolution and the turbulent mixing which would change the time-dependence of the solution.

Secondly, several approximations have been made regarding the surface-forcing. In particular, we have assumed that the forcing is spatially independent and applied to both top and bottom surfaces. The application of forcing to both the top and bottom boundaries was discussed at the end of Section 2 and more realistic forcing will likely give similar results. While spatial dependence of the forcing has been neglected throughout, it could be easily introduced since Equation (36) can be used with a spatially varying ${\tau }_0$ and $Q_1$. We expect this to lead to complicated dynamics in which different regions of the front evolve differently, potentially leading to along-front (y) variations.

Finally, we have assumed a scale-separation between the along-front and cross-front scales. While this assumption is valid for most observed fronts, it does prevent the formation of baroclinic instability, which we showed was significantly modified by vertical mixing in [25]. Baroclinic instability may act to modify small-scale turbulence (e.g., [37]) or stratification (e.g., [33]), which could, in turn, modify the dynamics considered here.