Potential Vorticity Generation in Breaking Gravity Waves

Abstract

Potential vorticity (PV) is an important quantity in stratified flows because it is conserved following the flow in the absence of forcing and viscous and diffusive effects. However, as shown by previous work for unstratified turbulence, viscosity and diffusion, when present, are not purely dissipative and can create potential vorticity even when none is present initially. In this work, we use direct numerical simulations to investigate the viscous and diffusive generation of potential vorticity and potential enstrophy (integrated square PV) in stratified turbulence. Simulations are initialized with a two-dimensional standing internal gravity wave, which has no potential vorticity apart from some low-level random noise; as a result, all potential vorticity and enstrophy comes from viscous and diffusive effects. Significant potential enstrophy is found when the standing wave breaks, and the maximum potential enstrophy increases with increasing Reynolds number. The mechanism for the initial PV generation is spanwise diffusion of buoyancy perturbations, which grow as the standing wave three-dimensionalizes, into the direction of spanwise vorticity. The viscous and diffusive terms responsible are small-scale and are sensitive to under-resolution, so high resolution is required to obtain robust results.

1. Introduction

The Ertel potential vorticity (PV) is an important quantity in geophysical fluid dynamics. PV characterizes the component of vorticity normal to the isosurfaces of a conserved scalar, and it is materially conserved in the absence of forcing and molecular viscosity and diffusion. In quasi-geostrophic (QG) dynamics, where the PV can be inverted to find the velocity and buoyancy, conservation of PV is everything. Even for the primitive dynamics of the atmosphere and ocean, where large-scale QG turbulence interacts with inertia–gravity waves, weakly rotating stratified turbulence at intermediate scales, and isotropic turbulence at small scales, PV is a useful quantity for analyzing the vortical part of the flow. Indeed, PV is the basis of various decomposition methods for separating slow, quasi-horizontal vortical motion, which has PV, from fast, propagating inertia–gravity waves, which have no PV [1,2,3,4].

QG dynamics break down at intermediate length scales in the atmosphere and ocean, where stratification is strong but rotation is weaker than at large scales. Stratified turbulence is often used as an idealization of turbulent dynamics at these scales [5,6]. In stratified turbulence, PV is a useful diagnostic for quantifying the role of vortices and gravity waves. In the Craya–Herring decomposition [7] and the closely related vortical/wave mode decomposition [2], the vortical component of the flow accounts for the linear PV. Lilly [5] proposed that conservation of linear PV might constrain stratified turbulence and lead to an inverse energy cascade, as in two-dimensional turbulence. However, an inverse cascade does not occur in stratified turbulence [6,8], due to a leakage of energy from vortices to gravity waves [9] and the importance of nonlinear terms in the PV [10,11] due to the development of small vertical scales [6,12]. Nevertheless, the question of whether the downscale cascade of energy through such scales in the atmosphere and ocean is dominated by vortices or waves has not yet been fully resolved [13,14].

Since PV is used to diagnose and separate vortices and gravity waves in geophysical flows and stratified turbulence, it is important to know where PV comes from. In QG turbulence, the integrated square PV, called potential enstrophy, cascades conservatively to smaller scales, analogously to the enstrophy cascade in two-dimensional turbulence. As a result, the QG enstrophy cascade from the large-scale circulation is one source of PV at intermediate and small scales in the atmosphere and ocean, where QG breaks down and stratified turbulence may be present.

However, viscosity and diffusion can also create PV. Because PV depends nonlinearly on vorticity and buoyancy, viscous and diffusive effects are not purely dissipative or necessarily restricted to small length scales, as they are for regular vorticity. This effect was explored in unstratified turbulence with a passive scalar by Herring, Kerr, and Rotunno [15]. By analyzing the evolution of PV in decaying turbulence, they found that viscous and diffusive effects on PV may occur at large scales. By considering a Taylor–Green vortex oriented with zero initial PV, they demonstrated that viscosity and diffusion can create large-scale PV even when none is initially present. These findings raised doubts about whether the notion of a potential enstrophy cascade is meaningful in the absence of QG turbulence. Subsequent work suggested that a potential enstrophy cascade may be possible in certain parameter regimes of geophysical turbulence [16], but the role of viscous and diffusive generation of potential vorticity in stratified turbulence has not been explored.

In this paper, inspired by [15], we use direct numerical simulation (DNS) to investigate the role of viscosity and diffusion on the evolution of PV in stratified turbulence. Compared to [15], we explore the effects of stratification and the impact of larger Reynolds numbers. A two-dimensional standing internal gravity wave is used to initialize simulations of decaying stratified turbulence. A standing wave is a linear solution to the equations of motion, but instabilities and nonlinearities inevitably cause the wave to break down into smaller-scale turbulence. Standing waves have been used in laboratory [17] and numerical [18,19,20] studies of stratified turbulence. More generally, numerical simulation is an important tool for studying breaking gravity waves, which are a key source of turbulence in the atmosphere [21,22,23]. A breaking standing wave is a useful flow configuration for quantifying the viscous generation of PV because, like the Taylor–Green case in [15], the initial wave has no PV. As a result, any PV in the resulting turbulence must come from viscous generation. We will explore the magnitude and mechanisms of PV generation, the associated length scales, and the dependence on viscosity and stratification.

2. Materials and Methods

2.1. Equations

We use the uniformly stratified, non-rotating Boussinesq equations to describe stratified turbulence. The equations of motion are

$$\begin{array}{cc}\hfill & \hfill \frac{\mathrm{D}\mathit{u}}{\mathrm{D}t}=-\nabla p+b\widehat{\mathbf{z}}+\nu {\nabla }^2\mathit{u},\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & \hfill \nabla ·\mathit{u}=0,\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill & \hfill \frac{\mathrm{D}b}{\mathrm{D}t}+N^{2}w=\kappa {\nabla }^{2}b,\hfill \end{array}$$

(3)

where $\mathrm{D}/\mathrm{D}t=\partial /\partial t+\mathit{u}·\nabla$ is the material derivative, $\mathit{u}=(u,v,w)$ is the velocity, b is the buoyancy, p is the pressure scaled by a reference density, N is the constant buoyancy frequency, $\nu$ is the kinematic viscosity, $\kappa$ is the buoyancy diffusivity, and $\widehat{z}$ is the vertical unit vector. There is no forcing. The buoyancy b is proportional to the negative fluctuation of density from the linear basic state, and $N^2$ is proportional to the background density gradient; the total buoyancy $b_{\mathrm{tot}}=N^{2}z+b$ is proportional to the (negative) density.

The Ertel PV for (1)–(3) is

$$q=\mathbf{\omega }·\left(N^2\widehat{\mathbf{z}}+\nabla b\right)=N^2{\omega }_z+\mathbf{\omega }·\nabla b,$$

(4)

where $\mathbf{\omega }=\nabla \times \mathit{u}=({\omega }_x,{\omega }_y,{\omega }_z)$ is the vorticity [24]. The PV evolves according to

$$\begin{array}{cc}\frac{\mathrm{D}q}{\mathrm{D}t}\hfill & =\nu \left(N^2\widehat{\mathbf{z}}+\nabla b\right)·{\nabla }^2\mathbf{\omega }+\kappa \mathbf{\omega }·{\nabla }^2\left(\nabla b\right)\hfill \\ \hfill & =\nu N^2{\nabla }^2{\omega }_z+\nu \nabla b·{\nabla }^2\mathbf{\omega }+\kappa \mathbf{\omega }·{\nabla }^2\left(\nabla b\right).\hfill \end{array}$$

(5)

In the absence of viscosity and diffusion, $\mathrm{D}q/\mathrm{D}t=0$ and the PV is materially conserved. When $\nu$, $\kappa >0$, the evolution of q is more complicated. Due to the nonlinear dependence of q on $\mathbf{\omega }$ and b, the viscous and diffusive terms in the q Equation (5) are not simply diffusive. Viscous effects can modify q by viscous diffusion of the vorticity into the direction of the buoyancy gradient, either by diffusion of ${\omega }_z$ along the background buoyancy gradient $N^2$ or by diffusion of $\mathbf{\omega }$ in the direction of $\nabla b$. In addition, buoyancy diffusion can modify q by diffusion of the buoyancy gradient into the direction of the vorticity. As a result, as indicated by [15], viscosity and diffusion may create PV even when none is initially present.

The potential enstrophy is the integrated square PV, and is given by

$$\begin{array}{cc}\hfill V& \hfill =\frac{1}{2}\langle q^2\rangle ,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & \hfill =\frac{1}{2}{N}^4\langle {\omega }_z^2\rangle +N^2\langle {\omega }_z\mathbf{\omega }·\nabla b\rangle +\frac{1}{2}\langle {\left(\mathbf{\omega }·\nabla b\right)}^2\rangle ,\hfill \end{array}$$

(7)

where $\langle ·\rangle$ denotes domain average. Assuming no flux through the boundary, the potential enstrophy satisfies

$$\begin{array}{cc}\frac{\mathrm{d}V}{\mathrm{d}t}\hfill & =\langle \nu N^{2}q{\nabla }^2{\omega }_z+\nu q\nabla b·{\nabla }^2\mathbf{\omega }+\kappa q\mathbf{\omega }·{\nabla }^2\left(\nabla b\right)\rangle \hfill \\ \hfill & =D_{\nu }+D_{\kappa },\hfill \end{array}$$

(8)

where $D_{\nu }$ is the term with viscosity $\nu$ and $D_{\kappa }$ is the term with buoyancy diffusivity $\kappa$. The viscous and diffusive terms are not necessarily dissipative, and may be a source or sink of potential enstrophy.

2.2. Scale Analysis

We can estimate the magnitude of the potential enstrophy in stratified turbulence using scale analysis. In stratified turbulence with large buoyancy Reynolds number $R{e}_b=ϵ/\left(\nu N^2\right)$, where $ϵ$ is the kinetic energy dissipation rate, there should be a Kolmogorov inertial range at small scales between the Ozmidov and Kolmogorov scales

$$L_O={\left(\frac{ϵ}{N^3}\right)}^{1/2},\eta ={\left(\frac{{\nu }^3}{ϵ}\right)}^{1/4},$$

(9)

respectively [25,26,27], since $L_O/\eta \sim {\left(R{e}_b\right)}^{3/4}$ [28,29]. Since the PV is a function of the velocity and buoyancy gradients, it should be dominated by the smallest scales when such an inertial range is present.

Consider the third term in (7), which is likely to be the largest term when the Reynolds number is large since it depends on the squared vorticity and squared buoyancy gradient. Using Kolmogorov scaling gives $\mathbf{\omega }\sim \upsilon /\eta$, where $\upsilon ={\left(\nu ϵ\right)}^{1/4}$ is the Kolmogorov velocity scale. We assume that the small-scale buoyancy fluctuations are $b\sim \upsilon N$. As a result, estimating V using the third term in (7) gives

$$\begin{array}{cc}\hfill V& \hfill \sim \frac{N^2{ϵ}^2}{{\nu }^2}\sim \frac{N^4{U}^2}{L^2}R{e}^{2}F{r}^2\sim \frac{N^4{U}^2}{L^2}R{e}_{b}Re,\hfill \end{array}$$

(10)

where the second and third ∼ use Taylor’s hypothesis $ϵ\sim U^3/L$ [6], U and L are the large-scale velocity and length scales, $Re=UL/\nu$ is the Reynolds number, and $Fr=U/\left(NL\right)$ is the Froude number. Note that $N^4{U}^2/L^2$ is the magnitude of the potential enstrophy from the large energy-containing scales. This estimate suggests that stratified turbulence may have significant potential enstrophy at small scales, and that potential enstrophy increases with increasing $Re$ and $R{e}_b$.

2.3. Numerical Approach

Simulations are performed in a cubic domain with size $2\pi$ and periodic boundary conditions. The initial conditions for a two-dimensional standing wave with unit wavenumber, following [19], are

$$\begin{array}{cc}\hfill u& \hfill =asinxcosz,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill w& \hfill =-acosxsinz,\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill b& \hfill =0,\hfill \end{array}$$

(13)

where $a>0$ is the amplitude and $v=0$. In the linear equations, these initial conditions excite a standing wave

$$\begin{array}{cc}\hfill u& \hfill =asinxcoszcos\left(\omega t\right),\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill w& \hfill =-acosxsinzcos\left(\omega t\right),\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill b& \hfill =\frac{a{N}^2}{\omega }cosxsinzsin\left(\omega t\right),\hfill \end{array}$$

(16)

with frequency $\omega =N/\sqrt{2}$. Since the initial buoyancy and vertical vorticity are zero, the initial standing wave has no PV. Standing waves eventually break and generate small-scale turbulence [19,20]. In particular, the standing wave (14)–(16) is overturning and subject to convective instability when $\sqrt{2}a/N>1$. In our simulations, the wave amplitude is fixed at $a=1$. Stratification and the overturning of the standing wave are controlled by varying N.

Equations (1)–(3) are simulated with a Fourier-based spectral transform model. Time-stepping is third-order Adams–Bashforth and the viscous and diffusion terms are handled with a Crank–Nicolson approach [30]. Aliasing is removed using the two-thirds rule, so with $n^3$ grid points, wave vectors $k=(k_x,k_y,k_z)$ are truncated at $|k_x|$, $|k_y|$, $|k_z|=k_T=n/3$, i.e., we keep $-n/3<k_x,k_y,k_z<n/3$. The resolution n is chosen to resolve the Kolmogorov scale with $\eta k_T\ge 2.5$, where $\eta$ is calculated using the maximum value of $ϵ=\nu \langle |\mathbf{\omega }{|}^2\rangle$. As discussed below, high resolution is required to resolve the viscous terms in the q equation properly. While recent studies have shown that stratified turbulence may be sensitive to increasing Prandtl number $Pr=\nu /\kappa$ at $Re$-values that are attainable in DNS [31,32], in this work we set $Pr=1$ for computational efficiency, so $\kappa =\nu$. A small amount of random noise is added to the initial velocity, with uniform distribution and magnitude $a/100$, to facilitate three-dimensionalization. This noise introduces a small amount of initial q and V that is quickly removed by viscous and diffusion terms, as shown below.

Since the initial standing wave has unit velocity and unit wavenumber, Equations (1)–(3) can be interpreted as dimensionless with $\nu =1/Re$, $N=1/Fr$, and unit turnover time. We use $\nu$, N and $Re$, $Fr$ interchangeably in what follows to distinguish stratification and viscosity. Several simulations at different $Re$ and $Fr$ are discussed below and summarized in

Table 1. Summary of simulations.

Run	$\mathit{Fr}$	$\sqrt{2}\mathit{a}/\mathit{N}$	$\mathit{Re}$	n	Max $\mathit{ϵ}$	$\mathit{\eta }{\mathit{k}}_{\mathit{T}}$	Max ${\mathit{Re}}_{\mathit{b}}$
Main	1	$\sqrt{2}$	3333	1024	0.0037	3.2	12.2
High Visc	1	$\sqrt{2}$	1667	512	0.0047	2.5	7.9
Low Visc	1	$\sqrt{2}$	4500	1728	0.0060	3.8	26.9
High Strat	$1/\sqrt{2}$	1	3333	1024	0.0052	2.9	8.6
Higher Strat	1/2	$1/\sqrt{2}$	3333	1034	0.0030	3.3	2.5
Low Res	1	$\sqrt{2}$	3333	512	0.0041	1.5	13.8
High Res	1	$\sqrt{2}$	3333	1536	0.0040	4.6	13.2

. The control simulation (Main) has $n=1024$, $N=1$ ($Fr=1$), $\nu =0.0003$ ($Re=3333$), and $\sqrt{2}a/N=\sqrt{2}$, so the standing wave overturns and quickly breaks down into stratified turbulence. Sensitivity to Reynolds number is explored in simulations High Visc and Low Visc, in which $Re$ is decreased to 1667 and increased to 4500. Sensitivity to stratification is explored in simulations High Strat and Higher Strat, in which $Fr$ is reduced to $1/\sqrt{2}$ and $1/2$. Finally, sensitivity to resolution is investigated in simulations Low Res and High Res by reducing and increasing the resolution to 512 and 1536 grid points.

3. Results

3.1. Main Simulation

We begin by examining the potential enstrophy evolution in the Main simulation, which has $Fr=1$, $Re=3333$, and a large, overturning wave amplitude. Time series of kinetic and potential energy, energy dissipation, turbulent Froude number, and potential enstrophy are shown in Figure 1. The wave period for this simulation is $T=2\pi \sqrt{2}\approx 8.9$. The wave makes two periods with negligible dissipation. By $t=20$, the energy dissipation begins increasing and the energy begins to decrease more significantly. The potential energy dissipation reaches a maximum at $t=26$, followed by the maximum kinetic energy dissipation at $t=30$. Table 1 gives $\eta k_T=3.2$ and $R{e}_b=12.2$, both of which are based on the maximum value of $ϵ$ in the simulation, which is also given in Table 1. Figure 1c shows that the turbulent Froude number $F{r}_t=ϵ/\left(N{E}_K\right)\lesssim 0.1$, where $F{r}_t$ is defined using velocity scale $\sqrt{E_K}$, where $E_k$ is the domain-averaged kinetic energy, and length scale from Taylor’s hypothesis [29].

The potential enstrophy in this simulation is shown in Figure 1d. The potential enstrophy is negligible until $t\approx 24$. The inset shows V at early times with a log scale; the small potential enstrophy introduced by the random noise decays rapidly until $t\approx 15$, and then begins to increase. A maximum value of $V=1.9$ is reached at $t=29$. To put this number in perspective, we can calculate the potential enstrophy that would be associated with the initial velocity field if it was rotated by 90${}^{\circ }$:

$$V_0=\frac{1}{2}{N}^4\langle {\left({\omega }_y^0\right)}^2\rangle =\frac{1}{2},$$

(17)

where ${\omega }_y^0$ is the vorticity associated with (11)–(13). Therefore, we can see that significant potential enstrophy is generated in this simulation, about four times larger than what one would expect from the initial velocity field and stratification alone.

Vertical $(x,z)$ slices of the total buoyancy $b_{\mathrm{tot}}$, normal component of vorticity ${\omega }_y$, and potential vorticity q are plotted in Figure 2 for a selection of times from $t=16$, before the onset of significant dissipation, to $t=36$, when the turbulence is decaying. First, we discuss the evolution of the buoyancy and normal vorticity (left and centre columns of Figure 2). At $t=16$, the buoyancy (Figure 2a) shows the structure of the standing wave with some overturning isopycnals, as expected. At $t=20$, the overturning isopycnals have collapsed (Figure 2b) with significant spanwise vorticity in the overturning regions (Figure 2g). By $t=24$, the buoyancy field of the standing wave has been significantly altered by breaking (Figure 2c), and smaller spanwise vortices are present (Figure 2h). At $t=28$ and 36, the wave has largely broken down into turbulence, with the vorticity dominated by small scales (Figure 2d–e,i–j).

The potential vorticity is shown in the right column of Figure 2. There is no significant PV visible until $t=24$. At $t=24$, thin bands of positive and negative PV appear, mainly in the lower overturning region (similar bands are apparent in the upper overturning region at other y). By $t=28$, the PV structure has become more rich and complex, with both small-scale bands and vortex structures visible. At this time, the PV is still mainly located in the upper and lower breaking regions, avoiding $z=0$, $\pi$, and $2\pi$. This time corresponds to the maximum V in Figure 1d. The PV occupies more of the domain at $t=36$, but the amplitude is lower, consistent with the decay of V seen in Figure 1d.

To investigate the increase in PV around $t\approx 24$, time series of the terms in the potential enstrophy Equation (8) are shown in Figure 3. These terms are negligible for $t\lesssim 24$, but begin to deviate from zero for $t\gtrsim 24$. From $24<t<26$, the buoyancy diffusivity term $D_{\kappa }$ dominates, and is responsible for the initial creation of V as the wave begins to break. At these times, the viscous terms are slightly dissipative. V grows rapidly from $24<t<26$, but the growth slows around $t=26$ (Figure 1d). At $t=26$, the growth in V increases again, this time, as seen in Figure 3, due to a sharp increase in the viscous term $D_{\nu }$. After around $t=28$, the combined viscous and diffusion term is negative and V decreases, apart from a short increase in V around $t=36$.

Figure 3 indicates that the initial generation of PV, as seen in the rise in V, is a result of the $\kappa$ term in the V equation, which corresponds to diffusion of the buoyancy gradient into the direction of vorticity. The $\kappa$ term can be decomposed as

$$D_{\kappa }=\kappa q\mathbf{\omega }·{\nabla }^2(\nabla b)=\kappa q{\omega }_x{\nabla }^2\frac{\partial b}{\partial x}+\kappa q{\omega }_y{\nabla }^2\frac{\partial b}{\partial y}+\kappa q{\omega }_z{\nabla }^2\frac{\partial b}{\partial z}.$$

(18)

For $24<t<26$, $D_{\kappa }$ is dominated by the second term, which involves ${\omega }_y$ and $\partial b/\partial y$ (not shown). The mechanism of PV generation by this term is illustrated in Figure 4, which shows a vertical slice of q and two horizontal slices of b at $t=24$. The flow is still mainly two-dimensional at this time, but three-dimensionalization of the buoyancy field is visible along a few lines parallel to the y-axis in Figure 4b,c. The planes $z=z_1$ and $z=z_2$ cut through bands of PV in the lower breaking region. Along $z_1$, the PV bands occur precisely where the buoyancy field is three-dimensionalizing, as seen in Figure 4b, where the ticks mark the locations of the PV along $z_1$ in Figure 4a. Similarly, the PV bands along $z=z_2$ occur at the same location as the buoyancy three-dimensionalization at that height. Three-dimensionalization of b promotes growth of $\partial b/\partial y$, which diffuses into the direction of ${\omega }_y$, which at this time is still the dominant component of the vorticity. Therefore, PV first emerges at the locations where three-dimensionalization of the buoyancy develops, via the second term in (18).

Kinetic energy spectra $E\left(k\right)$ and potential enstrophy spectra $V\left(k\right)$, plotted as functions of the total wavenumber $k=\left|k\right|$, are shown in Figure 5. Spectra are calculated at $t=28$, which is the time of maximum potential enstrophy. As expected, the energy spectrum is peaked at the primary standing wave wavenumber $k=1$, oscillates at the lowest few wavenumbers, and then falls off like a power law until the Kolmogorov wavenumber, which is marked in the figure. By contrast, the potential enstrophy spectrum is peaked at large wavenumbers. It has a positive slope of around 3.3 for $1\le k\le 20$. Despite the fact that V is generated by viscous and diffusive effects, $V\left(k\right)$ has a broad peak around $k\approx 20$ to 40, which is just beyond the Ozmidov wavenumber, which is marked in the figure. However, we will see below that this peak is not set by the Ozmidov wavenumber.

3.2. Sensitivity to Reynolds and Froude Numbers

We investigate the dependence of viscous and diffusive PV generation on Reynolds number with two additional simulations with lower and higher Reynolds numbers. The Reynolds number is varied by changing $\nu$, and the resolution is adjusted to resolve the Kolmogorov scale. Time series of the potential enstrophy for these simulations, along with the Main simulation, are shown in Figure 6a. With smaller $Re$, the maximum potential enstrophy is reduced from around 1.9 to 0.4. The time of maximum potential enstrophy is also delayed, from around $t=28$ in the Main simulation to $t=37$ in the High Visc simulation. On the other hand, for the Low Visc simulation with larger $Re$, the maximum potential enstrophy is much larger, reaching a value of around 64 at $t=27$, which is slightly earlier than the Main simulation. Potential enstrophy spectra from these simulations are shown in Figure 7a. The potential enstrophy spectra for the Main and High Visc experiments are similar at small k but increase with increasing $Re$ at larger k, indicating that the increased potential enstrophy with higher Reynolds number is primarily at small scales going from High Visc to Main. However, the potential enstrophy in the Low Visc case is increased at all scales, not just small scales.

Time series of potential enstrophy from simulations with increased stratification are shown in Figure 6b. Increasing the stratification delays the onset of wave breaking and, as a result, the growth in potential enstrophy is also delayed. The High Strat simulation has $\sqrt{2}a/N=1$ and is therefore on the threshold of convective instability; this simulation reaches a maximum potential enstrophy of 3.1 at $t=49$, which is larger and later than in the Main simulation. The Higher Strat simulation has $\sqrt{2}a/N=1/\sqrt{2}$ and reaches a maximum V of $2.5$ at $t=82$. This is a smaller maximum V than in the High Strat simulation, and it also occurs later. Potential enstrophy spectra from these simulations are shown in Figure 7b. The peak in potential enstrophy remains around $k\approx 20$ to 40, even as the Ozmidov wavenumber, which is shown in the figure, increases from 16 to 37 with increasing stratification. The peak in the potential enstrophy spectrum therefore does not scale with the Ozmidov wavenumber. In addition, the large-scale potential enstrophy, with $k<20$, increases with increasing stratification.

To make sense of the dependence of V on viscosity and stratification, we can compare our simulation results with the prediction (10). The dependence on stratification in particular is difficult to understand, because the maximum V seems to increase, then decrease, with increasing stratification. However, analyzing the dependence of V on $Fr$ or N alone presents an incomplete picture, because the strength of the turbulence, as captured by the maximum $ϵ$, also changes with stratification (Table 1). Equation (10) suggests that the maximum V should depend on $\nu$, N, and $ϵ$. Figure 8 shows the maximum V plotted against $N^2{ϵ}^2/{\nu }^2$, where the maximum value of $ϵ$ is used (maximum $ϵ$ generally occurs around the same time as maximum V). While the agreement with (10) is not particularly good, we do find that maximum V increases with increasing $N^2{ϵ}^2/{\nu }^2$, which accounts for the non-monotonic dependence on stratification. In all cases, the maximum potential enstrophy is less than what was predicted by (10), but the Low Visc case, which has the highest $Re$ and $R{e}_b$, is closest to (10).

3.3. Sensitivity to Numerical Resolution

High spatial resolution is required to accurately resolve the viscous and diffusive terms in the PV equation. Under-resolution can lead to spurious generation of PV. The dependence on resolution is illustrated in Figure 9, which shows time series of V for the Main simulation and simulations with higher and lower resolution. The Low Res simulation, which has $\eta k_T=1.5$, has a maximum V of around 4.6. This value is more than double the maximum V of 1.9 in the Main simulation, which has $\eta k_T=3.2$. By contrast, the High Res simulation, which has $\eta k_T=4.6$, has a maximum V of 1.8, which is just slightly smaller than the value in the Main simulation. The timing of maximum V is slightly accelerated in the Main simulation compared to the High Res simulation, but these simulations are much closer to one another than to the Low Res case. Inadequate resolution results in accelerated and artificially high potential enstrophy growth.

4. Discussion

A breaking standing gravity wave is a useful case for studying viscous and diffusive generation of PV in stratified turbulence because such flows have no initial PV. Any PV that develops must come from viscous and diffusive effects. Our simulations indicate that significant PV and potential enstrophy are created as these waves break, much more than would be expected from the velocity and length scales of the breaking wave. Indeed, more potential enstrophy is created at higher Reynolds numbers, because even though the viscosity and diffusion coefficients are smaller, the velocity and buoyancy gradients in the turbulence are correspondingly larger. The Reynolds numbers in this work, while large compared to what was possible in [15], are still small compared to geophysical fluid values. It therefore seems likely that viscous and diffusive generation of PV from breaking gravity waves may be even more significant in the atmosphere and ocean.

The mechanism of the initial PV generation is diffusion of spanwise buoyancy fluctuations, which develop as the waves three-dimensionalize, into the direction of the spanwise vorticity, which is initially dominant. As a result, the initial PV generation occurs at small length scales associated with three-dimensionalization of the flow. Once significant potential enstrophy is present, spectra are peaked at intermediate wavenumbers 20–40, corresponding to scales smaller than the wavelength of the breaking wave, but much larger than the Kolmogorov scale. Interestingly, this intermediate scale does not vary like the Ozmidov scale, and indeed does not change much with Reynolds or Froude number, indicating that it might be related to the scale of the breaking wave. The potential enstrophy at large scales, while relatively small, increases with increasing Reynolds numbers and decreasing Froude numbers.

The mechanism of PV generation in breaking standing waves is somewhat different from the Taylor–Green vortices considered by [15]. The Taylor–Green vortices were oriented so that they had no initial PV, like the standing wave. Potential enstrophy in their case was initially created at large scales, followed by more intense generation at small scales; they confirmed that the small-scale generation was due to viscous and diffusive effects rather than transfer from large scales. The standing wave does not have significant large-scale potential enstrophy at early times. However, for both the Taylor–Green vortex in [15] and our standing wave, viscous and diffusive effects are primarily a source of potential enstrophy; dissipation takes over only at late times. By contrast, for the decaying turbulence with Gaussian initial conditions considered in [15], viscous and diffusive effects were primarily dissipative. More work is required to determine if the differences between our results and the pioneering work of [15] are due to the different Reynolds numbers, initial configuration, or the presence of stratification. In addition, it would be interesting to examine the role of viscous and diffusive effects in stratified turbulence with large-scale PV from random initial conditions or forcing.

Viscous and diffusive generation of potential enstrophy is sensitive to numerical resolution. DNS requires resolution of the Kolmogorov scale. Resolutions with $\eta k_T\approx 1$ are typical for DNS of breaking gravity waves [23] and decaying stratified turbulence [33]. However, our findings show that under-resolution can lead to spurious generation of PV when $\eta k_T=1.5$. Instead, larger values of $\eta k_T\ge 2.5$ are needed for robust simulations of potential enstrophy generation. Such high resolution is not routinely employed for DNS of stratified turbulence. This finding, while perhaps not surprising given the small-scale mechanisms at play, suggests that some DNS studies of stratified turbulence and gravity wave breaking may have unrealistically high potential enstrophy due to under-resolution of the viscous generation processes. More work is required to understand the implications on the dynamics of turbulence with artificially enhanced potential enstrophy.