3.1. Flow visualization

Figure 3 shows an instantaneous density profile on an $x$–$z$ plane at $t=60$ and 200 in run 2. Vortices generated by the Kelvin–Helmholtz instability can be seen in figure 3(a). After these vortices collapse, the flow becomes turbulent with small-scale density fluctuations in figure 3(b). This temporal development of the stably stratified shear layer agrees well with the previous DNS with the same Reynolds number and Richardson number (Watanabe et al. Reference Watanabe, Riley, Nagata, Matsuda and Onishi2019a).

Figure 3. Density $\bar {\rho }$ on an $x$–$z$ plane at (a) $t=60$ and (b) $t=200$ in run 2. Only a small part of the computational domain is shown here.

Figure 4 visualizes streamwise velocity fluctuations $\bar {u}'$ in run 2 on the horizontal plane at $z=0$ from $t=2$ to 320. The characteristic length scale of the velocity fluctuations increases with time. At late time, the velocity distribution is highly anisotropic, and regions with positive and negative $\bar {u}'$ are elongated in the $x$ direction. These are related to the ELSS, which were also found in previous DNS (Watanabe et al. Reference Watanabe, Riley, Nagata, Matsuda and Onishi2019a). Figure 4(f) depicts the length of the shear layer thickness $\delta _u$, which is defined with the mean streamwise velocity $\langle \bar {u}\rangle$ as the vertical distance between $\langle \bar {u}\rangle =-0.49$ and 0.49. The vertical profile of $\langle \bar {u}\rangle$ is presented in the next subsection. The streamwise length of the ELSS is as large as $10\delta _u$ while the spanwise length is about $\delta _u$. The streamwise and spanwise sizes of the ELSS are smaller than $L_x$ and $L_y$ even at $t=320$, respectively, and the computational domain contains multiple ELSS identified by $\bar {u}'$. The longitudinal auto-correlation functions of $u$ and $v$ are defined as

(3.1)$$\begin{gather} f_u(r_x,z,t)=\langle u'(x,y,z,t)u'(x+r_x,y,z,t)\rangle / u_{rms}^2(z,t), \end{gather}$$

(3.2)$$\begin{gather}f_v(r_y,z,t)=\langle v'(x,y,z,t)v'(x,y+r_y,z,t)\rangle / v_{rms}^2(z,t). \end{gather}$$

Here, $r_i$ is the distance between two points separated in the $i$ direction. If the domain is small compared with the length scale of the flow, $f_u$ and $f_v$ do not reach zero at large scales because the periodic boundary conditions generate artificial correlation between two points (Davidson Reference Davidson2004). In run 2, $f_u(r_x)$ and $f_v(r_y)$ at $z=0$ and $t=320$ reach 0 at $r_x=41.1$ and $r_y=8.2$, which are about $L_x/10$ and $L_y/10$, respectively. Thus, the domain is large enough to prevent the finite domain size from affecting the development of turbulence. Velocity fluctuations at $t=2$ are related to the initial velocity fluctuations. The roller vortices of the Kelvin–Helmholtz instability exist at $t=60$. The diameter of the roller vortices at $t=60$ is ${{O}}(h_0)$ in figure 3(a), and this length is much smaller than the streamwise length of the ELSS, which is as large as $10^2h_0$. Thus, the velocity fluctuations after $t=120$ are hardly correlated with those at $t=2$ and 60.

Figure 4. Streamwise velocity fluctuations $\bar {u}'$ on the horizontal plane at the centre of the shear layer in run 2: (a) $t =2$; (b) $t =60$; (c) $t =120$; (d) $t =180$; (e) $t =240$; (f) $t =320$. Figure (f) depicts the size of $\delta _u$ and $10\delta _u$, where $\delta _u$ is the shear layer thickness defined with the mean velocity profile.

3.2. Fundamental statistics of stably stratified shear layers

The LES is validated by comparison between runs 1 (DNS) and 2 (LES). Figure 5(a) compares vertical profiles of mean streamwise velocity and density. These profiles show that the shear layer grows from $t=100$ to 240 because both mean streamwise velocity (momentum) and density are transferred in the vertical direction by turbulence. Figure 5(b) shows the temporal evolution of r.m.s. velocity fluctuations at $z=0$ and shear layer thickness $\delta _u$. The figure also includes $u_{rms}$ in five LES before the ensemble average is taken. The difference among the five simulations is small, and the domain size is large enough to obtain well-converged statistics for r.m.s. quantities. When the turbulence is generated by the instability, the r.m.s. velocity fluctuations increase with time. The fluctuations begin to decay after they reach a peak. The shear layer thickness also increases at early time. However, the growth of the shear layer thickness becomes slow after $t\approx 150$ because of the stable stratification. These statistics are almost identical between the LES and DNS.

Figure 5. Comparisons between DNS (run 1) and LES (run 2). (a) Vertical profiles of mean velocity and density at $t=100$ and 240. (b) Temporal variations of root-mean-squared velocity fluctuations at $z=0$ and shear layer thickness $\delta _u$. Grey broken lines represent $u_{rms}$ in five simulations of LES.

Figure 6(a) compares the energy spectrum of streamwise velocity, $E_u$, between runs 1 and 2. The LES and DNS results agree well with each other for a low wavenumber range while $E_u$ in the LES sharply drops for $k_x\gtrsim 10$. The spectrum at $t=60$ has a peak at $k_x\approx 0.6$, which is associated with the roller vortices. The wavelength of this peak is $\lambda _x\approx 10$, which is much smaller than the streamwise length of the ELSS at $t=320$. Therefore, the length scales of large-scale fluctuations are different between the transitional and fully developed states. As also found in figure 4, the velocity fluctuations during the transition are not directly related to those at late time. The figure also shows $E_u$ obtained in each run of the LES. The shape of $E_u$ at $t=60$ is qualitatively similar among nine simulations, and the size and spacing of the roller vortices are statistically identical among all simulations. Quantitative variations among different simulations are noticeable for a low wavenumber range, for which the convergence of spectra is improved by taking ensemble averages of 35 simulations.

Figure 6. (a) Energy spectrum of streamwise velocity $E_u$ at the centre of the shear layer at $t=60$, 140 and 320 in LES (run 2). Direct numerical simulation results at $t=140$ and 320 are also shown for comparison. Grey lines are $E_u$ in each simulation of LES (nine simulations). (b) Vertical profiles of turbulent kinetic energy dissipation rate at $t=140$ in DNS ($\varepsilon ^{(k)}$ in run 1) and LES ($\varepsilon _{G}^{(k)}+\varepsilon _{S}^{(k)}$ in run 2). The dissipation terms of grid scales ($\varepsilon _{G}^{(k)}$) and subgrid scales ($\varepsilon _{S}^{(k)}$) are also shown in the figure.

In implicit LES the dissipation terms of turbulent kinetic energy by the grid scales and the SGS model are given respectively by

(3.3a,b)\begin{equation} \varepsilon_{G}^{(k)}=\frac{1}{{Re}} \left\langle\frac{\partial \bar{u}_{i}'}{\partial x_{j}}\frac{\partial \bar{u}_{i}'}{\partial x_{j}}\right\rangle, \quad \varepsilon_{S}^{(k)}={-}\left\langle \bar{u}_{i}'F(\bar{u}_{i}')\right\rangle. \end{equation}

Figure 6(b) shows the turbulent kinetic energy dissipation rate at $t=140$ in runs 1 and 2. The total dissipation rate $\varepsilon _{G}^{(k)}+\varepsilon _{S}^{(k)}$ in the LES agrees well with the dissipation rate in the DNS, which is calculated as $\varepsilon ^{(k)}=\varepsilon _{G}^{(k)}$ by replacing $\bar {u}_{i}'$ with $u_i'$ in (3.3a,b). Therefore, the amount of dissipation is accurately controlled by the implicit SGS model. In the LES, $\varepsilon _{S}^{(k)}$ is as large as $\varepsilon _{G}^{(k)}$, and about one half of the energy dissipation is treated with the model at $t=140$. The implicit LES does not use the dynamic adjustment of the filter. The averaged energy dissipation rate in turbulence is generally controlled by the averaged interscale energy flux from large to small scales. The filter dissipates the energy at resolved small scales before the energy is transferred to further small scales. The good agreement in the dissipation rate between LES and DNS may be explained by the spatial resolution of LES, which is good enough to resolve large-scale fluctuations that determine the averaged interscale energy flux.

Figure 7(a) shows the temporal evolution of $\delta _u$ at ${Re}=40\,000$ in runs 3–6. In the non-stratified turbulent shear layer ($Ri=0$), $\delta _u$ almost linearly increases with time for $t\gtrsim 100$. As also found at ${Re}=2000$, the growth of $\delta _u$ at ${Re}=40\,000$ is very slow for $t\gtrsim 150$ at $Ri=0.06$ and 0.1. However, $\delta _u$ at $Ri=0.02$ gradually increases even for $t\gtrsim 150$. Figure 7(b) shows $u_{rms}$ on the shear layer centreline at ${Re}=40\,000$. In a self-similar turbulent mixing layer, $u_{rms}$ at $z=0$ does not vary with time (Rogers & Moser Reference Rogers and Moser1994). Indeed, $u_{rms}\approx 0.17$ at $Ri=0$ hardly varies after the turbulent shear layer has fully developed, and this value of 0.17 agrees well with previous studies on self-similar turbulent mixing layers (Rogers & Moser Reference Rogers and Moser1994). The time of the peak in $u_{rms}$ depends on $Ri$. Furthermore, comparison of $u_{rms}$ between figures 5(b) and 7(b) also confirms that $u_{rms}$ at $Ri=0.06$ reaches a peak slightly earlier for ${Re}=40\,000$ than for ${Re}=2000$. These differences in the initial growth of $u_{rms}$ imply that the detailed transition process depends on ${Re}$ and $Ri$.

Figure 7. Temporal variations of (a) shear layer thickness $\delta _u$ and (b) r.m.s. streamwise velocity fluctuations at $z=0$ for ${Re}=40\,000$ (runs 3–6).

The comparisons between the DNS and LES confirm that the LES accurately captures large-scale velocity and density fluctuations of the flow. The LES results will be used to assess the statistical properties of the stably stratified shear layer. In the rest of the paper most results are presented for $t=140$, 240 and 320. Table 2 summarizes important parameters of the flow obtained by LES at these time instances, where $L_O=\sqrt {\varepsilon ^{(k)}/N^3}$ is the Ozmidov scale defined with the buoyancy frequency $\tilde {N}=\sqrt {-(g/\rho _a) \langle \partial \tilde {\rho }/\partial \tilde {z} \rangle }$ $(N=\tilde {N}t_r)$, $\eta =({Re}^{3}\varepsilon ^{(k)})^{-1/4}$ is the Kolmogorov scale, ${Re}_b=(L_O/\eta )^{4/3}$ is the buoyancy Reynolds number and $Ri_g=N^2/\langle \partial u /\partial z\rangle ^{2}$ is the gradient Richardson number. The turbulent kinetic energy dissipation rate in LES is calculated as $\varepsilon ^{(k)}=\varepsilon _{G}^{(k)}+\varepsilon _{S}^{(k)}$. The table shows $L_{O}$, $\eta$, ${Re}_b$ and $Ri_g$ obtained at $z=0$. Here, $L_{O}$ and $\eta$ are normalized by $\delta _u$ because the energy spectra will be presented as functions of a wavelength normalized by $\delta _u$. In the rest of the paper, $z$ is also normalized by $\delta _u$ because $z/\delta _u$ is useful to indicate the location with respect to the shear layer thickness. In run 2 the spatial resolution $\varDelta _x$ is $2.2\eta$-$8.0\eta$, where $\varDelta _x/\eta$ decreases with time. Thus, most length scales are resolved in run 2, which will be used to evaluate the transport equations of energy spectra. Runs 3–6 have $\varDelta _x/\eta ={{O}}(10^1)$ and large-scale velocity and density fluctuations are well resolved in LES. The Ozmidov scale is often interpreted as the smallest length scale that is affected by the buoyancy, and it is considered that vertical turbulent motions at scales greater than $L_O$ are significantly damped. The buoyancy Reynolds number is defined as the ratio between $L_O$ and $\eta$. As $\eta$ represents the smallest scale of turbulence, ${Re}_b\gg 1$ indicates that small-scale turbulent motions are not directly influenced by the stable stratification. This is comfortably the case at $t=140$ with ${Re}_b=93$ for $({Re},Ri)=(2000,0.06)$. However, smaller values of ${Re}_b$ at the later time indicate that turbulent motions at the smallest scale are strongly suppressed by the buoyancy. At ${Re}=40\,000$, ${Re}_b$ is larger than ${{O}}(10^1)$ except at late time with $Ri=0.1$. At $Ri=0.06$ and 0.1, $Ri_g=0.42 - 0.56$ weakly depends on time after $t=140$. On the other hand, $Ri_g$ increases from 0.16 to 0.24 at $({Re}, Ri)=(40\,000, 0.02)$ because the shear layer thickness also grows with time.

Table 2. Statistical properties of the stably stratified shear layer in LES. Values of $L_{O}$, $\eta$, ${Re}_b$ and $Ri_g$ are taken at the centre of the shear layers.

3.3. Energy spectra of velocity and density

The energy spectra of velocity and density fluctuations are calculated with Fourier transform in the $x$ direction. These spectra, premultiplied by the streamwise wavenumber $k_x$, are plotted against the vertical location $z$ and the streamwise wavelength $\lambda _x=2{\rm \pi} /k_x$, where $\lambda _x$ is shown on a logarithmic scale. When $E_u(k_x)$ at a given vertical position is shown in a linear plot, the visible area beneath the curve between wavenumbers $k_{x1}$ and $k_{x2}$ can be interpreted as the contribution to $\langle u'^2\rangle$ from $k_{x1}\leqslant k_{x}\leqslant k_{x2}$ because of $\langle u'^2\rangle =\int _0^{\infty } E_u(k_x)\,\textrm {d}k_x$. The equivalent relation for the logarithmic coordinate of $k_{x}$ is $\langle u'^2\rangle =\int _{-\infty }^{\infty } k_{x}E_ud(\log _{10}{k_x})$, which is obtained with $dk_x=k_xd(\log _{10}{k_x})$. Therefore, when the spectrum is plotted on the logarithmic coordinate of $k_x$, $k_xE_u$ is useful because it is the contribution to $\langle u'^2\rangle$ on the coordinate of $\log _{10}{k_x}$. This meaning of the premultiplied spectrum is held even if $\lambda _x=2{\rm \pi} /k_x$ is used instead of $k_x$ because $d\log _{10}{k_x}$ is proportional to $-d\log _{10}{\lambda _x}$. As the turbulent kinetic energy is distributed over a wide range of scales, the premultiplied spectrum is widely used to assess the scale dependence of turbulent kinetic energy (Kim & Adrian Reference Kim and Adrian1999; Hutchins & Marusic Reference Hutchins and Marusic2007; Takamure et al. Reference Takamure, Ito, Sakai, Iwano and Hayase2018). Similarly, premultiplied cospectra related to turbulent fluxes of momentum and density have been used in previous studies on stratified turbulent flows (Lienhard & van Atta Reference Lienhard and van Atta1990; Gerz & Schumann Reference Gerz and Schumann1996; Komori & Nagata Reference Komori and Nagata1996).

Figure 8(a–c) shows the energy spectrum of streamwise velocity, $k_{x}E_u$, at $t=140$, 240 and 320 in run 2. As also observed in previous DNS (Watanabe et al. Reference Watanabe, Riley, Nagata, Matsuda and Onishi2019a), the spectral shape significantly changes with time. At $t=140$, $k_{x}E_{u}$ at each vertical location has a single peak around $\lambda _x/\delta _u\approx 2$. At late time, the middle of the shear layer has a spectral shape with two peaks at $\lambda _x/\delta _u\approx 1$ and 10 while a single peak can be identified for $z/\delta _u\gtrsim 0.3$. The peak in $k_{x}E_{u}$ at $\lambda _x/\delta _u \approx 10$ is related to the ELSS visualized in figure 4(f), where regions with positive and negative fluctuations extend over about $10\delta _u$ in the $x$ direction. A similar spectral evolution can be found for $E_v$ and $E_w$ shown in figure 8(d–i). Figure 8(j-l) shows $k_{x}E_{\rho }$ at the same time instances. As found for the kinetic energy spectra, $k_{x}E_{\rho }$ in the middle of the shear layer at $t=320$ has two peaks at $\lambda _x/\delta _u \approx 0.6$ and $8$. However, $k_{x}E_{\rho }$ decreases with $z$ while the velocity spectra have a large peak at $z/\delta _u\approx 0.4$. These spectral shapes in the LES generally agree with the DNS results in Watanabe et al. (Reference Watanabe, Riley, Nagata, Matsuda and Onishi2019a). However, the peak locations are more clearly identified in the LES since ensemble averages are taken in this study.

Figure 8. Premultiplied energy spectra of velocity and density, $k_{x}E_{\alpha }$ ($\alpha =u, v, w$ or $\rho$) in run 2 (${Re}=2000$, $Ri=0.06$), plotted against the streamwise wavelength $\lambda _x=2{\rm \pi} /k_x$ and the vertical location $z$: (a,d,g,j) $t=140$; (b,e,h,k) $t=240$; (c,f,i,l) $t=320$. (a–c) Streamwise velocity $E_u$; (d–f) spanwise velocity $E_v$; (g–i) vertical velocity $E_w$; (j–l) density $E_\rho$. The colour contours are normalized by the peak values at specific times.

The double-peak spectra at late time indicate that the flow in the middle of the shear layer has two energy-containing length scales, which can be identified as the peak wavelengths. One is the short peak-wavelength of the spectra, which is close to the shear layer thickness $\delta _u$. In turbulent free shear flows, such as jets, wakes and mixing layers, a characteristic length of LSM is the width of the flow (Pope Reference Pope2000). Therefore, the peak at $\lambda _x\sim \delta _u$ is associated with usual LSM, which appear in a turbulent shear layer even without stable stratification. Fluid motions at this length scale are simply referred to as ‘LSM’. The other energy-containing length scale is obtained as the long peak-wavelength, which is much larger than $\delta _u$. This length scale is related to the streamwise length of elongated patterns of velocity fluctuations visualized in figure 4(f). Therefore, fluid motions associated with fluctuations at the long peak-wavelength are called ‘ELSS’.

The peak in $k_{x}E_{u}$ associated with the ELSS is found only for $|z|/\delta _u \lesssim 0.3$. The ELSS hardly reach the interfacial region between the turbulent shear layer and the non-turbulent region, which are separated by a thin interfacial layer (turbulent/non-turbulent interface). Visualization of the turbulent/non-turbulent interface in a stably stratified shear layer showed that the interface geometry does not exhibit an imprint of the ELSS (Watanabe, Riley & Nagata Reference Watanabe, Riley and Nagata2016a). This is because the ELSS appear in the middle of the shear layer while the interface is formed at the outer edge of the shear layer. The ELSS and the interface are separated in space, and it is expected that the ELSS do not have significant influences on the turbulent/non-turbulent interface.

The existence of LSM and ELSS for ${Re}<2000$ was reported in Watanabe et al. (Reference Watanabe, Riley, Nagata, Matsuda and Onishi2019a). Runs 3–6 are used to examine them at high ${Re}$. Figure 9 shows $k_x E_{u}$ at $z=0$ for ${Re}=40\,000$. The spectrum $k_x E_{u}$ is normalized by its maximum value $[k_x E_{u}]_{max}$ at $z=0$ in order to compare the plots at different time. The formation of the double-peak spectrum occurs for $Ri=0.06$ and 0.1. For low $Ri$, $k_x E_{u}$ has a single peak and the flow has only one characteristic length scale at large scales. This large single peak is expected from the presence of large-scale roller vortices in a turbulent mixing layer without density stratification (Balaras, Piomelli & Wallace Reference Balaras, Piomelli and Wallace2001). The peak wavelength at $Ri=0$ and 0.02 is about $\lambda _x/\delta _u=2 - 3$, which is much smaller than the size of ELSS. The ELSS corresponding to the spectral peak at $\lambda _x/\delta _u\approx 10^{1}$ appear only at moderately large $Ri$. This peak at $\lambda _x/\delta _u\approx 10^{1}$ appears earlier at $Ri=0.1$ than at $Ri=0.06$. Thus, the ELSS develop faster for higher $Ri$.

Figure 9. Energy spectra of streamwise velocity at $z=0$ at $t=140, 240$ and 320 in runs 3–6 $({Re}=40\,000)$: (a) $Ri=0$ (run 3); (b) $Ri=0.02$ (run 4); (c) $Ri=0.06$ (run 5); (d) $Ri=0.1$ (run 6). For comparison of the spectra at different $t$, $k_x E_{u}$ is normalized by its maximum value $[k_x E_{u}]_{max}$ at each time instance.

3.4. Characteristics of large-scale velocity and density fluctuations

The scale dependence of turbulent statistics is investigated with the energy spectra in run 2. The characteristic length scales of LSM and ELSS are estimated with the spectral peaks. Figure 10(a) shows the turbulent kinetic energy spectrum $E_{k}=(E_u+E_v+E_w)/2$ and density spectrum $E_\rho$ at $t=260$ on the centreline, where two peaks can be identified for both spectra. Here, the short peak-wavelength of $E_{\rho }$ is denoted by $\varLambda _{S\rho }$ while the long peak-wavelength is denoted by $\varLambda _{L\rho }$. Similarly, $\varLambda _{Sk}$ and $\varLambda _{Lk}$ represent the short and long peak-wavelengths of $E_k$, respectively. The premultiplied spectra at $z=0$ are plotted on the discrete coordinate of $\lambda _x$, where one can find the peaks in $E_{k}$ and $E_\rho$. A third-order polynomial curve fitting with five discrete points centred at the peak of $E_k$ or $E_\rho$ is applied to determine the peak wavelength from the fitted curve. These peak wavelengths are the characteristic length scales of LSM and ELSS. The spectra before $t=220$ have only a single peak at $z=0$ and the ELSS have not formed yet. Figure 10(b) shows temporal variations of the peak wavelengths obtained at $z=0$ after $t=220$. Here, these wavelengths are normalized by $\delta _u$, which hardly increases with time for $t\geqslant 220$ as shown in figure 5(b). The short peak-wavelength weakly depends on time: $\varLambda _{S\rho }$ and $\varLambda _{Sk}$ stay between $0.5\delta _u$ and $0.7\delta _u$. On the other hand, the long peak-wavelengths increase with time approximately from $5\delta _u$ to $9\delta _u$. Therefore, the length scale of LSM hardly grows with time while the length of ELSS monotonically increases. The increase of $\varLambda _{Lk}$ and $\varLambda _{L\rho }$ can be caused by the deformation due to mean shear. In the stratified shear layer a time scale of turbulence increases as the turbulence decays, while a time scale of the mean shear $\sim U_0/\delta _u$ is almost independent of time because $\delta _u$ hardly changes with time once the stably stratified turbulent shear layer has fully developed. Therefore, the turbulence time scale becomes much longer than the shear time scale at late time. In the stably stratified shear layer investigated in this study, the integral shear parameter defined as a turbulence-to-shear time scale ratio exceeds 15 (Watanabe et al. Reference Watanabe, Riley, Nagata, Matsuda and Onishi2019a), which is larger than typical values of 3–6 in non-stratified turbulent free shear flows (Pope Reference Pope2000). In general, turbulence with larger scales evolves more slowly. Therefore, the turbulent motions, especially at scales of the ELSS, evolve much more slowly than the shear time scale, and its evolution is possibly approximated by rapid distortion theory (Savill Reference Savill1987). The theory predicts that the longitudinal integral length scale in the mean flow direction monotonically increases with time in turbulence with a simple shear (Lee, Kim & Moin Reference Lee, Kim and Moin1990). This behaviour is consistent with the growth of $\varLambda _{Lk}$ and $\varLambda _{L\rho }$ found in figure 10(b).

Figure 10. (a) Energy spectra of turbulent kinetic energy ($E_k$) and density ($E_\rho$) at $t=260$ with the definition of short and long peak-wavelengths. The spectra multiplied by $k_x$ are normalized by turbulent kinetic energy $k_T= (\langle {\bar {u}^\prime}{^2}\rangle +\langle {\bar {v}^\prime}{^2}\rangle +\langle {\bar {w}^\prime}{^2}\rangle )/2$ or density variance $\langle {\bar {\rho }^\prime}{^2}\rangle$. (b) Temporal evolution of short and long peak-wavelengths in $E_k$ and $E_\rho$. These results are taken at the centre of the shear layer in run 2.

The degree of anisotropy of the Reynolds stress tensor $R_{ij}=\langle \bar {u}_{i}'\bar {u}_{j}'\rangle$ is evaluated with the normalized anisotropy tensor defined as

(3.4)\begin{equation} b_{ij}=\frac{R_{ij}}{R_{kk}}-\frac{1}{3}\delta_{ij}. \end{equation}

The state of anisotropy of $R_{ij}$ is often expressed by the Lumley triangle, which is defined with $\eta _b=(b_{ij}b_{ji}/6)^{1/2}$ and $\xi _b=(b_{ij}b_{jk}b_{ki}/6)^{1/3}$ (Pope Reference Pope2000). The isotropic Reynolds stress has $(\xi _b,\eta _b)=(0,0)$. Although the Lumley triangle was also investigated in DNS of stratified shear layers (Smyth & Moum Reference Smyth and Moum2000a; Wingstedt et al. Reference Wingstedt, Fossum, Pettersson Reif and Werne2015), small computational domains used in previous DNS did not allow the growth of ELSS. However, the spectral energy density for ELSS is not negligible in figure 8, and the ELSS might have significant influences on the anisotropy of Reynolds stress. The Lumley triangle is often used for turbulence modelling but is used as a diagnostic tool for the isotropy of Reynolds stress in this study. Because $b_{ij}$ is defined with $R_{ij}$ normalized by $R_{kk}$, it is useful to compare $(\xi _b,\eta _b)$ among different flows even if these flows have different magnitudes of the turbulent kinetic energy. For example, Gargett & Wells (Reference Gargett and Wells2007) provided the discussion on $\eta _b$ and $\xi _b$ in observations of ocean turbulence by comparing them with canonical turbulent flows, such as turbulent boundary layers.

Figure 11(a) shows variations of $(\xi _b,\eta _b)$ from $z=0$ to $0.5\delta _u$ at $t=140$ and $320$. The LES results (run 2) are compared with the results obtained at different heights of a turbulent boundary layer (Watanabe et al. Reference Watanabe, Zhang and Nagata2019b) and at the centre of a non-stratified, fully developed turbulent mixing layer (Takamure et al. Reference Takamure, Ito, Sakai, Iwano and Hayase2018). Here, the vertical location in the turbulent boundary layer, $y$, is normalized by the viscous length $\delta _\nu =\nu /u_\tau$ defined with the friction velocity $u_\tau =\sqrt {\tau /\rho }$ ($\tau$: wall shear stress) as $y^{+}=y/\delta _\nu$. Three locations, $y^{+}=5$, $30$ and $800$, represent the approximate locations that separate the viscous sublayer ($y^{+}\lesssim 5$), the buffer layer ($5\lesssim y^{+}\lesssim 30$), the logarithmic layer ($30\lesssim y^{+}\lesssim 800$) and the outer region ($y^{+}\gtrsim 800$). At $t=140$, the centre of the stably stratified shear layer has $(\xi _b,\eta _b)\approx (0.08,0.1)$ (a green circle in figure 11a), which is close to the outer region of the turbulent boundary layer. At $t=320$, both $\xi _b$ and $\eta _b$ become large compared with early time. Here, the centre of the shear layer has $(\xi _b, \eta _b)=(0.16,0.17)$ (a red circle in figure 11a), and the Reynolds stress is more anisotropic than at the early time. In addition, $(\xi _b, \eta _b)\approx (0.16,0.17)$ is within the logarithmic layer in case of the turbulent boundary layer. On the other hand, the region away from the centreline still attains $(\xi _b, \eta _b)$ close to the outer region of the turbulent boundary layer.

Figure 11. (a) Variation of $\xi _b$ and $\eta _b$ between $z=0$ and $0.5\delta _u$ at $t=140$ and 320 in run 2. An enlarged figure is also shown for clarity. (b) Plots of $\xi _b$ and $\eta _b$ associated with large and small peak-wavelengths of the turbulent kinetic energy spectrum at $z=0$ and $t=320$ in run 2. Three thin broken lines correspond to an axisymmetric state with one large eigenvalue, $\eta _b=\xi _b$, or one small eigenvalue, $\eta _b=-\xi _b$, and a two-dimensional state $\eta _b=(1/27+2\xi _b^3)^{1/2}$. Figures also include $\xi _b$ and $\eta _b$ obtained in the turbulent boundary layer at the momentum thickness Reynolds number of ${Re}_{\theta }=13\,000$ in Watanabe et al. (Reference Watanabe, Zhang and Nagata2019b) and at the centre of the non-stratified, fully developed turbulent mixing layer (Takamure et al. Reference Takamure, Ito, Sakai, Iwano and Hayase2018).

The flow has two characteristic energy-containing length scales represented by $\varLambda _{Sk}$ and $\varLambda _{Lk}$. We introduce $k_{min}$, which represents the wavenumber corresponding to the minimum value of $k_{x}E_k$ at $z=0$ between the two peaks at $\varLambda _{Sk}$ and $\varLambda _{Lk}$ (see figure 10a). The contributions to the Reynolds stress from fluctuations associated with $\varLambda _{Sk}$ (LSM) and $\varLambda _{Lk}$ (ELSS) are respectively evaluated by integrating a corresponding cospectrum as

(3.5a,b) \begin{equation} R_{ij}^{(S)}=\int_{k_{min}}^{\infty}{{\rm Re}}[\langle \widehat{\bar{u}_{i}'}\widehat{\bar{u}_{j}'}^{*}\rangle]\,\textrm{d}k_x,\quad R_{ij}^{(L)}=\int_{0}^{k_{min}}{{\rm Re}}[\langle \widehat{\bar{u}_{i}'}\widehat{\bar{u}_{j}'}^{*}\rangle]\,\textrm{d}k_x. \end{equation}

The anisotropy of velocity fluctuations associated with the LSM and ELSS is evaluated with $\eta _b$ and $\xi _b$ calculated separately from $R_{ij}^{(S)}$ and $R_{ij}^{(L)}$. Figure 11(b) compares $\eta _b$ and $\xi _b$ obtained from $R_{ij}$, $R_{ij}^{(S)}$ and $R_{ij}^{(L)}$ at $z=0$ and $t=320$. For fluctuations associated with the large peak-wavelength, $(\eta _b,\xi _b)\approx (0.2,0.2)$ is obtained from $R_{ij}^{(L)}$. The turbulent boundary layer also has $(\eta _b,\xi _b)\approx (0.2,0.2)$ between the buffer layer and the logarithmic layer, and the anisotropy of the Reynolds stress for $\varLambda _{Lk}$ (ELSS) is similar to the near-wall region. However, $\eta _b$ and $\xi _b$ obtained from $R_{ij}^{(S)}$ are close to the values in the outer region. Therefore, the highly anisotropic Reynolds stress $R_{ij}$ near the centre of the shear layer is caused by the ELSS with an approximate length of $\varLambda _{Lk}$. In figure 11(a), $(\eta _b,\xi _b)$ near the edge of the shear layer at $t=320$ (a red cross in figure 11a) are still close to the values in the outer region of the turbulent boundary layer. This is explained by the absence of the ELSS near the edge of the stably stratified shear layer since the long-wavelength peak at $\lambda _x\approx \varLambda _{Lk}$ in the energy spectra does not exist at large $z$ in figure 8(c,f,i).

A fractional contribution to the total Reynolds stress from $R_{ij}^{(L)}$ and $R_{ij}^{(S)}$ can be defined as $F_S = R_{ij}^{(S)}/R_{ij}$ and $F_L = R_{ij}^{(L)}/R_{ij}$ for each component of $i$ and $j$. Similarly, $F_S$ and $F_L$ for density variance $\langle {\bar {\rho }^\prime}{^2}\rangle$ and vertical density flux $\langle \bar {\rho }'\bar {w}'\rangle$ can be defined with the integral of the corresponding spectrum and cospectrum. Figure 12(a) shows temporal variations of $F_S$ and $F_L$ at $z=0$ for velocity and density variances. For all components, $F_S$ decreases and $F_L$ increases with time, and the contribution from the ELSS relatively becomes more important with time for velocity and density fluctuations. Figure 12(b) shows $F_S$ and $F_L$ obtained for the vertical turbulent fluxes of $u$ and $\rho$. Unlike velocity and density variances, $F_L$ stays around 0 even at late time, and the LSM with scales close to $\varLambda _{Sk}$ dominate the vertical transfer of momentum and density. The Reynolds stress and the vertical fluxes become small at late time, and $F_S$ and $F_L$ only provide the relative contribution at a given time. Therefore, large $F_L$ does not mean that the ELSS have a dominant contribution to the energetic property over the shear layer development. Figure 12(c) plots $R_{ij}^{(L)}$ and $R_{ij}^{(S)}$ for $(i,j)=(1,1)$ and $(1,3)$. Although $R_{11}^{(S)}$ decays with time, $R_{11}^{(L)} \approx 1$ slowly varies and the velocity fluctuations of the ELSS hardly decay. Similar results are obtained for $R_{22}$ and $R_{33}$ (not shown here). On the other hand, $R_{13}^{(L)}$ is close to 0 while $R_{13}^{(S)}$ also decreases to 0 with time. Thus, the ELSS hardly contribute to the vertical transport of momentum. Large scatters are found for $F_L$ and $F_S$ at late time in figure 12(b). This is caused by small $\langle \bar {u}'\bar {w}'\rangle$ and $\langle \bar {\rho }'\bar {w}'\rangle$, which are the denominators of $F_S$ and $F_L$. Figure 12(a,b) confirm that the ELSS carry an important fraction of the turbulent kinetic energy while the vertical transport of momentum due to the ELSS is not significant. This is due to the weak correlation between $u$ and $w$. In turbulent shear flows the decay of turbulent kinetic energy is caused by the viscous dissipation while the destruction of the Reynolds stress $\langle u'w'\rangle$ is caused by the pressure–strain correlation (Pope Reference Pope2000). The viscous dissipation is active at small scales. However, the destruction of the Reynolds stress actively occurs at larger scales than the viscous dissipation (Takahashi et al. Reference Takahashi, Iwano, Sakai and Ito2019). The different length scales of the decay process between the turbulent kinetic energy and the Reynolds stress is a possible reason for the weak vertical momentum transport by the ELSS with large turbulent kinetic energy.

Figure 12. Fractional contributions of fluctuations associated with long and short peak-wavelengths ($F_L$ and $F_S$) of $E_k$ to (a) velocity and density variances, (b) vertical turbulent fluxes of streamwise velocity and density. (c) Reynolds stress associated with long and short peak-wavelengths, $R_{ij}^{(L)}$ and $R_{ij}^{(S)}$, for $(i,j)=(1,1)$ and $(1,3)$. Figure (d) shows $F_L$ and $F_S$ for the dissipation rates of turbulent kinetic energy and density variance. These results are taken at $z=0$ in run 2.

For the turbulent kinetic energy dissipation rate, $F_S$ and $F_L$ are obtained by integrating $\hat {\varepsilon }^{(k)}= \hat {\varepsilon }_{G}^{(u)} +\hat {\varepsilon }_{G}^{(v)} +\hat {\varepsilon }_{G}^{(w)} +\hat {\varepsilon }_{S}^{(u)} +\hat {\varepsilon }_{S}^{(v)} +\hat {\varepsilon }_{S}^{(w)}$, i.e.

(3.6a,b)\begin{equation} F_S=\frac{\displaystyle\int_{k_{min}}^{\infty} \hat{\varepsilon}^{(k)}\,\textrm{d}k_x} {\displaystyle\int_{0}^{\infty} \hat{\varepsilon}^{(k)} \,\textrm{d}k_x}, \quad F_L={=}\frac{\displaystyle\int_{0}^{k_{min}} \hat{\varepsilon}^{(k)} \,\textrm{d}k_x} {\displaystyle\int_{0}^{\infty} \hat{\varepsilon}^{(k)} \,\textrm{d}k_x}. \end{equation}

Similarly, $F_S$ and $F_L$ for the dissipation rate of density variance are defined with the integrals of the dissipation term $\hat {\varepsilon }^{(\rho )} = \hat {\varepsilon }_{G}^{(\rho )} + \hat {\varepsilon }_{S}^{(\rho )}$. Figure 12 shows $F_S$ and $F_L$ for $\hat {\varepsilon }^{(k)}$ and $\hat {\varepsilon }^{(\rho )}$. Even at late time, $F_L$ is close to 0, and the ELSS have only a small contribution to the dissipation. It should be noted that $F_S\approx 1$ does not mean that the LSM causes the dissipation because $F_S$ includes the contribution from velocity and density fluctuations at smaller scales than the LSM.

3.5. Energy budget in wavenumber space

The spectral budget equations are analysed to investigate temporal evolutions of the spectra. Here, each term in the budget equations is premultiplied by $k_x$, and is plotted as a function of $\lambda _x/\delta _u$ at a given location of $z/\delta _u$, where a logarithmic scale is used for $\lambda _x$. The LES results (run 2) are presented for $z/\delta _u=0$ and $0.4$ and for $t=140$ and 240 because the spectral shapes are very different between these heights and times. The diffusive terms $\hat {D}^{(u_\alpha )}$ and $\hat {D}^{(\rho )}$ in (2.14) and (2.16) are about $10^{1} - 10^{2}$ times smaller than other spatial transport terms, and these terms are not shown in this paper.

3.5.1. Spectral energy budget at early time

Figure 13(a–f) presents the spectral energy budget (2.14) at $t=140$ for $E_u$, $E_v$ and $E_w$. The left and right columns of the figure show the results at $z=0$ and $0.4\delta _u$, respectively. In figure 13(a,b) the production term $\hat {P}^{(u)}$ is mostly positive for a wide range of $\lambda _x$. At both vertical locations, $\hat {P}^{(u)}$ peaks at $\lambda _x/\delta _u\approx 1$, and the energy production is the most active for the length scale close to $\delta _u$, i.e.  the length scale of the LSM. This is consistent with canonical turbulent free shear flows, where the energy production occurs at large scales, which can be characterized by the width of the flow in the transverse (cross-flow) direction. Positive production can also be seen for very large scales $\lambda _x/\delta _u\approx 10^1$ at $z=0$, and the energy production at very large scales already occurs at the early time. This energy produced at $\lambda _x/\delta _u\approx 10^1$ is removed by the interscale transport term $\hat {T}_{S}^{(u)}$.

Figure 13. Spectral energy budgets of (a,b) streamwise velocity, (c,d) spanwise velocity, (e,f) vertical velocity and (g,h) density at $t=140$ in run 2. Results are shown for (a,c,e,g) $z=0$ and (b,d,f,h) $z=0.4\delta _u$. Each term in the energy budget equations is plotted against the streamwise wavelength $\lambda _x=2{\rm \pi} /k_x$.

In figure 13(a,b) the pressure–strain correlation $\hat {\varPi }^{(u)}$ is mostly negative, and large negative $\hat {\varPi }^{(u)}$ appears for the range of $\lambda _x$ with large positive $\hat {P}^{(u)}$. In figure 13(c–f), $\hat {\varPi }^{(v)}$ and $\hat {\varPi }^{(w)}$ tend to be positive for the wavelength with negative $\hat {\varPi }^{(u)}$. The energy of $\bar {u}'$, increased by the production, is partially redistributed to those of $\bar {v}'$ and $\bar {w}'$ at large scales close to $\lambda _x\approx \delta _u$. In figure 13(c), $\hat {\varPi }^{(v)}$ is negative for $\lambda _x/\delta _u\gtrsim 2$, for which $\hat {\varPi }^{(w)}$ in figure 13(e) is positive, and the energy is redistributed from $\bar {v}'$ mainly to $\bar {w}'$. The wavelength of this energy redistribution is larger than the peak wavelength of the positive energy production for $E_u$. Interestingly, similar behaviour of the pressure–strain correlation was reported in channel flows (Lee & Moser Reference Lee and Moser2019): the energy of the streamwise velocity is redistributed to the vertical (wall-normal) and spanwise velocity components at large scales, while for even larger scales, the energy of the spanwise velocity is redistributed to that of the vertical velocity.

The dissipation terms $\varepsilon ^{(u_\alpha )}=\varepsilon _{G}^{(u_\alpha )}+\varepsilon _{S}^{(u_\alpha )}$ have a negative peak at small scales as expected from the scale dependence of the energy dissipation in turbulence. The turbulent diffusion terms $\hat {T}_{T}^{(u_\alpha )}$ are relatively small compared with other important terms at $z=0$ in figure 13(a,c,e). However, $\hat {T}_{T}^{(u)}$ at $z/\delta _u=0.4$ is the main source of the energy growth at large scales in figure 13(b). Pham & Sarkar (Reference Pham and Sarkar2010) also found that the turbulent kinetic energy is transferred outward in the stably stratified shear layer. On the other hand, $\hat {T}_{T}^{(v)}$ and $\hat {T}_{T}^{(w)}$ are smaller than other important terms even at $z/\delta _u=0.4$, where the pressure–strain correlation is the source of the energy.

The interscale transport term $\hat {T}_{S}^{(u)}$ in figure 13(a,b) is positive for $\lambda _x/\delta _u \lesssim 0.8$ and is negative for larger scales. Similarly, positive and negative $\hat {T}_{S}^{(w)}$ in figure 13(e,f) can be found for small and large scales, respectively. Therefore, a forward energy transfer from large to small scales occurs for the energy of $\bar {u}'$ and $\bar {w}'$. However, $\hat {T}_{S}^{(v)}$ in figure 13(c,d) exhibits negative values for $0.2\lesssim \lambda _x/\delta _u \lesssim 2$, which are surrounded by positive values for larger and smaller scales. This profile indicates that the energy at scales of $0.2\lesssim \lambda _x/\delta _u \lesssim 2$ is transferred to both large and small scales. Positive $\hat {T}_{S}^{(v)}$ for $\lambda _x/\delta _u \gtrsim 2$ appears for the wavelength with negative $\hat {\varPi }^{(v)}$. Therefore, the interscale energy transfer of $\bar {v}'$ feeds the energy at large scales, which is further redistributed by the pressure–strain correlation. At $t=140$, the Ozmidov scale is $L_O/\delta _u=6.9\times 10^{-2}$ at the centre of the shear layer. Turbulent motions with scales greater than the Ozmidov scale are expected to be strongly suppressed by the buoyancy (Riley & Lindborg Reference Riley and Lindborg2008). However, $\hat {T}_{S}^{(u_\alpha )}$ confirms that the forward energy transfer occurs for length scales larger than the Ozmidov scale, indicating that the active interscale transport in the turbulent shear layer is possible even under the strong influence of stable stratification. It is interesting to note that the scale dependence of $\hat {T}_{S}^{(u_\alpha )}$ in the stably stratified shear layer is consistent with DNS results of channel flows (Lee & Moser Reference Lee and Moser2019). Above the buffer layer of the channel flow, an interscale transfer from large to small scales occurs for the energy of streamwise and vertical velocities while the energy of spanwise velocity is transferred to both large and small scales.

The pressure diffusion and buoyancy flux terms appear only in the equation for $E_w$ shown in figure 13(e,f). The pressure diffusion term for $\lambda _x/\delta _u\approx 10^0$ is positive at $z=0$ and negative at $z/\delta _u=0.4$, and transfers the energy toward the middle of the shear layer at large scales. The profile of the buoyancy flux is similar to the pressure–strain correlation with opposite signs, and the energy converted to $w'$ is partially removed by the buoyancy force.

Figure 13(g,h) shows the spectral energy budget of density fluctuations, (2.16), at $t=140$. Both streamwise velocity and density fluctuations are produced by the vertical gradients of $\langle u \rangle$ and $\langle \rho \rangle$, whose initial profiles are given by a hyperbolic tangent function. Therefore, the spectral energy budget is similar between $E_u$ and $E_\rho$ at $t=140$.

3.5.2. Spectral energy budget at late time

Figure 14(a–f) presents the spectral energy budget at $t=240$ for three velocity components. At $t=240$, $E_u$ and $E_\rho$ have a second peak associated with the ELSS while the long-wavelength peak is forming in $E_v$ and $E_w$ (figure 8). For $E_{u}$, the dissipation, turbulent transport and interscale transfer terms are qualitatively similar at $t=140$ and 240 while $\hat {P}^{(u)}$ and $\hat {\varPi }^{(u)}$ have different profiles. The production $\hat {P}^{(u)}$ at $t=240$ has a large peak at $\lambda _x/\delta _u\approx 0.6$ at both heights. At $z=0$, $\hat {P}^{(u)}$ has a second peak for $\lambda _x/\delta _u\approx 10$, which does not exist at $t=140$. The energy production for $\lambda _x/\delta _u\approx 10$ does not occur at $z/\delta _u=0.4$, where the ELSS do not appear. Negative $\hat {\varPi }^{(u)}$ appears for $0.1\lesssim \lambda _x/\delta _u \lesssim 1$, where the production term has large positive values. As also found at early time, the energy produced by the mean velocity gradient is redistributed to other velocity components for this wavelength range. However, this energy redistribution does not occur at $\lambda _x/\delta _u\approx 10$, where the energy production is active. Comparisons of all terms shown in figure 14(a) imply that the growth of the ELSS at $\lambda _x/\delta _x\approx 10^1$ is mainly caused by the energy production at very large scales.

Figure 14. Same as figure 13, but at $t=240$.

The most terms for the budget of $E_v$ are qualitatively similar between $t=140$ and 240. However, negative $\hat {\varPi }^{(v)}$ at $\lambda _x/\delta _u \gtrsim 5$ found at $t=140$ is not observed at $t=240$, and $\hat {\varPi }^{(v)}$ does not cause the decay of $E_v$ at very large scales. For $E_{w}$, time dependence is significant for $\hat {T}_{P}^{(w)}$ and $\hat {\varPi }^{(w)}$. The pressure diffusion $\hat {T}_{P}^{(w)}$ at $t=240$ is positive at $\lambda _x/\delta _u\approx 10$ while $\hat {\varPi }^{(w)}$ is close to 0 for $\lambda _x/\delta _u\gtrsim 2$ at $t=240$. Thus, the pressure effect contributes to the energy growth at very large scales at $z=0$.

Figure 14(g,h) shows the budget of $E_{\rho }$ at $t=240$. Difference between $t=140$ and 240 is similar for $E_u$ and $E_\rho$. A long-wavelength peak of the production term appears at $\lambda _x/\delta _u\approx 10$ at $t=240$ in figure 14(g), and the production at the very large scales generates the density fluctuations of ELSS.

One of the noteworthy features for the budgets of $E_{u}$ and $E_{\rho }$ is small negative production at small scales of $\lambda _{x}/\delta _{u}\approx 0.1$ in figure 14(a,b,g,h). Even at $t=140$, the production terms have small negative values for $\lambda _x/\delta _u\approx 0.1$ in figure 13(a). These negative values of the production terms are also found in DNS (not shown here) and are not artificially caused by the implicit SGS model. The negative production requires ${{Re}}[\langle \widehat {\bar {u}'}^{*}\widehat {\bar {w}'}\rangle ] > 0$ and ${{Re}}[\langle \hat {\bar {\rho }'}^{*}\widehat {\bar {w}'}\rangle ]<0$. These cospectra represent the spectral density of the turbulent fluxes. The signs of the cospectra corresponding to the negative production are associated with the counter-gradient transports of streamwise momentum and density. This is contrary to the spectral analysis of Reynolds stress of a fully developed turbulent mixing layer without density stratification, where the turbulent momentum transport is in the down-gradient direction for all length scales (Takamure et al. Reference Takamure, Ito, Sakai, Iwano and Hayase2018). Interestingly, negative values of the production term of the turbulent kinetic energy spectrum were also found in channel flows (Lee & Moser Reference Lee and Moser2019). Around vertical locations between the buffer layer and logarithmic layer of the channel flow, the production term has large positive values at large scales while negative production appears for slightly smaller wavelengths than the positive production. As in the stably stratified shear layer considered in this study, the negative production in the channel flow is small compared with their positive peak at large scales.

3.5.3. Formation of double-peak energy spectra on the shear layer centreline

From $t=140$ to $240$, the spectral shape at $z=0$ changes from a single-peak profile to a double-peak profile. The formation of the double-peak spectrum can be clearly appreciated by $T_\alpha ^{-1} \equiv (\partial E_{\alpha }/\partial t)/E_{\alpha }$, which can be interpreted as the rate (inverse of time scale) of the spectral evolution. Figure 15 shows $T_\alpha ^{-1}$ for three velocity components and density at $t=140$ and 240 at the centre of the shear layer. For $\lambda _x/\delta _u\lesssim 0.1$, $T_\alpha ^{-1}$ becomes large as $\lambda _x$ decreases. This is because the energy spectra become extremely small as the implicit SGS model dissipates the energy at resolved small scales. The energy spectra decay with time as also confirmed from $T_\alpha ^{-1}<0$ for most wavelengths. At $t=240$, the decay rate is small at $\lambda _x/\delta _u\approx 10$, which is the length scale of the ELSS. Therefore, the velocity fluctuations of the ELSS decay more slowly than those of the LSM. The scale dependence of $T_\alpha ^{-1}$ suggests that the spectrum with a single peak rapidly decays with time except for very large scales, and eventually the double-peak spectrum is formed once the spectrum decays sufficiently except for $\lambda _x/\delta _u\approx 10$. At $t=140$ in figure 15(a), the slow decay of energy spectra at $\lambda _x/\delta _u\approx 10$ is found only for $E_u$ and $E_\rho$. Thus, the decay of $E_v$ and $E_w$ is faster at larger scales. Consistently, the formation of the double-peak spectra occurs at earlier time for $E_u$ and $E_\rho$ than for $E_v$ and $E_w$ in figure 8.

Figure 15. The rate of the change in energy spectra defined as $T_{\alpha }^{-1}=(\partial E_{\alpha }/\partial t)/E_{\alpha }$ for $\alpha =u, v, w$ and $\rho$ at $z=0$ in run 2: (a) $t=140$; (b) $t=240$.

Figure 16(a) shows the production term $k_{x}\hat {P}^{(u)}$ normalized by its maximum value $[k_{x}\hat {P}^{(u)}]_{max}$ at $z=0$ at four time instances. At $t=140$, the energy is produced for a wide range of length scales. However, the shape of $k_{x}\hat {P}^{(u)}$ changes from $t=140$ to 230, and the production term at late time has two peaks corresponding to the peaks in $E_u$. The energy production at $\lambda _x/\delta _u\approx 10$ becomes more important at late time. The production is almost zero or can be negative for $\lambda _x/\delta _u\approx 6$. This negligible production is consistent with $E_u$ in figure 8(c), where $E_u$ is locally small for the wavelength with small $k_{x}\hat {P}^{(u)}$. The wavelength of positive production corresponding to the long-wavelength peak is shifted to larger scales with time as also expected from the increase of the long peak-wavelength in figure 10(b). Figure 16(b) shows the sum of the pressure terms, $\hat {P}_{r}^{(w)}=\hat {T}_{P}^{(w)}+\hat {\varPi }^{(w)}$, at $z=0$, where $\hat {P}_{r}^{(w)}$ is premultiplied by $k_x$ and normalized by the maximum value of $k_x\hat {P}_{r}^{(w)}$. The pressure terms of $E_w$ play a similar role to the production term of $E_u$: the energy growth at $\lambda _x/\delta _u\approx 0.5$ and 10 is caused by the pressure effects for $E_{w}$ and by the production for $E_{u}$. The pressure effects are responsible for the change in the shape of $E_w$. Similar time dependence was found for the production term of density (not shown here). For $E_v$, we cannot identify a single term that is responsible for the formation of the double-peak spectrum. The growth of $E_v$ at large $\lambda _x$ is dominated by the interscale energy transfer $\hat {T}_{S}^{(v)}$ in figure 14(c), where the energy transfer from small to large scales occurs. The decrease of $E_v$ by the pressure–strain correlation for $\lambda _x/\delta _u\approx 10$ occurs at $t=140$ but not at $t=240$ as shown in figures 13(c) and 14(c). The balance of these terms causes the slow decay rate of very large-scale fluctuations in figure 15(b).

Figure 16. Temporal variations of (a) production term $\hat {P}^{(u)}$ of $E_u$ and (b) pressure term $\hat {P}_{r}^{(w)}=\hat {\varPi }^{(w)}+\hat {T}_{P}^{(w)}$ of $E_w$ on the centreline ($z=0$) in run 2. The terms $k_{x}\hat {P}^{(u)}$ and $k_{x}\hat {P}_{r}^{(w)}$ are normalized by their maximum values, $[k_{x}\hat {P}^{(u)}]_{max}$ and $[k_{x}\hat {P}_{r}^{(w)}]_{max}$.

3.6. Scaling of energy spectrum and structure function

Another interesting feature in the energy spectra in figure 8 is that the premultiplied spectrum at a fixed height is roughly constant for a certain wavenumber range at $t=320$, implying a $k_x^{-1}$ power law. For example, $k_xE_{u}$ in figure 8(c) weakly depends on the wavenumber for $2\lesssim \lambda _x/\delta _u\lesssim 4$ near the centreline, while the range with constant $k_xE_{u}$ extends to smaller scales at $z/\delta _u\approx 0.25$, where large peaks at $\lambda _x/\delta _u\approx 1$ and 10 do not exist. As discussed above, some terms of the spectral energy equations have similar wavenumber dependence in the stably stratified shear layer and wall-bounded shear flows. Furthermore, turbulent structures that have been observed in the wall-bounded shear flows were also observed in the stably stratified shear layer, e.g.  hairpin vortices and ELSS (Watanabe et al. Reference Watanabe, Riley, Nagata, Matsuda and Onishi2019a). Both stable stratification and a wall cause the suppression of vertical turbulent motions, which may lead similar turbulent characteristics in these flows. Interestingly, the scaling $E_{u}\sim k_x^{-1}$ has also been predicted and experimentally confirmed in wall turbulence (Perry & Chong Reference Perry and Chong1982; Nickels et al. Reference Nickels, Marusic, Hafez and Chong2005). The equivalent scaling in real space can be written with the velocity structure function, which is defined as $\langle (\Delta u)^2\rangle =\langle [\bar {u}(x+r,y,z)-\bar {u}(x,y,z)]^2\rangle$ with a separation distance $r$ in the $x$ direction. The scaling law of $\langle (\Delta u)^2\rangle$ corresponding to $E_{u}\sim k_{x}^{-1}$ is the logarithmic law given by $\langle (\Delta u)^2\rangle \sim \ln (r)$ (Davidson, Nickels & Krogstad Reference Davidson, Nickels and Krogstad2006). The physical basis underlying these scaling laws is the existence of turbulent motions with a range of scales which are characterized by a common velocity scale. Davidson et al. (Reference Davidson, Nickels and Krogstad2006) also pointed out that the logarithmic law is more easily observed than the $k_x^{-1}$ law because of the aliasing problem of one-dimensional spectra.

Figure 17(a–c) shows the structure functions of $u$, $v$ and $w$ at $t=320$ at three vertical locations. Here, the structure functions of $v$ and $w$ are $\langle (\Delta v)^2\rangle =\langle [\bar {v}(x+r,y,z)-\bar {v}(x,y,z)]^2\rangle$ and $\langle (\Delta w)^2\rangle =\langle [\bar {w}(x+r,y,z)-\bar {w}(x,y,z)]^2\rangle$. Except near the edge of the shear layer ($z/\delta _u=0.46$), the logarithmic law is valid for a certain range of large scales. The broken lines in (a) represent $\langle (\Delta u)^2\rangle = A + B\ln (r)$ with $A$ and $B$ obtained by the least square method, which is applied for $1.2 \leqslant r/\delta _u \leqslant 2.8$ at $z=0$ and $0.4\leqslant r/\delta _u \leqslant 2.8$ at $z/\delta _u=0.23$. In figure 17(b,c) the logarithmic law is also shown with the broken lines, which are obtained with the least square method applied for $0.97 \leqslant r/\delta _u \leqslant 2.4$ at $z=0$ and $0.81 \leqslant r/\delta _u \leqslant 1.6$ at $z/\delta _u=0.23$. The logarithmic law holds for all velocity components. However, $A$ and $B$ are not universal and are different for each velocity component and vertical location. Figure 17(d) shows the premultiplied spectrum $k_{x}E_{u}$ at the same locations. Although there is a range of $k_x$ where $k_{x}E_{u}$ weakly depends on $k_x$, this wavenumber range of constant $k_{x}E_{u}$ at $z=0$ is not as clear as the $\ln (r)$ law in figure 17(a), as also expected from the comparison between the $k_{x}^{-1}$ and logarithmic laws in turbulent boundary layers (Davidson et al. Reference Davidson, Nickels and Krogstad2006). Near the centre of the shear layer, $E_{\alpha }\sim k_{x}^{-1}$ tends to appear for the scales between two peaks in figure 8(c,f,i). In figure 8(c), $E_u$ does not have large peaks at $z/\delta _u\approx 0.25$, where the logarithmic law is valid for a wide range of scales in figure 17(a). Therefore, the logarithmic law is clearly observed for large scales when the scales are not affected by the spectral peaks of LSM and ELSS. As the scale separation between LSM and ELSS increases with time, the $\ln (r)$ law and $k^{-1}$ law may be more clearly observed at further later time.

Figure 17. Structure functions of (a) streamwise velocity $\langle (\Delta u)^2\rangle$, (b) spanwise velocity $\langle (\Delta v)^2\rangle$ and (c) vertical velocity $\langle (\Delta w)^2\rangle$ at $t=320$. The broken lines represent $A \ln (r/\delta _u)+B$ with $A$ and $B$ obtained with the least square method applied for $z=0$ and $0.23\delta _u$. (d) Premultiplied energy spectrum of streamwise velocity, $k_xE_u$, at $t=320$. These results are taken at $z/\delta _u=0$, $0.23$ and $0.46$ in run 2.

3.7. Comparison between the stably stratified shear layer and wall turbulence

The ELSS develop at late time in the stably stratified shear layer. Similar structures, known as superstructures or very LSM, also exist in wall-bounded turbulent shear flows, e.g.  channel, pipe and boundary-layer flows (Hutchins & Marusic Reference Hutchins and Marusic2007; Monty et al. Reference Monty, Hutchins, Ng, Marusic and Chong2009). It was also shown that a large number of hairpin vortices exist in the stably stratified shear layer, where the hairpin vortices induce velocity fluctuations at the streamwise length scale of LSM (Watanabe et al. Reference Watanabe, Riley, Nagata, Matsuda and Onishi2019a). Thus, the stably stratified shear layer at late time resembles the wall turbulence in terms of the turbulent structures. Furthermore, the LES results presented above also confirm that some statistics are also qualitatively similar in these flows, e.g.  the energy budget in wavenumber space, the anisotropy of Reynolds stress tensor and the logarithmic law of the velocity structure functions at large scales.

As confirmed from the energy spectra, the ELSS emerge at late time for ${Re}=2000$ and 40 000 when $Ri$ is not too small. The detailed process of turbulent transition depends on the Reynolds number and Richardson number (Mashayek & Peltier Reference Mashayek and Peltier2012; Thorpe Reference Thorpe2012). Therefore, the transition process is not directly related to the formation of ELSS. This is also supported by the flow visualization in figure 4: the velocity fluctuations of the ELSS are hardly correlated with the velocity profiles associated with the Kelvin–Helmholtz instability and the initial velocity fluctuations. The spectral energy budget in § 3.5 also confirms that the flow dynamics is similar between the stably stratified shear layer and channel flow. Therefore, it is reasonable to assume that the ELSS in these flows are generated by the same mechanism. However, there is also a possibility that the ELSS is generated and sustained by long internal gravity waves, which may dominate the velocity fluctuations at the late time. Further investigation is required to clarify the generation mechanism of the ELSS. The ELSS develop even if ${Re}_b$ is larger than 100 (${Re}=40\,000$, $Ri=0.06$). The present LES suggests that $Ri$ or $Ri_g$ is an important parameter in the development of the ELSS because the ELSS are not observed for $({Re},Ri)=(40\,000,0.02)$, for which $Ri_g$ in table 2 is below the critical Richardson number 1/4 (Miles Reference Miles1961). On the other hand, the ELSS are observed for $Ri=0.06$ and 0.1, for which $Ri_g$ exceeds the critical value. At $Ri=0$ and 0.02, $\delta _u$ grows with time even at late time in figure 7(a). However, $\delta _u$ hardly increases with time for $Ri=0.06$ and 0.1 after the turbulent transition, and the vertical length scale of the flow does not vary with time. In wall turbulence the largest vertical length scale of turbulence at a given height is also fixed as the distance from the wall. The confinement of the growth of the vertical length scale is possibly important in the formation of the ELSS in these flows.

Negative production of $E_u$ and $E_\rho$ at small scales is observed in figure 14(a,b,g,h). As explained above, the negative production is due to the counter-gradient transport of streamwise momentum and density. Previous studies on stratified turbulence also observed the counter-gradient transport regardless of the presence of mean shear (Lienhard & van Atta Reference Lienhard and van Atta1990; Holt, Koseff & Ferziger Reference Holt, Koseff and Ferziger1992; Komori & Nagata Reference Komori and Nagata1996). However, the physical mechanism of the counter-gradient transport in stratified turbulence may depend on the mean shear as discussed in Gerz & Schumann (Reference Gerz and Schumann1996). They explained the counter-gradient transport at small scales in stratified turbulence with mean shear based on the collision of fluid parcels which are vertically transferred by inclined hairpin vortices generated by the mean shear. Indeed, a large number of hairpin vortices develop in the stably stratified shear layer (Watanabe et al. Reference Watanabe, Riley, Nagata, Matsuda and Onishi2019a). Gerz & Schumann (Reference Gerz and Schumann1996) also suggested that the same mechanism can also cause the counter-gradient transport of momentum in wall-bounded shear flows, where the hairpin vortices have also been found in numerous studies (e.g. Wu & Moin Reference Wu and Moin2009). The wavenumber dependence of the production terms of $E_u$ and $E_\rho$ is also qualitatively similar between the middle of the shear layer (figure 14a) and the buffer and logarithmic layers of channel flows (Lee & Moser Reference Lee and Moser2019). Thus, the mechanism of the counter-gradient momentum transport proposed by Gerz & Schumann (Reference Gerz and Schumann1996) seems to be valid in these flows.

The integral shear parameter $S^{*}=Sq^2/\varepsilon ^{(k)}$ ($q^2=\langle u'_i u'_i\rangle$) represents the time scale ratio of the eddy turnover time $q^2/\varepsilon ^{(k)}$ to the mean-shear time scale $S^{-1}=\langle \partial u / \partial z \rangle ^{-1}$. The evolution of turbulence is strongly influenced by the mean shear when $S^{*}$ is much larger than $1$. Non-stratified turbulent free shear flows have $S^{*}\approx 3$–$6$ (Pope Reference Pope2000) while the stably stratified turbulent shear layer at $z=0$ has $S^{*}\geqslant 10$, which continuously increases with time. At $({Re},Ri)=(2000,0.06)$, $S^{*}$ reaches about 30 at $t=320$ (Watanabe et al. Reference Watanabe, Riley, Nagata, Matsuda and Onishi2019a) and $S^{*}\approx 30$ is as large as in the buffer layer of channel flows (Jiménez Reference Jiménez2013). The counter-gradient transport by the mechanism proposed in Gerz & Schumann (Reference Gerz and Schumann1996) is simply due to the strong influence of mean shear and is not directly related to the stable stratification. Therefore, the counter-gradient momentum transport at small scales occurs in homogeneous shear flow with $S^{*}\approx 30$ even without stable stratification (Iida & Nagano Reference Iida and Nagano2007), where the shape of the cospectrum of the vertical momentum-flux is qualitatively similar to $\hat {P}^{(u)}$ in figure 16(a) except the scales of the ELSS. The mean shear $S$ in a temporally evolving, turbulent shear layer decreases with time in the absence of stable stratification because the shear layer thickness increases. Therefore, although the eddy turnover time $q^2/\varepsilon ^{(k)}$ becomes large with time, $S^{*}$ stays small and the counter-gradient transport does not occur without stable stratification (Takamure et al. Reference Takamure, Ito, Sakai, Iwano and Hayase2018). The stable stratification suppresses the growth of the shear layer thickness. Thus, $S\sim U_0/\delta _u$ hardly varies with time while $S^{*}$ increases with time because of the increase of the eddy turnover time. Therefore, the stably stratified shear layer at late time has large $S^{*}$, and the strong influence of the mean shear can result in the counter-gradient transport at small scales.