2.1. Spectral model formulation and solution method

We employ a Boussinesq, Fourier pseudo-spectral model for a DNS of multiple, interacting KH billows in a domain that allows varying phases and wavelengths along their axes and regions where they are misaligned as they arise. The model formulation, solution method, boundary and initial conditions and analysis and visualization methods are described below.

The spectral model is triply periodic and solves the Boussinesq Navier–Stokes equations expressed as

(2.1)\begin{gather} \frac{\partial u_j}{\partial x_j} = 0, \end{gather}

(2.2)\begin{gather} \frac{\partial u_i}{\partial t} + \frac{\partial u_iu_j}{\partial x_j} + \frac{1}{\rho_0}\frac{\partial p^\prime}{\partial x_i} = \frac{\theta^\prime}{\theta_0}g\delta_{i3} + \nu\frac{\partial^2 u_i}{\partial x_j\partial x_j} \end{gather}

and

(2.3)\begin{equation} \frac{\partial \theta}{\partial t} + \frac{\partial \theta u_j}{\partial x_j} = \kappa\frac{\partial^2 \theta}{\partial x_j\partial x_j}. \end{equation}

Here, $u_i=(u_1,u_2,u_3)$ are the velocity components along $(x,y,z)$, $p$, $\rho$ and $\theta$ are the pressure, density and potential temperature (with $\theta =T$ for a Boussinesq flow), $\nu$ and $\kappa$ are kinematic viscosity and thermal diffusivity, primes and ‘0’ subscripts denote perturbation and mean quantities and successive indices imply summation.

The Patterson & Orszag (Reference Patterson and Orszag1971) pseudo-spectral algorithm for computing nonlinear products, removal of aliasing errors using the 2/3rds truncation rule and a low-storage, third-order Runge–Kutta method Williamson (Reference Williamson1980) for time advancement were employed for accuracy and efficiency. The code is highly optimized for massively parallel computers, and we use a maximum of 3240 cores.

2.2. Boundary conditions, initial conditions and domain

We specify a shear flow representative of the atmosphere and oceans, where inertia-GWs (IGWs) dominate horizontal velocity and shear variances and IGW steepening yields localized shear and stability maxima, with mean wind and stability profiles that we approximate by

(2.4)\begin{equation} U(z) = U_0\cos\left(\frac{{\rm \pi} z}{Z}\right)\tanh\left(\frac{z}{h}\right)\end{equation}

and

(2.5)\begin{equation} N^2(z) = N_0^2 + (N_m^2-N_0^2) {\rm sech}^2\left(\frac{z}{h}\right). \end{equation}

Here, $Z$ is the domain depth, $z=0$ at the domain centre, $U_0$ is the half-shear velocity difference and we assume $N_0^2=10^{-4}\,{\rm s}^{-2}$ and $N_m^2=8N_0^2$, as large enhancements in thin layers are seen to arise due to GW interactions in multi-scale environments (Fritts, Wang & Werne Reference Fritts, Wang and Werne2013). The choice for $N_m$ yields a buoyancy period at the peak stratification of $T_b=2{\rm \pi} /N_m=222\,{\rm s}$. The remaining parameters are related by defining a minimum Richardson number, $Ri=N_m^2/({\rm d}U/{\rm d}z)^2=0.1$, a Reynolds number, $Re=U_0h/\nu =5000$, for kinematic viscosity, $\nu$, and a dominant KHI wavelength ${\lambda _h}_0=L$. We then specify a shear depth and $h=0.07L$, and a streamwise domain of length $X=3L$. This choice for $h$ is somewhat smaller than the linear stability estimate, $\lambda _h/(4{\rm \pi} )$, in order to enable primarily 3, but occasionally 4, initial KHI along $X$ arising from the initial noise seed. A large domain along the KHI billows of length $Y=9L$ enables several regions where emerging KH billows are misaligned or exhibit significant phase variability along $y$. A domain depth $Z=X=3L$ enables strong dynamics extending well away from the initial shear layer. The initial fields are shown in figure 1.

Figure 1. Initial fields in $U(z)$, $T(z)$, $N^2(z)$ and $Ri$ employed for the DNS discussed below.

An initial white-noise seed (zero spectral slope) with $u_{rms}=U_0\times 10^{-5}$ enables emergence of weak, initial KHI (defined as $t=0$) having varying wavelengths and phases along $y$ that lead to misaligned cores and significant phase variations that persist to finite amplitudes. A second noise seed having $u_{rms}=U_0\times 10^{-2}$ is added to the emerging multi-billow field (accompanying emergence of small-scale structures due to initial KH billow interactions at $t=3.1T_b$) to approximate the evolution of a multi-scale KHI field in the presence of background turbulence and ensure secondary KHI, as was inferred in previous observational and modelling studies by Hecht et al. (Reference Hecht, Wan, Gelinas, Fritts, Walterscheid and Franke2014) and Fritts et al. (Reference Fritts, Wan, Werne, Lund and Hecht2014) based on the observed spatial scales of secondary CIs. This acts to ensure the presence of secondary KHI in the braids between adjacent billows and their potential for competition with the dynamics of tubes and knots.

We note that the specified $Re=5000$ enables performance of multi-billow DNS describing the evolutions of tubes and knots. This $Re$ is appropriate for KHI $\lambda _h\sim 2\text {--}3\,{\rm km}$ at PMC altitudes of $\sim$82 km, but is orders of magnitude smaller than is appropriate for observed KHI scales in the middle stratosphere ($\lambda _h\sim 1\text {--}3\,{\rm km}$) or at lower altitudes. It is sufficiently large, however, to compare the relative contributions of KHI tubes and knots with those of secondary instabilities of individual billows in driving KHI turbulence and billow breakdown. For atmospheric reference, we dimensionalize our results by choosing a representative large-scale KHI $\lambda _h=2\,{\rm km}$ for the stratosphere, but which can be much smaller or up to $\sim$3 times larger at these and other altitudes. This choice, and the $N_m$ and $h$ specified above imply $U_0=12.5\,{\rm m}\,{\rm s}^{-1}$.

2.3. Spectral model resolution

DNS places severe constraints on resolution and major requirements for computational resources where high spatial resolution is required to describe turbulence dynamics approaching the Kolmogorov scale,

(2.6)\begin{equation} \eta=(\nu^3/\bar{\epsilon})^{1/4}. \end{equation}

Here, the local $\bar {\epsilon }$ is defined as the largest average over individual $x$–$y$, $x$–$z$ and $y$–$z$ planes spanning the computational domain (of which the $x$–$y$ means are typically the largest)

(2.7)\begin{equation} \bar{\epsilon}=2\nu\langle S_{ij}S_{ij}\rangle, \end{equation}

where brackets denote the spatial averaging and

(2.8)\begin{equation} S_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right). \end{equation}

We assume a Prandtl number, $Pr=1$, so as to require the same resolution for velocity and temperature fields. Additionally, we employ isotropic resolution and ensure that

(2.9)\begin{equation} \Delta x<1.8\eta \end{equation}

is satisfied throughout to ensure true DNS (Moin & Mahesh Reference Moin and Mahesh1998; Pope Reference Pope2000). This constraint is based on the assumption of isotropic homogeneous turbulence. However, emerging instabilities and localized small-scale turbulence imply a very much larger volume-averaged $\bar {\eta }$ than occurs locally and implied insufficient resolution where it is most needed accompanying initial turbulence sources.

In order to define a more realistic (but computationally costly) required resolution, we instead calculate local $\bar {\eta }$ separately on all two-dimensional (2-D) model planes (normal to $x$, $y$ and $z$) and use the smallest obtained $\bar {\eta }$ (most often occurring on the $x$–$y$ plane at $z=0$) to specify the required resolution. This significantly increases computational resource needs, as these vary for multi-radix spectral codes as ${\sim }{N_xN_yN_zN_t}\ln (N_x)\ln (N_y)\ln (N_z)$ for an interval $\Delta {t}=N_t\delta t$ and resolution $\Delta x=X/N_x$ for domain width $X$. Our highest spatial resolution of $(N_x,N_y,N_z)=(3240,9720,3240)$ requires $(1080,3240,1080)$ de-aliased complex spectral modes, thus costs $\sim$20 times more resources than half that resolution. Comparisons of local 1.8$\bar {\eta }/L$ and model resolution $\Delta x/L$, and of $N_x$ and local $\bar \epsilon$, demonstrating compliance with (2.9) throughout the simulation are shown in figure 2.

Figure 2. Temporal variations of 1.8$\bar {\eta }/L$ and $\Delta x/L$ (lines and bars in a), with dashed (solid) lines showing averaging centred at $z=0$ (over $|z|/L\le 0.3$). Expansion of the turbulence layer along $x$ and in $z$ by 5.5$T_b$ and further thereafter (see below) justifies the broader averaging after 5.5$T_b$ employed to define 1.8$\bar {\eta }/L$ and $\Delta x/L$ and demonstrates compliance with (2.9) extending to the later times for which the mean energy dissipation rate, $\bar {\epsilon }$, is evaluated by Fritts et al. (Reference Fritts, Wang, Thorpe and Lund2022). For reference, we assume $\nu =0.35\,{\rm m}^{2}\,{\rm s}^{-1}$, which corresponds to an $\sim$65 km altitude at mid-latitudes in winter.

2.4. Describing sources and evolutions of tubes, knots and twist waves

We employ volume visualization to reveal the vorticity dynamics of KHI tube and knot dynamics following earlier modelling studies of GW breaking by Andreassen et al. (Reference Andreassen, Hvidsten, Fritts and Arendt1998) and Fritts, Arendt & Andreassen (Reference Fritts, Arendt and Andreassen1998). These authors used volume visualization of vorticity sources and dynamics to reveal the evolutions from initial vorticity sheets to vortex tubes, interactions among vortex tubes exciting twist waves and their cascade to smaller scales via successive twist wave interactions thereafter. For these purposes, viscous diffusion can be neglected and the vorticity equation obtained from our Boussinesq equations reduces to

(2.10)\begin{equation} \frac{{\rm D}\zeta_i}{{\rm D}t} \simeq \zeta_jS_{ij} -\left[\left(\frac{\boldsymbol{\nabla}\theta}{\theta}\right)\times \left(\frac{\boldsymbol{\nabla} p}{\rho}\right)\right]_i, \end{equation}

where $\zeta _j=(\boldsymbol {\nabla }\times V)_j$ and repeated subscripts in (2.10) imply summation. The first source term includes the tilting, twisting and stretching components, the second includes the baroclinic components. Of these, the tilting and twisting terms alter vorticity orientations, but do not increase magnitudes. Hence, the stretching terms contribute most to vorticity intensification

(2.11)\begin{equation} \frac{\partial \zeta_i}{\partial t} \simeq \zeta_jS_{ij}, \end{equation}

and these influences will be discussed further in § 4.3.

We also employ the negative intermediate eigenvalue $\lambda _2$ of the tensor defined as $\boldsymbol {H}=\boldsymbol {S}^2+\boldsymbol {R}^2$ to identify flow features having strong rotational character, as described by Jeong & Hussain (Reference Jeong and Hussain1995), where $\boldsymbol {S}$ is the strain tensor defined in (2.8) and $\boldsymbol {R}$ is the rotation tensor with components

(2.12)\begin{equation} R_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}-\frac{\partial u_j}{\partial x_i}\right). \end{equation}

Straining motions without rotation make no contributions to $\lambda _2<0$, $\lambda _2$ is dimensionless and $\lambda _2$ magnitudes are a measure of only vortex rotation; $\lambda _2$ thus allows us to follow the transition from emerging vortex tubes, their interactions driving knot formation, their excitation of twist waves and the evolution to fully developed turbulence and its subsequent decay. Understanding the complex forms of KHI tubes and knots is significantly enhanced by viewing the evolving vortex structures in three-dimensional (3-D) volumes, especially from two viewing perspectives, as done in §§ 4.3.1, 4.3.2 and 4.4 below. Fritts et al. (Reference Fritts, Wang, Thorpe and Lund2022) also used $\lambda _2$ to reveal the evolution of tubes and knots due to misaligned KH billows employing a deep compressible model that was not a DNS.

4.1. Emergence of misaligned KH billows and vortex tubes

As noted above, the simulation was specified to have a large domain along the KH billow axes to allow exploration of the dynamics of interacting KH billows, the formation of tubes and knots, their breakdown to turbulence and the implications for KHI energetics, energy dissipation and mixing relative to instabilities at locations exhibiting relatively uniform KH billows. Emerging KH billow scales and phases arising from the weak initial noise fields are shown in figure 6 with $(x,y)$ cross-sections at $z=0$ of $\zeta$, $\lambda _2$ and $T^\prime /T_0$ spanning 3 $T_b$. Both $\zeta$ and $\lambda _2$ are employed because $\zeta$ displays all vorticity features whereas $\lambda _2$ shows only those vortical structures having rotational character. We also use $\zeta$ to denote $|\zeta |$ for convenience except where subscripted.

Figure 6. The $(x,y)$ cross-sections of $\zeta$, $\lambda _2$ and $T^\prime /T_0$ at $z=0$ from 0–3 $T_b$. Colour scales for $\zeta$ and $T^\prime /T_0$ are shown in each panel and are saturated for $\zeta$ by $\sim$2 and 3 times at 2 and 3 $T_b$. Here, $\lambda _2$ varies from 0 to the domain maximum at each time. Domain dimensions are in units of the nominal KH wavelength, ${\lambda _h}_0=L$.

The $\zeta$ fields at the first time shown, at which the maximum $|T^\prime /T_0|=0.003$ (denoted 0$T_b$), reveal weak, structured, small-scale features that are consistent with local optimal perturbations in sheared flows at early times that typically are oriented at $\sim$45$^{\circ }$ angles to the large-scale shear (Farrell & Ioannou Reference Farrell and Ioannou1996; Achatz Reference Achatz2005; Fritts et al. Reference Fritts, Wang, Werne, Lund and Wan2009). The initial $\lambda _2$ and $T^\prime /T_0$ fields do not exhibit the same structures because they are larger-scale integrated responses to global optimal perturbations that initiate the finite-amplitude KHI features at later times. These larger-scale responses correspond to the two KH wavelengths anticipated by linear stability theory (e.g. 3 or 4 KH billow $\lambda _h$ along $x$, denoted wavenumbers 3 and 4) for the chosen shear depth $h$, as described above.

The three fields in figure 6 exhibit similar larger-scale responses at 1$T_b$, but $\zeta$ and $\lambda _2$ still reveal smaller-scale responses from $y/L\sim 1.5\text {--}8.5$. Initially disjoint phases in $x$ (denoted $\phi$) begin to exhibit increasingly uniform $\phi$ variations along $y$ by 1$T_b$, though maintaining a larger-scale ${\rm d}\phi /{{\rm d}y}>0$ (due to the initial noise in a periodic $y$ domain). All fields show decreasing influences by wavenumber 4 KH billows at this time, but there remain multiple sites at which phases along $y$ are misaligned or discontinuous.

The KH billow phases exhibit increasingly uniform ${\rm d}\phi /{{\rm d}y}$ at most $y$ by 2$T_b$. See the third warm-to-cold phase transition in $T^\prime /T_0$ (the red/blue gradient with slowly varying $\phi$ along $y$ in figure 6) at $x/L\sim 1.5\text {--}2.5$, especially over $y/L\sim 0\text {--}4$. Corresponding features in $\zeta$ undulate somewhat, but are continuous along $y$ at this stage. The other KH billow phases seen in $T^\prime /T_0$ and $\zeta$ are locally more uniform as well, but exhibit rapid phase transitions or are discontinuous along $y$ at $z=0$; $\lambda _2$ continues to exhibit a lack of vortex structures at the sites of the strongest phase variations seen in the other fields.

All fields show evidence of narrowing phase transitions by 3$T_b$ at the sites exhibiting the largest initial, and most persistent, phase variations at earlier times; $\zeta$ has evolved thinner, more intense, boundaries between adjacent KH billows that correspond to the vortex sheets between them along $x$ at $z=0$. The regions exhibiting larger ${\rm d}\phi /{{\rm d}y}$ at 2$T_b$ are seen to evolve coherent vortices having ${\zeta }_z<0$ with large $\lambda _2$ at their outer edges. The $T^\prime /T_0$ fields also reveal vortices passing through $z=0$ having downward motions (red edges) at their positive $x$ edges that confirm vortex tubes also having $\zeta _x<0$ and $\zeta _y>0$ due to their confinement to the vortex sheets between KH billows.

The KH billow cores also exhibit significant structures at 3$T_b$ at sites adjacent to the vortex tubes noted above. The most obvious are seen in $\zeta$ and exhibit small-scale vorticity having expected links to the adjacent vortex tubes by analogy with the laboratory observations of emerging knots (see the dashed ovals at bottom and upper left in panels (b) and (c) of figure 3, respectively). Corresponding structures are also seen in $\lambda _2$ and more weakly in $T^\prime /T_0$ (see the kinks in the blue/green transitions). These sites will be seen below to be the major drivers of large-scale tube and knot dynamics.

4.2. Instability and turbulence evolutions of misaligned KH billows

We now explore the dynamics of KH billows that have variable phases along their axes, given the evidence for such dynamics from laboratory and atmospheric studies noted above. Our specific interest is the degree to which these dynamics accelerate KH billow breakdown and increase the impacts of KHI tube and knot dynamics relative to those without such influences. Because $\zeta$ at later times varies by orders of magnitude, we show cross-sections of these fields with a $\log _{10}\zeta$ colour scale hereafter.

Cross-sections of $\zeta (x,y)$ at $z=0$ and $\zeta (x,z)$ at $y/L=2$, 3, 4, 5.5 and 8.5 from 3-5$T_b$ are shown at top and bottom in figure 7, respectively. A subset of the $\zeta (x,y)$ domain including the strongest dynamics is shown in figure 8 from 2.75–4.25$T_b$. Corresponding $T^\prime /T_0(x,y)$ fields at $z=0$ and $T/T_0(x,z)$ fields at $y/L=2$, 3, 4, 5.5 and 8.5, with the latter extending to 10$T_b$ to show the relative tendencies for restratification, are shown in figure 9.

Figure 7. The $\zeta (x,y)$ and $\zeta (x,z)$ fields in (a) and (b), respectively, with those in (b) at $z/L=2$, 3, 4, 5.5 and 8.5 from 3–5$T_b$ at intervals of 0.5$T_b$ (left to right). Colour scales (upper right in panel (a) where they begin) are logarithmic and saturated to better reveal the evolving dynamics.

Figure 8. As in figure 7 at top for sub-sections of the $\zeta (x,y)$ fields from 2.5–4.25$T_b$ at 0.25$T_b$ intervals. See text for discussion of the highlighted features. Turbulence onset is considered to accompany the rapid emergence of resolved, but much smaller, vortices. Also see the discussion of figures 14 and 15 below.

Figure 9. As in figure 7 for $T^\prime /T_0(x,y)$ and $T/T_0(x,z)$ (a–c) showing the advective influences of the KH billows, their secondary CIs and KHIs and the tube and knot dynamics at $z=0$ and $y/L=2$, 3, 4, 5.5 and 8.5, respectively. The $T/T_0(x,z)$ panels are extended to 10$T_b$ at bottom in order to show the approach to a turbulent mixing layer at the various $y$ at later times. Colour scales (upper right in panel a) are shown where they begin.

Figure 7 provides an overview of the evolution from the onset of KH billow interactions and initial tube formation to widespread turbulence $\sim$2$T_b$ later. These dynamics are rapid and strong, yield ‘turbulence onset’ at $\sim$3.25$T_b$ (i.e. emergence of very small scales that are not resolved in the image cadence shown in figure 8) and require only $\sim$0.75$T_b$ to achieve peak $\zeta$. Importantly, the strongest turbulence transitions occur in regions of misaligned KH billows driving tube and knot dynamics. Secondary KHI and CI also exhibit stronger and more rapid interactions at these sites. In order of emerging turbulence intensities, the major contributors to KH billow breakdown in these regions are, (i) tube and knot dynamics where KH billows have large ${\rm d}\phi /{{\rm d}y}$ and vortex tubes link adjacent billows, (ii) secondary KHI on vortex sheets distorted by tube and knot dynamics and (iii) secondary KHI and CI interactions also influenced by tube and knot flow distortions.

Intensities of tube and knot dynamics driving strong vortex interactions and the rapid cascade to turbulence are not reflected in the colours of the images in figure 7 because the colour scales are saturated in order to reveal dynamics having intermediate $\zeta$, and the maxima increase by $\sim$20 times from 3 to 4$T_b$. The true maxima of the $\zeta$ fields at 3, 3.5, 4, 4.5, 5 and 5.5$T_b$ are $\sim$0.8, 7.5, 7.7, 7.2, 7.7 and 6.9, but the colour scale remains the same at 3.5$T_b$ and thereafter in order illustrate the decreasing peak turbulence intensities even as additional secondary KHI and CI continue to form and intensify. The influences of tube and knot dynamics are also revealed in the $\zeta (x,z)$ cross-sections in figure 7(b), where the cross-sections at $y/L=4$ and 5.5 reveal more rapid transitions to billow core turbulence and breakdown thereafter.

To further quantify the evolution seen in figure 7(a), expanded views of the most dynamic $\zeta (x,y)$ domain from 2.5–4.25$T_b$ at 0.25$T_b$ are shown in figure 8. Inspection of these panels reveals the following additional details of features seen in figure 7:

1. (i) billow core interactions where initial ${\rm d}\phi /{{\rm d}y}$ increases and drives differential rotation and kinking along their axes in time (see the billow cores at $y/L\sim 4.5\text {--}6$ at $2.5\text {--}3T_b$);

2. (ii) large-scale roll-up of the vorticity sheets into tilted vortex tubes where they arise between misaligned KH billow cores (see white arrows at 2.5–3.25$T_b$ in figure 8);

3. (iii) interactions of vortex tubes and billow cores initiate knots (see red rectangles at 2.75 and 3$T_b$ in figure 8;

4. (iv) a systematic intensification of the narrow vorticity sheets between adjacent billows (see pink arrows at 3$T_b$ in figure 8);

5. (v) small-scale billow pairing at $y/L\sim 2$ (see orange arrows in figure 8a);

6. (vi) initial turbulent breakdown of vortex tubes and knots (see pink rectangles from 3.25-3.75$T_b$ in figure 8);

7. (vii) secondary KHI on intensifying vorticity sheets $\sim$0.75-1$T_b$ after initial vortex tube formation and breakdown of vortex knots (see yellow arrows at 3.5 and 3.75$T_b$ in figure 8);

8. (viii) secondary CI in the outer KH billow cores emerging after initial secondary KHI, and occurring first in regions where KH billows are most distorted (see red arrows at 3.5 and 3.75$T_b$ and elsewhere at later times in figure 8);

9. (ix) turbulence expansion in and between billow cores in the regions of most intense billow core distortions (see yellow rectangles at 4 and 4.25$T_b$ in figure 8); and

10. (x) delayed and weaker secondary instabilities and turbulence where tube and knot dynamics do not initiate turbulence transitions (see $y/L=8.5$ in figure 8).

Of these various dynamics, the most significant in driving rapid transitions to turbulence are the evolutions of misaligned KH billows, their stretching of intermediate vortex sheets that drive nearly orthogonal vortex tubes and subsequent interactions and entwining of the vortex tubes and billow cores forming knots. These dynamics account for the accelerated emergence of intense turbulence compared with the evolutions of relatively undisturbed KH billows (such as at $y/L\sim 0\text {--}1$ and 8–9 spanning the $y$ domain boundary and at $x$ and $y$ where tubes and knots do not occur). For reference, these billows have depths consistent with the laboratory experiments by Thorpe (Reference Thorpe1973b) for $Ri=0.1$.

At 3.5$T_b$, turbulence is already strong, but is entirely confined to very small regions within emerging tube and knot dynamics. Secondary KHI in the stratified braids and secondary CI in the billow exteriors exhibit enhanced amplitudes in the vicinity of tube and knot dynamics, but only contribute to turbulence where they are entrained into these dynamics. Coherent, larger-scale tube and knot dynamics yield intense, billow-scale turbulence patches, but no billow breakdown or merging by 4$T_b$ (see figures 7 and 8). As noted in the overview of figure 7 above, turbulence achieves the maximum local $\zeta$ at these times, but will be seen in Fritts et al. (Reference Fritts, Wang, Thorpe and Lund2022) to achieve the maximum local $\epsilon$ from $\sim$4–6.5$T_b$, depending on whether or not the tube and knot dynamics were the source. The evolutions, vorticity dynamics, transitions to turbulence and relative importance of these various initial tube and knot dynamics for two specific cases are addressed in much greater detail in our exploration via 3-D volumetric imaging of $\lambda _2$ spanning their transitions from initial laminar flows to fully developed turbulence in § 4.3.

The most turbulent KH billows at $y/L\sim 3$, 4 and 5.5 exhibiting tube and knot dynamics are first to exhibit loss of coherence and initial breakdown to a sheared mixing layer (see these sites at 4.5$T_b$ in figure 7). Those at $y/L=5.5$ having the strongest initial tube and knot dynamics are seen in figure 7(b) to undergo turbulent billow breakdown up to $\sim$0.5$T_b$ faster than at $y/L=3$ and 4 and $\sim$1$T_b$ faster than the KH billow evolutions that are relatively undisturbed by tube and knot dynamics at the extrema of the $y$ domain. Specifically, the three billow cores at $y/L=8.5$ and one billow core at $y/L=2$, where the tube and knot dynamics had small influences remain laminar at 5.5$T_b$. Regions exhibiting the strongest initial tube and knot dynamics, $y/L\sim 5.5$, 4 and 3 (in order of decreasing influence), also exhibit the most rapid restratification of the turbulent shear layers at these locations. Also seen at $y/L=2$ is a billow-pairing event where billows exhibit displaced vertical positions along $x$ accompanying large-scale vertical motions induced by the tube and knot dynamics.

A complementary perspective on the dynamics beyond 3$T_b$ is provided by $T^\prime /T_0(x,y)$ and $T/T_0(x,z)$ fields shown in figure 9 that correspond to the $\zeta$ fields in figure 7. As $T$ is constant following parcel motions in a Boussinesq fluid (apart from thermal diffusion), $T^\prime /T_0$ and $T/T_0$ are close tracers of larger- and smaller-scale motions and mixing accompanying these KHI dynamics to late times; $T^\prime /T_0(x,y)$ and $T/T_0(x,z)$ in figure 9 reveal emerging responses to billow roll-up, core structure, breakdown and re-stratification in the absence and presence of tube and knot dynamics that are not revealed as clearly by $\zeta$ in figure 7.

Specific examples of additional insights provided by the $T^\prime /T_0(x,y)$ and $T/T_0(x,z)$ fields relative to those provided in figures 7 and 8 include the following:

1. (i) $T^\prime /T_0(x,y)$ and $T/T_0(x,z)$ provide much clearer evidence of the vertical displacements, billow entrainment and layering causing static stability variations within the billow cores than seen in $\zeta$;

2. (ii) $T^\prime /T_0(x,y)$ reveal the inferred $\zeta _-$ orientations of secondary KHI between adjacent billows at $z=0$, as for the initial vortex tubes noted above;

3. (iii) $T^\prime /T_0(x,y)$ and $T/T_0(x,z)$ reveal the scales and character of secondary CI within the billow cores, and regions and extents of mixing events, none of which are apparent in the $\zeta$ fields; and

4. (iv) mean $T^\prime /T_0(x,y)$ along $x$ at $z=0$ varying in $y$ seen in figure 9 revealing large-scale induced flows extending to late times.

The more significant implications of our discussion of figures 7 and 9 are that the tube and knot dynamics, where they arise, have major influences on the character, intensities and time scales of turbulence transitions relative to those occurring in their absence. These dynamics also suggest more widespread influences in adjacent regions not exhibiting initial tube and knot dynamics due to induced larger-scale motions along $x$, $y$ and $z$.

4.3. Tube and knot dynamics revealed by 3-D visualization

The fields displayed in figures 6–9 reveal several examples of the tube and knot initiation represented in the cartoons in figure 5, approximating sites A−D and E−F, respectively, in figure 3. These are highlighted in figure 10, where Region 1 shows the primary region of initial billow misalignments and Region 2 shows adjacent billows exhibiting weaker phase variations. The early stages in the evolution of the vortex dynamics in Region 1 are shown for reference in figure 11. As the billow cores intensify (blue sheets and arrows in figure 11), they extend along $y$. Where these billow ‘ends’ are misaligned in $x$ at adjacent $y$, they cause differential stretching of the intermediate (weaker and unseen) vortex sheet along $x$ and $z$ yielding vortex tube formation exhibiting downward and rearward (upward and forward) advection toward smaller (larger) $x$ and $z$ and vortex tube intensification having $\zeta _\pm$ (denoting $\zeta _+$ and $\zeta _-$, respectively, see figure 5) for the KH billow end at larger (smaller) $x$ and $y$ given by

(4.1)\begin{equation} \frac{\partial \zeta_\pm}{\partial t}\simeq\zeta_\pm\frac{\partial \zeta_\pm}{\partial r_\pm}\end{equation}

for $r_\pm$ the distance along $\zeta _\pm$.

Figure 10. As in figure 8 for $y/L=2.5\text {--}6.5$ at 3.5$T_b$. Four regions in which 3-D visualization of tube and knot dynamics are presented below are shown with yellow, pink, white, and red rectangles labelled Regions 1, 2, 3 and 4, respectively.

Figure 11. Three-dimensional imaging of $\lambda _2$ in the yellow sub-domain in figure 10 (Region 1) and extending over $|z/L|\le 0.2$ from the emergence of misaligned KH billow cores to strong interactions among the billow cores and vortex tubes leading to knots. Times are shown in each panel. Peak vorticity increases slowly at early stages, but increases by $\sim$12 times from 3 to 3.5$T_b$. The more intense features at each time represent <0.1 % of the displayed volume. The colour scales vary from red/orange/yellow/green/blue from large to small negative $\lambda _2$ at each time. The $\lambda _2$ colour and opacity scales also vary between successive images in this and subsequent figures in order to span the range of features and scales appropriate for each local region and evolution.

Vortex tubes having $\zeta _+$ are shown with the green arrows at 2.5$T_b$. Those having $\zeta _-$ have roughly orthogonal orientations and attach to the trailing KH billow end below and rearward (yellow arrows at 2.5$T_b$). The initial tube and knot features (white arrows at 2.75$T_b$) exhibit roughly orthogonal vortices in close proximity that yield intense stretching (dashed arrows at 3$T_b$), or axial compression, of both vortices and intensifying interactions, enstrophy cascades and enhanced energy dissipation thereafter. Successive images in figure 11 show further intensification of the initial vortex tubes and knots, both knot and initial KH billow core breakup due to twist wave generation, and the emergence of secondary KHI and tube and knot dynamics on the intensified and distorted vortex sheets in close proximity to the initial tube and knot dynamics. Also see the more detailed discussions of these tube and knot dynamics in §§ 4.3.1–4.3.3 and 4.4.)

The results of these dynamics are the mutual intensification of the kinked billow cores and orthogonal vortex tubes yielding wrapping of these features at the initial billow core ‘ends’ and formation of features referred to as knots in the early laboratory shear-flow studies by Thorpe (Reference Thorpe1985, Reference Thorpe1987, see figure 3). These knots entwine the intensifying, orthogonal vortex tubes, and drive them increasingly into closer proximity. This results in dramatic accelerations of these interactions and induced differential axial and radial advection of the vortex cores that result in excitation of Kelvin twist waves at larger scales initially and successively smaller scales rapidly thereafter (see figure 11 at 3.25 and 3.5$T_b$). The latter panels reveal the evolution from a laminar initial flow to intense smaller-scale interactions within 0.5$T_b$ measured with respect to the large initial $N$ at $z=0$. These dynamics are explored in detail below for (i) the single vortex tube linking two adjacent KH billow cores in the pink sub-domain (Region 2) and (ii) the two vortex tubes linking to a single KH billow core in the white sub-domain (Region 3) in figure 10.

4.3.1. Tube and knot dynamics of one tube linking two KH billows

Tube and knot dynamics accompanying the evolution of the single vortex tube linking adjacent KH billow cores in Region 2 is shown with 3-D imaging viewed from above in figure 12 and from above, more positive $x$, and more negative $y$ in figure 13. These figures employ colour and opacity scales that vary in time as $\lambda _2$ magnitudes increase and between the two views at common times. This is because no colour and opacity scales allow inter-comparisons spanning the weak initial magnitudes, the very large range of magnitudes accompanying strong turbulence at later times, and only relative $\lambda _2$ magnitudes are relevant at each stage of the evolution. Supplementary Movie 1 and Supplementary Movie 2 provide animations of these figures at https://doi.org/10.1017/jfm.2021.1085.

Figure 12. As in figure 11 for $\lambda _2$ from 2.75–4$T_b$ viewed from above spanning the interval of strong initial KH billow and tube interactions and initial knot dynamics. Earlier times are shown in the larger $y$ portion of Region 2. Labels denote features described in § 4.3.1.

Figure 13. As in figure 12 from 3.1–4$T_b$ viewed from larger $x$, smaller $y$ and positive $z$.

The $\lambda _2$ fields at 2.75$T_b$ in figures 12 and 13 reveal that both KH billow cores aligned approximately along $y$ and the inclined vortex tube exhibit differential rotation along their axes due to axial stretching and compression where they are in close proximity. The vortex tube rotation, in particular, has increased dramatically relative to the billow cores by this time due to stretching of the intermediate vortex sheet described above and their increasing proximity thereafter. These initial large-scale dynamics drive an accelerating cascade of vortex interactions leading to intense, fully developed turbulence within 1$T_b$. Specific features of the subsequent dynamics include the following:

1. (i) the vortex tube is stretched and intensified where it is advected over and under the billow cores forming knots as the billow amplitudes increase (see sites labelled a from 2.75–3.1$T_b$ in figures 12 and 13);

2. (ii) vortex and billow core interactions in the knot regions begin to drive smaller-scale perturbations that manifest as twist waves at 3.2 and 3.3$T_b$ (sites b) that are much stronger than the larger-scale twist waves emerging in the billow cores in figure 11; they also induce vortex tube perturbations that intensify in time (site c);

3. (iii) seen emerging in the knot regions by 3.3 $T_b$ are smaller-scale vortex features in and external to the knots that had their roots in the initial distorted vortex sheet that is intensified due to stretching by the evolving knots (sites d from 3.4–3.6$T_b$);

4. (iv) a sheath of small-scale vortices is seen wrapping around the vortex tubes that are distinct from the induced twist waves and result from intensification of small-scale KHI on the evolving and entraining vortex sheet beginning at 3.3$T_b$ (site e) noted in the discussion of figure 11, but seen more clearly here;

5. (v) larger-scale twist waves arise due to unravelling of perturbed vortex tube ends that propagate inward toward the centre of the vortex tube extending from $\sim$3.3–3.6$T_b$ (sites d at larger scales);

6. (vi) the initial vortex tube weakens and begins to disappear due to fragmentation accompanying larger-scale mode-2 twist wave excitation by $\sim$3.5$T_b$;

7. (vii) secondary KHI emerge primarily along, but external to, the initial KH billow cores by 3.5$T_b$ (sites f);

8. (viii) weak secondary CI also arise having initial streamwise and vertical alignments in the billow exteriors by 3.5$T_b$ (sites g);

9. (ix) the knot regions remain the sites of most intense turbulence and exhibit rapid cascades of fully resolved, smaller-scale vortex dynamics that yield a broad turbulence inertial range with $\eta < L$/1000 by 4$T_b$;

10. (x) the turbulence fields exhibit discernible vortex structures extending down to scales of <0.01$L$ (thus $\sim$6$\Delta x$ at these times); and

11. (xi) small-scale secondary KHI exhibit interactions yielding secondary tube and knot dynamics in close proximity to the larger-scale knots by 3.6–3.8$T_b$ (sites h) and are explored in greater detail below.

Our discussion above reveals both a significant acceleration of KHI breakdown to turbulence due to tube and knot dynamics and apparent major increases in the scales and intensities of the dominant turbulence events in these regions relative to KH billows not undergoing tube and knot dynamics. A second region of tube and knot dynamics having different character is examined below to assess the universality of these dynamics within the larger-scale multiple KH billow environment.

4.3.2. Tube and knot dynamics of two tubes linking to one KH billow

Tube and knot dynamics arising where two tubes link to one KH billow core are shown in figures 14 and 15 employing similar viewing perspectives to those in § 4.3.1, but for Region 3 shown with the white rectangle in figure 10. The 3-D fields in figures 14 and 15 reveal that these dynamics emerge at very nearly the same time and span roughly the same interval as those described in § 4.3.1. They exhibit similar tube and knot evolutions and influences on the initial larger-scale vortex features in many respects, but there are also quantitative differences having implications for turbulence evolutions and intensities. For ease of comparisons of the two types of vortex dynamics, we employ the same feature labelling as in § 4.3.1. Supplementary Movie 3 and Supplementary Movie 4 provide animations of these figures.

Figure 14. As in figure 12 for $\lambda _2$ in Region 3 highlighted by the white rectangle in figure 10 viewed from above and spanning times from 2.75–4$T_b$ (see axes at 2.75$T_b$). Note the more rapid excitation of larger-scale twist waves that drive faster breakdown of the billow core and vortex tubes relative to those seen in Region 2 in figures 12 and 13.

Figure 15. As in figure 13 for $\lambda _2$ in Region 3 in figure 10 from 3–4$T_b$ viewed from larger x, smaller y and larger z, and spanning the interval of strong initial KH billow and tube interactions and initial knot dynamics. Earlier times are shown in the upper portion of the full subdomain.

Specific features of these dynamics, and differences from those in Region 2 discussed in § 4.3.1, include the following:

1. (i) vortex tubes at smaller (larger) $x$ (sites c) wrap over (under) the KH billow core because of their formation and intensification on the vortex sheet between adjacent billows (sites a in figures 14 and 15 which closely resemble sites a in Region 2);

2. (ii) the tubes and billow core exhibit mutual stretching (compression) along their axes due to differential advection by each other (sites a and b at 2.75 and 3$T_b$), but stretching by closely spaced tubes acts more strongly on the KH billow core in Region 3;

3. (iii) sheaths of smaller-scale vortices are seen wrapping around both vortex cores due to entraining of the initial vortex sheet having embedded and enhanced small-scale secondary KHI, as in Region 2 (see sites c from 3–3.7$T_b$);

4. (iv) KH billow core stretching between the vortex tubes in Region 3 increases its rotation and enables $\lambda _2$ revealing the emergence of billow core twist waves at 3.2 and 3.5$T_b$ not seen in Region 2;

5. (v) large-scale twist wave propagation along the vortex tubes away from their source regions (sites b and d) toward negative (positive) $x$ along the tube at smaller (larger) $y$, that are much more dramatic than seen in Region 2;

6. (vi) smaller-scale twist waves also arise on larger-scale twist waves where they are perturbed by other vortex features (sites d at $\sim$3.5$T_b$ and beyond);

7. (vii) as in Region 2, secondary KHI form and initiate interactions and secondary tube and knot dynamics on the vortex sheet deformed and intensified by the larger-scale tube and knot dynamics (sites f and g); and

8. (viii) as in Region 2, the knots in Region 3 remain the sites of the most intense turbulence, achieve a broad turbulence inertial range with $\eta < 3L/1000$ and exhibit clear vortex structures extending down to scales of <0.01$L$ by 4$T_b$.

In summary, the tube and knot dynamics in Region 3 exhibit several clear differences relative to those occurring in Region 2. The most significant arise from the stretching by two vortex tubes of the common KH billow core, which enables much stronger interactions among larger-scale vortices in close proximity having roughly orthogonal alignments. In Region 2, in contrast, the vortex tube ‘ends’ interact with relatively weaker orthogonal billow cores. These differences clearly account for the very different twist wave formation dynamics seen in these two regions, and the larger energy dissipation rates in Region 3 as described in the companion paper by Fritts et al. (Reference Fritts, Wang, Thorpe and Lund2022).

4.3.3. Secondary KHI and CI evolutions influenced by adjacent tube and knot dynamics

We noted the tendency for enhanced secondary KHI and CI in close proximity to emerging tube and knot dynamics in the discussions of figures 11–15 in Regions 1, 2 and 3. To explore these dynamics in greater detail, figure 16 provides an expanded time series of figure 12(f–h) from 3.6 and 4$T_b$ using optimized $\lambda _2$ and opacity scales in the Region 4 red rectangle in figure 10. These show the region between the two initial KH billows at smaller $x$ and $y$ relative to the large initial vortex tube at upper right linking the KH billows in Region 2. For reference, these images span only $\sim$5 % of the full horizontal domain seen at top in figure 7.

Figure 16. As in figure 12 for $\lambda _2$ in the lower left portion of Region 2 shown with the red rectangle in figure 10 (Region 4) viewed from above and spanning the interval from 3.6–4$T_b$. It includes the right portion of the large-scale KH billow that is relatively uniform along $y$ at left and the left edge of the KH billow at right having ${\rm d}\phi /{{\rm d}y}>0$.

As discussed in § 4.3.1, the upper left and lower right portions of this region exhibit intensifying to strongly interacting secondary KHI (sites a) where the vortex sheet is enhanced and deformed by the KH billows in close proximity. Weaker secondary KHI emerge on the less enhanced vortex sheet between adjacent KH billows (site b) at 3.6$T_b$. These vortex evolutions progress rapidly thereafter, and we employ different site labels with advancing subscripts to track the evolutions of individual features in time. Features labelled c–i in the panel at 3.7$T_b$ evolving to features labelled ${\rm c}_3-{\rm i}_3$ in the panel at 4$T_b$ reveal often significant evolutions of these small-scale dynamics over only 0.3$T_b$.

Key features of the secondary KHI, vortex tubes and secondary CI evolutions in this region include the following:

1. (i) initial secondary KHI have primary $\zeta _-$ alignments opposite that of the initial KH billow cores driving the large-scale tube and knot dynamics (see figures 11 and 12 at 2.75$T_b$ and sites a);

2. (ii) secondary CI in the outer billows are roughly orthogonal to the billow cores, see the blue/green features at sites b, i and j outside the tube and knot dynamics where the KH billows are more nearly aligned along $y$, but they fail to compete with the secondary tube and knot dynamics in these regions;

3. (iii) the secondary KHI and vortex tubes exhibit variations in their orientations along their axes at earlier times due to weak perturbations of the deforming initial vortex sheet (see sites c, d and e at 3.7$T_b$ and figure 11 at 3.5$T_b$);

4. (iv) the secondary KHI having $\zeta _-$ are more nearly continuous along their axes, whereas the vortex tubes having $\zeta _+$ alignments instead terminate at the orthogonal $\zeta _-$ tubes, as also seen in the initial large-scale dynamics, but rotated $\sim$90$^\circ$ clockwise;

5. (v) secondary KHI and vortex tubes with $\zeta _-$ and $\zeta _+$ alignments intensify due to vortex sheet stretching, lead to an extended ‘web’ of nearly orthogonal secondary billow cores and vortex tubes, and rapid evolutions of small-scale tube and knot dynamics mirroring the initial large-scale tube and knot alignments, but at scales $\sim$10 times smaller, and thus strongly constrained by a very small $Re$ (see sites $c_i-h_i$ for $i=1\text{--}3$ from 3.8−4$T_b$); and

6. (vi) emerging complexity in the outer KH billows not undergoing tube and knot dynamics where secondary CI and KHI interact and lead to mutual stretching and deformation, but at very small scales that exhibit relatively small intensities due to their very small $Re$ (see sites i$_i$ and j$_i$ from 3.7–4$T_b$).

These results reveal that KHI tube and knot dynamics are widespread across scales within individual events, and that they arise readily wherever initial vortex sheets intensify and form KHI for sufficiently small local $Ri$, even when their scales and intensities are constrained by low local $Re$. More importantly, secondary KHI tube and knot dynamics contribute weakly to turbulence intensities where these smaller KHI, by themselves, would fail to yield turbulence. This suggests that such dynamics may contribute to turbulence and mixing where traditional KHI is not expected to do so.

4.4. Secondary KHI and CI evolutions without strong tube and knot influences

As noted above, the absence of significant KHI ${\rm d}\phi /{{\rm d}y}$ at small and large $y$ enables KH billow breakdown via secondary KHI and CI largely without influences of large-scale tube and knot dynamics. These secondary instabilities arise from the same initial noise seed, but where it does not induce significant KHI phase variations leading to large-scale tube and knot dynamics. The primary KHI and secondary instability dynamics centred at $y/L=8.5$ are shown for the two central billow cores along $x$ with 3-D imaging of $\lambda _2$ viewed from above and from larger $x$, smaller $y$ and positive $z$, respectively, at left and right in figure 17 at the later times discussed above. The later emergence of these secondary instabilities, and their relatively slow progressions to initial turbulence compared with the tube and knot dynamics discussed in § 4.3, are our focus here. Supplementary Movie 5 and Supplementary Movie 6 provide animations of this figure.

Figure 17. As in figure 12 for $\lambda _2$ from $x/L=0.5\text {--}2.5$ and $y/L=8\text {--}9$ (a,c,e) showing emergence of initial secondary instabilities yielding initial transitions to turbulence. The same 3-D fields are shown at right viewed from larger $x$, smaller $y$, and positive $z$ (b,d,f). See text for details.

The 3-D imaging of $\lambda _2(x,y)$ viewed from above at left in figure 17 from 3.5–3.9$T_b$ exhibits two KH billows having coherent cores with small ${\rm d}\phi /{{\rm d}y}$ and ${\rm d}^2\phi /{{\rm d}y}^2<0$ and $\sim$0 at smaller and larger $x$ at 3.5$T_b$, and rapid decreases in billow core $\lambda _2$ relative to secondary KHI and CI thereafter. The imaging also reveals the emergence of secondary KHI and CI that drive secondary tube and knot events that dominate the weak turbulence transitions at later times. Specific features seen in figure 17 include the following:

1. (i) evidence of initial weak secondary KHI and CI arises in the vortex sheets and the KH billow interiors (see sites a and b, respectively, at 3.5$T_b$);

2. (ii) increasing secondary KHI amplitudes exhibit distinct $\zeta _-$ and $\zeta _+$ and undergo initial interactions at 3.7$T_b$ having the same forms seen in regions adjacent to the larger-scale tube and knot dynamics (see sites c, and compare with figure 16 at earlier times);

3. (iii) small-scale secondary CI in the KH billow exteriors intensify dramatically by 3.7$T_b$, comprise counter-rotating vortices that wrap around the billow cores and induce intensifying vortex sheets where vortex pairs stretch and distort adjacent, small-scale vortex sheets having primarily $\zeta _y$ (see sites d and e);

4. (iv) interactions among secondary KHI and vortex tubes on the vortex sheets between and around the KH billows induce smaller-scale knots, twist wave generation, and initial small-scale twist waves driving turbulence (see sites f in figure 17 at 3.9$T_b$);

5. (v) interactions between secondary KHI and CI at the outer edges of the KH billows cause intensification of $\lambda _2$ where CI vortices cross secondary KHI on entraining vortex sheets (see site g and others that are similar at 3.9$T_b$); and

6. (vi) regions where interacting secondary KHI and CI induce further tertiary interactions comprising largely orthogonal vortices in close proximity that are very small-scale versions (scales of $\sim$0.03$L$) of primary and secondary tube and knot interactions discussed above, but now at dramatically smaller $Re$ (see sites h).

This discussion of the interactions and dynamics of secondary instabilities of individual KH billows reveals features of these interactions that could not be inferred in the analyses of the larger-scale tube and knot dynamics in § 4.3. This is due to the inability of these larger-scale, more intense dynamics to also reveal the much weaker dynamics at much smaller scales, despite the expectation that they will also occur. These various instability events are compared in spectral assessments spanning their evolutions in § 5. Their potential implications for enhanced turbulence and mixing are discussed in § 6.