3.1. Evolution of the mixing layer

Figure 1.

Visualisations of the mass fraction fields $Y$ at dimensionless time $\widehat {t}=t \sqrt {\mathrm {At} g/L_x}=4.5$ for (a) low-At case at $\mathrm {At}=0.15$ , and (b) mid-At case at $\mathrm {At}=0.5$ .

Figure 1 shows the visualisations of the mass fraction fields $Y$ for the low-At and mid-At cases at dimensionless time $\widehat {t} = t\sqrt {\mathrm {At} g/L_x}=4.5$ , where $\sqrt {L_x/(\mathrm {At} g)}$ is the characteristic time scale of RT flows. The flow encompasses a wide range of length scales, and the mixing zone heights are similar for the two cases. A quantitative measure of the mixing zone growth is the mixing width $h(t)$ , defined as (Cabot & Cook Reference Cabot and Cook2006; Zhang et al. Reference Zhang, Betti, Yan, Zhao, Shvarts and Aluie2018b )

(3.1) \begin{align} h(t)=\frac {2}{\rho _h-\rho _l} \int _{-\infty }^\infty \mathrm {min}\left (\langle \rho \rangle _{xy}-\rho _l, \rho _h-\langle \rho \rangle _{xy}\right )\ {\rm d}z, \end{align}

where $\langle \cdot \rangle _{xy}=( {1}/{L_x L_y}) \iint (\cdot ){\rm d}x {\rm d}y$ is a horizontally averaged quantity (a function of $z$ ) and $\rho _l, \rho _h$ are densities of the light and heavy fluids. In figure 2(a) the mixing width exhibits quadratic growth with time for both cases after $\widehat {t}\gt 2$ , suggesting a self-similar expansion of the RT mixing zone. The fitted growth rate coefficient $\alpha$ , appearing in the quadratic growth expression $h(t) = \alpha \mathrm {At}\,gt^2$ , is 0.0232 for the low-At case and 0.0209 for the mid-At case, both of which fall within the range reported in the literature (Zhou 2017a).

A related metric is the mixedness parameter, $\Theta$ (see introduction), which measures the degree of molecular mixing along horizontal planes and is defined as

(3.2) \begin{align} \Theta = \frac {\int _{-\infty }^\infty \langle Y(1-Y)\rangle _{xy} {\rm d}z}{\int _{-\infty }^\infty \langle Y\rangle _{xy} \langle 1-Y\rangle _{xy} {\rm d}z}. \end{align}

The parameter $\Theta$ is intended to measure the rate of reaction in a mixture, in which the molar fraction shall be adopted in the above definition instead of mass fraction (Youngs Reference Youngs1991; Wilson & Andrews Reference Wilson and Andrews2002). Here in our RT configuration, the molecular weights of the two fluids are assumed to be equal, so the mass fraction is equal to the molar fraction, making the use of the above definition appropriate.

Figure 2(b) shows that initially, the presence of a small transition region near the fluid interface results in well-mixed flow within this region, leading to $\Theta \approx 1$ . The parameter stays around 1 for $\widehat {t}\lt 1$ during which the diffusion effect is dominant and the instability has yet to grow. After this transient period, the instability grows and the fluid is dominated by advection. The fluid interface stretches and folds, causing large regions of unmixed pure fluids to be entrained into the mixing zone, which reduces the mixedness parameter $\Theta$ . As the entrainment process progresses, the interfacial area increases, enhancing diffusion at these interfaces, which is reflected by a subsequent rise in $\Theta$ between $1.5 \leqslant \widehat {t} \leqslant 2$ . After $\widehat {t} \gt 2$ , the effect of advection and entrainment of pure fluids (which decreases $\Theta$ ) balances with the effect of interface diffusion (which increases $\Theta$ ), leading to a saturation of $\Theta$ around 0.8, consistent with results in Cook et al. (Reference Cook, Cabot and Miller2004); Banerjee & Andrews (Reference Banerjee and Andrews2009).

Figure 2.

Evolution of (a) the total mixing width $h$ and the area $A/(L_xL_y)$ , and (b) the mixedness parameter $\Theta$ of the isosurface $Y=0.5$ for the low-At and mid-At simulation cases. In the inset of panel (a), the square root of the mixing width is depicted to measure the growth coefficient $\alpha$ corresponding to the quadratic growth rate $h(t)=\alpha \mathrm {At}gt^2$ marked by the black solid lines. The low-At case in the inset is shifted upward by 0.1 to distinguish it from the mid-At case, and the corresponding $\alpha$ values are 0.0232 and 0.0209 for the low-At and mid-At cases, respectively. Throughout this work, solid and dashed lines represent the low-At and mid-At cases, respectively, unless otherwise specified.

To further reveal the roles of advection and diffusion, in figure 2(a) we also show the evolution of the normalised area $A/(L_xL_y)$ of the $Y=0.5$ isosurface, obtained using the marching cube algorithm (Lorensen & Cline Reference Lorensen and Cline1998). The algorithm divides space into voxels (small cubes), interpolates isosurface intersections along their edges and connects these points based on vertex positions to form a piecewise-linear isosurface approximation. In accordance with the evolution of the mixedness $\Theta$ , the isosurface areas are close to flat with a value of $L_xL_y$ during $\widehat {t}\lt 1$ , while they increase rapidly within $1\lt \widehat {t}\lt 1.5$ and the growths are close to exponential, similar to the area increase of a material surface in isotropic turbulence (Batchelor Reference Batchelor1952). After $\widehat {t}\gt 1.5$ , the enhanced diffusional effect hinders the growth of the isosurface area, and the exponential growth rate shifts towards a nearly linear growth in time after $\widehat {t}\gt 2$ , during which the system is self-similar and the mixedness parameter remains roughly constant.

Figure 3.

(a) Visualisation of the isosurface of the low-At case at $\widehat {t}=5$ . (b) The temporal evolution of the fractal dimension $D_3$ of the isosurfaces, measured with the box-counting algorithm (Mandelbrot Reference Mandelbrot1982). (c) Visualisation of isolines on three 2-D slices along the $x$ -, $y$ -, and $z$ -normal directions for the low-At case at $\widehat {t}=5$ . (d) The fractal dimension $D_2$ of the isolines at different slice locations $d$ , where $d$ denotes the $x$ , $y$ , or $z$ location of the slices relative to the domain centre. The low-At cases are shown in solid lines, while the mid-At cases are in dashed lines.

Unlike the mixing width, the interfacial area of the mid-At case is consistently larger than that of the low-At case, suggesting that the isosurface is more wrinkled with a greater heavy--light density contrast, and therefore, more potential energy is injected into the system. A quantitative measure of these interfacial topologies, shown in figure 3(a) for isosurfaces or figure 3(c) for cross-sectional isolines, is the fractal dimension, $D_3$ and $D_2$ , embedded in 3-D and two-dimensional (2-D) space, respectively. The fractal dimension is calculated using the box-counting algorithm (Mandelbrot Reference Mandelbrot1982; Sreenivasan & Meneveau Reference Sreenivasan and Meneveau1986), as detailed in Appendix B. Here, $D_3 = 1 + D_2$ (Sreenivasan & Meneveau Reference Sreenivasan and Meneveau1986), and a larger fractal dimension indicates a more complex interfacial manifold.

Temporal evolution of the fractal dimension is depicted in figure 3(b), which tracks the evolution of $D_3$ for the $Y=0.5$ isosurfaces. Between $1 \leqslant \widehat {t} \leqslant 1.5$ , the isoscalar surfaces rapidly evolve from a flat surface with a dimension of 2.0 to a highly distorted level set. After this transient phase, in the self-similar regime for $\widehat {t} \gt 2$ , the fractal dimension stabilises around 2.41 for the low-At case and 2.46 for the mid-At case. Compared with reported values of $D_2 = 1.7{-}1.8$ for 2-D RT turbulence (Hasegawa, Nishihara & Sakagami Reference Hasegawa, Nishihara and Sakagami1996) and $D_2 \leqslant 1.3$ for 3-D RT turbulence (Linden, Redondo & Youngs Reference Linden, Redondo and Youngs1994), both obtained from simulations without viscosity or mass diffusion, the RT isosurfaces in our 3-D simulations display different fractal behaviour when diffusion effects are included. Furthermore, the anisotropy in RT turbulence is manifested in the fractal dimension of isolines on different cross-sections. In figure 3(d) the fractal dimension of the $Y=0.5$ isolines on cross-sections in the $x$ -, $y$ -, and $z$ -normal planes is shown for both low-At and mid-At cases at $\widehat {t}=5$ . The variable $d$ represents the location of the embedded plane relative to the domain centre. In both cases, isolines on $x$ - and $y$ -normal planes exhibit similar fractal dimensions, with mean values of 1.4 for the low-At case and 1.45 for the mid-At case. In comparison, isolines on the $z$ -normal planes exhibit a fractal dimension of approximately 1 at the edges of the mixing zone (large $d$ magnitude), while in the middle region ( $d\approx 0$ ), the peaks of $D_2$ exceed the values of the $x$ - and $y$ -normal isolines. This anisotropic behaviour of $D_2$ reflects the anisotropic transport of the scalar field.

3.2. Evolution of scalar variance and gradient magnitude

In RT turbulence the evolution of the scalar field variation is crucial in determining the mixing dynamics. Similar to the equations for turbulent KE and strain/vorticity dynamics, we analyse the governing equations for scalar field variance and its gradient, which are

(3.3) \begin{gather} \rho \frac {D}{Dt}\frac {1}{2} Y^{\prime2}= -\underbrace {\rho D|\nabla Y^{\prime}|^2}_{\mathcal {\varepsilon }^Y} + \nabla \cdot \left (\rho D \nabla \frac {Y^{\prime2}}{2}\right ), \end{gather}

(3.4) \begin{gather} \frac {D}{Dt}\frac {1}{2}|\nabla Y|^2= \underbrace {-\nabla Y \cdot \mathbf {S}\cdot \nabla Y }_{\mathcal {S}^Y} + \underbrace {\nabla Y \cdot \nabla \left [ \frac {1}{\rho } \nabla \cdot (\rho D\nabla Y)\right ]}_{\mathcal {\varepsilon }^{\nabla Y}}, \end{gather}

where $D/Dt$ is the material derivative, $Y^{\prime}=Y-Y_0$ is the fluctuating part of mass fraction and $Y_0=\langle \rho Y\rangle /\langle \rho \rangle$ is the density-weighted mean concentration. Note that the equation for $Y^2$ has the same expression as (3.3) for $Y^{\prime2}$ , since $Y_0$ remains constant over time due to the conservations of both mass and mass fraction. Here, we define new terms $\mathcal {\varepsilon }^Y$ , $\mathcal {S}^Y$ and $\mathcal {\varepsilon }^{\nabla Y}$ , which denote the scalar dissipation rate, the scalar element stretching rate and the scalar gradient dissipation rate, respectively.

Figure 4.

( $a$ ) Temporal evolution of the the mean scalar dissipation $\langle \mathcal {\varepsilon }^Y\rangle$ , the mean scalar energy $\langle ({1}/{2})\rho Y^2 \rangle$ , as well as the three components of the mean squared mass fraction gradient. ( $b$ ) The evolution of the mean scalar element stretching $\langle \mathcal {S}^Y\rangle$ and the mean scalar gradient dissipation $\langle \mathcal {\varepsilon }^{\nabla Y}\rangle$ . Solid and dashed lines represent the low-At and mid-At cases, respectively.

Equation (3.3) indicates that the spatial mean concentration variance, when scaled by density, decreases monotonically over time, with a decreasing rate of $\langle \mathcal {\varepsilon }^Y \rangle$ . In figure 4(a) we present the temporal evolution of the mean mass fraction energy, $\langle ({1}/{2}) \rho Y^2 \rangle \equiv \langle ({1}/{2}) \rho Y^{\prime2} \rangle + ({1}/{2})\langle \rho \rangle Y_0^2$ , instead of the variance $\langle ({1}/{2}) \rho Y^{\prime2} \rangle$ , to have the same initial value between the low-At and mid-At cases. The results confirm the monotonic decay of concentration energy, and consequently, the mass fraction variance. Initially, the concentration energy is slightly less than 0.25, corresponding to fully segregated heavy and light fluids ( $Y=1, \rho _h=1$ on the upper half-domain so that $\langle ({1}/{2}) \rho Y^2 \rangle |_{t=0}=0.25$ ), except in a narrow region near the fluid interface. The decay rate is slow for $\widehat {t} \lt 1.5$ , before the onset of turbulence, but after $\widehat {t} \gt 2$ , during the self-similar phase, the decay rates for both cases increase, accompanied by a growing scalar gradient magnitude $\langle |\nabla Y|^2\rangle$ . Notably, the decay rate in the low-At case is faster than in the mid-At case, despite the reverse trend in the relative magnitude of the scalar gradient. This trend indicates that the density difference of the light fluids between the two cases outweighs the differences in scalar gradient magnitudes, and the constant mass diffusivity adopted in the simulations leads to a higher scalar dissipation rate $\mathcal {\varepsilon }^Y$ for the low-At case.

The scalar gradient magnitude budget equation (3.4) indicates that $ ({1}/{2})|\nabla Y|^2$ is not a conserved invariant due to surface stretching, unlike the scalar field itself. The second term on the right-hand side, $\mathcal {\varepsilon }^{\nabla Y}$ , represents the dissipation of the scalar gradient, with its mean value being consistently negative and approximately equal to $\langle -D (\nabla ^2Y)^2 \rangle$ in the $\mathrm {At}\to 0$ limit. The blue lines in figure 4(b) confirm that the spatial mean value of $\mathcal {\varepsilon }^{\nabla Y}$ is always negative. In contrast, the first term $\mathcal {S}^Y$ on the right-hand side is a primary source term in the growth of $({1}/{2})|\nabla Y|^2$ for RT flows. To confirm that $\mathcal {S}^Y$ indeed represents the stretching of scalar elements, we express it as

(3.5) \begin{align} \mathcal {S}^Y = -\nabla Y\nabla Y:\mathbf {S}= -|\nabla Y|^2\boldsymbol {n}\boldsymbol {n}:\nabla \boldsymbol {u} \approx |\nabla Y|^2 (\mathbf {I}-\boldsymbol {n}\boldsymbol {n}):\nabla \boldsymbol {u}, \end{align}

where $\boldsymbol {n}=\nabla Y/|\nabla Y|$ is the normal vector of the mass fraction isosurface pointing from the light fluid towards the heavy fluid. The incompressibility condition $\nabla \cdot \boldsymbol {u}=0$ approximately holds in the current compressible RT simulations, given that the initial conditions are isopycnal and the Mach numbers are low (Zhao et al. Reference Zhao, Betti and Aluie2022). The fractional rate of change of a surface element $A$ is given by $({1}/{A})({DA}/{Dt}) \equiv (\mathbf {I} - \boldsymbol {n}\boldsymbol {n}):\nabla \boldsymbol {u}$ when surface propagation due to diffusion is negligible (Poinsot & Veynante Reference Poinsot and Veynante2005). Therefore, the right-hand side of (3.5) represents the stretching rate of the mass fraction isosurface along the tangential plane, confirming that $\mathcal {S}^Y$ is a weighted stretching term.

Due to the initial stratification, the vertical scalar gradient is greater than the horizontal terms at early times, as is indicated in figure 4(a). However, when $\widehat {t} \gt 1$ , the horizontal derivatives of the mass fraction exceed the vertical component, consistent with the formation of coherent bubbles and spikes that are primarily elongated in the vertical direction. Between $1 \lt \widehat {t} \lt 2$ , all components of the scalar gradient undergo a rapid increase followed by a decrease. This behaviour results from the competing effects of advection-induced surface stretching and diffusion-induced scalar gradient dissipation, as described by the right-hand side of (3.4), and corresponds with the non-monotonic evolution of the mixedness parameter $\Theta$ during this period, as seen in figure 2(b). After $\widehat {t} \gt 2$ , the mean values of the mass fraction gradients increase over time, with directional anisotropy between horizontal and vertical components persisting, especially in the mid-At case with greater density contrast. These results indicate that the horizontal scalar gradient is larger than the vertical one, implying a smaller horizontal gradient length scale, $L_{\varDelta , x}$ , compared with the vertical length scale, $L_{\varDelta , z}$ :

(3.6) \begin{align} \frac {1}{L_{{\varDelta}, x}} \equiv \left (\frac {\langle |\partial _x Y|^2 \rangle }{\langle Y^2\rangle }\right )^{1/2} \; \geqslant \; \frac {1}{L_{\varDelta , z}} \equiv \left (\frac {\langle |\partial _z Y|^2 \rangle }{\langle Y^2\rangle }\right )^{1/2}. \end{align}

Physically, this means that at intermediate to small scales (associated with the $\nabla$ operator), the ‘eddies’ formed by the mass fraction field are elongated in the vertical direction relative to the horizontal direction (Zhao & Aluie Reference Zhao and Aluie2023).

The budget equations for scalar variance (3.3) and scalar gradient (3.4) respectively describe, at large and small scales, the annihilation and creation of the spatial variation of the transported scalar, suggesting an intrinsic connection between scalar dissipation, $\mathcal {\varepsilon }^Y$ , and scalar element stretching, $\mathcal {S}^Y$ . The time evolution of their spatial mean values shown in figures 4(a) and 4(b) supports this speculation, while the joint PDFs of the two quantities further validate their pointwise correlation (see figure 18 in Appendix C). The results suggest a shared underlying physical mechanism between $\mathcal {\varepsilon }^Y$ and $\mathcal {S}^Y$ , which is analogous to the connection between the pseudo-dissipation rate of KE, $\varepsilon = \nu |\boldsymbol {\omega }|^2$ , and the vortex stretching, $\boldsymbol {\omega }\cdot \mathbf {S}\cdot \boldsymbol {\omega }$ , where $\mathbf {S}=({1}/{2})(\nabla \boldsymbol {u}+\nabla \boldsymbol {u}^T)$ and $\boldsymbol {\omega }=\nabla \times \boldsymbol {u}$ are the strain rate tensor and vorticity. Vortex stretching increases local vorticity by stretching vortex tubes, leading to a higher energy dissipation rate. Similarly, in the case of mass fraction, scalar element stretching $\mathcal {S}^Y$ amplifies the local scalar gradient, thereby increasing the scalar dissipation rate $\mathcal {\varepsilon }^Y$ .

To further elucidate the connection between $\mathcal {S}^Y$ and $\mathcal {\varepsilon }^Y$ , we diagonalise the strain rate tensor via eigendecomposition, yielding $\mathbf {S}= \lambda _\alpha \boldsymbol {e}_\alpha \boldsymbol {e}_\alpha + \lambda _\beta \boldsymbol {e}_\beta \boldsymbol {e}_\beta + \lambda _\gamma \boldsymbol {e}_\gamma \boldsymbol {e}_\gamma$ , where $\lambda _\alpha \geqslant \lambda _\beta \geqslant \lambda _\gamma$ are the eigenvalues and $\boldsymbol {e}_\alpha , \boldsymbol {e}_\beta , \boldsymbol {e}_\gamma$ are the corresponding eigenvectors. This decomposition allows us to express the scalar element stretching as

(3.7) \begin{align} \begin{split} \mathcal {S}^Y &= -\lambda _\alpha |\nabla Y\cdot \boldsymbol {e}_\alpha |^2-\lambda _\beta |\nabla Y\cdot \boldsymbol {e}_\beta |^2-\lambda _\gamma |\nabla Y\cdot \boldsymbol {e}_\gamma |^2\\ &\approx -\lambda _\gamma |\nabla Y|^2 =\frac {-\lambda _\gamma }{\rho D}\mathcal {\varepsilon }^Y, \end{split} \end{align}

where the approximation holds because the scalar gradient $\nabla Y$ is primarily aligned with the eigenvector associated with the most compressive eigenvalue, $\lambda _\gamma \lt 0$ , as will be verified in figure 5. In addition, the coefficient of variation (standard deviation divide by mean) of both $-\lambda _\gamma$ and $\rho D$ are much smaller than that of $|\nabla Y|^2$ , indicating that the relative spatial variation of the field $\mathcal {\varepsilon }^Y$ is minimally influenced by the coefficient $-\lambda _\gamma /\rho D$ . Similarly, the connection between vortex stretching and the pseudo-dissipation rate $\varepsilon =\nu |\boldsymbol {\omega }|^2$ is

(3.8) \begin{align} \boldsymbol {\omega }\cdot \boldsymbol {S} \cdot \boldsymbol {\omega } \approx \lambda _\beta |\boldsymbol {\omega }|^2 = \frac {\lambda _\beta }{\nu } \varepsilon \end{align}

since vorticity is preferentially aligned with $\boldsymbol {e}_\beta$ , and that the coefficient of variation of $\lambda _\beta$ is smaller than that of $|\boldsymbol {\omega }|^2$ .

Figure 5.

( $a$ ) The PDF of the absolute cosine of the angle between mass fraction gradient $\nabla Y$ and the three strain eigenvectors. ( $b$ ) The PDF of the absolute cosine of the angle between vorticity $\boldsymbol {\omega }$ and the strain eigenvectors. ( $c$ ) The PDF of the ratios of the largest and intermediate eigenvalues, $\lambda _\alpha$ and $\lambda _\beta$ , to the magnitude of the most compressive eigenvalue, $|\lambda _\gamma |$ . The two vertical dashed lines correspond to the values of 0.285 and 0.710. ( $d$ ) The PDF of the cosine value between vorticity and $\nabla Y$ . The solid, dashed and dotted lines represent the cases of low-At, mid-At, and an incompressible passive scalar turbulence, respectively.

The preferential alignment of the scalar gradient with the eigenvector corresponding to the most compressive eigenvalue is well established for passive scalars in isotropic turbulence (Batchelor Reference Batchelor1959; Ashurst et al. Reference Ashurst, Kerstein, Kerr and Gibson1987). Similarly, in RT turbulence, figure 5(a) shows that the scalar gradient $\nabla Y$ indeed aligns with the eigenvector corresponding to the most compressive eigenvalue, $\boldsymbol {e}_\gamma$ , consistent with the findings of Cabot & Zhou (Reference Cabot and Zhou2013) in incompressible variable-density RT simulations. It indicates that $\nabla Y$ most likely aligns with $\boldsymbol {e}_\gamma$ and is perpendicular to $\boldsymbol {e}_\beta$ , while its alignment with $\boldsymbol {e}_\alpha$ appears more random. This alignment pattern is observed in both the low-At and mid-At cases, as well as the passive scalar in isotropic turbulence, shown in solid, dashed and dotted lines, respectively. The isotropic turbulence case corresponds to a stationary forced turbulence with passive scalar using a $512^3$ grid. The alignment of vorticity with strain eigenvectors, which are shown in figure 5(b), also resembles that in incompressible isotropic turbulence, where vorticity preferentially aligns with the eigenvector $\boldsymbol {e}_\beta$ corresponding to the intermediate eigenvalue. In figure 5(c) the PDFs of the eigenvalue ratios $\lambda _\alpha /|\lambda _\gamma |$ and $\lambda _\beta /|\lambda _\gamma |$ reveal that the most probable eigenvalue triplet $\lambda _\alpha :\lambda _\beta :\lambda _\gamma$ follows the ratio 2.84 : 1.14 : –4, close to the result of 3 : 1: –4 reported by Ashurst et al. (Reference Ashurst, Kerstein, Kerr and Gibson1987). Finally, figure 5(d) shows that the vorticity field is most likely perpendicular to the scalar gradient, with $\cos (\nabla Y, \boldsymbol {\omega })$ peaking at 0. This trend is applicable to the cases of low-At, mid-At and incompressible turbulence.

By combining the results from figure 5 with (3.7), we observe that in developing RT turbulence, the specific alignment between the strain field and the scalar gradient links the scalar dissipation rate $\mathcal {\varepsilon }^Y$ to the scalar element stretching rate $\mathcal {S}^Y$ . This connection bridges the large-scale scalar variance ${1}/{2}\rho |Y^{\prime}|^2$ with the small-scale scalar gradient magnitude $({1}/{2})|\nabla Y|^2$ . In the following § 4, we present scale decomposition to directly demonstrate that the transfer of scalar variance across length scales is equivalent to the scalar element stretching at a coarse-grained level.

3.3. Evolution of the direction of the scalar gradient

In addition to the scalar gradient magnitude, the evolution of the normal direction plays a crucial role in the alignment between $\nabla Y$ and the strain eigenframe, thereby influencing the magnitude of surface stretching term $\mathcal {S}^Y$ . The governing equation for the mass fraction normal vector, $\boldsymbol {n}=\nabla Y/|\nabla Y|$ , is given by (Dopazo et al. Reference Dopazo, Martin, Cifuentes and Hierro2018)

(3.9) \begin{align} \frac {D}{Dt}\boldsymbol {n} = \underbrace {-(\mathbf {I}-\boldsymbol {n}\boldsymbol {n}) \cdot \boldsymbol {S} \cdot \boldsymbol {n}}_{\mathcal {N}^S} + \underbrace {\frac {1}{2}\boldsymbol {\omega }\times \boldsymbol {n}}_{\mathcal {N}^\omega } + \underbrace {\frac {1}{|\nabla Y|} (\mathbf {I} -\boldsymbol {n}\boldsymbol {n})\cdot \nabla \left (\frac {1}{\rho }\nabla \cdot (\rho D\nabla Y)\right )}_{\mathcal {N}^D}. \end{align}

Unlike the evolution of the scalar gradient magnitude described in (3.4), which depends solely on the strain rate tensor in the scalar element stretching term, the evolution of the normal vector depends on both the symmetric and anti-symmetric parts of the velocity gradient tensor. The first term on the right-hand side of (3.9), $\mathcal {N}^S$ , represents the strain-induced rotation of the normal vector on the scalar isosurface, causing it to deviate from the direction of the hydrodynamic force acting on the surface. The second term, $\mathcal {N}^\omega$ , accounts for the rotation of the normal vector due to vorticity, while the third term, $\mathcal {N}^D$ , describes the reorientation of the normal vector related to scalar diffusion.

Figure 6.

( $a$ ) The PDFs of the dot product between the rate of change of the scalar gradient direction and the strain eigenvectors, where the strain eigenvectors are aligned to have $\cos (\boldsymbol {e}, \nabla Y)\geqslant 0$ everywhere. (b)–(d) The individual contributions from the strain term, the vorticity term and the diffusion term on the right-hand side of (3.9). Inset of panel (c) shows the PDF of $\mathcal {N}^\omega \cdot \boldsymbol {e}$ conditioned on $Q\gt 0, R\gt 0$ , where $Q,R$ are the second and third invariants of the velocity gradient tensor. All panels show both the low-At (solid lines) and the mid-At (dashed lines) cases at $\widehat {t}=5$ .

Figure 6 presents the contributions to the rate of change of the scalar gradient direction, illustrating the roles of different components in the preferential alignment between $\nabla Y$ and the strain eigenvectors. The total contribution is shown in panel (a), while panels (b)–(d) show the individual contributions from the terms $\mathcal {N}^S$ , $\mathcal {N}^\omega$ and $\mathcal {N}^D$ , respectively, at $\widehat {t}=5$ (similar results are observed at other time instants). We have aligned the eigenvectors such that $\cos (\nabla Y, \boldsymbol {e})\geqslant 0$ everywhere, i.e. $\cos (\boldsymbol {n}, \boldsymbol {e})=|\cos (\boldsymbol {n}, \boldsymbol {e})|$ . Thus, positive values of $({D}/{Dt})\boldsymbol {n} \cdot \boldsymbol {e}$ indicate a tendency for $\nabla Y$ to align more strongly with the respective eigenvector, while negative values suggest a weakening of the alignment. In figure 6(a) the PDFs of $({D}/{Dt})\boldsymbol {n} \cdot \boldsymbol {e}$ are roughly symmetric about zero, with a slight bias towards $({D}/{Dt})\boldsymbol {n} \cdot \boldsymbol {e}_\alpha \lt 0$ and $({D}/{Dt})\boldsymbol {n} \cdot \boldsymbol {e}_\gamma \gt 0$ . This suggests that, overall, the alignment between $\nabla Y$ and the strain eigenvectors remains temporally stable, with a minor tendency to enhance the alignment between $\nabla Y$ and $e_\gamma$ , assuming that the rate of rotation of the strain eigenbasis is a slow process (Nomura & Post Reference Nomura and Post1998) so that $({D}/{Dt})\boldsymbol {n} \cdot \boldsymbol {e} \approx ({D}/{Dt})(\boldsymbol {n} \cdot \boldsymbol {e})$ . This assumption has been verified numerically for both low- and mid-At cases at $\widehat {t}=5$ (not shown). For individual contributions, panel (b) reveals that $\mathcal {N}^S \cdot \boldsymbol {e}_\alpha \lt 0$ and $\mathcal {N}^S \cdot \boldsymbol {e}_\gamma \gt 0$ , indicating a strong reinforcement of the preferential alignment between $\nabla Y$ and $\boldsymbol {e}_\gamma$ . In contrast, panel (c) shows that $\mathcal {N}^\omega$ exhibits an opposite tendency, weakening the alignment between $\nabla Y$ and $\boldsymbol {e}_\gamma$ while promoting alignment with $\boldsymbol {e}_\alpha$ . This effect is more pronounced in vorticity-dominated regions (with positive velocity gradient invariants $Q$ and $R$ ), as highlighted in the inset of panel (c). Panel (d) shows that the diffusion term, $\mathcal {N}^D$ , has no significant preferential contribution to any of the three eigenvectors, likely due to its passive nature. Therefore, while the strain term tends to enhance alignment between the scalar gradient and the most compressive eigenvector, the vorticity term weakens this alignment. Furthermore, figures 6(b) and 6(c) show that the PDF tails exhibit an approximate exponential decay, with fatter tails compared with the Gaussian distribution, underscoring the intermittency effect associated with RT velocity gradients.

4.1. Filtering spectra of the velocity and scalar fields

To understand the scalar transfer across different length scales, we first analyse the distribution of KE and scalar variance among different length scales and present their spectra at various times. Since the presence of walls at the top and bottom boundaries complicates the use of a straightforward Fourier transform to decompose scales along the vertical direction, we employ the filtering spectrum proposed by Sadek & Aluie (Reference Sadek and Aluie2018). This approach defines the spectra for the KE and mass fraction fields as

(4.1) \begin{align} E_{\mathrm {KE}}(k_\ell ) \equiv \frac {\rm d}{{\rm d} k_\ell } \langle \overline {\rho }_\ell |\widetilde {\boldsymbol {u}}_\ell (\boldsymbol {x})|^2\rangle /2, \qquad E_Y(k_\ell )\equiv \frac {\rm d}{{\rm d} k_\ell } \langle |\overline {Y}_\ell (\boldsymbol {x})|^2 \rangle /2, \end{align}

respectively, where $k_\ell = L_z/\ell$ is the filtering wavenumber, $\overline {f}_\ell \equiv \int G_\ell ({\boldsymbol {r}}) f(\boldsymbol {x}-{\boldsymbol {r}}) d{\boldsymbol {r}}$ is the spatial coarse graining and $\widetilde {f}_\ell \equiv \overline {\rho f}/\overline {\rho }$ is the Favre density-weighted filtering. The filtering spectrum has been adopted in analysing 2-D and 3-D incompressible turbulence (Sadek & Aluie Reference Sadek and Aluie2018), RT turbulence (Zhao et al. Reference Zhao, Betti and Aluie2022), oceanic flows (Storer et al. Reference Storer, Buzzicotti, Khatri, Griffies and Aluie2022) and laser induced plasma dynamics (Yin et al. Reference Yin, Shang, Blackman, Collins and Aluie2023).

Figure 7 shows the filtering spectra for KE and mass fraction at different time instants. In panel (a) the KE spectra increase over time across all scales, with their peaks gradually shifting toward lower $k_\ell$ , corresponding to larger scales. This behaviour, when interpreted through the Favre scale decomposition (Zhao & Aluie Reference Zhao and Aluie2018), reflects the conversion of potential energy to KE at the largest scale, followed by a cascade of KE to progressively smaller scales driven by energy fluxes in RT turbulence (Zhao et al. Reference Zhao, Betti and Aluie2022).

In contrast, the mass fraction field, whose variance is an ideal invariant with no input at any scale, behaves differently compared with KE. Figure 7 (b) shows that the spectra $E_Y$ increase at small scales while decreasing at large scales. After $\widehat {t} \geqslant 4$ , the small-scale growth of $E_Y$ reaches saturation, but the decrease at large scales continues. We also include $E_Y$ at the initial instant ( $\widehat {t} = 0$ ) to illustrate the persistent trend of decreasing large-scale and increasing small-scale contents. The behaviour of $E_Y$ is closely linked to the cross-scale transfer of mass fraction variance, as we shall discuss shortly.

Figure 7.

The filtering spectra of the ( $a$ ) KE normalised by $(\rho _h-\rho _l) g L_x$ and ( $b$ ) mass fraction fields for the low-At (solid lines) and mid-At (dashed lines) cases at different times. The mid-At spectra of mass fraction are shifted downwards by 1 unit to distinguish from the low-At cases. The black dashed lines represent the $k_\ell ^{-5/3}$ scaling.

4.2. Budgets of coarse-grained scalar variance and scalar gradient

Based on Favre decomposition, the budget equations for large-scale scalar variance and the scalar gradient are

(4.2) \begin{align} &\overline {\rho } \widetilde {\frac {D}{Dt}} \frac {1}{2}\widetilde {Y}_\ell ^2 = -\nabla \cdot \left [\overline {\rho } \widetilde {Y}\widetilde {\tau }(\boldsymbol {u}, Y) - \widetilde {Y} \overline {\rho D\nabla Y}\right ] - \widetilde {\Pi }_\ell ^Y - \overline {\mathcal {\varepsilon }}^Y_\ell ,\\\notag \end{align}

(4.3) \begin{align} &\overline {\rho } \widetilde {\frac {D}{Dt}}\frac {1}{2}|\nabla \widetilde {Y}_\ell |^2 =- \nabla \cdot \left [\nabla \widetilde {Y} \nabla \cdot \left (\overline {\rho }\, \widetilde {\tau }(\boldsymbol {u}, Y)\right )\right ]+\overline {\rho } \widetilde {\mathcal {S}}^Y_\ell - \widetilde {\Pi }_\ell ^{\nabla Y} - \overline {\mathcal {\varepsilon }}^{\nabla Y}_\ell . \end{align}

Here, $\widetilde {({D}/{Dt}})\equiv ({\partial }/{\partial t}) + \widetilde {\boldsymbol {u}}_\ell \cdot \nabla$ is the coarse-grained material derivative, $\widetilde {\Pi }^Y_\ell$ and $\widetilde {\Pi }_\ell ^{\nabla Y}$ are the fluxes of scalar variance and scalar gradient, $\overline {\mathcal {\varepsilon }}^Y_\ell$ and $\overline {\mathcal {\varepsilon }}^{\nabla Y}_\ell$ are the coarse-grained dissipation rate of scalar and scalar gradient, $\widetilde {\mathcal {S}}^Y_\ell$ is the filtered scalar element stretching and $\widetilde {\boldsymbol \tau }(\boldsymbol {u},Y)$ is the generalised second-order moment. These terms are defined as

(4.4) \begin{align} \begin{split} \widetilde {\Pi }^Y_\ell &= -\overline {\rho }\nabla \widetilde {Y}\cdot \widetilde {\boldsymbol \tau }(\boldsymbol {u}, Y), \;\; \widetilde {\Pi }_\ell ^{\nabla Y}=-\frac {1}{\overline {\rho }}\nabla \cdot \left (\overline {\rho }\, \widetilde {\boldsymbol \tau }(\boldsymbol {u}, Y)\right ) \nabla \cdot \left (\overline {\rho }\nabla \widetilde {Y}\right ), \; \\ \overline {\mathcal {\varepsilon }}^Y_\ell &= \nabla \widetilde {Y}\cdot \overline {\rho D \nabla Y}, \;\; \overline {\mathcal {\varepsilon }}^{\nabla Y}_\ell = -\overline {\rho } \nabla \widetilde {Y} \cdot \nabla \left (\frac {1}{\overline {\rho }}\nabla \cdot \overline {\rho D\nabla Y}\right ), \; \\ \widetilde {\mathcal {S}}^Y_\ell &= -\nabla \widetilde {Y}\cdot \nabla \widetilde {\boldsymbol {u}}\cdot \nabla \widetilde {Y}, \;\; \widetilde {\boldsymbol \tau }(\boldsymbol {u}, Y) = \widetilde {Y \boldsymbol {u}} - \widetilde {Y}\widetilde {\boldsymbol {u}}. \end{split} \end{align}

The coarse-grained scalar stretching $\widetilde {\mathcal {S}}_\ell ^Y$ approaches the fine-grained scalar stretching $\mathcal {S}^Y$ defined in (3.4) as $\ell$ approaches zero, i.e. $\lim _{\ell \to 0} \widetilde {\mathcal {S}}_\ell ^Y = \mathcal {S}^Y$ .

Figure 8.

(a) The spatial mean values of the scalar variance flux $\widetilde {\Pi }_\ell ^Y$ with respect to the filtering width $\ell$ , normalised by the Kolmogorov scale $\eta$ for both the low-At (solid lines) and mid-At (dashed lines) cases over the time interval $\widehat {t} = 2$ -- $5$ . (b) The mean coarse-grained scalar dissipation over scale, $\overline {\mathcal {\varepsilon }}^Y_\ell$ .

Figure 8(a) illustrates the mean values of the scalar variance flux as functions of the filtering width $\ell$ at different time instants. All flux terms are positive on average, indicating an overall forward transfer of scalar variance from large to small scales. The transfer shows a steady increase over time ( $\widehat {t} = 2-5$ ) as RT turbulence develops, and the length scale associated with the $\widetilde {\Pi }_\ell ^Y$ peak shifts towards larger scales. Since $\widetilde {\Pi }^Y_\ell \equiv -\overline {\rho }\nabla \widetilde {Y}\cdot \widetilde {\boldsymbol \tau }(\boldsymbol {u},Y)$ depends on both the concentration and the velocity fields, the spectra in figure 7 demonstrate that the shift of the length scale associated with the peak $\widetilde {\Pi }^Y_\ell$ is mainly influenced by the velocity field. Figure 8(b) shows that the dissipation of scalar variance occurs at small length scales.

The influence of velocity fluctuations on the scalar variance flux can be manifested in its connection with the coarse-grained velocity gradient. We now explore the relationship of $\widetilde {\Pi }_\ell ^Y$ with the filtered velocity gradient, specifically focusing on both the strain rate and vorticity fields. The joint PDF between $\widetilde {\Pi }_\ell ^Y$ and the magnitude of the filtered strain rate $|\overline {\mathbf {S}}_\ell |$ at $\ell = 64 \eta$ , as shown in figure 9(a), reveals a positive correlation between these two quantities. In contrast, the joint PDF between $\widetilde {\Pi }_\ell ^Y$ and the magnitude of the filtered vorticity $|\overline {\boldsymbol {\omega }}_\ell |$ in figure 9(b) shows only a very weak correlation. These trends are consistent across different filtering widths, as illustrated in the conditional mean plots in figures 9(c) and 9(d). The average $\widetilde {\Pi }_\ell ^Y$ conditioned on $|\overline {\mathbf {S}}_\ell |$ is a monotonically increasing function across all cases, whereas when conditioned on $|\overline {\boldsymbol \omega }_\ell |$ , the results are not monotonic and display significantly lower magnitudes. Overall, the findings in figure 9 suggest that the scalar variance flux $\widetilde {\Pi }_\ell ^Y$ is predominantly influenced by the symmetric part of the filtered velocity gradient (strain rate field) rather than the anti-symmetric part (vorticity field).

Figure 9.

The correlation between scalar variance flux $\widetilde {\Pi }_\ell ^Y$ and the magnitude of the filtered strain rate and vorticity fields. Panels (a) and (b) show the joint PDF of $\widetilde {\Pi }_\ell ^Y$ with strain rate magnitude $|\overline {\boldsymbol {S}}|$ and vorticity magnitude $|\overline {\boldsymbol \omega }|$ , respectively, for the low-At case at filtering width $\ell =64\eta$ and at $\widehat {t}=5$ . The black dashed lines are the conditional averaged scalar variance flux. Panels (c) and (d) are the conditional averaged $\widetilde {\Pi }_\ell ^Y$ at three filtering widths over $|\overline {\boldsymbol {S}}|$ and $|\overline {\boldsymbol \omega }|$ , respectively. Both the low-At (solid lines) and the mid-At (dashed lines) cases at $\widehat {t}=5$ are included.

This behaviour can be well explained by the nonlinear model (Borue & Orszag Reference Borue and Orszag1998; Meneveau & Katz Reference Meneveau and Katz2000b ; Lees & Aluie Reference Lees and Aluie2019). For the moment, we equate the Favre-averaged quantities with plain spatial-averaged quantities, i.e. $\widetilde {f}_\ell \approx \overline {f}_\ell$ , which incurs an error of the order $O(\delta \rho /\rho )\propto O(\mathrm {At})$ (Eyink & Drivas Reference Eyink and Drivas2018). Based on this approximation, we have $\widetilde {\boldsymbol \tau }(\boldsymbol u, Y) \approx \overline {\boldsymbol \tau }(\boldsymbol u, Y)$ . By further invoking the scale-locality property (Eyink Reference Eyink2005) for scalar fields, $\overline {\boldsymbol \tau }_\ell (\boldsymbol u, Y)\approx \overline {\boldsymbol \tau }_\ell (\overline {\boldsymbol u}_\ell , \overline {Y}_\ell )$ , and the exact identity connecting the subscale flux term for any two fields $f(\boldsymbol {x}), g(\boldsymbol {x})$ and their increments (Constantin, Weinan & Titi Reference Constantin, Weinan and Titi1994)

(4.5) \begin{align} \overline {\tau }_\ell (f,g) = \int d^3{\boldsymbol {r}} G_\ell ({\boldsymbol {r}}) \delta f(\boldsymbol {x}; {\boldsymbol {r}}) \delta g(\boldsymbol {x};{\boldsymbol {r}}) - \int d^3{\boldsymbol {r}} G_\ell ({\boldsymbol {r}})\delta f(\boldsymbol {x};{\boldsymbol {r}})\int d^3{\boldsymbol {r}} G_\ell ({\boldsymbol {r}}) \delta g(\boldsymbol {x};{\boldsymbol {r}}), \end{align}

where $\delta f(\boldsymbol {x};{\boldsymbol {r}}) = f(\boldsymbol {x}+{\boldsymbol {r}})-f(\boldsymbol {x})$ and $\delta g(\boldsymbol {x};{\boldsymbol {r}}) = g(\boldsymbol {x}+{\boldsymbol {r}})-g(\boldsymbol {x})$ , we have

(4.6) \begin{align} \begin{split} \widetilde {\Pi }_\ell ^Y &= -\overline {\rho }\nabla \widetilde {Y}\cdot \widetilde {\tau }(\boldsymbol u, Y) \approx -\overline {\rho }\nabla \overline {Y}\cdot \overline {\tau }(\overline {Y}, \overline {\boldsymbol {u}})\\ &= -\overline {\rho }\nabla \overline {Y}\cdot \left (\int d^3{\boldsymbol {r}} G_\ell ({\boldsymbol {r}}) \delta \overline {Y} \delta \overline {\boldsymbol {u}} - \int d^3{\boldsymbol {r}} G_\ell ({\boldsymbol {r}}) \delta \overline {Y} \int d^3{\boldsymbol {r}} G_\ell ({\boldsymbol {r}})\delta \overline {\boldsymbol {u}}\right ) \\ &\approx -\overline {\rho }\nabla \overline {Y}\cdot \left (\int d^3{\boldsymbol {r}} G_\ell ({\boldsymbol {r}}) {\boldsymbol {r}} \cdot \nabla \overline {Y}\ {\boldsymbol {r}}\cdot \nabla \overline {\boldsymbol {u}} - \int d^3{\boldsymbol {r}} G_\ell ({\boldsymbol {r}}) {\boldsymbol {r}} \cdot \nabla \overline {Y} \int d^3{\boldsymbol {r}} G_\ell ({\boldsymbol {r}}) {\boldsymbol {r}} \cdot \nabla \overline {\boldsymbol {u}} \right ) \\ &= -\overline {\rho }\nabla \overline {Y}\cdot \ell ^2 \partial _i \overline {Y}\partial _k \overline {\boldsymbol {u}} \int d^3\left (\frac {{\boldsymbol {r}}}{\ell }\right ) G\left (\frac {{\boldsymbol {r}}}{\ell }\right ) \frac {r_i}{\ell }\frac {r_k}{\ell } \\ &= -C_0\ell ^2 \overline {\rho } \nabla \overline {Y} \cdot \nabla \overline {\boldsymbol {u}}\cdot \nabla \overline {Y} \equiv \widetilde {\Pi }_{\mathrm {NL},\ell }^Y, \end{split} \end{align}

where between the second the third line we have adopted the Taylor expansion $\delta \overline {Y} \approx {\boldsymbol {r}} \cdot \nabla \overline {Y}, \delta \overline {\boldsymbol {u}} \approx {\boldsymbol {r}}\cdot \nabla \overline {\boldsymbol {u}}$ to further approximate the expression, and the coefficient $C_0=\int d^3{\boldsymbol {r}} G({\boldsymbol {r}})|{\boldsymbol {r}}|^2$ is the second moment of the filtering kernel and, hence, is a constant. For the Gaussian kernel adopted here in (2.6), we have $C_0=1/12$ .

The nonlinear model suggests that the scalar variance flux $\widetilde {\Pi }_\ell ^Y$ across scale $\ell$ explicitly depends only on the strain rate tensor filtered at $\ell$ , with no direct dependence on the vorticity field. This is consistent with the conditional averaged results presented in figure 9. In deriving this model, several simplifications were made that could affect its performance. We now evaluate the model by comparing its prediction with the actual data.

Figure 10.

The joint PDF between $\widetilde {\Pi }_\ell ^Y$ and its nonlinear model $\widetilde {\Pi }^Y_{\ell ,\mathrm {NL}}$ for the low-At case at $\widehat {t}=5$ , evaluated at two widths ( $a$ ) $\ell =16\eta$ and ( $b$ ) $\ell =64\eta$ . An auxiliary diagonal dashed line is included for reference. ( $c$ ) The correlation coefficients between $\widetilde {\Pi }_\ell ^Y$ and its nonlinear model for both the low-At (solid lines) and mid-At (dashed lines) cases at $\widehat {t}=2-5$ . ( $d$ ) The correlation coefficients between $\widetilde {\Pi }_\ell ^Y$ and $\overline {\mathcal {\varepsilon }}^Y_\ell$ .

Figures 10(a) and 10(b) show the joint PDF of $\widetilde {\Pi }_\ell ^Y$ and its nonlinear model $\widetilde {\Pi }^Y_{\text {NL}} = -({1}/{12})\ell ^2\overline {\rho } \nabla \overline {Y} \cdot \nabla \overline {\boldsymbol {u}} \cdot \nabla \overline {Y}$ for the low-At case at $\widehat {t}=5$ . The joint PDFs for both small and intermediate length scales, $\ell = 16\eta$ and $\ell = 64\eta$ , align closely with the diagonal line, with correlation coefficients of 0.98 and 0.94, respectively. This indicates that the nonlinear model provides a strong pointwise approximation of the actual scalar variance flux.

Figure 10(c) further illustrates the correlation coefficient between the model and actual flux at different times ( $\widehat {t} = 2-5$ ) for both the low-At and mid-At cases. The results show that the correlation improves over time as the RT turbulence develops, suggesting that the model performs better at higher Reynolds numbers where the scale-locality assumption is applicable to a wider scale range. In addition, the correlation coefficient is higher at smaller scales, where the Taylor expansion adopted in deriving the nonlinear model is more accurate. Additionally, because of the approximation that replaces Favre-filtered quantities with plain filtered quantities, $\widetilde {f}(\boldsymbol {x}) \to \overline {f}(\boldsymbol {x})$ , the low-At case shows slightly better correlation than the mid-At case at small to intermediate scales, due to the smaller error of the order $O(\mathrm {At})$ . Interestingly, the correlation coefficients for the mid-At case at the largest filtering widths are slightly better than those in the low-At case, despite the overall relatively weak correlation. This trend is likely due to the higher Reynolds number in the mid-At case than the low-At case, as indicated in table 1, leading to a slightly better performance of the model. In addition, the time evolution of the mean values of $\widetilde {\Pi }_\ell ^Y$ and $\widetilde {\Pi }^Y_{\text {NL}}$ , along with their correlation coefficients, are illustrated in Appendix D. These results again demonstrate the good approximation of the nonlinear model within the self-similar regime of RT turbulence.

Figure 11.

Schematic of the pathway of scalar variance and its gradient transfer across scales.

The nonlinear model $\widetilde {\Pi }^Y_{\ell ,\mathrm {NL}}$ bears a strong resemblance to the filtered scalar stretching term $\widetilde {\mathcal {S}}^Y_\ell = -\partial _i \widetilde {Y}\partial _j \widetilde {Y}\partial _i \widetilde {u}_j \approx -\partial _i \overline {Y} \partial _j\overline {Y}\partial _i \overline {u}_j = \widetilde {\Pi }_{\ell ,\mathrm {NL}}^Y/(C_0\ell ^2 \overline {\rho })$ . This similarity indicates a shared physical mechanism underlying both the scalar variance flux, a sink term for large-scale scalar variance in (4.2), and the filtered scalar stretching, which acts as a source term in the coarse-grained scalar gradient equation (4.3). Given that the preferential alignment between the scalar gradient and the strain eigenbasis also holds for coarse-grained fields (see figure 13 a), the estimation in (3.7) that connects the fine-grained stretching term and the scalar dissipation rate is valid for the filtered results as well, i.e.

(4.7) \begin{align} \widetilde {\Pi }^Y_\ell \approx \widetilde {\Pi }^Y_{\ell , \mathrm {NL}} \propto \widetilde {\mathcal {S}}^Y_\ell \approx -\overline {\lambda }_\gamma |\nabla \overline {Y}_\ell |^2 \approx \frac {-\overline {\lambda }_\gamma }{\overline {\rho }_\ell D}\overline {\mathcal {\varepsilon }}_\ell ^Y \quad \mathrm {given }\; \nabla \overline {Y}_\ell \parallel \boldsymbol {e}_{\gamma }, \end{align}

where $\overline {\lambda }_\gamma$ is the eigenvalue corresponding to the most compressive eigenvector of the filtered strain field. The correlation between $\widetilde {\Pi }^Y_\ell$ and the large-scale scalar dissipation $\overline {\mathcal {\varepsilon }}_\ell ^Y$ is illustrated in figure 10(d), where the correlation coefficients are above 0.7 at small scales, confirming the validity of (4.7). This relation reflects a form of turbulent eddy diffusivity commonly adopted in RANS or LES. Note that in figure 10(d) the unexpectedly high correlation between $\widetilde {\Pi }^Y_\ell$ and $\overline {\mathcal {\varepsilon }}_\ell ^Y$ at large filtering scales arises from the contribution of long-wavelength modes in RT turbulence, generated by the bubble merger process. The vertical velocity and mass fraction field of these modes combine to produce a second-order moment $\widetilde {\tau }_\ell (\boldsymbol {u}, \boldsymbol {Y})$ that resembles the filtered mass fraction gradient, resulting in a strong correlation $C(\widetilde {\Pi }^Y_\ell , \overline {\mathcal {\varepsilon }}_\ell ^Y) \approx 1$ .

The pathways of scalar variance and scalar gradient transfer across scales and the underlying physical mechanism is depicted in figure 11. At intermediate to large scales ( $\gtrapprox \ell$ ), turbulent stirring by straining motion disperses the bulk scalar field at large scale into sheet structures associated with smaller length scales, thus transferring scalar variance from large to small scales. Therefore, the strain-induced surface stretching $\widetilde {\mathcal {S}}_\ell ^Y$ is closely related to the transfer term $\widetilde {\Pi }_\ell ^Y$ . On the other hand, the surface stretching $\widetilde {\mathcal {S}}_\ell ^Y$ acts on the tangential plane of scalar isosurfaces, accompanied by a compression in the normal direction, which brings adjacent isoscalar elements closer together and increases the coarse-grained scalar gradient $\nabla \overline {Y}_\ell$ . This effect further facilitates the dissipation rate $\overline {\mathcal {\varepsilon }}^Y_\ell$ associated with the filtered concentration field. The close proximity among $\widetilde {\Pi }_\ell ^Y \sim \widetilde {\mathcal {S}}_\ell ^Y \sim \overline {\mathcal {\varepsilon }}^Y_\ell$ is thus established. Additionally, since $\lim _{\ell \to 0} \widetilde {\mathcal {S}}_\ell ^Y = \mathcal {S}^Y$ and $\lim _{\ell \to 0}\overline {\mathcal {\varepsilon }}^Y_\ell = \mathcal {\varepsilon }^Y$ , the fine-grained quantities $\mathcal {S}^Y$ and $\mathcal {\varepsilon }^Y$ also bear similarities, as shown earlier,

4.3. Backscatter of the scalar variance flux

As observed in figures 9(a) and 10(a,b), while the mean value of $\widetilde {\Pi }_\ell ^Y$ is positive, causing a globally forward transfer of the scalar variance flux from large to small scales, locally negative values of $\widetilde {\Pi }_\ell ^Y$ do exist, particularly at small filtering widths $\ell$ , representing a backward transfer from small to large scales. Figure 12(a) illustrates the mean values of the positive and negative components of $\widetilde {\Pi }_\ell ^Y$ , showing that their magnitudes increase over time for both the low-At and mid-At cases. The positive component, $\widetilde {\Pi }^Y_{\ell ,+}$ , is significant across almost all scale ranges, whereas the negative component, $\widetilde {\Pi }^Y_{\ell ,-}$ , is only prominent at relatively small length scales ( $\ell \lt 160\eta$ ) and peaks around $16 \eta$ . These results indicate that inverse scalar variance transfer predominantly occurs at small scales. Although the mean value of the negative subset of $\widetilde {\Pi }_\ell ^Y$ is small compared with the positive subset, the proportion of the number of $\widetilde {\Pi }^Y_{\ell ,-}$ samples depicted in figure 12(a) shows that the negative region could exceed 40 % at intermediate to small scales. This suggests that scalar variance backscatter has a widespread influence.

Figure 12.

( $a$ ) The mean positive and negative components of the scalar variance flux, with the subscripts $\langle \cdot \rangle _+$ and $\langle \cdot \rangle _-$ indicating averages over its positive and negative subsets, respectively. The black lines represent the proportion of the negative volume of $\widetilde {\Pi }_\ell ^Y$ within the domain at $\widehat {t}=5$ . ( $b$ ) The conditional mean value $\langle \widetilde {\Pi }_\ell ^Y | |\overline {\mathbf {S}}_\ell |, |\overline {\boldsymbol {\omega }}_\ell | \rangle$ on the joint PDF of filtered strain and vorticity magnitudes of the low-At case at $\widehat {t}=5$ , with a filtering width $\ell =16 \eta$ . ( $c, d$ ) The conditional mean $\langle \widetilde {\Pi }_\ell ^Y | \overline {Q}^*_\ell , \overline {R}^*_\ell \rangle$ for the low-At and mid-At cases at $\widehat {t}=5$ with the filtering width $\ell =16\eta$ . The second and third velocity gradient invariants $\overline {Q}_\ell$ and $\overline {R}_\ell$ are normalised by $\langle |\overline {\boldsymbol {\omega }}_\ell |^2\rangle$ and $\langle |\overline {\boldsymbol {\omega }}_\ell |^2\rangle ^{3/2}$ , respectively.

To gain insight into the spatial distribution and topology of the backscatter of scalar variance flux, we obtain the mean value of $\widetilde {\Pi }_\ell ^Y$ conditional on the $|\overline {\mathbf {S}}_\ell |$ − $|\overline {\boldsymbol {\omega }}_\ell |$ space and the $\overline {Q}_\ell$ − $\overline {R}_\ell$ space in figure 12(b)–(d), where $\overline {Q}_\ell \equiv -({1}/{2})\partial _k \overline {u}_i\partial _i \overline {u}_k$ and $\overline {R}_\ell \equiv -({1}/{3})\partial _j \overline {u}_i\partial _k\overline {u}_j\partial _i\overline {u}_k$ are the second and third invariants of the filtered velocity gradients. In the $|\overline {\mathbf {S}}_\ell |$ − $|\overline {\boldsymbol {\omega }}_\ell |$ space, the negative $\widetilde {\Pi }_\ell ^Y$ occurs at regions with high-vorticity magnitude and low to intermediate strain magnitude, implying that vorticity plays an active role in the scalar variance backscatter. Figure 12(b) shows the result for the low-At case, with similar results observed for the mid-At case (not shown). In the $\overline {Q}_\ell$ − $\overline {R}_\ell$ space shown in figures 12(c) and 12(d), the conditional joint PDFs clearly show that, for both the low-At and mid-At cases, negative mean values of $\widetilde {\Pi }_\ell ^Y$ primarily occur in the first quadrant of the $\overline {Q}_\ell{-}\overline {R}_\ell$ plane, where $\overline {Q}_\ell \gt 0$ and $\overline {R}_\ell \gt 0$ . This region, located above the dashed curve on the right (the Vieillefosse tail, ${27}/{4} \overline {R}_\ell ^2 + \overline {Q}_\ell ^3 = 0$ ), corresponds to the local flow pattern of vortex compression (Meneveau Reference Meneveau2011), and that local vorticity strength dominates over strain ( $\overline {Q}_\ell \gt 0$ ). Thus, the vorticity vector influences the backscatter region with $\widetilde {\Pi }_\ell ^Y\lt 0$ , although its effect on the overall flux term $\widetilde {\Pi }_\ell$ is diminished, as indicated in the conditional mean of figure 9(d).

The negative scalar variance transfer is due to a modified preferential alignment between the filtered scalar gradient and the filtered strain rate eigenvectors. Due to the proximity between scalar variance flux and surface stretching, we have

(4.8) \begin{align} \begin{split} \widetilde {\Pi }_\ell ^Y \propto \widetilde {\mathcal {S}}^Y &=-\partial _i \widetilde {Y}\partial _j \widetilde {Y}\widetilde {\mathrm S}_{ij} \approx -\partial _i \overline {Y}\partial _j \overline {Y}\overline {\mathrm S}_{ij}\\ &= -|\nabla \overline {Y}|^2\left (\overline {\lambda }_\alpha \cos ^2(\nabla \overline {Y}, \boldsymbol {e}_{\alpha })+ \overline {\lambda }_\beta \cos ^2(\nabla \overline {Y}, \boldsymbol {e}_{\beta })+ \overline {\lambda }_\gamma \cos ^2(\nabla \overline {Y}, \boldsymbol {e}_{\gamma })\right ), \end{split} \end{align}

where $\overline {\lambda }_\alpha , \overline {\lambda }_\beta , \overline {\lambda }_\gamma$ and $\boldsymbol {e}_{\alpha }, \boldsymbol {e}_{\beta }, \boldsymbol {e}_{\gamma }$ correspond to the eigenvalues and eigenvectors of the filtered strain rate tensor $\overline {\boldsymbol {S}}$ , with $\overline {\lambda }_\alpha \gt 0$ and $\overline {\lambda }_\gamma \lt 0$ . Thus, if $\nabla \overline {Y}$ predominantly aligns with the most compressive eigenvector $\boldsymbol {e}_{\gamma }$ then $\widetilde {\Pi }_\ell ^Y \propto -\overline {\lambda }_\gamma |\nabla \overline {Y}|^2 \cos (\nabla \overline {Y}, \boldsymbol {e}_{\gamma })^2 \gt 0$ . On the other hand, if $\nabla \overline {Y}$ predominantly aligns with the most expansive eigenvector $\boldsymbol {e}_{\alpha }$ then $\widetilde {\Pi }_\ell ^Y \propto -\overline {\lambda }_\alpha |\nabla \overline {Y}|^2 \cos (\nabla \overline {Y}, \boldsymbol {e}_{\alpha })^2 \lt 0$ .

We have already shown in the scalar normal budgets (figure 6) that, for the unfiltered fields, the strain rate tensor tends to reinforce the preferential alignment between $\nabla Y$ and $\boldsymbol {e}_\gamma$ , while the vorticity vector tends to oppose this trend and promote the alignment between $\nabla Y$ and $\boldsymbol {e}_\alpha$ . This trend persists to the filtered fields, especially at small filtering scales where the scalar variance backscatter is most intense. Therefore, the vorticity-induced rotation of the scalar normal generates regions where $\nabla \overline {Y}$ is parallel to $\boldsymbol {e}_{\alpha }$ , and hence, $\widetilde {\Pi }_\ell ^Y\lt 0$ . Figures 13(a) and 13(b) show the PDFs of cosine values between $\nabla \overline {Y}_\ell$ and strain eigenvectors $\boldsymbol {e}$ filtered at $\ell =64\eta$ and conditioned on both positive and negative $\widetilde {\Pi }^Y_\ell$ , respectively (comparable results are observed with other filtering widths). The PDF conditional on the forward flux in panel (a) closely resembles the unconditioned results shown in figure 5(a), where $\nabla \overline {Y}$ aligns with the most compressive eigenvector $\boldsymbol {e}_\gamma$ . However, when conditioned on the backscatter with $\widetilde {\Pi }^Y_\ell \lt 0$ , the best alignment occurs between $\nabla \overline {Y}$ and the most expansive eigenvector $\boldsymbol {e}_\alpha$ . In summary, the vorticity vector rotates the scalar isosurface, generating subregions with $\nabla \overline {Y} \parallel \boldsymbol {e}_{\alpha }$ , and therefore, leads to the contraction of concentration isosurface and a negative scalar variance transfer.

Figure 13.

The alignment between the scalar gradient and the strain rate eigenvectors at $\widehat {t}=5$ , both filtered at $\ell =64\eta$ . Panel (a) is conditional on the positive subset of $\widetilde {\Pi }_\ell ^Y$ , while panel (b) is conditional on the negative subset. The results for both the low-At (solid lines) and mid-At (dashed lines) cases are similar.