3.1. Transition mechanism

We begin by examining the flow evolution for the case where $({\textit{Re}}, {\textit{Fr}}) = (1600, 1.0)$ , with a focus on the onset of turbulence and its driving mechanisms. Figure 1 shows the evolution of the vertical vorticity, highlighting key stages in the flow’s development. As the initial, weak disturbances to the base flow grow via zigzag instability, the columnar vortices become visibly distorted by $t \leqslant 44$ , introducing variations in the vertical direction. By $t = 48$ , localised, intense patches of prolonged vortical structures, which we will consider as the finite-amplitude manifestation of the primary zigzag instability, appear between the vortex columns and quickly trigger turbulence. The turbulence reaches its peak intensity by $t = 56$ , and subsequently decays over time as the initial columnar vortices disappear completely.

Figure 1. Volume-rendering images of the vertical vorticity field $\omega _z=\partial v/ \partial x- \partial u/ \partial y$ , from the DNS with ${\textit{Re}}=1600$ and ${\textit{Fr}}=1.0$ . The first image in the sequence was captured at time $t=36$ . Subsequent images are spaced at regular intervals of 4.0 units of time, arranged from left to right and top to bottom. An animation of the flow evolution can be found as movie 1 in the supplementary movies available at https://doi.org/10.1017/jfm.2025.10420.

In figure 2, we analyse the evolution of the vertically averaged and fluctuating components of kinetic energy. The flow field $\boldsymbol{u}=(u,v,w)$ is decomposed into the vertically averaged component $\overline {\boldsymbol{U}} = (\overline {U}, \overline {V}, 0)$ and the fluctuation $\boldsymbol{u}' = (u', v', w')$ . The vertically invariant and fluctuating kinetic energy are defined as

(3.1) \begin{equation} \overline {E}_k = \tfrac {1}{2} \langle \overline {\mathbf U} \boldsymbol{\boldsymbol{\cdot}} \overline {\mathbf U} \rangle , \quad E_k^{\prime} = \tfrac {1}{2} \langle \mathbf u' \boldsymbol{\boldsymbol{\cdot}} \mathbf u' \rangle , \end{equation}

where $\langle {\boldsymbol{\cdot}} \rangle \equiv (2\unicode{x03C0} )^{-3} \iiint ({\boldsymbol{\cdot}} )\, \mathrm{d}x\, \mathrm{d}y\, \mathrm{d}z$ denotes the volume average over the 3-D domain. Figure 2(a) shows that the vertically invariant component $\overline {E}_k$ dominates the kinetic energy for $t \lesssim 40$ . After this point, the vertically varying component grows rapidly, extracting energy from the vertically uniform background flow. In figure 2(b), where $E_k^{\prime}$ is replotted on a logarithmic scale, its linear growth can be observed approximately between $t = 15$ and $t=45$ . Figure 2(c) plots the volume-averaged kinetic energy pseudo-dissipation rate $\varepsilon _k = \langle \epsilon _k \rangle$ , where

(3.2) \begin{equation} \epsilon _k = \frac {1}{\textit{Re}} \frac {\partial {u_i}}{\partial {x_j}} \frac {\partial {u_i}}{\partial {x_j}}. \end{equation}

As is common in similar numerical studies (e.g. Lucas, Caulfield & Kerswell Reference Lucas, Caulfield and Kerswell2019; Howland, Taylor & Caulfield Reference Howland, Taylor and Caulfield2020), pseudo-dissipation is used in place of the full dissipation rate from the turbulent kinetic energy budget (see e.g. equation 1 in Caulfield Reference Caulfield2021) for ease of computation. The value of $\varepsilon _k$ reaches its peak at $t \approx 56$ , shortly after the fluctuating kinetic energy $E_k^{\prime}$ reaches its maximum at $t \approx 54$ .

Figure 2. Time series of (a,b) kinetic energy components and (c) dissipation rate for the case with ${\textit{Re}}=1600$ and ${\textit{Fr}} = 1.0$ . Circles in (b,c) correspond to $t=40$ , 44 and 48.

To understand the mechanisms behind the onset of turbulence, we focus on three specific time instances: $t=40$ , $44$ and $48$ , which are marked with circles in figure 2(b,c). During these times, the primary instability is approaching the end of its linear growth phase, as indicated in figure 2(b). At the same time, the dissipation is beginning to increase significantly, as shown in figure 2(c). Additional visualisations of these instances are provided in figure 3, where the isosurfaces of vertical vorticity are plotted (in red and blue) in conjunction with regions that are convectively unstable, i.e. $\partial \rho _{\textit{tot}} / \partial z \gt 0$ , highlighted in green. In this set-up, the vertical gradient of total density is given by

(3.3) \begin{equation} \dfrac {\partial \rho _{\textit{tot}}}{\partial z} = \dfrac {\partial \rho }{\partial z} -1. \end{equation}

This expression consists of the vertical gradient of density perturbation $\rho$ , which evolves following (2.1c ), as well as the background density gradient, $-1$ .

Figure 3. A 3-D visualisation of the onset of turbulence in a columnar TG vortex array ( ${\textit{Re}}=1600$ , ${\textit{Fr}} = 1.0$ ). The red and blue surfaces represent isosurfaces of vertical vorticity $\omega _z = 0.75$ (red) and $\omega _z = -0.75$ (blue). Convectively unstable regions, where $\partial \rho _{\textit{tot}} / \partial z \gt 0$ , are highlighted in dark green. The snapshots are shown at times (a) $t = 4$ , (b) $t = 40$ , (c) $t = 44$ and (d) $t = 48$ .

Figure 3 highlights key features of the zigzag instability after reaching finite amplitude. By $t=40$ , the isosurfaces of $\omega _z$ (in red and blue) show corrugations, and small linearly shaped patches of convectively unstable regions (in green) have appeared in the gaps between adjacent vortex columns. At $t=44$ , these corrugations have deepened, and the convectively unstable regions have expanded, particularly near the ‘creases’ of the isosurfaces of $\omega _z$ . By $t=48$ , more regions have become convectively unstable, and a closer look at the visualisation reveals structures resembling overturning ‘billows’ associated with shear-induced instabilities, as indicated by the two white arrows in figure 3(d).

The detailed flow structure during the onset of turbulence is further examined in figures 4 and 5. Figure 4 shows $x$ – $z$ cross-sections of the flow field at $y=\unicode{x03C0} /2$ , passing through the centres of two adjacent columnar vortices, indicated by the red dashed lines in figure 3(a). Distortion of the vortex columns is evident in figure 4(a), which visualises the vertical vorticity $\omega _z$ . The complex structures in the $\omega _z$ field between the vortex columns resemble instabilities found in horizontally sheared flows, such as those illustrated in figure 14 of Basak & Sarkar (Reference Basak and Sarkar2006) and figure 5 of Lewin & Caulfield (Reference Lewin and Caulfield2024). Pairs of counter-rotating vortices form between the vortex columns. Figure 4(b), which shows $\partial \rho _{\textit{tot}}/\partial z$ , reveals that these vortices are associated with density overturns, similar to the convective instability mechanism that arises from interactions between oppositely signed vortex columns, as illustrated by Waite & Smolarkiewicz (Reference Waite and Smolarkiewicz2008) in their figure 10. Figure 4(c,d) present the velocity field around the overturns. A detailed examination of the $w$ field indicates that isopycnal displacement is correlated with vertical velocity. Furthermore, significant vertical shear in the horizontal velocity $u$ also develops in association with the overturns.

Figure 4. Visualisation of $x$ – $z$ transects showing (a) the vertical vorticity $\omega _z$ , (b) the vertical gradient of $\rho _{\textit{tot}}$ , (c) the vertical velocity $w$ , and (d) the horizontal velocity $u$ for the case with ${\textit{Re}}=1600, \textit{Fr} = 1.0$ . The transects are taken along $y=\unicode{x03C0} /2$ , as indicated by the red dashed lines in figure 3(a). Black solid lines represent isopycnals, spaced at intervals $\Delta \rho _{\textit{tot}}=1/3$ . Snapshots are shown for times $t=40$ , $44$ and $48$ .

Figure 5. (a) Streamlines (represented by solid lines with arrows) at $y=\unicode{x03C0} /2$ for the case with ${\textit{Re}}=1600$ and ${\textit{Fr}}=1.0$ . (b) A zoomed-in view of the region outlined by the blue box in (a). Red solid lines represent isopycnals, spaced at intervals $\Delta \rho _{\textit{tot}}=0.16$ .

The detailed flow structure around the overturns is visualised in figure 5, where streamlines are plotted over the transects shown in figure 4. A schematic of this mechanism is presented in figure 6. In G24, we found that the zigzag instability generates vertical velocity at the edges of the vortex columns (see e.g. figures 4 and 5 of G24). This vertical velocity is small at a neutral, horizontal plane, such as $z\approx 1.5\unicode{x03C0}$ , as shown in the zoomed-in view in figure 5(b). At this neutral plane, the horizontal velocity reaches its maximum, separating two counter-rotating circulation cells, and causing distortion to the vortex columns (see figure 4 a). As shown in figure 6, the non-propagating zigzag instability grows over time, with the increasing vertical velocity causing larger isopycnal displacements. The system’s stability is further weakened by vertical shear in the horizontal velocity, which steepens the isopycnals and leads to their eventual overturning, as illustrated in the third sketch in figure 6. This process is also visible in figure 5(b), where isopycnals (in red) are plotted over the streamlines. Between $t=40$ and $t=44$ , the isopycnal displacements reach a sufficient amplitude to overturn due to the combined effects of vertical velocity $w$ , which displaces the isopycnals, and horizontal velocity $u$ , which generates shear and contributes to the roll-up of the ‘billows’.

Figure 6. Schematic illustrating the overturning mechanism. Arrows represent the fluid velocity, while the lines correspond to two isopycnal surfaces located above and below the neutral plane where $w=0$ . The heights of neutral planes are marked by grey dashed lines. Time progresses from left to right.

Figure 7. Volume-rendering images of (a) the vorticity component $\omega _y$ , and (b) the normal strain rate $\partial v/\partial y$ , both along the $y$ -axis, for the case with ${\textit{Re}} = 1600$ and ${\textit{Fr}} = 1.0$ at $t = 40$ . The green arrows indicate the overturning billows highlighted in figure 3(d).

Figure 8. Sample scatter plots of squared vertical shear $S^2$ versus vertical density gradient $\partial \rho _{\textit{tot}}/\partial z$ for times corresponding to the snapshots in figure 3. The colour scale is adjusted to display only points where $\widehat {{\textit{Re}}_b} \geqslant 1$ . The dotted vertical line indicates $\partial \rho _{\textit{tot}}/\partial z = -1$ , and the dashed horizontal line marks $S^2 = 4/{\textit{Fr}}^2$ , corresponding to the shear strength required for the gradient Richardson number to reach $1/4$ under background uniform stratification. Regions I, II and III, as labelled in (a), delineate different potential instability mechanisms: regions I and II correspond to strong shear ( $S^2 \gt 4/{\textit{Fr}}^2$ ), while regions II and III correspond to weakened (or even reversed) vertical density gradient from the background value ( $\partial \rho _{\textit{tot}}/\partial z \gt -1$ ).

Figure 7 provides perspective on how such billows, which are essentially vortex tubes aligned along horizontal axes, amplify within the TG base flow. As previously noted, the zigzag mode generates vertical velocity at the edges of the vortex columns, as shown in figures 4(c) and 5(d) of Guo et al. (Reference Guo, Taylor and Zhou2024). This vertical velocity is sheared horizontally, producing horizontal vorticity. In the case shown in figure 7(a), we focus on the $\omega _y$ component, the vorticity aligned along the $y$ -axis, generated as a result of the zigzag deformation. Figure 7(b) shows that the locations of these horizontal vortex tubes largely overlap with regions of positive normal strain along the $y$ -axis, i.e. where $\partial v / \partial y \gt 0$ , induced primarily by the TG base flow. This spatial coincidence suggests that vortex stretching occurs: the normal strain in the base flow amplifies the horizontal vorticity structures generated by the zigzag instability. The vortex-stretching-driven amplification strengthens the horizontal vortex tubes, eventually leading to the density overturns and contributing to the breakdown of the base flow into turbulence.

To quantitatively examine the onset of turbulence caused by the zigzag instability, we present scatter plots of the local vertical shear and density gradient at representative times, shown in figure 8. The squared vertical shear of horizontal velocities is given by

(3.4) \begin{equation} S^2 = \left ( \frac {\partial u}{\partial z} \right )^2 + \left ( \frac {\partial v}{\partial z} \right )^2. \end{equation}

From a stability perspective, large values of both $S^2$ and the density gradient $\partial \rho _{\textit{tot}} / \partial z$ indicate destabilising effects on the local dynamics. The scatter plots between $S^2$ and $\partial \rho _{\textit{tot}} / \partial z$ are colour-coded to highlight the most dissipative regions in the flow, measured by the local (pointwise) buoyancy Reynolds number

(3.5) \begin{equation} \widehat {{\textit{Re}}_b} = \epsilon _k\, {\textit{Re}}\, {\textit{Fr}}^2. \end{equation}

A dotted vertical line in figure 8 marks the background value $\partial \rho _{\textit{tot}} / \partial z = -1$ . At $t=40$ , the $\widehat {{\textit{Re}}_b}$ values are $O(10)$ or smaller, and the most dissipative regions (marked by yellow points) have $\partial \rho _{\textit{tot}} / \partial z$ values around $-1$ , as shown in figure 8(b), indicating that the isopycnals are only weakly disturbed from the background stratification.

As the perturbations strengthen at $t=44$ and $48$ , shown in figure 8(c,d), a larger fraction of the flow exhibits both significant vertical shear ( $S^2 \gt 4/{\textit{Fr}}^2$ ) and weakened stratification from the uniform background ( $\partial \rho _{\textit{tot}} / \partial z \gt -1$ ). These regions are associated with enhanced dissipation, with pointwise $\widehat {{\textit{Re}}_b}$ reaching $O(10^2{-}10^3)$ . In terms of the regions defined in figure 8(a), the most intense turbulence develops predominantly within region II, where both shear and convective mechanisms could coexist.

At $t=48$ , the highest dissipation regions (dark red points, $\widehat {{\textit{Re}}_b} \gt 100$ ) cluster near areas of strong vertical shear ( $S^2 \gtrsim 4/Fr^2$ ) and weak or reversed density gradients ( $\partial \rho _{\textit{tot}} / \partial z$ approximately zero). This suggests that the transition to turbulence involves a combined mechanism: vertical shear and convective overturning are not isolated but often collocated in space, reinforcing each other to trigger vigorous turbulence.

Thus far, our discussion has primarily focused on the representative case $({\textit{Re}}, {\textit{Fr}}) = (3200, 1.0)$ , where the instability mechanism characterised by overturning motions and vertical displacement of isopycnals closely resembles the convective-type instability described by Augier & Billant (Reference Augier and Billant2011), specifically in runs denoted as asterisks in their parameter space plot (figure 4). We observe this convective mechanism consistently across simulations with ${\textit{Fr}} \geqslant 0.5$ , suggesting that it is a robust feature of the moderately stratified regime examined here. In contrast, we do not find clear evidence of the separate shear instability identified by Augier & Billant (Reference Augier and Billant2011) in our runs with ${\textit{Fr}} \geqslant 0.5$ . At lower Froude numbers ( ${\textit{Fr}} \leqslant 0.25$ ), our simulations begin to show signs of inclined shear layers that may serve as precursors to shear-induced instabilities. A more detailed investigation of shear-driven transitions is deferred to future work, once additional simulations in the low- ${\textit{Fr}}$ , high- ${\textit{Re}}$ regime become available.

3.2. Scaling of buoyancy Reynolds number and turbulent Froude number

Our next goal is to identify the regions in the ${\textit{Re}}{-}{\textit{Fr}}$ parameter space where the TG vortex array can generate vigorous turbulence. To quantify the turbulence intensity in a given TG flow, we use the volume-averaged buoyancy Reynolds number ${\textit{Re}}_b = \langle \widehat {{\textit{Re}}_b} \rangle$ . The maximum value of ${\textit{Re}}_b$ observed in each simulation, denoted as ${\textit{Re}}_b^{\star }$ , is listed in table 1. These maximum values cover a reasonably broad range, from $O(0.01)$ to $O(100)$ .

To understand the driving mechanisms of the turbulence, we examine how ${\textit{Re}}_b^{\star }$ depends on the volume of ‘unstable’ regions during the period when the volume-averaged dissipation $\langle \epsilon _k \rangle$ reaches its maximum (i.e. when ${\textit{Re}}_b = {\textit{Re}}_b^{\star }$ ). The local stability with respect to vertical shear is measured using the commonly used gradient Richardson number,

(3.6) \begin{equation} Ri= - \dfrac {1}{{\textit{Fr}}^2} \dfrac {\partial \rho _{\textit{tot}}}{\partial z} \dfrac {1}{S^2}, \end{equation}

while the potential for convective instability is characterised by the density gradient $\partial \rho _{\textit{tot}} / \partial z$ . We define regions as shear unstable where $Ri \lt 0.25$ , and convectively unstable where $\partial \rho _{\textit{tot}} / \partial z \gt 0$ . In figure 9, ${\textit{Re}}_b^{\star }$ is plotted against the volume fraction of the 3-D domain occupied by these unstable regions at the time of maximum dissipation. The results show a strong increase in ${\textit{Re}}_b^{\star }$ with the volume fraction of unstable regions. When the volume fraction drops to zero, ${\textit{Re}}_b^{\star }$ falls below unity, indicating that either shear or convective instability must be present to drive ${\textit{Re}}_b^{\star }$ to values of $O(1)$ or higher. The dependence of ${\textit{Re}}_b^{\star }$ on the unstable volume fraction appears consistent across different values of ${\textit{Re}}$ and ${\textit{Fr}}$ . The similarity between the two criteria shown in figure 9 confirms our earlier observation (in figure 8) that regions of high shear and convective instability are, to a large extent, spatially collocated.

Figure 9. Volume fraction of (a) convectively unstable regions ( $\partial \rho _{\textit{tot}} /\partial z \gt 0$ ) and (b) shear unstable regions ( $Ri \lt 1/4$ ), sampled at the time of peak volume-averaged dissipation rate, as functions of the maximum buoyancy Reynolds number ${\textit{Re}}_b^\star$ . The dotted lines represent a linear fit between $\log {\textit{Re}}_b^\star$ and the volume fractions.

Figure 10. Maximum buoyancy Reynolds number ${\textit{Re}}_b^\star$ versus ${\textit{Re}}\, {\textit{Fr}}^2$ for all cases simulated.

Figure 11. Maximum turbulent Froude number ${\textit{Fr}}_t^\star$ versus (a) ${\textit{Fr}}$ and (b) ${\textit{Re}}_b^\star$ for all cases simulated.

To predict ${\textit{Re}}_b^{\star }$ as a function of ${\textit{Re}}$ and ${\textit{Fr}}$ , we propose the following scaling arguments. Let the characteristic velocity and length scale of the turbulent eddy generated by the breakdown of zigzag instability be denoted $U^*$ and $L^*$ , respectively. We hypothesise that $U^*$ and $L^*$ scale with the characteristic velocity and length scales of the initial vortex array, $U_0$ and $L$ , respectively. The kinetic energy of the eddy, $U^{*2}$ , therefore scales as $U_0^2$ , and the eddy turnover time $L^*/U^*$ scales with $L/U_0$ . The dissipation rate is $\varepsilon ^*\propto U^{*2}/(L^*/U^*) = U^{*3}/L^*$ , thus it scales with $U_0^3/L$ . This hypothesis that $\varepsilon ^* \propto U^{*3}/L^*$ was first proposed by Taylor (Reference Taylor1935) for isotropic turbulence, and has been applied to stratified turbulence by many authors (e.g. as reviewed by Riley & Lindborg Reference Riley, Lindborg, Davidson, Kaneda and Sreenivasan2012). In dimensionless terms, the scaling $\varepsilon ^* \propto U_0^3/L$ leads to the relation

(3.7) \begin{equation} {\textit{Re}}_b^{\star }\propto {\textit{Re}}\, {\textit{Fr}}^2, \end{equation}

which is tested in figure 10. This scaling appears to collapse the data across various ${\textit{Re}}$ and ${\textit{Fr}}$ reasonably well. The empirical relation ${\textit{Re}}_b^{\star } = 0.0088\, {\textit{Re}}\, {\textit{Fr}}^2$ , derived from the DNS data, suggests that ${\textit{Re}}\, {\textit{Fr}}^2$ must reach at least $O(100)$ for ${\textit{Re}}_b^{\star }$ to reach or exceed $O(1)$ such that the flow resulting from the zigzag instability can be considered turbulent.

To assess the influence of stratification on turbulence, a useful metric is the turbulent Froude number, e.g. similar to the definition used by Lewin et al. (Reference Lewin, de Bruyn Kops, Caulfield and Portwood2023) in characterising decaying stratified turbulence from an initially isotropic state:

(3.8) \begin{equation} {\textit{Fr}}_t=\frac {\varepsilon _k\, Fr}{\langle u^2+v^2\rangle }. \end{equation}

As the horizontal kinetic energy $\langle u^2+v^2\rangle$ decays monotonically over time, and $\varepsilon _k$ peaks during the burst of turbulence (see e.g. figure 2 c), $Fr_t$ reaches its maximum value approximately when $\varepsilon _k$ is at its peak. The peak value of $Fr_t$ , denoted $Fr_t^\star$ , is tabulated in table 1 and plotted against ${\textit{Fr}}$ in figure 11(a). An approximately linear relationship $Fr_t\propto Fr$ is observed, consistent with the findings that neither $\varepsilon _k$ nor $\langle u^2+v^2\rangle$ depends significantly on ${\textit{Fr}}$ . Similar to the DNS cases analysed by Lewin et al. (Reference Lewin, de Bruyn Kops, Caulfield and Portwood2023), the scatter plot in figure 11(b) between $Fr_t^\star$ and ${\textit{Re}}_b^\star$ suggests a strong correlation between these parameters. For all cases in which ${\textit{Re}}_b^\star$ exceeds unity, the value of $Fr_t^\star$ is greater than or equal to 0.016 (see table 1). This suggests that the turbulence resulting from the breakdown of TG vortices in the current study approaches or perhaps marginally fall within the ‘strongly stratified turbulence’ regime proposed by Brethouwer et al. (Reference Brethouwer, Billant, Lindborg and Chomaz2007), which requires ${\textit{Re}}_b^\star \gtrsim 1$ and $Fr_t^\star \ll 1$ simultaneously.

4.1. Formation of shear layers

In their seminal work on stratified turbulence generated by 3-D TG flows, Riley & de Bruyn Kops (Reference Riley and de Bruyn Kops2003) noted the spontaneous layering of horizontal momentum as the TG flow transitions into turbulence. Interestingly, the number of shear layers that they observed initially match the number of layers present in the initial condition (see their figure 3). Now, as we examine a base flow with no $z$ dependence, our focus is on whether layers will form, and if so, what determines the characteristic depth of these layers.

In figure 12, the 3-D field of the horizontal velocity $u$ is visualised for the case $({\textit{Re}},Fr)=(1600,1.0)$ previously analysed in § 3.1. The first image, taken at $t=48$ , captures a moment when the flow is transitioning to turbulence (see figure 2). At this point, the signature of the base flow’s columnar structure is still visible. In the next image, at $t=56$ , the flow has fully developed into turbulence, with the dissipation rate reaching its peak. The bottom two images, taken at $t=64$ and $72$ , reveal the clear layering of the horizontal momentum. At $t=64$ , the shear layers are still overlaid with a significant level of turbulence. By $t=72$ , the layers become more distinct, with the sign of $u$ alternating along the vertical axis, and each layer extending across the entire horizontal plane. This layering is further illustrated by the $x$ – $z$ transects of $u$ shown in figure 13(a).

Figure 12. Volume-rendering images of the horizontal velocity $u$ from the DNS with ${\textit{Re}}=1600$ and $Fr=1.0$ . The images are shown at (a) $t=48$ , (b) $t=56$ , (c) $t=64$ and (d) $t=72$ . An animation of the flow evolution can be found as movie 2 in the supplementary material.

Figure 13. Visualisation of $x$ – $z$ transects showing (a) the horizontal velocity $u$ , and (b) the density perturbation $\rho$ , for the case with ${\textit{Re}}=1600, \textit{Fr} = 1.0$ . The transects are taken along $y=\unicode{x03C0} /2$ , as indicated by the red dashed lines in figure 3(a). Snapshots are shown for times $t=48$ , $56$ , $64$ and $72$ .

Such layering is observed across the various $({\textit{Re}},Fr)$ combinations simulated. In figure 14, we examine the vertical structure of $\langle u \rangle _h$ as it evolves over time, where $\langle {\boldsymbol{\cdot}} \rangle _h \equiv (2\unicode{x03C0} )^{-2} \iint ({\boldsymbol{\cdot}} )\, \mathrm{d}x\, \mathrm{d}y$ represents the spatial average over the $x$ – $y$ plane. The vertical structure of $\langle u \rangle _h$ shows layers, typically appearing around or after $t \approx 50$ when the instability has broken down the TG base flow into turbulence. The magnitude of $\langle u \rangle _h$ is of $O(0.1)$ , and the horizontally averaged mean flow appears to persist after $t \gtrsim 50$ . The number of layers over the domain height $2\unicode{x03C0}$ seems to depend strongly on ${\textit{Fr}}$ , with layer depth increasing as ${\textit{Fr}}$ increases, while the dependence on ${\textit{Re}}$ is more subtle.

Figure 14. Horizontally averaged velocity $\langle u \rangle _h$ as a function of time $t$ and vertical position $z$ . Each column corresponds to a different Froude number, from left to right: $Fr = 0.25$ , $0.5$ , $1.0$ and $2.0$ . Each row represents a different Reynolds number: (a) ${\textit{Re}} = 800$ , (b) $Re=1600$ , and (c) $Re=3200$ .

Figure 15. The average vertical wavenumber $\tilde {k}$ representing the horizontally averaged mean flow, $\langle u \rangle _h$ , as a function of (a) ${\textit{Re}}$ , and (b) ${\textit{Fr}}$ .

The dependence of shear layer depth on $({\textit{Re}}, Fr)$ is further illustrated in figure 15. The energy spectrum $E(k)$ is computed for the vertical profile of the horizontally averaged mean flow $\langle u \rangle _h(z)$ , where $k$ represents the vertical wavenumber. As such, $ E(k)$ represents the vertical wavenumber spectrum of the horizontally averaged velocity $\langle u \rangle _h(z)$ , previously shown in figure 14, rather than the full 3-D velocity field. This choice is motivated by our focus on the persistent layered flow structure and contrasts with previous studies (e.g. Brethouwer et al. Reference Brethouwer, Billant, Lindborg and Chomaz2007; Zhou & Diamessis Reference Zhou and Diamessis2019), which compute the spectrum from the 3-D local velocity field to characterise the vertical scales of turbulent motions. An energy-weighted mean wavenumber $\tilde {k}$ is defined as

(4.1) \begin{equation} \tilde {k} = \frac {\sum _k k\,E(k) }{\sum _k E(k)}, \end{equation}

where the summation over $k$ spans from $k=1$ to the Nyquist wavenumber (which depends on the grid resolution tabulated in table 1). The value of $\tilde {k}$ is estimated every two time units and subsequently averaged over the interval $50 \leqslant t \leqslant 100$ , during which the shear layers are observed to persist in the flow across all cases (see figure 14).

As shown in figure 15(a), the average vertical wavenumber $\tilde {k}$ , which describes the inverse shear layer depth, exhibits a weak dependence on ${\textit{Re}}$ at a given ${\textit{Fr}}$ . Generally, $\tilde {k}$ increases slightly as ${\textit{Re}}$ becomes larger for a fixed ${\textit{Fr}}$ . In figure 15(b), $\tilde {k}$ decreases with ${\textit{Fr}}$ for $Fr \leqslant 1$ and becomes approximately constant as ${\textit{Fr}}$ increases from 1 to 2. According to the self-similar scaling of the vertical length scale proposed by Billant & Chomaz (Reference Billant and Chomaz2001), $\tilde {k}$ is expected to scale with $Fr^{-1}$ in the strongly stratified limit $Fr \ll 1$ , which seems to hold approximately for $Fr \leqslant 0.25$ , as shown in figure 15(b). In dimensional terms, this scaling corresponds to a vertical length scale proportional to $ U_0/N$ .

Figure 16. The time-averaged mean vertical wavenumber $\tilde {k}$ representing the horizontally averaged velocity $\langle u \rangle _h$ observed in the DNS, versus the fastest-growing wavenumber $k_{\textit{fgm}}$ , according to the LSA.

To examine the role of linear instability in determining the length scale of the mean shear layers, we identify the wavenumber corresponding to the fastest-growing mode of the primary zigzag instability, $k_{\textit{fgm}}$ , from the LSA described in G24. While G24 focused on ${\textit{Re}}=1600$ , we extend this analysis to two additional Reynolds numbers. The wavenumber predicted by the LSA, $k_{\textit{fgm}}$ , is compared to the observed wavenumber $\tilde {k}$ in figure 16. Since the LSA is performed for integer values of $k$ , the value of $k_{\textit{fgm}}$ is quantised, i.e. integers only, as shown in figure 16. Despite this, the wavenumbers associated with the mean shear layers, $\tilde {k}$ , show reasonable agreement with the LSA-predicted $k_{\textit{fgm}}$ . This suggests that linear instability plays a role in shaping the structure of the momentum layers that result from the turbulent breakdown of columnar TG vortices.

4.2. Formation of buoyancy layers

The phenomenon of buoyancy layering appears more subtle than the layering of horizontal momentum (§ 4.1). As shown in figure 13(b), while weak signs of layering in $\rho$ are visible at $t=56$ , the magnitude of the density perturbation decreases significantly over time. Unlike the $u$ velocity field, which shows bands of constant $u$ extending across the entire horizontal plane by $t=72$ , the layering of $\rho$ is more transient and weaker in magnitude. The vertical structure of the horizontally averaged density, $\langle \rho \rangle _h(z)$ , is displayed in figure 17. Although some persistent layering of the horizontal averaged density is present, the signal in $\langle \rho \rangle _h$ is considerably weaker in magnitude than in $\langle u \rangle _h$ (shown in figure 14). This difference is evident, as $\langle \rho \rangle _h$ only becomes visible for the low- ${\textit{Fr}}$ cases (e.g. $Fr=0.25$ ) when the colour bar range is adjusted to $\pm 0.1$ in figure 17, as compared to $\pm 0.3$ for $\langle u \rangle _h$ as in figure 14. The layers of $\langle \rho \rangle _h$ are finer than those of $\langle u \rangle _h$ for the same set of $({\textit{Re}},Fr)$ when the layers initially form at $35\lesssim t \lesssim 40$ . The number of layers decreases over time, presumably due to the merging of layers.

Figure 17. Same layout as in figure 14, but showing the horizontally averaged density perturbation $\langle \rho \rangle _h$ .

To better understand the dynamics of these persistent buoyancy layers, we apply the framework proposed by Winters & D’Asaro (Reference Winters and D’Asaro1996) and Nakamura (Reference Nakamura1996), as employed by Zhou et al. (Reference Zhou, Taylor, Caulfield and Linden2017) in their analysis of density layers in stratified shear flows. In this approach, we introduce a hypothetical ‘sorted’ vertical coordinate $z_*$ , which is related to the ‘background’ buoyancy profile obtained by sorting all fluid parcels adiabatically to minimise the system’s potential energy. For buoyancy $b\equiv -\rho _{\textit{tot}}$ , the sorted profile $b(z_*,t)$ can be obtained by evaluating the probability density function of $b$ (Tseng & Ferziger Reference Tseng and Ferziger2001), which is equivalent to sorting the density parcels over the entire 3-D volume. Within this coordinate system, the buoyancy field $b(x,y,z,t)$ is averaged over each isopycnal surface, thereby reducing the full 3-D field to a one-dimensional, ‘sorted’ profile $b(z_*,t)$ . The $b(z_*,t)$ profile characterises any structural changes to the background stratification. A useful measure for diagnosing how the background stratification has been altered by the flow is the vertical gradient of the background buoyancy profile:

(4.2) \begin{equation} N_*^2(z_*,t)\equiv \frac {\partial b(z_*,t)}{\partial z_*}. \end{equation}

In the undisturbed state where $\partial \rho _{\textit{tot}}/\partial z = -1$ , we have $N_*^2 = 1$ . Any ‘sharpening’ of the background buoyancy profile due to the flow would result in $N_*^2 \gt 1$ .

Zhou et al. (Reference Zhou, Taylor, Caulfield and Linden2017) derived the following equation governing the evolution of $N^2_*$ ,

(4.3) \begin{equation} \frac {\partial N^2_*}{\partial t}= \underbrace {\frac {\partial ^{2}\kappa _{e}}{\partial z_{*}^{2}}N^2_*}_{\text{Source}\ \mathbb{S}(t)}+ \underbrace {2\frac {\partial \kappa _{e}}{\partial z_{*}}\frac {\partial N^2_*}{\partial z_{*}}}_{\text{Advection}\ \mathbb{A}(t)}+ \underbrace {\kappa _{e}\frac {\partial ^{2}N^2_*}{\partial z_{*}^{2}}}_{\text{Diffusion}\ \mathbb{D}(t)}, \end{equation}

and highlighted the relevance of the first term on the right-hand side which acts as a ‘source’ term for $N^2_*$ . For sharpening of the buoyancy profile to occur, the curvature of $\kappa _e$ , i.e. $\partial ^2\kappa _e/\partial z_*^2$ , must be positive, ensuring that the source term contributes positively to the $N_*^2$ budget and potentially overcomes the ‘diffusion’ of $N^2_*$ , i.e. the third term on the right-hand side of (4.3). The effective diffusivity of buoyancy $\kappa _e$ is defined as

(4.4) \begin{equation} \kappa _e=\kappa \left (\frac {\partial z_{*}}{\partial b}\right )^{2}\langle |\boldsymbol{\nabla} b|^{2}\rangle _{z_{*}}, \end{equation}

where $\kappa =1/({\textit{Re}}\, Pr)$ is the molecular diffusivity, and $\langle {\boldsymbol{\cdot}} \rangle _{z_*}$ denotes an average over the isopycnal surface that corresponding to $b=b(z_*)$ . Winters & D’Asaro (Reference Winters and D’Asaro1996) provided a geometrical interpretation of (4.4), rewriting $\kappa _e$ as

(4.5) \begin{equation} \kappa _e=\kappa \left (\frac {A_s}{A}\right )^{2}. \end{equation}

Here, $A_s$ is the isopycnal surface area corresponding to $b(z_*)$ , and $A$ is the domain’s horizontal area. In a turbulent flow, $\kappa _e$ is typically greater than $\kappa$ , as the distortion of isopycnal surfaces acts to increase $\kappa _e$ according to (4.5).

Figure 18. (a) The sorted buoyancy gradient $N^2_*$ , and (b) the scaled effective diffusivity $\kappa _e/\kappa$ as they vary in the tracer-based, sorted coordinate $z_*$ and time for the case $({\textit{Re}},Fr)=(1600,1.0)$ .

Figure 19. Plots of (a) $N^2_*$ and (b) $\kappa _e\, Pr$ as functions of $z_*$ at various times for the case $({\textit{Re}},Fr)=(1600,1.0)$ .

The evolution of $N^2_*$ and $\kappa _e$ for the case $({\textit{Re}}, Fr) = (1600, 1.0)$ is shown in figures 18 and 19. In figure 18(a), red bands where $N^2_*\gt 1$ appear periodically along the $z_*$ -axis, indicating the formation of buoyancy ‘interfaces’ at these locations for $t \gtrsim 45$ , coinciding with the onset of turbulence in the vortex array (§ 3.1). Between these interfaces, blue bands emerge, which correspond to considerably persistent ‘layers’ where the buoyancy gradient is reduced from the undisturbed value, i.e. $N^2_* \lt 1$ . These modifications in the $N^2_*$ structure due to the turbulent breakdown of the vortex array are also evident in figure 19(a), where spikes of $N^2_*\gt 1$ appear at the interfaces, with layers in between. The magnitude of $N^2_*$ reaches its peak at $t \approx 60$ , as the turbulence begins to decay substantially (see figure 2). As the turbulence continues to dissipate, the contrast between the layers and interfaces gradually fades, presumably due to diffusion, i.e. the third term on the right-hand side of (4.3).

Further insights into the layering dynamics are provided in figures 18(b) and 19(b), which display the effective diffusivity ( $\kappa _e$ ) scaled by its molecular value ( $\kappa$ ), i.e. $\kappa _e/\kappa$ . Figure 18(b) reveals patches of elevated $\kappa _e$ at specific heights (in terms of $z_*$ ) during the time interval $50 \lesssim t \lesssim 65$ . The occurrence of these high $\kappa _e$ values at particular heights could be linked to the instability mechanism discussed in § 3.1. As shown in figure 5(b) and sketched in figure 6, the isopycnals are distorted by the $w$ velocity only above and below certain heights (i.e. the neutral planes where $w=0$ ), forming a pair of circulation cells (or vortex tubes under stretching, as shown in figure 7) that eventually overturn the isopycnals. This process leads to elevated $\kappa _e$ values between the neutral planes and near-molecular values at the neutral plane itself, which could explain the pattern of $\kappa _e$ seen in figure 18(b). The overturning motions that trigger turbulence reduce the vertical density gradient at specific heights, corresponding to the local minima in the $N_*^2(z_*)$ profile shown in figure 18(a). It is at these heights of reduced $N_*^2$ that large values of $\kappa _e$ emerge in figure 18(b), providing the link between zigzag-instability-induced overturns and the observed structure of $\kappa _e$ .

A visual inspection of figures 19(a,b), showing $N^2_*$ and $\kappa _e/\kappa$ , respectively, suggests that the source term in (4.3) plays a key role in the layering dynamics of buoyancy. The local minima of $\kappa _e$ in figure 19(b), where the curvature is positive (i.e. $\partial ^2\kappa _e/\partial z_*^2 \gt 0$ ), align with the locations of spikes in $N^2_*$ (representing interfaces) in figure 19(a). Between these interfaces, the overturning of isopycnals generates intense turbulence and enlarges the isopycnal surfaces, leading to large values of $\kappa _e$ . Conversely, at the interfaces, where $N^2_*$ is large, turbulence is suppressed, resulting in low values of $\kappa _e$ . This contrast creates the positive curvature at the interface, which locally sharpens the buoyancy gradient through the source term in (4.3). However, since the turbulence is not sustained, $\kappa _e$ decreases as the turbulence decays, reducing the curvature and diminishing the interface’s ability to sharpen itself. This could explain, at least in part, why buoyancy layering is more transient and weaker than momentum layering (§ 4.1).