3.1. Derivation of the Kármán–Howarth–Monin equations

Consider two measured points located at $\boldsymbol {x}$ and $\boldsymbol {x}'=\boldsymbol {x}+\boldsymbol {r}$, where $\boldsymbol {r}$ is the displacement of these two points. We denote the quantities evaluated at $\boldsymbol {x}'$ with a prime, e.g. $u' = u(\boldsymbol {x}') = u(\boldsymbol {x}+\boldsymbol {r})$. Homogeneity implies that for two-point statistical quantities, the spatial derivatives follow

(3.1)\begin{equation} \boldsymbol{\nabla} ={-}\boldsymbol{\nabla}' ={-}\boldsymbol{\nabla}_{\boldsymbol{r}}, \end{equation}

where $\boldsymbol {\nabla }'$ and $\boldsymbol {\nabla }_{\boldsymbol {r}}$ denote the gradients taken with respect to $\boldsymbol {x}'$ and $\boldsymbol {r}$, respectively.

For statistically steady states, by multiplying $\psi '$, $\sigma '$ and $\theta '$ to (2.5a), (2.5b) and (2.5c), and then adding the corresponding conjugate equations, respectively, we obtain the Kármán–Howarth–Monin (KHM) equations (cf. Monin & Yaglom Reference Monin and Yaglom1975; Frisch Reference Frisch1995) for (2.5a)–(2.5c):

(3.2a)$$\begin{gather} {-}\frac{1}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}\overline{\delta \boldsymbol{u} \left( \delta u^2 + \delta w^2 \right)} + \left( \overline{u\sigma'} + \overline{u'\sigma} - \overline{w\theta'} - \overline{w'\theta} \right) = D_\psi, \end{gather}$$

(3.2b)$$\begin{gather}{-}\frac{1}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}\overline{\delta \boldsymbol{u}\,\delta \sigma^2 } + \left( \overline{u\sigma'} + \overline{u'\sigma} \right) = D_\sigma, \end{gather}$$

(3.2c)$$\begin{gather}{-}\frac{1}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}\overline{\delta \boldsymbol{u}\,\delta \theta^2 } - \left( \overline{w\theta'} + \overline{w'\theta} \right) = D_\theta, \end{gather}$$

where

(3.3a)$$\begin{gather} D_\psi ={-}2\nu\,\overline{\nabla^4\psi\,\nabla'^4\psi'} - 2\alpha \overline{\psi\psi'}, \end{gather}$$

(3.3b)$$\begin{gather}D_\sigma ={-}2\kappa\,\overline{\nabla^3\sigma\,\nabla'^3\sigma'} - 2\mu \overline{\varDelta^{1/2}\sigma\varDelta'^{1/2}\sigma'}, \end{gather}$$

(3.3c)$$\begin{gather}D_\theta ={-}2\kappa\,\overline{\nabla^3\theta\,\nabla'^3\theta'} - 2\mu \overline{\varDelta^{1/2}\theta\varDelta'^{1/2}\theta'}, \end{gather}$$

are the effects of dissipation and diffusivity. Note that $D_\psi |_{\boldsymbol {r}=\boldsymbol {0}}$, $D_\sigma |_{\boldsymbol {r}=\boldsymbol {0}}$ and $D_\theta |_{\boldsymbol {r}=\boldsymbol {0}}$ are all negative.

For the present isotropic convection system, we know qualitatively that the instability brings about energy injection through the buoyancy terms, i.e. the quadratic terms on the left-hand sides of (3.2a)–(3.2c), then the nonlinear advection transfers energy across scales, and finally, the energy is dissipated due to dissipation and diffusivity. Since an important feature of the BO theory is the interaction between the potential energy and kinetic energy, this physical picture inspires us to argue the BO scaling from the energy transfer processes.

3.2. Argument based on energy flux

Taking the Fourier transform of the KHM equations (3.2a)–(3.2c), we obtain the spectral energy equation, which can be written symbolically as

(3.4a)$$\begin{gather} \partial_{K} F_K + B = D_K, \end{gather}$$

(3.4b)$$\begin{gather}\partial_{K} F_P + B = D_P, \end{gather}$$

where $F$ is the energy flux across scales, $B$ is the buoyancy effect that injects energy into the system, $D$ is the dissipation, and the subindices $K$ and $P$ denote the kinetic and potential energy, respectively.

We now argue the existence of BO scaling based on the transfer directions of the kinetic and potential energy. To be specific, we assume that for the present 2-D isotropic system, the kinetic energy transfers upscale, while the potential energy transfers downscale. This assumption is based on the 2-D turbulence's upscale energy transfer (Kraichnan Reference Kraichnan1967) and the forward cascade of temperature field in a turbulent velocity field (Obukhov Reference Obukhov1949; Corrsin Reference Corrsin1951). Although in some cases, the existence of potential parts may lead to a downscale flux of the kinetic energy in two dimensions (Xie & Bühler Reference Xie and Bühler2019b), in § 4 we will demonstrate numerically the validity of this assumption for the present system.

The potential energy transfers downscale, so if the buoyancy generation is a power function of wavenumber (i.e. a simple cascade picture), say $B\sim K^{-a}$ ($a>1$), then the potential energy flux can be expressed as

(3.5)\begin{equation} F_P \propto \int_{K_0}^{K} \tilde K^{{-}a} \,\mathrm{d} \tilde K = K_0^{{-}a+1} - K^{{-}a+1}, \end{equation}

where $K_0$ is the large-scale dissipation wavenumber. This formula immediately implies that a constant potential energy flux is seen at small scales when $K/K_0\gg 1$. We can also obtain the relative importance of the constant flux and the energy injection by comparing $F_P$ and $K\,\partial _{K}F_P$:

(3.6)\begin{equation} \left| \frac{F_P}{K\,\partial_{K}F_P} \right| = \frac{1}{a-1}\left( \left( \frac{K_0}{K} \right)^{1-a} - 1 \right), \end{equation}

which tends to infinity, implying that the flux is dominant as $K\to \infty$. Therefore, the BO scaling induced by constant potential energy flux should be observed at small scales.

For the kinetic energy, which transfers upscale, we have

(3.7)\begin{equation} F_K \propto \int_{-\infty}^{K} \tilde K^{{-}a} \,\mathrm{d} \tilde K ={-}\frac{1}{1-a}\,K^{{-}a+1}. \end{equation}

Note that because of the opposite directions of energy transfer, the integration directions in (3.7) and (3.5) are opposite. Similarly, we have

(3.8)\begin{equation} \left| \frac{F_K}{K\,\partial_{K}F_K} \right| = \frac{1}{a-1}, \end{equation}

implying a balance between the transfer and injection of kinetic energy, which completes the BO scenario.

It is noteworthy that because of the inverse cascade of kinetic energy, the Bolgiano scale is not important in the 2-D case. To be specific, we can estimate the Bolgiano scale using the expression (cf. Lohse & Xia Reference Lohse and Xia2010)

(3.9)\begin{equation} l_B = \epsilon_K^{5/4}\epsilon_P^{{-}3/4}(\beta g)^{{-}3/2}, \end{equation}

where $\epsilon _K$ and $\epsilon _P$ are the small-scale kinetic and potential energy dissipation rates, respectively. Because $\epsilon _K$ is approximately zero as the downscale energy flux is negligibly small in 2-D turbulence, we have $l_B=0$ theoretically, which in numerical simulations is approximately equal to the smallest resolved length scale (Celani et al. Reference Celani, Matsumoto, Mazzino and Vergassola2002). The negligibly small $l_B$ also implies the non-existence of the KO scaling in the present 2-D system.

3.3. Scalings of structure functions

From the perspective of the KHM equations (3.2a)–(3.2c), the BO scaling is consistent with the knowledge of 2-D and 3-D turbulence, which transfers energy upscale and downscale, respectively. The key is the link between the direction of energy flux and the anomaly's influence on the procedure of obtaining the expressions for the third-order structure functions. Specifically, in 3-D turbulence, where the kinetic energy transfers downscale and dissipates, we can ignore the effect of energy injection, and take the leading-order approximation of the dissipation as the energy dissipation in the KHM equations to obtain the Kolmogorov $4/5$-law (Kolmogorov Reference Kolmogorov1941; Frisch Reference Frisch1995). While for 2-D turbulence, where the kinetic energy transfers upscale, we need to subtract the leading-order constant energy dissipation from the dissipation term in the KHM equations (Lindborg Reference Lindborg1999; Xie & Bühler Reference Xie and Bühler2018).

Therefore, assuming that the buoyancy potential energy transfers downscale, following the argument of 3-D turbulence, (3.2b) implies

(3.10)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\overline{\delta \boldsymbol{u}\,\delta \sigma^2 } = {\rm constant}. \end{equation}

Similarly, as the kinetic energy transfers upscale in our 2-D system, after subtracting the leading-order energy injection by buoyancy terms and the constant energy dissipation, (3.2a) results in

(3.11)\begin{equation} -\frac{1}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}\overline{\delta \boldsymbol{u} (\delta u^2 + \delta w^2 )} + f(-\overline{\delta u\,\delta \sigma} + \overline{\delta w\,\delta\theta}) = 0. \end{equation}

Combining (3.10) and (3.11), and invoking the dimensional analysis, we obtain the third- and second-order structure functions

(3.12a)$$\begin{gather} V_K^{(1)} \equiv \overline{\delta u( \delta u^2 + \delta w^2 )} \sim r^{9/5}, \end{gather}$$

(3.12b)$$\begin{gather}V_P^{(1)} \equiv \overline{\delta u(\delta \sigma^2 + \delta \theta^2)} \sim{-}r, \end{gather}$$

(3.12c)$$\begin{gather}\overline{\delta u^2} \sim \overline{\delta w^2} \sim r^{6/5}, \end{gather}$$

(3.12d)$$\begin{gather}\overline{\delta \sigma^2} \sim \overline{\delta \theta^2} \sim r^{2/5}, \end{gather}$$

which are exactly the predictions of the BO theory. Here, the upper index $^{(1)}$ denotes the horizontal component of the structure function vector, and correspondingly, we define $V_K^{(2)} \equiv \overline {\delta w\left ( \delta u^2 + \delta w^2 \right )}$ and $V_P^{(2)} \equiv \overline {\delta w\left ( \delta \sigma ^2 + \delta \theta ^2 \right )}$ for the vertical components. The signs of the third-order structure functions correspond to the direction of energy flux. The positive and negative signs correspond to the downscale and upscale fluxes, respectively (Cho & Lindborg Reference Cho and Lindborg2001; Xie & Bühler Reference Xie and Bühler2019a). Note that to obtain (3.12b), we do not need to invoke the dimensional analysis, so it is more robust compared with other relations. Also, in above derivation we make use of the inverse cascade of kinetic energy to obtain (3.11), therefore we claim that the BO scaling coexists with the inverse kinetic energy cascade, which is supported by our numerical results shown in § 4.

The high-order structure functions may also be power functions of the two-point distance, i.e.

(3.13a,b)\begin{equation} \overline{\delta u^n} \sim r^{\zeta_u} \quad \mathrm{and} \quad \overline{\delta \theta^n} \sim r^{\zeta_\theta}, \end{equation}

where $\zeta _u$ and $\zeta _\theta$ are both functions of $n$. In a non-intermittent situation, the BO scenario implies

(3.14a,b)\begin{equation} \zeta_u = \tfrac{3}{5}n \quad \mathrm{and} \quad \zeta_\theta = \tfrac{1}{5}n. \end{equation}

4.1. Turbulent fields

We first present some important features of the turbulent fields using the numerical results with $\nu =10^{-16}$. Figure 1 shows snapshots of vorticity $\nabla ^2 \psi$ and two buoyancy fields. It is seen that the vertical buoyancy field $\theta$, which contains up- and down-moving plumes, just resembles the temperature field in the 2-D homogeneous RBC system (Celani et al. Reference Celani, Lanotte, Mazzino and Vergassola2000). However, it differs from the horizontal buoyancy field $\sigma$, where left- and right-moving plumes are present. Their difference can be seen from the asymmetry of sharp interfaces: in the $\theta$ field, the horizontal sharp interfaces prefer negative gradients of $\theta$, while the vertical interfaces consist of nearly symmetric positive and negative gradient of $\theta$; and vice versa for the $\sigma$ field. Thus the combination of $\theta$ and $\sigma$ leads to the present isotropic convection.

Figure 1.

Snapshots of $\nabla ^2\psi$, $\sigma$ and $\theta$ at a statistically steady state.

The isotropic feature of the turbulent fields is illustrated further in figure 2, where the second- and third-order isotropic fields, including $\overline {\delta u^2}+\overline {\delta w^2}$, $\overline {\delta \sigma ^2}+\overline {\delta \theta ^2}$, $\boldsymbol {\nabla }\boldsymbol {\cdot } V_{K}$ and $\boldsymbol {\nabla }\boldsymbol {\cdot } V_{P}$, are plotted separately. As expected, all these fields show isotropic behaviours, which will be confirmed quantitatively later. Although anisotropy is present at large scales, owing to the effect of the square periodic box, the flows at these scales are dominated by large-scale damping and therefore do not affect the statistics at scales where BO scaling presents.

Figure 2.

Isotropic fields of $\overline {\delta u^2}+\overline {\delta w^2}$, $\overline {\delta \sigma ^2}+\overline {\delta \theta ^2}$, $\boldsymbol {\nabla }\boldsymbol {\cdot } V_{K}$ and $\boldsymbol {\nabla }\boldsymbol {\cdot } V_{P}$.

4.2. Low-order structure functions

In § 3.2, we discussed that the key for the BO scenario is associated with the transfer directions of the potential and kinetic energy. Here, we present the numerical evidence. Figure 3 shows the kinetic and potential energy fluxes, and their difference is also presented because the buoyancy terms in the kinetic and potential energy equations cancel. It is seen that the potential energy transfers downscale and the kinetic energy transfers upscale, thus supporting the assumption proposed in § 3.2. And the wide plateau appearing in their difference justifies not only the negligible effects of dissipation and diffusivity, but also the correctness of our simulations. We also compare the potential energy fluxes in simulations with different $\nu$. It is seen that decreasing the hyperviscosity results in a wider range for the approximately constant potential energy flux, and therefore extends the possible range for detecting the BO scaling. In other words, the present isotropic convection system would exhibit the BO scaling more evidently with smaller hyperviscosities (equivalently, at larger ${{Ra}}$ values).

Figure 3.

Kinetic ($F_{K}$) and potential ($F_{P}$) energy fluxes, as well as their difference. Potential energy fluxes with different hyperviscosities are also plotted for comparison.

To verify the above discussion quantitatively, and most importantly the scaling relations (3.12) derived in § 3.3, we now check in figure 4 the third-order structure functions corresponding to the kinetic and potential energy fluxes. There exists a wide range of length scales in which the data closely follow the BO scaling. Moreover, the results in figure 4(b) show that the BO scaling range extends much further as the hyperviscosity decreases, confirming the discussion above. Here, the isotropic features of the present system are verified quantitatively by comparing the results in different directions. Similarly, the second-order structure functions shown in figure 5 also follow the BO scalings. In addition, the well-collapsed data for different directions and components further demonstrate the isotropy of the flow fields. Thus combining the information of figures 4 and 5, we can conclude safely that the BO scaling is well observed in the low-order structure functions of the present isotropic convection system.

Figure 4.

The third-order structure functions corresponding to the kinetic and potential energy fluxes at different directions. The BO scalings are plotted for reference. Again, the data with different hyperviscosities are plotted in panel (b) for comparison.

Figure 5.

The second-order structure functions of (a) velocity and (b) temperature fields for different components and directions. The BO scalings are plotted for reference.

4.3. High-order structure functions

In this section, we study the high-order structure functions. Because the flow fields well satisfy the isotropic condition, we focus on the statistics of the horizontal velocity $u$ and the temperature field $\theta$. In figure 6, we plot the structure functions of $\delta u$ and $\delta \theta$ from order $2$ to order $9$, and the corresponding scaling exponents $\zeta _u$ and $\zeta _\theta$ (cf. (3.13a,b)) in the inertial range are also shown. Note that only the absolute values of $S_n^\theta$ are exhibited, as the high-order structure functions of $\theta$ can be negative. Figures 6(c,d) show the compensated plots of the ninth-order structure functions of $u$ and $\theta$, to show the quality of scaling behaviour. The compensated scaling $r^{27.3/5}$ in the velocity structure function is very close to the BO scaling $r^{27/5}$. A bottleneck effect (cf. Frisch et al. Reference Frisch, Kurien, Pandit, Pauls, Ray, Wirth and Zhu2008) is seen in the structure function of $\theta$, but the scaling recalls similar behaviour observed in 3-D isotropic turbulence by Haugen & Brandenburg (Reference Haugen and Brandenburg2004). The bottleneck effect is not seen in the velocity structure function; this is because the kinetic cascades upscale while the potential energy cascades downscale. The high-order structure functions do have power-function dependence on two-point distance. The non-intermittent BO scaling $\zeta _u = 3n/5$ (cf. (3.14a,b)) well captures the behaviour of the high-order velocity structure functions, but the temperature structure functions show strong intermittency. It is seen that the scaling exponent $\zeta _\theta$ deviates from the BO prediction (i.e. $\zeta _\theta = n/5$) and saturates to value 0.75 as $n\geq 5$, which is similar to the results in 2-D homogeneous RBC (Celani et al. Reference Celani, Matsumoto, Mazzino and Vergassola2002). However, an interesting observation is that $\overline {\delta \theta ^3} \sim r^{3/5}$, which was not shown in previous studies.

Figure 6.

(a,b) Structure functions $S_n^u$ and $S_n^\theta$ with $n=2,3,4,5,6,7,8,9$. The curves are shifted vertically. (c,d) Compensated ninth-order structure functions of $u$ and $\theta$, where the horizontal dashed lines are presented for reference. The power-function fittings to the data in the BO scaling range are also presented. (e, f) Corresponding scaling exponents $\zeta _u$ and $\zeta _\theta$.

The intermittent effects can be also observed in the p.d.f.s of, for example, the normalized longitudinal velocity increment $\delta u_L$ and temperature increment $\delta \theta _L$, as shown in figure 7. Here, $\delta u_L$ is evaluated with the $x$-direction displacement, and $\delta \theta _L$ is evaluated with the $z$-direction displacement. It is seen that the p.d.f.s of $\delta u_L$ for different scales are well collapsed and all Gaussian-like, though the upscale energy transfer brings about positive skewness, which resembles that in the inverse cascade range of 2-D turbulence (cf. Boffetta, Celani & Vergassola Reference Boffetta, Celani and Vergassola2000; Boffetta & Ecke Reference Boffetta and Ecke2012). However, the p.d.f.s of $\delta \theta _L$ for different scales do not fall on top of each other, and their departure from the Gaussian distribution becomes much stronger at smaller scales, manifesting the effects of intermittency (Celani et al. Reference Celani, Mazzino and Vergassola2001).

Figure 7.

P.d.f.s of normalized $\delta u_L$ and $\delta \theta _L$; $r$ is the distance between two measured points. The black dashed line is the normal distribution for reference.

Figure 8 shows the integration kernels for the sixth- and eighth-order structure functions of $u$ and $\theta$, which justify the convergence of high-order structure functions. Figures 8(a,b) show a larger probability of the positive velocity difference, which is in consistent with the positive third-order structure function $\overline {\delta u_L^3} > 0$ corresponding to an upscale kinetic energy flux shown in figure 3. This contrasts with the integration kernels observed in 3-D RBC with a forward kinetic energy cascade (cf. Sun et al. Reference Sun, Zhou and Xia2006). The collapse of integration kernels with different displacements again shows the non-intermittency of velocity differences. The integration kernels of temperature structure functions shown in figures 8(c,d) show a preference of negative values, which is in accord with the plumes structures (cf. figure 1). The hot plumes travel upwards and cold plumes travel downwards, so collisions of plumes give strong negative and weak positive $\theta$ gradients in the vertical direction. This negative-value preference is also observed in 3-D RBC (Sun et al. Reference Sun, Zhou and Xia2006). Different from the integration kernels of velocity difference, those of temperature difference with different displacements do not collapse, indicating strong intermittency, which is consistent with figure 7(b).

Figure 8.

Integration kernels for the sixth- and eighth-order longitudinal structure functions of $u$ and $\theta$ with different displacements of two measured points.

4.4. RBC in a periodic domain

To gain a better understanding of the present isotropic convection system, we compare it with the homogeneous RBC system, i.e. the canonical anisotropic RBC in a periodic domain. The simulation uses the same algorithm as for the isotropic system with $\nu =10^{-16}$, but without the horizontal component in the buoyancy field. To obtain trustable statistics, we run the simulation to a statistically steady state and maintain for $8\times 10^4$ eddy turnover time, which is calculated as the inverse of the mean square vorticity. The convergence of the statistics has also been checked by calculating the integral kernels of structure functions (not shown here) as in figure 8.

Figure 9 shows some statistical quantities of the homogeneous RBC, which are the counterparts of those shown in figure 2. It is seen that even the second-order statistics have exhibited obvious anisotropic features, which elongate along the vertical direction as expected. When it comes to the third-order statistics, the anisotropy is more dramatic, resulting in the longitudinal and transversal structure functions being distinctive.

Figure 9.

Isotropic fields of $\overline {\delta u^2}+\overline {\delta w^2}$, $\overline {\delta \theta ^2}$, $\boldsymbol {\nabla }\boldsymbol {\cdot } V_{K}$ and $\boldsymbol {\nabla }\boldsymbol {\cdot } V_{P}$ in 2-D homogeneous (but anisotropic) RBC.

To further illustrate the impact of anisotropy, we show in figure 10 the third-order structure functions corresponding to kinetic and potential energy fluxes in different directions in the anisotropic RBC, which are the counterparts of those shown in figure 4. Figure 10(a) shows that both the horizontal and vertical components of the kinetic energy structure function deviate from the BO scaling, with a steeper scaling for the horizontal component and a shallower scaling for the vertical component. For the potential energy structure function, the vertical component is close to the BO scaling, while the horizontal component is shallower. The anisotropic behaviour of both the kinetic and potential energy third-order structure functions, which differ from the isotropic structure functions of our isotropic system (cf. figure 4), shows the anisotropic energy flux phenomenon in the homogeneous RBC.

Figure 10.

The third-order structure functions corresponding to the kinetic and potential energy fluxes at different directions in anisotropic RBC. The BO scalings are plotted for reference.

The anisotropy traces back to instability, which drives the turbulent dynamics. In homogeneous RBC, the linear growth rate $\lambda _{RBC}$ with respect to the zero state of (2.1a)–(2.1b) reads

(4.1)\begin{equation} \lambda_{RBC} = \frac{1}{2} \left( -(\nu_0+\kappa_0) K^2 + \sqrt{ (\nu_0+\kappa_0)^2 K^4 - 4\nu_0\kappa_0K^4 + \frac{4 g\beta k ^2}{K^2} } \right), \end{equation}

which is anisotropic. For horizontal modes with $\boldsymbol {k} = (0,m)$, the growth rate $\lambda _{RBC}$ is negative, therefore these modes are linearly stable and do not inject energy into the system, while the most unstable mode is a so-called elevator mode with $m=0$. Therefore, it is expected to observe more flow structures – plumes – propagating vertically. So in figures 9(c,d), the divergence of third-order structure functions, which correspond to energy sources and sinks, concentrate around locations with zero vertical displacement. Considering that in the spectral space the nonlinear advection transfers energy from sources, where the modes are unstable, to sinks, the anisotropic growth rate (cf. (4.1)) leads to anisotropic energy transfer (cf. figure 10), which further leads to anisotropic structure functions in inertial ranges. It is this anisotropic energy transfer that contaminates the scaling behaviour in the canonical RBC. Note that an alternative way to analyse the anisotropy in turbulent flows is by examining structure functions with different symmetry (cf. Biferale & Procaccia Reference Biferale and Procaccia2005), which is beyond the scope of this paper.