2D and 3D Properties of Stably Stratified Turbulence

Abstract

Through interactions between the modes of “waves” and “vortices”, stably stratified turbulence exhibits characteristic features both in 2D and 3D. Using DNS of the Navier–Stokes equations coupled to an equation for temperature (the Bousinessq system), we investigate dynamical properties of stably stratified turbulence. After reviewing the importance of three characteristic length scales and their relations in stably stratified turbulence, we present numerical results showing how structures and spectra develop with the growth of the length scales. It is demonstrated that the temperature fluctuations make sharp fronts vertically which result in non-symmetric PDFs of the vertical derivative ${\partial }_z\theta$.

1. Introduction

The first time YK met Jack was at the IUTAM Symposium on “Turbulence and Chaotic Phenomena in Fluids” which professor T. Tatsumi organized in Kyoto in 1983. That symposium was really a big event, and YK also meet many researchers whom he would re-encounter later in the US and elsewhere, for example, Annick Pouquet also attended the symposium. The next time he crossed paths with Jack was at the APS meeting at SUNY, Buffalo in 1988 where Jack organized a specialty conference on Multipoint Turbulent Closure. The speakers of the conference were Bill George, Marcel Lesieur, Bill Dannevik, Chuck Leith, Akira Yoshizawa, Jean-Pierre Bertoglio, and Jack. At that time YK had just begun his career as a postdoc at UCSD and subsequently at Los Alamos National Laboratory (LANL). At Los Alamos, he was affiliated with the Center for Nonlinear Studies (CNLS), and worked with Bob Kraichnan. (His official host at CNLS was Darryl Holm since Bob Kraichnan was by then a consultant at LANL). After Los Alamos, YK was awarded an Advanced Study Program (ASP) postdoc scholarship at NCAR and came to Boulder in 1991, and started full time collaborations with Jack.

The first project we conducted was “diffusion in stably stratified turbulence”. Using decaying turbulence simulations without external forcing in a $2\pi$-periodic box, we showed that randomly scattered pancake vortices are characteristic vortical structures in homogeneous turbulence [1]. This paper also considered the effect of stratification on eddy diffusion, which was stimulated by previous papers [2,3,4,5,6]. Then, we included a narrow band forcing term in the momentum equation in Fourier space to produce stationary turbulence, and investigated the power transition of the horizontal kinetic energy spectra in stably stratified turbulence [7]. In that paper, we confirmed the atmospheric observations of Nastrom and Gage (1985) that the horizontal spectrum shows a transition from $k_{\perp }^{-3}$ at the large scales to $k_{\perp }^{-5/3}$ at small scales where $k_{\perp }$ is the horizontal wave number [8].

Adding a narrow band forcing to the momentum equation is a simple but an artificial method to produce stationary turbulence, and we need to have a specific experimental design to choose a proper type of forcing. So far, various kinds of forcing are used for different purposes by many researchers [9,10,11,12,13,14]. The purpose of the present paper is to investigate how two and three-dimensional properties of stably stratified turbulence are developed by the interaction of “waves” and “vortices” in the turbulence. For this purpose, we first decouple the temperature equation from the momentum equation, and develop two-dimensional turbulence in the $x,y$ coordinates which is uniform in the z coordinate in a 3D periodic cubic box. After two-dimensional turbulence is fully developed, the temperature equation is turned on, and we observe how the two-dimensional turbulence adjusts to the structures in stably stratified turbulence.

Here we explain the system of equations of the velocity and the temperature fluctuations and the method of the simulation [7]. Under the Boussinesq approximation, the equation of motion of stably stratified flow for the velocity vector, $\mathbf{u}=(u,v,w)$ and the temperature fluctuation from the linear (stable) mean, $\theta$ are written

$$\begin{array}{cc}\hfill \frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}·\mathbf{\nabla })\mathbf{u}& =-\mathbf{\nabla }p+\nu {\nabla }^2\mathbf{u}+\theta \widehat{\mathbf{z}}+\mathbf{f}\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \frac{\partial \theta }{\partial t}+(\mathbf{u}·\mathbf{\nabla })\theta & =\kappa {\nabla }^2\theta -N^{2}w\hfill \end{array}$$

(2)

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

(3)

where $p$ is pressure, $\widehat{\mathbf{z}}$ is the unit vector in the vertical direction, $\mathbf{f}$ is an external force, and $\nu$ and $\kappa$ are kinematic viscosity and temperature diffusivity, respectively. N is the Brunt–Väisälä frequency, defined as $N=\sqrt{g\alpha (\partial \overline{T}/\partial z)/T_0}$, where $\alpha$ is the coefficient of thermal expansion. Here we assume $\overline{T}\left(z\right)$ satisfies $d\overline{T}/dz\equiv -N^2$. Periodic boundary conditions are used in the calculation of $\mathbf{u}$ and $\theta$ by assuming a linear ambient density profile. We set the Prandtl number $\nu /\kappa =1$. For the time advancement, the low-storage third-order Runge–Kutta method is used [15]. The simulation code is decomposed into horizontal slabs and parallelized using the Message Passing Interface (MPI) system.

2. Vortical Structures in Stratified Turbulence

In stably stratified flows, waves produced by buoyancy and vortices coexist with turbulence, and their interaction produces many length scales. In general, at large scales the effect of buoyancy is assumed to be dominant, and by regarding the Brunt–Väisälä frequency as a frequency of stable oscillation, we define the buoyancy scale as $L_B(=U/N)$, where U is a characteristic horizontal scale.

At smaller scales, the turbulence characteristics become clearer. If we focus on the energy cascade, the energy dissipation rate $\epsilon$ is given by the velocity and the length scales, u and l, of the Kolmogorov inertial range as $\epsilon \sim u^3/l$. If $l\sim L_B$ and $u\sim U$ are substituted in the dissipation expression, and U is eliminated, we obtain the Ozmidov scale $L_O(=\sqrt{\epsilon /N^3})$ [16,17]. The Ozmidov scale is at the boundary between the buoyancy and the inertial range; scales smaller than $L_O$ are unstable and tend to overturn. At even smaller scales, the Kolmogorov scale $L_K(={({\nu }^3/\epsilon )}^{1/4})$ appears, and eddies smaller than $L_K$ dissipate the kinetic energy by viscosity.

Thus we have three scales, $L_B,L_O,L_K$, as the characteristic length scales for stably stratified turbulence, and combinations of these describe different features of turbulence. The ratio of $L_O$ and $L_K$ is

$$\frac{L_O}{L_K}=\frac{{(\epsilon /N^3)}^{1/2}}{{\nu }^3{/\epsilon )}^{1/4}}=R_B^{3/4}$$

(4)

where $R_B\left(=\epsilon /\left(\nu N^2\right)\right)$ is called the buoyancy Reynolds number [7,14,18]. When $R_B\gg 1$, the effect of viscosity is confined to scales much smaller than $L_O$, with the scale separation expected at the boundary $L_O$.

Figure 1 is the celebrated mesoscale horizontal spectra of Nastrom and Gage collected from aircraft flights [8] of the zonal (east–west) wind, the meridional (north–south) wind, and the potential temperature. In this plot, the abscissa denotes the magnitude of the horizontal wave number $k_{\perp }=\sqrt{x_x^2+k_y^2}$.

Each spectrum has the functional form of $k_{\perp }^{\alpha }$, and the power $\alpha$ is between $-2$ and $-3$ at large scales and is close to $-5/3$ at small scales. Historically, this $-5/3$ range was first suspected to be the result of a 2D inverse energy cascade, but later was corroborated by analysis of the 3rd order horizontal velocity structure functions; the direction of energy cascade is regular, i.e., large to small scales [19,20]. The boundary wavenumber between the different powers corresponds to $k_O=2\pi /L_O$, and for small scales below $L_O$ the energy cascade is enhanced by eddy overturning.

When $R_B\gg 1$, since $L_O$ is the smallest scale in the buoyancy dominated region, it can be used as a metric to characterize other length scales. If we assume that the horizontal length scale $L_h$ is associated with the turbulent energy cascade process it can be expressed in terms of the characteristic horizontal velocity U and the rate of energy dissipation $\epsilon$ as $L_h\sim U^3/\epsilon$ [21]. We further assume that the vertical length scale $L_v$ corresponds to the buoyancy scale $L_B=U/N$. Then, by taking the ratio with $L_O$, and by expressing the results by the horizontal Froude number, $F{r}_h\equiv U/\left(N{L}_h\right)=\epsilon /\left(N{U}^2\right)$, we have the following relations:

$$\frac{L_h}{L_O}=F{r}_h^{-3/2},\frac{L_v}{L_O}=F{r}_h^{-1/2},$$

(5)

which indicates that $L_h/L_O\sim {(L_v/L_O)}^3$ holds in the buoyancy dominated region [14].

It was shown that randomly scattered horizontal pancake vortices are the characteristic vortical structures in stably stratified turbulence [7]. Based on the above relations, the ratio of horizontal and vertical length scales of the vortices is given as $L_h/L_O\sim {(L_v/L_O)}^3$, $L_h/L_v\sim F{r}_h^{-1}$.

3. Development of 2dim and 3dim Energy and Spectrum

3.1. Craya–Herring Decomposition

The basic idea of the Craya–Herring decomposition is the following [22]. The Fourier transform of the velocity incompressibility condition provides $\mathbf{k}·\tilde{\mathbf{u}}\left(\mathbf{k}\right)=0$, where $\mathbf{k}={(k_x,k_y,k_z)}^T$ is the wavenumber vector, directly shows that the Fourier transform of the velocity vector, $\tilde{\mathbf{u}}\left(\mathbf{k}\right)$ spans two independent basis vectors ${\mathbf{e}}_1\left(\mathbf{k}\right)$, ${\mathbf{e}}_2\left(\mathbf{k}\right)$ in the plane perpendicular to $\mathbf{k}$ as $\tilde{\mathbf{u}}\left(\mathbf{k}\right)={\varphi }_1{\mathbf{e}}_1\left(\mathbf{k}\right)+{\varphi }_2{\mathbf{e}}_2\left(\mathbf{k}\right)$. As the orthonomal basis vectors, the following expressions follow:

$$\begin{array}{cc}\hfill {\mathbf{e}}_1\left(\mathbf{k}\right)& =\frac{\mathbf{k}\times \widehat{\mathbf{z}}}{\left|\right|\mathbf{k}\times \widehat{\mathbf{z}}\left|\right|}=\frac{1}{\sqrt{k_x^2+k_y^2}}\left(\begin{array}{c}{k}_y\\ -k_x\\ 0\end{array}\right),\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {\mathbf{e}}_2\left(\mathbf{k}\right)& =\frac{\mathbf{k}\times (\mathbf{k}\times \widehat{\mathbf{z}})}{\left|\right|\mathbf{k}\times (\mathbf{k}\times \widehat{\mathbf{z}})\left|\right|}=\frac{1}{\sqrt{k_x^2+k_y^2+k_z^2}\sqrt{k_x^2+k_y^2}}\left(\begin{array}{c}{k}_z{k}_x\\ k_z{k}_x\\ -(k_x^2+k_y^2)\end{array}\right),\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\mathbf{e}}_3\left(\mathbf{k}\right)& =\frac{\mathbf{k}}{\left|\right|\mathbf{k}\left|\right|}=\frac{1}{\sqrt{k_x^2+k_y^2+k_z^2}}\left(\begin{array}{c}{k}_x\\ k_y\\ k_z\end{array}\right).\hfill \end{array}$$

(8)

Based on these expressions, ${\varphi }_1$ and ${\varphi }_2$ can be calculated explicitly as

$$\begin{array}{cc}\hfill {\varphi }_1& =\tilde{\mathbf{u}}·{\mathbf{e}}_1=\frac{1}{\sqrt{k_x^2+k_y^2}}(k_y\tilde{u}-k_x\tilde{v})=\frac{i}{\sqrt{k_x^2+k_y^2}}{\tilde{\omega }}_z\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\varphi }_2& =\tilde{\mathbf{u}}·{\mathbf{e}}_2=-\frac{\sqrt{k_x^2+k_y^2+k_z^2}}{\sqrt{k_x^2+k_y^2}}\tilde{w},\hfill \end{array}$$

(10)

where ${\tilde{\omega }}_z$ is the Fourier transform of the z-component of vorticity. These expressions define the energy densities ${\Phi }_1\left(\mathbf{k}\right)$ and ${\Phi }_2\left(\mathbf{k}\right)$ respectively as

$${\Phi }_1\left(\mathbf{k}\right)=\frac{1}{2}|{\varphi }_1{|}^2=\frac{1}{2}\frac{|k_y\tilde{u}-k_x\tilde{v}{|}^2}{k_x^2+k_y^2},{\Phi }_2\left(\mathbf{k}\right)=\frac{1}{2}{\left|{\varphi }_2\right|}^2=\frac{1}{2}\frac{(k_x^2+k_y^2+k_z^2){\left|\tilde{w}\right|}^2}{k_x^2+k_y^2}.$$

(11)

From the above equations; we see that ${\varphi }_1$ is related to the horizontal vorticity while ${\varphi }_2$ is related to the vertical velocity. Notice that the vertical velocity is coupled with the temperature fluctuation in an oscillatory manner, such as waves, we see that ${\varphi }_1$ is horizontal and vortical and ${\varphi }_2$ is vertical and wavy [7].

3.2. Numerical Results

We Fourier transform Equations (1)–(3) and integrate them pseudo-spectrally in a 2$\pi$ periodic box in the $x,y,z$ dimensions. Here, the results with ${1024}^3$ grid points are presented. As the external force, we add an integrated white-noise to the horizontal velocity components, $\tilde{u}(k_x,k_y,k_z)$ and $\tilde{v}(k_x,k_y,k_z)$ in a narrow waveband at the horizontal wave number $k_f=\sqrt{k_x^2+k_y^2}=4,10$.

$$\tilde{\mathbf{f}}={(f_x(k_x,k_y,0),f_y(k_x,k_y,0),0)}^T$$

(12)

Starting with zero kinetic energy, we first decouple the equation for the temperature fluctuation from the momentum equations for some time from $t=0$, thus allowing the forcing to build two-dimensional turbulence uniformly in all horizontal layers in the z direction; ${\varphi }_2=0$ is imposed in this first stage. After two-dimensional turbulence is developed, the temperature equation is turned on to see how the energy in ${\varphi }_2$ grows.

Figure 2 shows the time development of the energy in ${\varphi }_1$ and ${\varphi }_2$, ${\displaystyle \frac{1}{2}}{\left|{\varphi }_1\right|}^2$, and ${\displaystyle \frac{1}{2}}{\left|{\varphi }_2\right|}^2$, for four different cases, $(N^2,k_f)=(10,4),(100,4),(10,10),(100,10)$, respectively. For all cases, the development of the features are summarized as follows:

• The ${\varphi }_1$ energy grows rapidly in a fraction of time after $t=0$, and continues to increase gradually. Then the equation for the temperature fluctuations is turned on. After a short period of time, the ${\varphi }_2$ energy begins to grow dramatically at an exponential growth, $e^{\alpha t}$. The growth rate $\alpha$ is larger for $k_f=10$ than $k_f=4$. As the energy in ${\varphi }_2$ grows comparable to ${\varphi }_1$, they begin to interact and the speed and the growth of ${\varphi }_2$ decreases (and slightly that of ${\varphi }_1$ as well).
• As the interaction of ${\varphi }_1$ and ${\varphi }_2$ grows, and around the time when ${\varphi }_2$ energy becomes almost 10% of that of ${\varphi }_1$, the ${\varphi }_1$ energy reaches its maximum value, and then begins to decrease. Slightly after that time, ${\varphi }_2$ reaches its maximum value, followed by a gradual decrease in ${\varphi }_1$ and ${\varphi }_2$.
• At long times, the interaction between ${\varphi }_1$ and ${\varphi }_2$ becomes almost stationary and produces small oscillations. The fluctuations in the forcing seems to be the cause of these oscillations. ${\varphi }_1$ and ${\varphi }_2$ oscillate almost in phase, which is presumably due to the incompressibility that the Boussinesq approximation assumes.

Next, we examine the development of statistical quantities in Figure 2d; the case $N^2=100,k_f=10$ is an example. Figure 3 spans three time slots marked by orange solid lines, (i) $(0\le t\le 22.79)$, (ii) $(22.84\le t\le 45.04)$, and (iii) $(45.12\le t\le 76.43)$, corresponding to the above three characteristic features (1)∼(3). Also in Figure 3, the dashed lines show the times $(t_1\sim t_5)$ at which the flow features are analyzed in the later analysis.

Figure 4i shows the development of the ${\varphi }_1\left(k_{\perp }\right)$ spectra during the (i) time slot in Figure 3. This time interval corresponds to the process of forced two-dimensional turbulence with $-5/3$ spectra at scales smaller than the forcing wavenumbers; this corresponds to an inverse cascade of energy. At larger scales, spectra with $-3\sim -4$ slopes appear as typical spectra indicating a forward cascade of kinetic energy.

After the equation of the temperature fluctuations is turned on and ${\varphi }_2$ grows rapidly, the large scale part (small wavenumber part) of the ${\varphi }_1$ spectra tends to be flat as the small scale part (large wavenumber part) is elevated to be shallow spectra with an exponent approaching $-1$ (Figure 4ii). It seems obvious that the large scale energy is transferred significantly because of the interaction between ${\varphi }_1$ and ${\varphi }_2$.

When the interaction between ${\varphi }_1$ and ${\varphi }_2$ becomes almost stationary, the ${\varphi }_1$ spectra shows small changes, but tends to show a clear $-5/3$ slope at the high $k_{\perp }$ regions; at the same time the low $k_{\perp }$ spectra become flatter (Figure 4iii).

To understand the change in behavior in the ${\varphi }_1$ spectra, Figure 5 shows the change in the buoyancy, the Ozmidov and the Kolmogorov scales, explained in the previous section. A noteworthy phenomenon is that the Ozmidov, and the Kolmogorov scales cross around $t=23$ at which the buoyancy scale reaches its maximum. In the previous figures, this time corresponds to the time when the 2D state begins to change to 3D by the interaction between waves and vortices. Another interesting point is that the development of the buoyancy scale follows a similar path as ${\varphi }_1$ in Figure 2 and Figure 3. Throughout the simulation, the buoyancy scale is the largest scale in the flow. The similarity between the buoyancy scale and ${\varphi }_1$ implies that the horizontal flow continues to dominant, but we have to be a little careful to reach this conclusion—an explanation follows.

Figure 6 shows the development of the horizontal average of the variances, $\langle u^2\rangle /2$, $\langle v^2\rangle /2$, $\langle w^2\rangle /2$, and the total energy $(\langle u^2\rangle +\langle v^2\rangle +\langle w^2\rangle )/2$ for varying z. Results are shown at $t_1\sim t_4$ in Figure 3 (upper row), and the horizontal average of the horizontal energy, $(\langle u^2\rangle +\langle v^2\rangle )/2$, and monitored also at $t_1\sim t_4$ are the horizontal average of the square of the vertical derivative of temperature fluctuations, $\langle {(\partial \theta /\partial z)}^2\rangle$ (lower row). The following observations can be made:

(1)

At $t_1$, the energy in all directions is uniform in z. The vertical gradient of temperature fluctuations $\langle {(\partial \theta /\partial z)}^2\rangle$ shows small oscillation due to the initial small growth of ${\varphi }_2$. The total flow is mainly two-dimensional at this time.

(2)

At $t_2$, fluctuations in the horizontal energy develop, but the vertical energy is still very close to zero. The oscillation of the energy is not uniform, and small deviations in frequency and amplitude are seen. However, the structures in the x and y directions are similar. More striking is the vertical gradient of the temperature fluctuations that shows a few periodic humps with the tops showing stronger oscillations.

(3)

At $t=t_3$, the energy in the x and y directions, $\langle u^2\rangle /2,\langle v^2\rangle /2$, are in a similar range of magnitude with sharp peaks often at different locations in z. This implies that the velocity field is directionally intermittent, with the kinetic energy sometimes highly concentrated in one direction. Without looking at the coherence of the energy, we cannot readily identify the flow structure, but one possibility is there may at times be uni-directional wind flows.

(4)

At $t=t_4$, the intermittency in both the energy and the temperature derivative distributions grows. In particular, the peaks of the latter become sharply outstanding from the signals of low but rapid oscillations.

Figure 6. Upper row: The horizontal average of the energy, $\langle u^2\rangle /2$ (violet), $\langle v^2\rangle /2$ (green), $\langle w^2\rangle /2$ (blue) and the total energy $(\langle u^2\rangle +\langle v^2\rangle +\langle w^2\rangle )/2$ (orange) as a function of z. Lower row: The horizontal average of the horizontal energy, $(\langle u^2\rangle +\langle v^2\rangle )/2$ (violet) and the horizontal average of the square of the vertical derivative of temperature fluctuations, $\langle {(\partial \theta /\partial z)}^2\rangle$ (green). All statistics are monitored at $t_1\sim t_4$ in Figure 3.

Figure 6. Upper row: The horizontal average of the energy, $\langle u^2\rangle /2$ (violet), $\langle v^2\rangle /2$ (green), $\langle w^2\rangle /2$ (blue) and the total energy $(\langle u^2\rangle +\langle v^2\rangle +\langle w^2\rangle )/2$ (orange) as a function of z. Lower row: The horizontal average of the horizontal energy, $(\langle u^2\rangle +\langle v^2\rangle )/2$ (violet) and the horizontal average of the square of the vertical derivative of temperature fluctuations, $\langle {(\partial \theta /\partial z)}^2\rangle$ (green). All statistics are monitored at $t_1\sim t_4$ in Figure 3.

The above monitoring times, $t_1\sim t_4$, are all in time slot (ii) in Figure 3 and Figure 4. The characteristic strong intermittent growth during this time slot is attributed to the phenomenon that the 2D energy in the large scales is destroyed with the energy transferred into small scales producing spectra with a nearly $-1$ slope, much shallower than the Kolmogorov’s $-5/3$ spectrum.

The dramatic growth in intermittency in time interval (ii) in Figure 3 and Figure 4 is moderated in time interval (iii). In latter time range, the change in the large scale spectrum and the subsequent energy outflow towards small scales weakens (Figure 4iii), and the small scale spectrum has fallen to the Kolmogorov $-5/3$ spectrum. The end part of the spectra at the smallest scales has become steeper than the $-5/3$ power, showing a hint of dissipation spectrum. This observation has support if we look at the growth of the Kolmogorov scale at late time in Figure 5. It is interesting that the Kolmogorov scale increases in a stepped manner at $t\sim 56$ and $t\sim 67$.

The horizontal average of the horizontal energy, $(\langle u^2\rangle +\langle v^2\rangle )/2$, and the horizontal average of the square of the vertical derivative of temperature fluctuations, $\langle {(\partial \theta /\partial z)}^2\rangle$, at the time $t_5$ in Figure 3 are shown in Figure 7. After a long time, the oscillatory interaction of ${\varphi }_1$ and ${\varphi }_2$ is presumably near equilibrium. The intervals of the peaks are almost equally periodic, which means that the high intermittency of the temperature gradient is softened. The characteristic frequency of the peaks can be compared with the value of $2\pi U/N$ in Figure 7. However, the frequency should be related with that of the forcing, and the true reason for the interval is still an open question.

Figure 7. The horizontal average of the horizontal energy, $(\langle u^2\rangle +\langle v^2\rangle )/2$, and the horizontal average of the square of the vertical derivative of temperature fluctuations, $\langle {(\partial \theta /\partial z)}^2\rangle$, at about the final time $t_5$ in Figure 3. The characteristic frequency of the peaks is compared with the value of $2\pi U/N$.

Figure 7. The horizontal average of the horizontal energy, $(\langle u^2\rangle +\langle v^2\rangle )/2$, and the horizontal average of the square of the vertical derivative of temperature fluctuations, $\langle {(\partial \theta /\partial z)}^2\rangle$, at about the final time $t_5$ in Figure 3. The characteristic frequency of the peaks is compared with the value of $2\pi U/N$.

Now we look at the contour plot of the temperature fluctuations in an $x-z$ plane at $y=\pi$ (grid number 512) for $N^2=10,k_f=4$ at $t=64.7340$ (Figure 8). (To see the structures clearly, a grey scale is used. The grid number is depicted in the figure as the x and z coordinates).

Many horizontally extended layers are observed which have sharp interfaces between black and white regions. The existence of such interfaces implies that the temperature fluctuations jump at the interfaces. To observe such jumps in temperature fluctuations, its value at $x=134$ is monitored along the white line in Figure 8.

Figure 9a shows the graph of the temperature fluctuation $\theta$ as a function of z (in terms of grid number) along $x=134$. Horizontally flat and straight lines are observed periodically which are bridged by finely oscillatory stairs.

In Figure 9b, the saw-tooth pattern of Figure 9a is superposed with the mean temperature, $\langle \theta \rangle$, corresponding to $N^2=10$. The overall saw-tooth structures in Figure 9a produce cliff–ramp structures in the total temperature field. (Major cliffs are marked with red circles). Temperature jumps at the height of the cliffs. In other words, the cliffs form the interface between cold and warm fluid masses. The most characteristic feature of the cliff–ramp structures is that ramps are mostly leftward (i.e., cooling as it goes up) with jitters. Similar temperature front structures are observed in simulations of homogeneous shear flows [23,24] as well as in large-eddy simulations of the stably stratified atmospheric boundary layer [25].

To compensate for the cooling ramps, cliffs should be warming as it goes up. Therefore, there is an asymmetry in the slope of the temperature distribution. Such asymmetry is clearly seen statistically if we look at the PDFs of temperature fluctuation and its derivatives in the whole domain (Figure 10).

In Figure 10, the PDF of temperature fluctuation is close to Gaussian, while those of its derivatives show exponential tails. Among the three derivatives, the derivatives in x- and y- are symmetric, while the z- derivative is skewed with the positive values being higher probability. This statistics supports the phenomenon that the cliffs are associated more with warming.

4. Discussion

Apparently it appears odd that the integration of the Navier–Stokes equations with periodic boundary conditions produces non-symmetric statistics of the derivative of the temperature fluctuations. However, we can check that the integration of ${\theta }_z$ in z is zero which is compatible with the periodic boundary condition. Even under this circumstance, we suspect that the intermittent velocity components (Figure 6) acting as a strong horizontal “wind” may produce a non-symmetric gradient of temperature fluctuations. Nevertheless, clarification of the true mechanism is still an open question.

The present research on the development of structures and spectra in stably stratified turbulence is largely related to the broad study of mixing and transport in turbulence. This old and important subject has been investigated by many researchers. The reader is referred to some of the basic reviews and the references therein for the historic and recent developments [18,26,27,28,29].