2.1. Energy cascade rate/flux in the scale-integrated local KH equation

The KH equation is a generalized Kármán–Howarth equation that is derived directly from the incompressible Navier–Stokes equations without any modelling. Before averaging, the instantaneous KH equation with no mean flow and neglecting the forcing term reads (Hill Reference Hill2001, Reference Hill2002)

(2.1)\begin{equation} \frac{\partial \delta u _i^2}{\partial t} + u^*_{j}\,\frac{\partial \delta u _i^2}{\partial x_j} =-\frac{\partial \delta u _j\,\delta u _i^2}{\partial r_j}-\frac{8}{\rho}\,\frac{\partial p^*\,\delta u _i}{\partial r_i} +\nu\,\frac{1}{2}\,\frac{\partial^2 \delta u _i^2}{\partial x_j\,\partial x_j }+ 2\nu\,\frac{\partial^2 \delta u _i^2}{\partial r_j\,\partial r_j} - 4\epsilon^*, \end{equation}

where $\delta u_i = \delta u_i({\boldsymbol {x}},{\boldsymbol {r}}) = u_i^+ - u_i^-$ is the velocity increment vector in the $i{\textrm {th}}$ Cartesian direction over displacement vector ${\boldsymbol {r}}$. The superscripts $+$ and $-$ represent two points, ${\boldsymbol {x}}+{\boldsymbol {r}}/2$ and ${\boldsymbol {x}}-{\boldsymbol {r}}/2$, in the physical domain that have a separation vector $r_i = x^+_i - x^-_i$ and middle point $x_i = (x^+_i + x^-_i)/2$ (see figure 1a). The superscript $*$ denotes the average value between two points, e.g. the two-point average dissipation is defined as $\epsilon ^*({\boldsymbol {x}},{\boldsymbol {r}}) = (\epsilon ^+ +\epsilon ^-)/2$, and $\epsilon ^\pm$ here is the ‘pseudo-dissipation’ defined at every point as $\epsilon =\nu ({\partial u_i}/{\partial x_j})^2$. (In Hill (Reference Hill2002), an alternate expression involving the real dissipation was introduced – his (2.13) – at the cost of including an additional pressure term.) Note that throughout this paper, when referring to ‘dissipation’ we will mean the pseudo-dissipation. Also, we will use ${\boldsymbol {r}}_s = {\boldsymbol {r}}/2$ to denote the radial coordinate vector from the local ‘origin’ ${\boldsymbol {x}}$.

Figure 1.

(a) Sketch showing local domain of integration over a ball of diameter $\ell$ used in the symmetric Hill (Reference Hill2002) structure-function approach in which pairs of points separated by distances $r=2r_s$ up to $\ell$ are used. (b) Integration up to a ball of radius $\ell$ in which pairs of points separated by distances $r$ up to $\ell$ are used as in the approach of Duchon & Robert (Reference Duchon and Robert2000). For volume averaging, in (a) 3-D integration over the vector ${\boldsymbol {r}}_s$ is performed at fixed ${\boldsymbol {x}}$, while in (b) 3-D integration over the vector ${\boldsymbol {r}}$ is performed at fixed ${\boldsymbol {x}}$. For surface integrations, in (a) integration is done over the spherical surface of radius $\ell /2$, while in (b) it is done over a spherical surface of radius $\ell$.

As remarked by Hill (Reference Hill2001, Reference Hill2002), it is then instructive to apply integration over a sphere in ${\boldsymbol {r}}_s$-space up to a radius $\ell /2$, i.e. over a sphere of diameter $\ell$. The resulting equation is divided by the sphere volume $V_\ell =\frac {4}{3}{\rm \pi} ( {\ell }/{2})^3$ and a factor 4, and Gauss’ theorem is used for the ${\boldsymbol {r}}$-divergence terms (recalling that $\partial {\boldsymbol {r}}=2\,\partial {\boldsymbol {r}}_s$), yielding

(2.2)\begin{align} &\frac{1}{2 V_\ell}\mathop{\iiint}\limits_{V_{\ell}} \left(\frac{\partial \delta u_i^2/2}{\partial t} + u^*_{j}\, \frac{\partial \delta u _i^2/2}{\partial x_j} \right) {\rm d}^3{\boldsymbol{r}}_s\nonumber\\ &\quad =-\frac{3}{4\ell}\,\frac{1}{S_\ell}\oint_{S_{\ell}} \delta u _i^2\,\delta u _j\,\hat{n}_j \,{\rm d}S - \frac{6}{\rho \ell}\,\frac{1}{S_\ell}\oint_{S_{\ell}}p^*\,\delta u _j\,\hat{n}_j \,{\rm d}S \nonumber\\ &\qquad + \frac{\nu}{4}\,\frac{1}{V_\ell}\iiint\limits_{V_{\ell}}\left(\frac{1}{2}\,\frac{\partial^2 \delta u_i^2}{\partial x_j\,\partial x_j } + 2\,\frac{\partial^2 \delta u_i^2}{\partial r_j\,\partial r_j} \right) {\rm d}^3{\boldsymbol{r}}_s - \frac{1}{V_\ell}\iiint\limits_{V_{\ell}} \epsilon^* \,{\rm d}^3{\boldsymbol{r}}_s, \end{align}

where $S_\ell$ represents the bounding sphere's surface of area $S_\ell = 4{\rm \pi} (\ell /2)^2$, and $\hat {n}_j$ is the radial unit vector normal to the sphere surface. Equation (2.2) suggests defining a structure-function-based kinetic energy at scale $\ell$ according to

(2.3)\begin{equation} k_{{sf},\ell}({\boldsymbol{x}},t) = \frac{1}{2 V_\ell}\iiint\limits_{V_{\ell}} \frac{1}{2}\,\delta u _i^2({\boldsymbol{x}},{\boldsymbol{r}}) \, {\rm d}^3{\boldsymbol{r}}_s, \end{equation}

so that the first term in (2.2) corresponds to $\partial k_{{sf},\ell }/\partial t$. The $1/2$ factor in front of the integral is justified since the volume integration over the entire sphere will double count the energy contained in $\delta u _i^2 = (u_i^+ - u_i^-)^2$. Equation (2.2) thus describes the transport of two-point, structure-function energy $k_{{sf},\ell }$, which represents energy within eddies with length scales up to $\ell$ (Davidson Reference Davidson2015) in both the length scale $\ell$ and physical position $\boldsymbol {x}$ spaces. The last term in (2.2) represents the $r$-averaged rate of dissipation, with the radius vector ${\boldsymbol {r}}_s = {\boldsymbol {r}}/2$ being integrated up to magnitude $\ell /2$:

(2.4)\begin{equation} \epsilon_\ell({\boldsymbol{x}},t) \equiv \frac{1}{V_\ell}\iiint\limits_{V_{\ell}} \epsilon^*({\boldsymbol{x}},{\boldsymbol{r}}) \, {\rm d}^3{\boldsymbol{r}}_s. \end{equation}

As remarked by Hill (Reference Hill2001, Reference Hill2002), this quantity corresponds directly to the spherical average of local dissipation at scale $\ell$, and plays a central role in the celebrated KRSH (Kolmogorov Reference Kolmogorov1962).

The local energy cascade rate in the inertial range at position ${\boldsymbol {x}}$ and time $t$ is defined as

(2.5)\begin{equation} \varPhi_\ell({\boldsymbol{x}},t) \equiv-\frac{3}{4\ell}\,\frac{1}{S_\ell}\oint_{S_{\ell}} \delta u _i^2\,\delta u _j\,\hat{n}_j \, {\rm d}S =-\frac{3}{4\ell}\,[\delta u_i^2\,\delta u_j\,\hat{n}_j ]_{S_\ell}, \end{equation}

where $[\cdot ]_{S_\ell }$ indicates area averaging over the sphere of diameter $\ell$. We note that in this definition, $\varPhi _\ell ({\boldsymbol {x}},t)$ represents the surface average of a flux that is defined positive if energy is flowing into the sphere in the ${\boldsymbol {r}}$-scale space. The position is fixed at ${\boldsymbol {x}}$, thus the quantity $\varPhi _\ell ({\boldsymbol {x}},t)$ does not contain possible confounding spatial transport effects.

In terms of the overall average of (2.2), under the assumptions of homogeneous isotropic flow and statistical steady-state conditions, and for $\ell$ in the inertial range of turbulence, the unsteady transport and viscous terms vanish. The pressure term is also zero due to isotropy and incompressibility. Therefore, (2.2) can be simplified and yields, as expected,

(2.6)\begin{equation} \langle \varPhi_\ell \rangle = \langle \epsilon_\ell \rangle=\langle \epsilon \rangle, \end{equation}

or equivalently, $[\delta u_i^2\,\delta u_j\,\hat {n}_j ]_{S_\ell }=-(4/3)\ell \langle \epsilon \rangle$, the $4/3$ law (Frisch Reference Frisch1995).

In this paper, the focus will be mainly on the flux term $\varPhi _\ell$, with some attention also on the dissipation term $\epsilon _\ell$. Analysis of the time derivative, spatial advection terms and pressure terms is left for other ongoing studies. The viscous flux terms (in both spatial and scale spaces) are also not considered, since our present interest concerns the inertial range.

2.2. Energy cascade rate/flux in the filtering approach

In this subsection, we review the transport equation of the SGS kinetic energy (Germano Reference Germano1992) for $k_{{sgs},\ell } \equiv \frac {1}{2} \tau _{ii}$, where $\tau _{ij}=\widetilde {u_i u_j} - \tilde {u}_i \tilde {u}_j$ is the SGS stress tensor, with the tilde symbol $(\sim )$ denoting spatial filtering of variables. The transport equation for $k_{{sgs},\ell }$ reads (Germano Reference Germano1992)

(2.7)\begin{align} \frac{\partial k_{{sgs},\ell}}{\partial t} + \tilde{u}_j\,\frac{\partial k_{{sgs},\ell} }{\partial x_j}&=-\frac{1}{2}\,\frac{\partial}{\partial x_j}\left(\widetilde{u_i u_i u_j}- 2\tilde{u}_i\,\widetilde{u_i u_j}-\tilde{u}_j\,\widetilde{u_i u_i}-\tilde{u}_i \tilde{u}_i \tilde{u}_j\right) -\frac{\partial}{\partial x_j}\left( \widetilde{p u_j} - \tilde{p}\tilde{u}_j \right)\nonumber\\ &\quad +\frac{\partial}{\partial x_j}\left( \nu\,\frac{\partial k_{{sgs},\ell} }{\partial x_j} \right) + \nu\,\frac{\partial \tilde{u}_i}{\partial x_j}\,\frac{\partial \tilde{u}_i}{\partial x_j} -\nu\,\widetilde{\frac{\partial u_i}{\partial x_j}\,\frac{\partial u_i}{\partial x_j}} -\tau_{ij}\tilde{S}_{ij} . \end{align}

The last term is called the SGS rate of dissipation at position $({\boldsymbol {x}})$, and is often denoted as

(2.8)\begin{equation} \varPi_\ell({\boldsymbol{x}},t) \equiv-\tau_{ij}\tilde{S}_{ij}. \end{equation}

For filtering, in the present work, we consider a spherical-shaped sharp top-hat filter in physical space with diameter $\ell$. Therefore, for any field variable $A({\boldsymbol {x}})$, we define the filtered variable as $\tilde {A}({\boldsymbol {x}}) = {V_\ell }^{-1}\iiint _{V_{\ell }} A({\boldsymbol {x}}+{\boldsymbol {r}}_s) \, {\rm d}^3{\boldsymbol {r}}_s$. Note that each term in (2.2) and (2.7) is thus evaluated at the same length scale. Terms in (2.7) can be compared directly to terms in (2.2); in particular, the local dissipation terms are exactly the same, i.e.

(2.9)\begin{equation} -\!\nu\,\widetilde{\frac{\partial u_i}{\partial x_j}\,\frac{\partial u_i}{\partial x_j}} =-\frac{1}{V_\ell}\iiint\limits_{V_{\ell}} \nu\,\frac{\partial u_i}{\partial x_j}\,\frac{\partial u_i}{\partial x_j} \,{\rm d}^3{\boldsymbol{r}}_s = \epsilon_\ell({\boldsymbol{x}},t). \end{equation}

Again, for homogeneous steady-state turbulence in the inertial range (neglecting viscous diffusion and resolved dissipation terms), upon averaging, (2.7) simplifies to

(2.10)\begin{equation} \langle \varPi_\ell \rangle = \langle \epsilon_\ell \rangle =\langle \epsilon \rangle, \end{equation}

which is similar to (2.6); thus on average, certainly the two definitions of energy cascade rate/flux agree with each other, i.e. $\langle \varPi _\ell \rangle = \langle \varPhi _\ell \rangle$.

It is also of interest to compare the average value of the two definitions of kinetic energy used in both definitions of energy cascade rate/flux. In the inertial range of high Reynolds number turbulence, both $\langle k_{{sf},\ell }\rangle$ and $\langle k_{{sgs},\ell }\rangle$ can be evaluated based on the Kolmogorov $r^{2/3}$ law and $k^{-5/3}$ spectrum, respectively. The result is (see Appendix B for details) $\langle k_{{sf},\ell } \rangle \approx 1.6 \, \langle \epsilon \rangle ^{2/3} \ell ^{2/3}$ and $\langle k_{{sgs},\ell }\rangle \approx 1.2 \, \langle \epsilon \rangle ^{2/3} \ell ^{2/3}$. In other words, they are of similar order of magnitude but the SGS kinetic energy is slightly smaller.

2.3. Other relationships between the structure-function and filtering approaches

In the present paper, we will perform the data analysis and comparisons using the two approaches mentioned above (scale-integrated local KH and filtering formulations). However, it is useful at this stage to include some remarks regarding other structure-function and energy definitions used in earlier works by Vreman, Geurts & Kuerten (Reference Vreman, Geurts and Kuerten1994), Constantin, Weinan E & Titi (Reference Constantin, Weinan and Titi1994), Duchon & Robert (Reference Duchon and Robert2000), Eyink (Reference Eyink2002) and Dubrulle (Reference Dubrulle2019). Those approaches focus typically on the structure function written at one of the endpoints instead of the midpoint. Duchon & Robert (Reference Duchon and Robert2000) and Dubrulle (Reference Dubrulle2019) focus on the two-point correlation quantity $C({\boldsymbol {x}},{\boldsymbol {r}}) = u_i({\boldsymbol {x}})\,u_i({\boldsymbol {x}} + {\boldsymbol {r}})$ (see figure 1b). Local averaging over all values of ${\boldsymbol {r}}$ from ${\boldsymbol {r}}=0$ up to scale $|{\boldsymbol {r}}|=\ell$ at any given ${\boldsymbol {x}}$ then corresponds to the ‘mixed’ energy quantity $u_i \tilde {u}_i/2$ (denoted as $E^{\ell }$ in Dubrulle Reference Dubrulle2019), where the filtering is over a sphere of diameter $2\ell$ so as to combine two points with separation distances up to $\ell$. The quantity $C({\boldsymbol {x}},{\boldsymbol {r}})$ combines filtered and unfiltered velocities, hence it is more difficult to interpret for comparisons of structure-function and LES filtering approaches. In their transport equation, Duchon & Robert (Reference Duchon and Robert2000) show that a term similar to the third-order structure-function term of (2.5) arises. However, in order for the structure function to correspond to scale $\ell$, one has to choose to integrate over a sphere of diameter $2 \ell$ (the locally integrated dissipation rate would then be $\epsilon _{2\ell }$). In a spherical integration over ${\boldsymbol {r}}$ of powers of the velocity difference $[u_i({\boldsymbol {x}}+{\boldsymbol {r}})-u_i({\boldsymbol {x}})]$, only the first term is affected by filtering or averaging over the spherical shell, while the centre velocity $u_i({\boldsymbol {x}})$ remains fully local. Note that in the scale-integrated local KH equation, the averaging affects both endpoint velocities in the same way, and both become averaged at scale $\ell$ in a formally symmetric way.

An early connection between structure-function and filtering approaches was developed by Vreman et al. (Reference Vreman, Geurts and Kuerten1994). In the Vreman analysis, the structure function is defined based on the difference of velocity $u_i({\boldsymbol {x}}+{\boldsymbol {r}})$ and the locally filtered velocity $\tilde {u}_i$ centred at ${\boldsymbol {x}}$. Spherical integration of $(u_i({\boldsymbol {x}}+{\boldsymbol {r}}) - \tilde {u}_i)^2$ over a sphere of radius $\ell /2$ then yields equivalence with the SGS kinetic energy at scale $\ell$. But $(u_i({\boldsymbol {x}}+{\boldsymbol {r}}) - \tilde {u}_i)^2$ does not equal the usual structure-function definition, now due to a mixture of filtered and unfiltered quantities at two points even before local filtering.

Another interesting approach was presented in Constantin et al. (Reference Constantin, Weinan and Titi1994) and connected to the LES filtering approach by Eyink (Reference Eyink1995, Reference Eyink2006) (equations (2.12)–(2.14) in the latter). In fact, as recounted in the review by Eyink & Sreenivasan (Reference Eyink and Sreenivasan2006), early unpublished work by Onsager anticipated such expressions half a century prior. Written in terms of the sharp spherical filter that we use here, the expression for the trace of the SGS stress reads

(2.11) \begin{align} \tau_{ii}({\boldsymbol{x}}) = \frac{1}{V_\ell} \iiint\limits_{V_\ell} [u_i({\boldsymbol{x}}+{\boldsymbol{r}})-u_i({\boldsymbol{x}})]^2 \, {\rm d}^3{\boldsymbol{r}} -\Biggl(\frac{1}{V_\ell}\iiint\limits_{V_\ell} [u_i({\boldsymbol{x}}+{\boldsymbol{r}})-u_i({\boldsymbol{x}})] \, {\rm d}^3{\boldsymbol{r}} \Biggr)^2. \end{align}

This equation represents an exact relationship between two-point structure functions and the SGS kinetic energy. But for the right-hand side to correspond to structure functions up to scale $\ell$, the integration must be done over a sphere of radius $\ell$ and thus a filtering scale of $2\ell$ for the stress tensor in the filtering formulation. The suggested relationship then appears to be between SGS stress kinetic energy at scale $2\ell$ and structure functions up to two-point separations $\ell$ but averaged over a local domain of size $2\ell$, similarly to the Duchon & Robert (Reference Duchon and Robert2000) approach. Note that while each of the terms in (2.11) is also a mixture of filtered and unfiltered velocities, the subtraction cancels the local term and restores the fully filtered property inherent in the definition of $\tau _{ii}$.

While not expecting qualitatively different results (except perhaps using the diameter instead of the radius as a name for ‘scale’), we here continue our focus on the more ‘symmetric’ formulation by Hill, with fixed position ${\boldsymbol {x}}$ specified at the midpoint between two points separated by vector ${\boldsymbol {r}}$ whose magnitude then spans up to scale $\ell$ (or integration radius ${\boldsymbol {r}}_s$ up to radius $\ell /2$).

5.1. Conditional statistics based on strain rate ($S^2_\ell$) and rotation rate ($\varOmega ^2_\ell$)

Figures 5(a,b) show distinct regions including both local forward (red area) and inverse (blue area) cascade rates. It is apparent visually that the presence of a strong local forward energy cascade is associated with increased local strain rate, as indicated by the solid black circles in both figures 5(a,b) and the corresponding black solid circles in figure 5(c). Similarly, a strong local inverse energy cascade is observed alongside a significant local rotation rate, depicted by the dashed red circles in figures 5(a,b) and the corresponding red dashed circles in figure 5(d). The strong correlation between forward cascade and local straining is consistent with multiple earlier observations and prior works in the literature (e.g. Borue & Orszag Reference Borue and Orszag1998, and recently Johnson Reference Johnson2021; Carbone & Bragg Reference Carbone and Bragg2020). We focus attention on the regions with negative energy cascade rates. Conditional averaging can elucidate the statistical significance of these regions. Specifically, we inquire whether there are large-scale flow local features as characterized by the filtered velocity gradient invariants that are systematically accompanied by inverse cascade, i.e. negative $\varPhi _\ell$. Thus we perform conditional averaging of $\varPhi _\ell$ based on the invariants $S_\ell ^2$ and $\varOmega _\ell ^2$, and repeat the analysis for the SGS energy flux quantity $\varPi _\ell$.

Figure 6 shows the joint conditionally averaged $\varPhi _\ell$ and $\varPi _\ell$ based on $S^2_\ell$ and $\varOmega ^2_\ell$, denoted as $\langle \varPhi _\ell | S^2_\ell, \varOmega ^2_\ell \rangle$ and $\langle \varPi _\ell | S^2_\ell, \varOmega ^2_\ell \rangle$, respectively. The analysis is performed by computing averages over two million randomly distributed points ${\boldsymbol {x}}$. In the presented results, $\varPhi _\ell$ is normalized by $\langle \epsilon \rangle$, while $S^2_\ell$ and $\varOmega ^2_\ell$ are normalized by $\langle Q_w \rangle = \frac {1}{2}\langle \varOmega ^2_\ell \rangle$. Figures 6(a–d) present the joint conditionally averaged $\langle \varPhi _\ell | S^2_\ell, \varOmega ^2_\ell \rangle$ at four different length scales, namely $\ell = \{30, 45, 60 , 75\} \eta$, highlighting the dominance of the forward cascade by the extensive red region. This magnitude is many times larger than the maximum magnitude observed in the blue region, representing the inverse cascade. The red region covers a wide range of $S^2_\ell$ and $\varOmega ^2_\ell$ values, consistent with the expectation that the global average would favour a forward cascade ($\langle \varPhi _\ell >0 \rangle$). The highest positive values of $\langle \varPhi _\ell | S^2_\ell, \varOmega ^2_\ell \rangle$ correspond to high strain rates and low rotation rates, and they decrease as the strain rate decreases. Interestingly, the inverse cascade appears explicitly in the lower-right corners of the plots, where the rotation rate is strong but the strain rate is weak. It is worth noting that the conditionally averaged values shown in figure 6 reflect the combined outcome of the forward and inverse cascades. Consequently, in specific regions characterized by distinct strain and rotation rates, events with forward and inverse cascades can cancel each other out. Only in the lower right corner is there an indication of net inverse cascade when the cascade rate is defined using the structure-function approach.

Figure 6.

Plots of (a–d) $\langle \varPhi _\ell | S^2_\ell, \varOmega ^2_\ell \rangle$ and (e–h) $\langle \varPi _\ell | S^2_\ell, \varOmega ^2_\ell \rangle$ at $\ell = \{30, 45, 60, 75\}\eta$. The black dashed lines in (d) represent isolines of $Q_\ell \equiv -\frac {1}{2}(S_\ell ^2 - \varOmega _\ell ^2)$; $\varPhi _\ell$ and $\varPi_\ell$ is normalized by $\langle \epsilon \rangle$; and $S^2_\ell$ and $\varOmega ^2_\ell$ is normalized by $\langle Q_w \rangle = \frac {1}{2}\langle \varOmega _\ell ^2 \rangle$.

In figure 6(d), we superimpose dashed lines representing the isolines of $Q_\ell$, with the $Q_\ell = 0$ line indicating the condition of equal strain and rotation rates. The parallel dashed lines correspond to $Q_\ell =-10$, $-5$, $0$, $5$, $10$ and $15$, respectively. The $Q_\ell =15$ line appears near the boundary separating the red and blue regions. However, the boundary of the blue region does not appear to align well with the $Q_\ell$ isoline. This observation suggests that $Q_\ell$ might be not enough to provide an adequate threshold for distinguishing the net forward and inverse cascade regions.

Figures 6(d–g) present results for the joint conditionally averaged $\varPi _\ell$ based on $S^2_\ell$ and $\varOmega ^2_\ell$, corresponding to the same filter scales as in figures 6(a–d). It is evident that trends for the positive cascade rate (red region) for $\varPi _\ell$ resemble closely those of $\varPhi _\ell$, with the peak of the forward cascade occurring at a high strain rate and low rotation rate. The magnitude of the maximum forward cascade rate for $\varPi _\ell$ is slightly weaker compared to that of $\varPhi _\ell$. The most significant difference is that only a few instances of blue squares are observed in regions characterized by strong rotation and weak strain, indicating that the overall predominance of the forward cascade persists regardless of the local values of $S^2_\ell$ and $\varOmega ^2_\ell$. These results highlight some important statistical differences between $\varPhi _\ell$ and $\varPi _\ell$. We further evaluate the conditional averages of the energy fluxes $\varPhi _\ell$ and $\varPi _\ell$, conditioned on either $S^2_\ell$ and $\varOmega ^2_\ell$ individually. This analysis is motivated by the work of Buaria & Pumir (Reference Buaria and Pumir2022), which highlighted different scalings of conditional averages with respect to strain rate compared to rotation rates even though one would expect similar results based on dimensional arguments. In figure 7(a), we demonstrate that both $\varPhi _\ell$ and $\varPi _\ell$ exhibit a power-law relationship when conditioned on $S^2_\ell$, with slope $3/2$. This scaling aligns with dimensional analysis and Kolmogorov scaling, $\langle \varPhi _\ell | S_\ell ^2 \rangle \sim [S_\ell ^2]^{3/2}$ and $\langle \varPi _\ell | S_\ell ^2 \rangle \sim [S_\ell ^2]^{3/2}$. In contrast, when conditioning $\varPhi _\ell$ and $\varPi _\ell$ on $\varOmega ^2_\ell$ (figure 7b), a different trend emerges, with much weaker dependence on rotation rate.

Figure 7.

(a) Plots of $\langle \varPhi _\ell | S^2_\ell \rangle$ (black triangles and lines) and $\langle \varPi _\ell | S^2_\ell \rangle$ (blue open triangles and lines). (b) Plots of $\langle \varPhi _\ell | \varOmega ^2_\ell \rangle$ (black triangles and lines) and $\langle \varPi _\ell | \varOmega ^2_\ell \rangle$ (blue open triangles and lines). The red dashed lines represent lines with slope $3/2$ in the log-log plot.

To develop a more detailed understanding of the inverse cascade region within $\langle \varPhi _\ell | S^2_\ell, \varOmega ^2_\ell \rangle$, we perform a further analysis by dividing the samples based on $\varPhi _\ell > 0$ and $\varPhi _\ell < 0$ for $\ell = 45 \eta$. We then calculate the conditional average of these separated samples considering $S^2_\ell$ and $\varOmega ^2_\ell$, denoted as $\langle \varPhi _\ell | S^2_\ell, \varOmega ^2_\ell, \varPhi _\ell > 0 \rangle$ and $\langle \varPhi _\ell | S^2_\ell, \varOmega ^2_\ell, \varPhi _\ell < 0 \rangle$. The results are presented in figure 8. From figure 8(a), it can be observed that the forward cascade clearly increases with $S^2_\ell$, with the highest values of $\varPhi _\ell$ concentrated in the upper left corner. It increases also with $\varOmega ^2_\ell$, but less rapidly. Combined, the trend seems to be an increase approximately proportional to ${\sim }S^2_\ell + 0.75 \, \varOmega ^2_\ell$. Differently, figure 8(b) illustrates that the inverse cascade is approximately proportional to ${\sim }S^2_\ell + 0.5 \, \varOmega ^2_\ell$, i.e. shallower isolines extending more in the horizontal direction than in the vertical compared to the forward cascade case shown in figure 8(a). This observation elucidates why the strongest red region in figure 6(b) emerges at the largest $S^2_\ell$, while below this threshold, the forward cascade events weaken progressively and are gradually cancelled out by the inverse cascade. Finally, in regions characterized by a weak strain rate and strong rotation rate, the inverse cascade becomes the dominant contribution.

Figure 8.

(a,b) The conditional averages of cascade rates obtained by sampling positive and negative signs, i.e. $\langle \varPhi _\ell | S^2_\ell, \varOmega ^2_\ell, \varPhi _\ell > 0 \rangle$ and $\langle \varPhi _\ell | S^2_\ell, \varOmega ^2_\ell, \varPhi _\ell < 0 \rangle$, respectively. (c,d) The logarithm base 10 of the number of samples on the $S^2_\ell, \varOmega ^2_\ell$ map (out of a total of $2\times 10^6$ samples). The isolines in (c,d) are the values corresponding to the contour.

Figures 8(c,d) display the distribution of the number of samples corresponding to positive and negative cascade rates in logarithmic scale (out of the 2 million samples (balls) considered). Our focus is directed specifically towards the bottom right corners of the plots, which correspond to the region where the inverse cascade is observed in figure 6(b). Interestingly, we observe that at $\varOmega ^2_\ell /\langle Q_w \rangle \approx 40$ and $S^2_\ell /\langle Q_w \rangle < 10$, the numbers of samples representing both inverse and forward cascade rates are approximately equivalent, falling within the range $10^{1}$–$10^{1.5}$. This implies that within this region, the magnitude of the inverse cascade must be significant to achieve net negative values for the conditional average. Still, for the conditions with net inverse cascade, the number of occurrences for both forward and inverse cascade rates is quite small, of the order of only $1/10^5$ of the total samples, indicating a very low frequency. We point out that in the extreme bins, the conditionally averaged fluxes are unlikely to be fully converged. Therefore, our observations are meant to be mostly qualitative in these regions. In particular, the inverse cascade regions depicted in figure 6 are attributed primarily to rare but intense events. In the next subsection, we will show that inverse cascade can be better characterized by conditioning on $Q_\ell$ and $R_\ell$ invariants.

5.2. Conditional statistics based on $Q_\ell$ and $R_\ell$ invariants

In the context of the $\varOmega ^2_\ell$ and $S^2_\ell$ map shown in figure 6, we observe the presence of a distinct inverse energy cascade in the region characterized by strong rotation but weak strain (corresponding to large $Q_{\ell}$ values) for $\varPhi _\ell$. However, such observations did not hold for $\varPi _\ell$. But these results do not preclude the possibility that net forward and inverse cascade may be associated with other invariants of the filtered velocity gradient tensor $\tilde {A}_{ij}$.

Figure 9 shows the joint conditionally averaged $\langle \varPhi _\ell | Q_\ell, R_\ell \rangle$ and $\langle \varPi _\ell | Q_\ell, R_\ell \rangle$ at four different scales, namely, $\ell = \{30, 45, 60, 75\} \eta$. Across all plots, we can observe the distinctive teardrop-shape pattern on the $Q$–$R$ map, as reported in previous studies (Chong et al. Reference Chong, Perry and Cantwell1990; Meneveau Reference Meneveau2011). The black solid lines are the boundaries, separating the four quadrants based on the signs of $Q_{\ell}$ and $R_{\ell}$. Notably, it becomes evident that both $\varPhi_{\ell}$ and $\varPi_{\ell}$ exhibit a strong and dominant inverse cascade in the quadrant characterized by $Q_{\ell} > 0$ and $R_{\ell} > 0$. It is useful to recall that the variable $R_\ell$ is associated with the rate of change of $Q_\ell$ (in fact, assuming restricted Euler dynamics, they are related by ${\rm d}Q_\ell /{\rm d}t = -3R_\ell$; Cantwell Reference Cantwell1992; Meneveau Reference Meneveau2011). Therefore $R_\ell >0$ corresponds to a decreasing trend of $Q_\ell$, i.e. vortex compression. We find that such events are associated with inverse energy cascade. This observation is particularly interesting considering the absence of an observable inverse cascade for $\varPi _\ell$ in the $S^2_\ell$ and $\varOmega ^2_\ell$ map. Hence these results indicate that the variables $Q_{\ell}$ and $R_{\ell}$ provide a more effective characterization to identify inverse cascade using conditional averaging.

Figure 9.

Plots of (a–d) $\langle \varPhi _\ell | Q_\ell, R_\ell \rangle$ and (e–h) $\langle \varPi _\ell | Q_\ell, R_\ell \rangle$ at $\ell = \{30, 45, 60, 75\}\eta$. The black solid lines separate the four quadrants, and the two lines $Q=-(\frac {27}{4} R^2)^{1/3}$ (Vieillefosse lines) are also shown. Here, $\varPhi _\ell$ and $\varPi _\ell$ are normalized by $\langle \epsilon \rangle$; $Q_\ell$ is normalized by $\langle Q_w \rangle$; and $R_\ell$ is normalized by $\langle Q_w \rangle ^{3/2}$. The data and editable analysis code that generated these joint p.d.f.s (for the case at $\ell =45\eta$) can be found at https://multimediafluids.uam.elogim.com/S0022112023010662/JFM-Notebooks/files/figure9.

The region characterized by $Q_{\ell}<0$ and $R_{\ell}>0$, which corresponds to the strain-dominated region, exhibits the most pronounced forward cascade. This observation aligns with the findings of many prior analyses in the literature (Borue & Orszag Reference Borue and Orszag1998; van der Bos et al. Reference van der Bos, Tao, Meneveau and Katz2002; Carbone & Bragg Reference Carbone and Bragg2020; Johnson Reference Johnson2021) as well as those depicted in figure 9, further emphasizing that the strong local strain rate plays a crucial role in driving the forward energy cascade. We note that Borue & Orszag (Reference Borue and Orszag1998) and van der Bos et al. (Reference van der Bos, Tao, Meneveau and Katz2002) display conditional averages weighted by the joint p.d.f. of $R_{\ell}$ and $Q_{\ell}$. In their results, there was hardly any indication of backscatter/inverse cascade in the $Q_{\ell}>0$ and $R_{\ell}>0$ ‘vortex compression’ region, because the overall probability density of that region is smaller than for the other regions. However, the unweighted conditional averaging represents the relevant values if the large-scale flow is in that particular state ($Q_{\ell}>0$ and $R_{\ell}>0$), and is therefore relevant to our analysis. We also notice a small blue region in the $Q_\ell < 0$, $R_\ell < 0$ quadrant of $\langle \varPhi _\ell | Q_\ell, R_\ell \rangle$ (but not seen for $\langle \varPi _\ell | Q_\ell, R_\ell \rangle$). The occurrence of inverse cascade in strain-dominated, strain self-stretching regions appears intriguing. However, the small number of samples (${\sim }O(10)$) in the bin showing inverse cascade precludes us from ascribing much significance to this observation for now.

In a manner similar to that in figure 8, we perform further conditional averaging, also distinguishing positive and negative cascade rates. Figures 10(a,b) present $\langle \varPhi _\ell | Q_\ell, R_\ell, \varPhi _\ell > 0 \rangle$ and $\langle \varPhi _\ell | Q_\ell, R_\ell, \varPhi _\ell < 0 \rangle$ at $\ell = 45 \eta$. In the case of the inverse cascade, it is observed to occur in all four quadrants (figure 10b), with a more evenly distributed and symmetric presence in the upper two quadrants associated with $Q_{\ell}>0$, i.e. the rotation-dominated regions. The characteristic teardrop shape is less prominent and exhibits a shorter tail compared to the forward cascade (figure 10a). Regarding the forward cascade, it is evident that it is most dominant in the $Q_{\ell}<0$, $R_{\ell}>0$ quadrant, consistent with figure 9. However, in the $Q_{\ell} > 0$, $R_{\ell} >0$ quadrant, the forward cascade is weaker and is overall cancelled out by the stronger inverse cascade in that particular region.

Figure 10.

(a,b) Plots of $\langle \varPhi _\ell | Q_\ell, R_\ell, \varPhi _\ell > 0 \rangle$ and $\langle \varPhi _\ell | Q_\ell, R_\ell, \varPhi _\ell < 0 \rangle$. (c,d) Plots of the logarithm base 10 of the number of samples on the $Q_{\ell}, R_{\ell}$ map. The isolines in (c,d) are the values corresponding to the contours.

Figures 10(c,d) display the distributions of number of samples of forward and inverse cascade rates, respectively. The shapes of the distributions align with figures 10(a,b), but a majority of the samples are concentrated at the centre, corresponding to small values of $Q_{\ell}$ and $R_{\ell}$. This observation confirms that the strong instances of inverse cascade and forward cascade observed in figure 9 are determined primarily by infrequent but extreme events (intermittency). Note that the joint p.d.f. in figure 10(d) is more left–right symmetric than that in figure 10(c), suggesting a less non-Gaussian behaviour of the flow in regions of inverse cascade than in the forward cascade regions.

Finally, to provide a visual impression of the spatial distribution of regions of negative $\varPhi _\ell$ in the flow, in figure 11 we provide a 3-D visualization of two instances in specific small subdomains of size $150^3$ grid points (i.e. $(225\eta )^3$ out of the overall $8192^3$ DNS domain. The selection of these subdomains was based on the condition that $Q_\ell /\langle Q_w \rangle > 15$ and $R_\ell /\langle Q_w \rangle ^{3/2} > 15$ at the centre of each subdomain, such that the centre is at a strong vortex compression region. We then calculate the values of $\varPhi _\ell$, $\varPi _\ell$, $Q_\ell$ and $R_\ell$ at every second grid point. In figures 11(a,b), the light blue regions correspond to the isosurface of a large negative value of $\varPhi _\ell /\langle \epsilon \rangle = - 60$, indicating the presence of an inverse cascade with significant magnitude. Clearly, we can see that the occurrence of a strong inverse cascade is associated closely with the presence of the vortices. Figure 11(a) depicts that large negative $\varPhi _\ell$ appears near the centre and not at the core of the yellow tube, although one should recall that $\varPhi _\ell ({\boldsymbol {x}},t)$ is defined locally as centred at ${\boldsymbol {x}}$ but represents the energy cascade into balls of diameter $45\eta$, i.e. comparable to the diameter of the vortex (yellow region) shown. The blue region is also largely connected with the black isosurface ($R_\ell =20 \langle Q_w\rangle ^{3/2}$), indicating a strong ‘vortex compression’ region within the yellow tube. Interactive 3-D versions of the figure that can be accessed following the links in the figure caption help to elucidate the spatial structure. Figure 11(b) is an entirely different instance of similar conditions, showing a more broken up vortex, and showing that $\varPhi _\ell$ can also peak near the sides, and appear more scattered within the vortex. Coupled with the results shown in figure 9, the visualizations suggest that the strong inverse cascade occurs along the large-scale vortices, in regions of these vortices in which $R_{\ell}>0$, i.e. the vortex compression regions. We can also observe some yellow tubes, within which inverse cascade and compression are both absent. This is consistent with the statistics such that when conditionally averaging in terms of $Q_\ell$ but irrespective of $R_{\ell}$, the inverse cascade becomes very weak and almost non-existent. However, once one considers only $R_{\ell}>0$ regions, inverse cascade can be observed clearly in high vortical regions.

Figure 11.

(a,b) Isosurfaces of local $\varPhi _\ell /\langle \epsilon \rangle = -60$ (light blue), $Q_\ell /\langle Q_w \rangle = 20$ (yellow) and $R_\ell /\langle Q_w \rangle ^{3/2} = 30$ (black) in two different 3-D subdomains. (c,d) Isosurfaces of local $\varPi _\ell /\langle \epsilon \rangle = -20$ (green-blue) in the same subdomains and isosurface of $Q_\ell /\langle Q_w \rangle$ and $R_\ell / \langle Q_w \rangle ^{3/2}$ as in (a,b). Interactive visualizations are available for each panel at Panel (a), Panel (b), Panel (c), Panel (d). The link to the directory containing the visualization code and the 3-D fields for these data can be found at https://multimediafluids.uam.elogim.com/S0022112023010662/JFM-Notebooks/files/figure11.

We also show $\varPi _\ell /\langle \epsilon \rangle = -20$ (green-blue isosurface) in the corresponding 3-D subdomains shown as figures 11(c,d). Clearly, we can see that the green-blue and black regions largely overlap within the yellow region in figure 11(c). In figure 11(d), the overlapping between green-blue, yellow and black regions occurs at the centre and top right region of the subdomain, indicating that strong negative $\varPi _\ell$ is also associated with strong vortex compression within a high vortical region, consistent with figure 9. However, the patterns of the green SGS flux regions are smoother, consistent with the two-dimensional visualisation in figure 5.

Caution must be expressed that visualizations provide only qualitative impressions, and more quantitative analysis requires structure-based conditional averaging, such as undertaken recently in Park & Lozano-Duran (Reference Park and Lozano-Duran2023). While such analysis is beyond the scope of the present paper, the conditional statistics presented in figure 9 already by themselves provide the strong statistically robust connection between cascade rate measures and features of the large-scale velocity gradient tensor.