3.1 Effects of initial perturbation amplitude

In all simulations shown here there is a strong transient adjustment associated with the initial perturbation, with a decrease in kinetic energy, growth of buoyancy fluctuations and enhanced dissipation. This adjustment takes place over the first few non-dimensional time units, on scales comparable to the advective timescale ( $t^{\ast }\sim h^{\ast }/U_{0}^{\ast }$ ) and buoyancy timescale ( $t^{\ast }\sim h^{\ast }/\sqrt{B_{0}}^{\ast }$ ), both of which are $O(1)$ quantities when non-dimensionalized.

Profiles of the stratification and the gradient Richardson number of the mean flow (horizontally averaged) at $t=1$ (i.e. after the flow has started to adjust to the initial perturbation) are shown in figure 1. It is clear that the initial turbulence has led to a more complicated mean flow, particularly for the intermediate and strongly turbulent cases ( ${\mathcal{A}}=0.05$ and ${\mathcal{A}}=0.1$ ). However, it should be emphasized that in spite of the initial turbulence, the large-scale stratification still varies over the same length scale, and the mean Richardson number still exhibits a minimum on the order of $Ri_{b}$ . As such, the initial perturbation does not appear to significantly change the bulk stability of the flow.

Figure 1.

Vertical profiles of (a) mean buoyancy gradient, $\text{d}B_{mean}/\text{d}z$ , and (b) gradient Richardson number of the mean flow, $Ri_{mean}$ at $t=1$ . $Ri=1/4$ is denoted by a thin dashed line in panel (b). Note that the line corresponding to ${\mathcal{A}}=0.0025$ lies almost directly beneath the line for ${\mathcal{A}}=0.01$ .

In addition to the characteristic amplitude of the velocity perturbation, we can further quantify the initial turbulence in terms of a buoyancy Reynolds number, $Re_{b}=\unicode[STIX]{x1D700}^{\ast }/\unicode[STIX]{x1D708}^{\ast }N^{\ast 2}$ , which describes the turbulence intensity. Characteristic values of $Re_{b}$ at $t=1$ corresponding to each individual simulation are listed in table 1. Values vary from $Re_{b}\lesssim 10$ for the laminar and weak initial turbulence cases ( ${\mathcal{A}}=0.0025,0.01$ ), at which viscous effects are expected to be important (Shih et al. Reference Shih, Koseff, Ivey and Ferziger2005) to $Re_{b}\sim O(10^{2})$ for the strong initial turbulence cases, corresponding to highly energetic turbulence. For context, oceanic observations suggest a broad range of $Re_{b}\sim O(10)-O(10^{4})$ , with $Re_{b}\gtrsim 10^{2}$ during vigorous mixing (Smyth, Moum & Caldwell Reference Smyth, Moum and Caldwell2001).

We now present a qualitative overview of the impact of perturbation amplitude on the evolution of the shear instability. Figure 2 shows vertical slices of buoyancy through $y=0$ at different times for cases with weak, intermediate and strong initial turbulence; the laminar case follows a very similar evolution to the case with weak initial turbulence and is therefore not shown here. For laminar and weakly turbulent initial conditions, the flow follows the usual lifecycle of a classical Kelvin–Helmholtz instability (Peltier & Caulfield Reference Peltier and Caulfield2003; Smyth & Moum Reference Smyth and Moum2012): two-dimensional perturbations grow and form finite-amplitude billows (figure 2 b,c), which are susceptible to three-dimensional secondary instabilities leading to a transition to turbulence (figure 2 d). As the flow considered here is unforced, it eventually relaminarizes after a period of turbulent mixing (figure 2 e).

As the initial perturbation amplitude is increased, the details of the flow evolution are modified. For the case of intermediate-amplitude initial turbulence ( ${\mathcal{A}}=0.05$ , figure 2 f–j), a growing perturbation still leads to billow formation, transition to turbulence, and subsequent relaminarization. However, the growing perturbation is no longer two-dimensional, and the billow becomes turbulent much earlier in its development, without rolling up to the extent seen in the laminar and weakly turbulent cases (compare figure 2 c,h). As the strength of the initial perturbation is further increased, the characteristic billows associated with Kelvin–Helmholtz instability are no longer observed (figure 2 k–o). Instead, the initial turbulence fills the shear layer and subsequently decays.

Figure 2.

Vertical slices of buoyancy through $y=0$ at various times for (a–e) weakly turbulent initial conditions ( ${\mathcal{A}}=0.01$ ), (f–j) intermediate initial conditions ( ${\mathcal{A}}=0.05$ ), and (k–o) strongly turbulent initial conditions ( ${\mathcal{A}}=0.1$ ). Non-dimensional times corresponding to each snapshot are shown in the bottom left corner of each panel.

Figure 3(a) compares the initial buoyancy given in (2.2) with the horizontally-averaged buoyancy $B_{mean}$ after the flow has relaminarized. The vertical extent of the relaminarized layer depends strongly on whether finite-amplitude billows (figure 2 c,h) form and move denser fluid upwards (and vice versa). As a result, the final state of the flow with strongly turbulent initial conditions (without billow formation) has a much smaller vertical extent than any of the other cases shown. A similar trend is observed for the mean streamwise velocity (not shown here).

In addition to $B_{mean}$ , the gradient Richardson number of the final mean flow, $Ri_{mean}$ , is shown in figure 3(b). The minimum Richardson number at the end of the mixing event is well above the critical Miles-Howard value of $1/4$ for the laminar, weak and intermediate initial turbulence cases, consistent with experimental results suggesting that the Richardson number after a Kelvin–Helmholtz mixing event should lie in the range $0.32\pm 0.06$ (Thorpe Reference Thorpe1973). In contrast, strong initial turbulence results in much lower Richardson numbers after relaminarization. The final state is only marginally stable at the centre of the shear layer (green curve).

Figure 3.

Initial (thin black lines) and final (thick coloured lines) vertical profiles of (a) mean buoyancy, $B_{mean}$ , and (b) gradient Richardson number of the mean flow, $Ri_{mean}$ . $Ri=1/4$ is denoted by a thin dashed line in panel (b).

The effect of varying the initial perturbation amplitude can be studied more quantitatively by considering the evolution of the total kinetic and potential energy of the flow. An equation governing the evolution of the kinetic energy, ${\mathcal{K}}\equiv \overline{u_{j}u_{j}/2}$ (where the overline denotes a volume average) is formed by taking the inner product of (2.3) with $\boldsymbol{u}$ and then averaging over the flow domain, i.e.

(3.1) $$\begin{eqnarray}\frac{\text{d}{\mathcal{K}}}{\text{d}t}={\mathcal{B}}-\unicode[STIX]{x1D700},\end{eqnarray}$$

where ${\mathcal{B}}\equiv Ri_{b}\overline{bw}$ is the vertical buoyancy flux and $\unicode[STIX]{x1D700}\equiv \overline{\unicode[STIX]{x2202}_{j}u_{i}\unicode[STIX]{x2202}_{j}u_{i}}/Re$ is the (positive-definite) kinetic energy dissipation. Similarly, an equation governing the total potential energy ${\mathcal{P}}=-Ri_{b}\overline{bz}$ is constructed by multiplying (2.4) by $Ri_{b}z$ and volume-averaging, i.e.

(3.2) $$\begin{eqnarray}\frac{\text{d}{\mathcal{P}}}{\text{d}t}=-{\mathcal{B}}+{\mathcal{D}}_{p},\end{eqnarray}$$

where ${\mathcal{D}}_{p}\equiv Ri_{b}(b_{top}-b_{bottom})/(RePrL_{z})$ denotes the change in potential energy due to diffusion of the mean laminar flow (i.e. the diffusive change in potential energy that would occur even in the absence of any fluid motion). In contrast, the buoyancy flux ${\mathcal{B}}$ describes the exchange of energy between kinetic and potential forms due to fluid motion, including both reversible stirring and irreversible mixing (Winters et al. Reference Winters, Lombard, Riley and D’Asaro1995; Peltier & Caulfield Reference Peltier and Caulfield2003; Davies Wykes, Hughes & Dalziel Reference Davies Wykes, Hughes and Dalziel2015).

To isolate the irreversible mixing from the total buoyancy flux, we decompose the total potential energy into available and background components, ${\mathcal{P}}={\mathcal{P}}_{a}+{\mathcal{P}}_{b}$ . The background potential energy refers to the minimum potential energy state of the system which would be achieved by adiabatically rearranging the buoyancy field into a gravitationally stable configuration,

(3.3) $$\begin{eqnarray}{\mathcal{P}}_{b}=-Ri_{b}\overline{b^{\star }z}.\end{eqnarray}$$

The sorted profile $b^{\star }(z)$ may be calculated from the full three-dimensional buoyancy field either by sorting or by computing its probability density function, as we do here (Winters et al. Reference Winters, Lombard, Riley and D’Asaro1995; Tseng & Ferziger Reference Tseng and Ferziger2001; Davies Wykes et al. Reference Davies Wykes, Hughes and Dalziel2015). The available potential energy ${\mathcal{P}}_{a}$ is then the energy that may be used to mix the flow and arises due to horizontal buoyancy gradients or from statically-unstable regions of the flow (Lorenz Reference Lorenz1955). The governing equation for potential energy (3.2) can then be split into equations governing the evolution of these two components,

(3.4a,b ) $$\begin{eqnarray}\frac{\text{d}{\mathcal{P}}_{a}}{\text{d}t}=-{\mathcal{B}}-{\mathcal{M}},\quad \frac{\text{d}{\mathcal{P}}_{b}}{\text{d}t}={\mathcal{M}}+{\mathcal{D}}_{p}.\end{eqnarray}$$

It should be noted that Winters et al. (Reference Winters, Lombard, Riley and D’Asaro1995) refers to the sum of the terms ${\mathcal{M}}+{\mathcal{D}}_{p}$ as $\unicode[STIX]{x1D6F7}_{d}$ , describing the total change in background potential energy. However, here we are interested primarily in the irreversible mixing due to the evolution of the shear instability and its transition to turbulence. As such, we follow the convention of Peltier & Caulfield (Reference Peltier and Caulfield2003) and distinguish between changes in ${\mathcal{P}}_{b}$ due to fluid motion (described by ${\mathcal{M}}$ ) and changes due to laminar diffusion of the mean flow (described by ${\mathcal{D}}_{p}$ ).

The evolution of the kinetic and potential energies are shown in figure 4. In the laminar and weakly turbulent cases ( ${\mathcal{A}}=0.0025$ and $0.01$ ; blue and red curves) there is little change in kinetic or potential energy until $t\sim 100$ when the linearly-growing perturbation reaches finite amplitude. The change in ${\mathcal{P}}$ is initially associated with the growth of available potential energy arising from the formation of the billow structure. Then, as the secondary instabilities grow and trigger turbulence, the background potential energy begins to increase irreversibly while the available potential energy decays. Finally, as the flow relaminarizes, the energy changes approach near-constant values.

For the case with strongest initial turbulence ( ${\mathcal{A}}=0.1$ ; green curves), the energy evolves in a very different manner, as might be expected based on the buoyancy evolution (figure 2 k–o). The kinetic and potential energies change quickly as the initial turbulence decays and mixes the mean flow (also seen to some extent for the case with intermediate initial turbulence). Changes in available potential energy are small, and changes in total potential energy largely correspond to changes in ${\mathcal{P}}_{b}$ from irreversible mixing. After the initial vigorous mixing, the small remaining changes in energy are largely due to the effects of viscosity and diffusion on the relaminarized mean flow, with a plateau in ${\mathcal{P}}_{b}-{\mathcal{D}}_{p}t$ (figure 4 c), and no further increases associated with the growth of the Kelvin–Helmholtz instability.

The case with intermediate-amplitude initial turbulence ( ${\mathcal{A}}=0.05$ ; yellow curves) shares features of both the weak and strong initial turbulence cases. There is an initial change in kinetic and potential energy owing to the turbulent initial conditions. Then, as the billow begins to form (figure 2 h), available potential energy rises, as for the classical Kelvin–Helmholtz instability. However, the flow becomes turbulent before the billow forms completely, and correspondingly the growth of ${\mathcal{P}}_{a}$ is arrested as ${\mathcal{P}}_{b}$ begins to rise due to turbulent mixing. Finally, as in the other three cases, the energies plateau as the turbulence decays and the flow relaminarizes. Overall potential energy changes depend non-monotonically on the perturbation amplitude: larger changes in total and background potential energy are observed for ${\mathcal{A}}=0.01$ when compared with ${\mathcal{A}}=0.0025$ and ${\mathcal{A}}=0.05$ .

Figure 4.

Time evolution of changes in (a) kinetic energy ${\mathcal{K}}$ , (b) total potential energy ${\mathcal{P}}$ (solid) and available potential energy ${\mathcal{P}}_{a}$ (dashed), and (c) background potential energy associated with fluid motion, ${\mathcal{P}}_{b}-{\mathcal{D}}_{p}t$ .

Using the above description of the flow energetics, the mixing efficiency can be defined instantaneously as

(3.5) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{i}=\frac{{\mathcal{M}}}{{\mathcal{M}}+\unicode[STIX]{x1D700}}.\end{eqnarray}$$

In the laminar and weak turbulence cases, the evolution of ${\mathcal{M}}$ , $\unicode[STIX]{x1D700}$ and $\unicode[STIX]{x1D702}_{i}$ (blue and red curves, figure 5 a,b) is similar to that found in previous studies of mixing by Kelvin–Helmholtz instability (Smyth et al. Reference Smyth, Moum and Caldwell2001; Peltier & Caulfield Reference Peltier and Caulfield2003; Mashayek et al. Reference Mashayek, Caulfield and Peltier2013). There is an initial phase during which the irreversible mixing rate is larger than the dissipation, with very high mixing efficiency, $\unicode[STIX]{x1D702}_{i}\sim 0.6$ . This occurs shortly after the peak in available potential energy (figure 4 b), and is associated with the rollup of the shear instability into the classical billow structure. Then, as the flow becomes turbulent, the mixing efficiency drops to a value $\unicode[STIX]{x1D702}_{i}\sim 0.3$ as $\unicode[STIX]{x1D700}$ peaks. As the flow relaminarizes, ${\mathcal{M}}\rightarrow 0$ while $\unicode[STIX]{x1D700}$ remains non-zero owing to the continued viscous dissipation of the mean flow, and so $\unicode[STIX]{x1D702}_{i}$ decreases further.

Figure 5.

(a) Time variation of mixing rate, ${\mathcal{M}}$ (solid), and dissipation, $\unicode[STIX]{x1D700}$ (dashed). (b) Time variation of instantaneous mixing efficiency, $\unicode[STIX]{x1D702}_{i}$ . (c) Time variation of turbulent diffusivity, $K_{\unicode[STIX]{x1D70C}}$ . (d) Cumulative mixing efficiencies $\unicode[STIX]{x1D702}_{c}$ calculated over an entire mixing event for three definitions of the initial time $t=t_{i}$ . The error bars correspond to mixing efficiencies calculated with final times $t_{f}={\mathcal{D}}_{p}/2$ and $t_{f}=2{\mathcal{D}}_{p}$ .

In cases with increased initial perturbation amplitude, the mixing efficiency evolves differently. As figure 2 shows, the billow structure breaks down into turbulence much earlier in the case of intermediate initial turbulence, without rolling up to the same extent as in the classical Kelvin–Helmholtz case. Correspondingly, the peak in $\unicode[STIX]{x1D702}_{i}$ associated with laminar rollup of the billow is absent, and the instantaneous irreversible mixing rate ${\mathcal{M}}$ is always less than the dissipation $\unicode[STIX]{x1D700}$ (yellow curves, figure 5 a,b). The maximum instantaneous mixing efficiency in this case is only approximately $0.35$ , as in the later turbulent mixing stage of the lower initial amplitude cases. As the initial perturbation amplitude is further increased and the billow structure no longer arises (figure 2 k–o), the mixing efficiency only achieves values of approximately $0.2$ during the initial turbulent mixing phase before decaying at later times.

As with the overall changes in potential energy shown in figure 4(b,c), the peak value of the irreversible mixing rate and the dissipation depend non-monotonically on the initial perturbation amplitude ${\mathcal{A}}$ , with maximum values occurring for ${\mathcal{A}}=0.01$ . Indeed, the maximum values of both ${\mathcal{M}}$ and $\unicode[STIX]{x1D700}$ are very similar for the laminar ( ${\mathcal{A}}=0.0025$ ) and intermediate initial turbulence ( ${\mathcal{A}}=0.05$ ) cases, despite the marked differences in both instantaneous and cumulative mixing efficiencies between these two cases. This further emphasizes the role of the highly-efficient pre-turbulent phase associated with billow rollup in determining the overall efficiency of a given mixing event.

As mentioned above, the viscous dissipation can be related to a turbulent density flux via a turbulent diffusivity, $K_{\unicode[STIX]{x1D70C}}$ , expressed in non-dimensional form as (Salehipour & Peltier Reference Salehipour and Peltier2015)

(3.6) $$\begin{eqnarray}K_{\unicode[STIX]{x1D70C}}=\unicode[STIX]{x1D6E4}PrRe_{b}=\frac{\unicode[STIX]{x1D702}}{1-\unicode[STIX]{x1D702}}PrRe_{b}.\end{eqnarray}$$

All cases show non-zero values of $K_{\unicode[STIX]{x1D70C}}$ at similar times to the changes in ${\mathcal{P}}_{b}$ (figure 4 c). For strongly turbulent initial conditions ( ${\mathcal{A}}=0.1$ , green line in figure 5 c), there is a small increase in $K_{\unicode[STIX]{x1D70C}}$ at early times and then subsequent decay. For the other three cases, much larger peak values of $K_{\unicode[STIX]{x1D70C}}$ are observed. In all three cases there is a peak in $K_{\unicode[STIX]{x1D70C}}$ associated with peak values of $\unicode[STIX]{x1D700}$ (and peak values of $Re_{b}$ , not shown), with the highest values associated with ${\mathcal{A}}=0.01$ . However, there is an additional peak for the two cases in which billow formation was not impeded ( ${\mathcal{A}}=0.0025$ and $0.01$ ). This peak is associated with the onset of turbulence after the rolling up of the billows and the associated peak in mixing efficiency. While the viscous dissipation is still relatively low at this point in the flow evolution, the highly efficient motions lead to a correspondingly high turbulent buoyancy flux. As a result, $K_{\unicode[STIX]{x1D70C}}$ remains elevated for a much longer period of time in these two cases relative to the case with intermediate initial turbulence ( ${\mathcal{A}}=0.05$ ), leading to larger changes overall in the background potential energy.

It is also of interest to describe an individual mixing event in terms of an overall mixing efficiency. A cumulative analogue to (3.5) is defined by integrating ${\mathcal{M}}$ and $\unicode[STIX]{x1D700}$ over the duration of a mixing event (from some initial time $t_{i}$ to some final time $t_{f}$ ) as

(3.7) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{c}=\frac{\displaystyle \int _{t_{i}}^{t_{f}}{\mathcal{M}}\,\text{d}t}{\displaystyle \int _{t_{i}}^{t_{f}}{\mathcal{M}}\,\text{d}t+\int _{t_{i}}^{t_{f}}\unicode[STIX]{x1D700}\,\text{d}t}.\end{eqnarray}$$

As previously discussed, and as figures 4 and 5(a) show, there is a strong initial transient adjustment for the cases with intermediate and strong initial turbulence. We find that the computed values of $\unicode[STIX]{x1D702}_{c}$ are sensitive to initial integration times $t_{i}$ within this adjustment period until approximately $t_{i}\sim 1$ ; the computed values of $\unicode[STIX]{x1D702}_{c}$ do not change significantly for $1\lesssim t_{i}\lesssim 10$ . We define $t_{f}$ as the time when ${\mathcal{M}}={\mathcal{D}}_{p}$ after the mixing event, i.e. the time at which changes in ${\mathcal{P}}_{b}$ due to turbulent mixing become comparable to changes due to laminar diffusion in the absence of fluid motion. Note, however, that the computed mixing efficiencies in figure 5(d) are not strongly sensitive to the choice of $t_{f}$ , as long as $t_{f}$ is after the shear layer begins to relaminarize.

The cases with laminar and weakly turbulent initial conditions (i.e. the cases that follow a typical Kelvin–Helmholtz lifecycle) yield similar cumulative mixing efficiencies, $\unicode[STIX]{x1D702}_{c}\approx 0.27$ , consistent with previous results for simulations of the Kelvin–Helmholtz instability (Peltier & Caulfield Reference Peltier and Caulfield2003; Mashayek et al. Reference Mashayek, Caulfield and Peltier2013). As the initial perturbation amplitude is increased and the large peak in $\unicode[STIX]{x1D702}_{i}$ associated with laminar rollup is lost, $\unicode[STIX]{x1D702}_{c}$ decreases to values of approximately $0.18$ for $t_{i}=0$ in the intermediate initial turbulence case and $0.025$ in the strong initial turbulence case. This decrease in $\unicode[STIX]{x1D702}_{c}$ is consistent with the evolution of ${\mathcal{M}}$ , $\unicode[STIX]{x1D700}$ and $\unicode[STIX]{x1D702}_{i}$ for these amplitudes: irreversible mixing is essentially restricted to times when the flow is turbulent, without the highly efficient pre-turbulent phase associated with billow structure and its large available potential energy. If the initial transient adjustment is excluded ( $t_{i}\geqslant 1$ ), $\unicode[STIX]{x1D702}_{c}$ increases to values of approximately $0.22$ and $0.13$ for ${\mathcal{A}}=0.05$ and ${\mathcal{A}}=0.1$ , respectively. This increase in $\unicode[STIX]{x1D702}_{c}$ is due to the fact that, at very early times, the buoyancy field has not yet adjusted to the initial perturbation, and so the flow is dissipating kinetic energy without mixing buoyancy.

3.2 $Re$ and $Ri_{b}$ effects

Thus far, we have focused on a single choice of the mean flow parameters $Re$ and $Ri_{b}$ , showing that turbulent initial conditions can impact the evolution of a stratified shear layer. To examine the effect of the mean flow parameters, we now consider simulations at a variety of initial $Re$ and $Ri_{b}$ for the same four values of ${\mathcal{A}}$ as considered in the previous section. While the initial perturbation varies slightly depending on the domain size and grid resolution (and, therefore, the corresponding values of $k_{h}$ in (2.8)), the kinetic energy of the initial perturbation and the turbulence intensity at early times are similar for the same value of ${\mathcal{A}}$ across the range of simulations presented in table 1.

Billow growth due to shear instability is associated with an increase in the available potential energy, ${\mathcal{P}}_{a}$ . As such, we look for local maxima in the growth rate of available potential energy, $\text{d}{\mathcal{P}}_{a}/\text{d}t$ , as a means of identifying billows in our suite of simulations. In figure 6 we show the billow structure corresponding to a range of different values of $Re$ and $Ri_{b}$ and all four values of ${\mathcal{A}}$ , identified using this metric. Note that the specific times corresponding to each snapshot in figure 6 differ depending on both the flow parameters and the initial perturbations. Finite-amplitude billows typically form at later non-dimensional times for higher values of $Ri_{b}$ (owing to the slower growth rate of the underlying shear instability) and earlier for higher ${\mathcal{A}}$ (owing to the larger amount of energy initially present in the unstable mode). However, despite the different times highlighted in figure 6, for the cases with clear billow structures the instability follows a similar evolution to that shown in figure 2: the billow becomes unstable to secondary instabilities, the shear becomes fully turbulent, and after some time the turbulence decays and the flow relaminarizes.

Figure 6.

Vertical slices of buoyancy through $y=0$ at various times corresponding to local maxima in $\text{d}{\mathcal{P}}_{a}/\text{d}t$ . (a–e) Laminar initial conditions, ${\mathcal{A}}=0.0025$ . (f–j) Weakly turbulent initial conditions, ${\mathcal{A}}=0.01$ . (k–o) Intermediate turbulence initial conditions, ${\mathcal{A}}=0.05$ . (p–t) Strongly turbulent initial conditions, ${\mathcal{A}}=0.1$ . Each column corresponds to a different set of simulations with the specified values of $Re$ and $Ri_{b}$ . Non-dimensional times corresponding to each snapshot are shown in the bottom left corner of each panel. Panels with thick black borders correspond to the value of ${\mathcal{A}}$ which gives the maximum change in ${\mathcal{P}}_{b}$ for a given $(Re,Ri_{b})$ .

Several trends become apparent upon consideration of figure 6. As the initial turbulence amplitude is increased, the development of the shear instability is accordingly modified: the billow structure becomes susceptible to secondary instabilities before fully rolling up, and is prevented from forming for strongly turbulent initial conditions (as shown in the bottom two rows of figure 6). This effect is consistent across the range of $Re$ and $Ri_{b}$ investigated here.

However, figure 6 also reveals additional effects related to the initial mean flow parameters. The initial Reynolds number of the shear layer appears to affect the sensitivity of the billow structure to modification by pre-existing turbulence. For example, with lower $Re$ and weak initial turbulence (figure 6 f,g,i), the resulting billow appears similar to the corresponding structure with laminar initial conditions (figure 6 a,b,d). However, as $Re$ is increased, the billow structure is modified; for example, in comparing figures 6(c) and 6(h), it is clear that the weak initial turbulence has led to asymmetry in the fully developed billow before transition to turbulence. Similarly, $Re$ effects are apparent for the partially developed billows in the case of intermediate initial turbulence. As figure 6(k,l,m) show, when the Reynolds number is increased for the same initial $Ri_{b}$ , the height of the overturn before breaking down into turbulence is decreased (figure 6 k,l,m).

Figure 6 also suggests a strong dependence on the initial bulk Richardson number, $Ri_{b}$ . The underlying shear instability is less sensitive to the initial turbulence at lower $Ri_{b}$ , with the billow able to grow to its full vertical height even in the presence of intermediate initial turbulence (figure 6 n). At the other extreme, for higher $Ri_{b}$ , intermediate initial turbulence is sufficiently strong to almost prevent billow formation altogether (figure 6 o).

Figure 7(a) shows the overall change in potential energy associated with turbulent irreversible mixing (i.e.  $\unicode[STIX]{x0394}{\mathcal{P}}_{b}-{\mathcal{D}}_{p}t$ ) as a function of the buoyancy Reynolds number at $t=1$ . As observed in figure 4(b), there is a non-monotonic variation of total change in ${\mathcal{P}}_{b}$ with $Re_{b}$ . The specific amplitude at which the change in background potential energy peaks depends on $Ri_{b}$ : while the cases with $Ri_{b}=0.16$ (all $Re$ ) and the cases with $Ri_{b}=0.20$ and $Re=4000$ exhibit maximum mixing for weak initial turbulence ( ${\mathcal{A}}=0.01$ ), the cases with $Ri_{b}=0.12$ instead show a clear peak for intermediate initial conditions ( ${\mathcal{A}}=0.05$ ). This can be understood by considering again the billow structures in figure 6 – in particular, those cases indicated by thick black borders which correspond to the large changes in ${\mathcal{P}}_{b}$ . The cases with the most mixing correspond to cases in which there is some initial turbulence present, but not so much as to prevent the billow from growing to its maximum vertical extent. That is, the overturn associated with the underlying instability is able to reach its maximum amplitude, but then transitions to turbulence before completely rolling up. This is especially noticeable in figure 6(n). As initial turbulence strength is further increased and the billow structure is not able to form completely (e.g. figure 6 k,l,m), the total change in background potential energy correspondingly begins to decrease. Finally, very small changes in background potential energy are observed for strong initial turbulence across the range of $Re$ and $Ri_{b}$ considered here. The primary exception to this trend is the set of simulations with $Ri_{b}=0.20$ and $Re=2000$ (yellow squares), in which diffusion of the mean flow plays a significant role before the instability can reach finite amplitude for ${\mathcal{A}}=0.0025$ and ${\mathcal{A}}=0.01$ ; for these cases, very little mixing occurs due to a lack of vigorous turbulence.

Figure 7.

(a) Overall potential energy change due to irreversible mixing, ${\mathcal{P}}_{b}-D_{p}t$ , as a function of $Re_{b}$ at $t=1$ . (b) Final centreline $Ri_{mean}$ as a function of ${\mathcal{P}}_{b}-D_{p}t$ . (c) Cumulative mixing efficiency, $\unicode[STIX]{x1D702}_{c}$ , calculated with $t_{i}=1$ and $t_{f}$ corresponding to ${\mathcal{M}}={\mathcal{D}}_{p}$ as a function of $Re_{b}$ at $t=1$ for a range of different $Re$ and $Ri_{b}$ . (d) Cumulative mixing efficiency, $\unicode[STIX]{x1D702}_{c}$ , as a function of the maximum available potential energy during a simulation, ${\mathcal{P}}_{a,max}$ .

As figure 3(b) shows, the stability of the final mean flow (quantified by its Richardson number) also depends on the initial turbulence amplitude. In figure 7(b) we show the final value of $Ri_{mean}$ at the centre of the shear layer after the mixing event as a function of total change in ${\mathcal{P}}_{b}$ . Unsurprisingly, there is a clear relationship between the total amount of irreversible mixing that occurs and the stability of the final flow, with large values of $\unicode[STIX]{x0394}{\mathcal{P}}_{b}-{\mathcal{D}}_{p}t_{f}$ leading to large values of centreline $Ri_{mean}$ . This trend is observed across the range of simulations presented here, suggesting that the degree to which a billow is able to form determines the degree to which the mean flow is stabilized.

Our analysis of the cases with $Re=2000$ and $Ri_{b}=0.16$ in the previous subsection showed a strong dependence of the cumulative mixing efficiency on the initial perturbation amplitude, with values dropping from $\unicode[STIX]{x1D702}_{c}\approx 0.27$ to ${\sim}0.13$ as the initial turbulence amplitude is increased (figure 5 d). This trend is observed for the majority of the cases considered here: individual mixing events are typically more efficient for lower values of ${\mathcal{A}}$ (figure 7 c). There is, however, a clear dependence of the mixing efficiency on $Ri_{b}$ , with much higher cumulative mixing efficiencies ( ${\sim}0.3$ ) observed for cases with $Ri_{b}=0.12$ and lower efficiencies ( ${\sim}0.2$ ) for cases with $Ri_{b}=0.20$ . As with the total change in background potential energy, the main exceptions to this trend are cases in which turbulence is weak and, therefore, the diffusion of the background flow plays an important role, namely, the cases with $Re=2000$ and $Ri_{b}=0.20$ , as well as the simulation with $Re=1000$ , $Ri_{b}=0.16$ and ${\mathcal{A}}=0.002$ . The largest overall change in potential energy is not necessarily associated with more efficient mixing (compare figure 7 a,c).

In the previous subsection we showed that a key factor in the efficiency of each mixing event is the degree to which the billow structure, with its associated pre-turbulent peak in instantaneous mixing efficiency, is able to form. Given that the available potential energy ${\mathcal{P}}_{a}$ increases with the size of the billow structure (see figure 4 b), we can use the maximum value of ${\mathcal{P}}_{a}$ as a measure of billow development. Figure 7(d) shows the cumulative mixing efficiency for each of our simulations as a function of ${\mathcal{P}}_{a,max}$ . There is a clear trend in the efficiency of our simulated mixing events across all values of $Re$ , $Ri_{b}$ and ${\mathcal{A}}$ considered here. The cases with strong initial turbulence in which there is little billow growth (and, accordingly, ${\mathcal{P}}_{a,max}$ is small) have lower mixing efficiencies of the order of $\unicode[STIX]{x1D702}_{c}\sim 0.1{-}0.15$ . When larger overturns are able to form and ${\mathcal{P}}_{a,max}$ increases, the mixing efficiency also increases before eventually reaching a plateau at values of $\unicode[STIX]{x1D702}_{c}\sim 0.3$ , consistent with values from previous studies of the Kelvin–Helmholtz instability in laminar initial flows (Peltier & Caulfield Reference Peltier and Caulfield2003; Mashayek et al. Reference Mashayek, Caulfield and Peltier2013; Mashayek & Peltier Reference Mashayek and Peltier2013). This trend even includes the previously identified outlier cases from figure 7(a,c) in which diffusion of the background flow played a key role and led to smaller maximum billow sizes (figure 6 a).

Altogether, figures 6 and 7 demonstrate that initial turbulence significantly affects the development of a stratified shear layer, both qualitatively and quantitatively, across a range of $Re$ and $Ri_{b}$ . These effects are consistent with the results presented in § 3.1, in which initial turbulence may impede the development of the classical Kelvin–Helmholtz billow and, thus, alter the details of the turbulent mixing. For similar initial turbulence amplitudes, these effects appear more pronounced for higher $Re$ and higher $Ri_{b}$ , suggesting a relationship between the rate of decay of the initial turbulence (which is slower for higher $Re$ ) and the growth rate of the underlying instability (which increases for lower $Ri_{b}$ ).