3.1. Problem set-up and methodology

We revisit the DNS dataset of homogeneous, forced, stably stratified turbulence studied by Yi & Koseff (Reference Yi and Koseff2022a ), which is governed by the incompressible, Navier–Stokes equations under the Boussinesq approximation with linear velocity forcing (Lundgren, T. S. Reference Lundgren2003; Rosales & Meneveau Reference Rosales and Meneveau2005):

(3.1) \begin{align}\partial _{j} u_{j} &= 0,\end{align}

(3.2) \begin{align} \partial _{t} u_{j} + u_{m} \partial _{m} u_{j} &= -\frac {1}{\rho _{0}} \partial _{j} p - \frac {g}{\rho _{0}} \rho \delta _{j3} + \nu \partial _{m} \partial _{m} u_{j} + A u_{j} , \end{align}

(3.3) \begin{align} \partial _{t} \rho + u_{m} \partial _{m} \rho &= -w {\textrm d}_{z} \overline {\rho }_{{B}} + D \partial _{m} \partial _{m} \rho , \end{align}

where $u_j$ , $p$ , $\rho$ represent the velocity, pressure and density fields, respectively; $\rho$ is defined as the perturbation from $\rho _{0} + \overline {\rho }_{{B}}(z)$ , with $\rho _{0}$ being the reference density and $\overline {\rho }_{{B}}(z)$ being the prescribed, stable, linearly varying background density field ( ${d}_{z} \overline {\rho }_{{B}} \lt 0$ and constant); $g$ is the gravitational acceleration; $\nu$ is the kinematic viscosity of the fluid; $D$ is the molecular diffusivity of density; and $A$ is the momentum forcing rate. A point to highlight is that this momentum forcing is active in all three directions and is proportional to the velocity perturbations in that direction. Imposing a stability that damps $w$ (via the buoyancy force), hence, also reduces $P_w = A w^{2}$ relative to $P_{u} = Au^{2}$ and $P_{v} = Av^{2}$ , further damping the vertical component. Tensor indices (1, 2, 3) correspond to spatial directions ( $x$ , $y$ , $z$ ) and velocity fields ( $u$ , $v$ , $w$ ) with gravity acting along the $z$ -axis. Repeated indices imply summation. Equations (3.1)–(3.3) were solved in a triply periodic, cubic domain of length $L = 2\pi$ using a Fourier pseudospectral solver with a fourth-order Runge–Kutta time-stepping scheme. Nonlinear terms were dealiased exactly by zero-padding (i.e. they are computed using zero-padded arrays of size $3N/2$ and retaining the first $N$ Fourier modes in each dimension). The code has been verified and used in several studies (Yi & Koseff Reference Yi and Koseff2022a , Reference Yi and Koseffb , Reference Yi and Koseff2023).

A key strength of using this dataset rather than the three shear-forced datasets studied by Yi & Koseff (Reference Yi and Koseff2023) is that the linear forcing ( $A u_j$ ) is divergence free, which allows for the pressure fluctuations to be decomposed into just two components associated with the non-divergent parts of the nonlinear and buoyancy terms in (3.2), which will be discussed further in § 4.1.

3.2. Second-moment budgets

Here, we overview the volume-averaged (indicated by overbars) budgets of the turbulent kinetic energy (TKE, $k= ({1}/{2}) \overline {u_{j} u_{j}}$ ), turbulent potential energy (TPE, $k_p = ({1}/{2}) \alpha ^{2} \overline {\rho \rho }$ ), Reynolds stress ( $\overline {u_{i} u_{j}}$ ) and buoyancy flux ( $B_j = ({g}/{\rho _{0}}) \overline {u_{j} \rho }$ ) associated with (3.1)–(3.3). Here, $\alpha = g/(\rho _{0} N)$ is a dimensional coefficient used to convert the dimensions of $\rho$ to match those of velocity (and thus convert the density variance into potential energy units). Because the system is axisymmetric about the vertical axis, we simply consider the horizontal and vertical Reynolds stresses (diagonal entries of $\overline {u_{i} u_{j}}$ ) and the vertical buoyancy flux ( $B = ({g}/{\rho _{0}}) \overline {w \rho }$ ). This reduced set of budgets is given as

(3.4) \begin{align} d_{t} k &= \underbrace {2Ak}_{\text{TKE production}} - \underbrace {\frac {g}{\rho _{0}} \overline {w \rho }}_{\substack {\text{vertical} \\ \text{buoyancy flux}}} - \underbrace {\nu \overline {\partial _{m} u_{j} \partial _{m} u_{j}}}_{\text{TKE dissipation}} = P_{k} - B - \epsilon _{k} , \end{align}

(3.5) \begin{align} d_{t} k_{p} &= \underbrace {\frac {g}{\rho _{0}} \overline {w \rho }}_{\substack {\text{vertical} \\ \text{buoyancy flux}}} - \underbrace {D \alpha ^{2} \overline {\partial _{m} \rho \partial _{m} \rho }}_{\text{TPE dissipation}} = B - \epsilon _{p} , \end{align}

(3.6) \begin{align} d_{t} k_{H} &= \underbrace {2Ak_{H}}_{\substack {\text{horizontal TKE} \\ \text{production}}} + \underbrace {\frac {1}{\rho _{0}} \overline {p s_{HH}}}_{\substack {\text{horizontal} \\ \text{pressure-strain}}} - \underbrace {\nu \overline {\partial _{m} u_{H} \partial _{m} u_{H}}}_{\substack {\text{horizontal TKE} \\ \text{dissipation}}} = P_{H} + R_{H} - \epsilon _{H} , \end{align}

(3.7) \begin{align} d_{t} k_{w} &= \underbrace {2Ak_{w}}_{\substack {\text{vertical TKE} \\ \text{production}}} + \underbrace {\frac {1}{\rho _{0}} \overline {p s_{33}}}_{\substack {\text{vertical} \\ \text{pressure-strain}}} - \underbrace {\frac {g}{\rho _{0}} \overline {w \rho }}_{\substack {\text{vertical} \\ \text{buoyancy flux}}} - \underbrace {\nu \overline {\partial _{m} w \partial _{m} w}}_{\substack {\text{vertical TKE} \\ \text{dissipation}}} = P_{w} + R_{w} - B - \epsilon _{w} , \end{align}

(3.8) \begin{align} d_{t} B &= \underbrace {A B}_{\substack {\text{vertical buoyancy flux} \\ \text{production by forcing}}} + \underbrace {\frac {g}{\rho _{0}^{2}} \overline {p \partial _{z} \rho }}_{\substack {\text{vertical} \\ \text{pressure scrambling}}} + \underbrace {2N^{2} (k_{w} - k_{p})}_{\substack {\text{source/sink due to}\\ {k_{w}\;\text{and}\;k_{p}}}} - \underbrace {(\nu + D)\left (\frac {g}{\rho _{0}} \right ) \overline {\partial _{m} w \partial _{m} \rho }}_{\substack {\text{vertical buoyancy flux} \\ \text{dissipation}}} \nonumber \\ &= P_{B} + R_{B} + 2 N^{2} (k_{w} - k_{p}) - \epsilon _{B} . \end{align}

Equations (3.4) and (3.5) are the volume-averaged kinetic and potential energy budgets, where $P_k$ is the production rate of TKE by the linear forcing term; $B$ is the vertical buoyancy flux that here converts TKE to TPE; and $\epsilon _k$ and $\epsilon _p$ are the dissipation rates of TKE and TPE, respectively. Equations (3.6) and (3.7) are the volume-averaged budgets for the horizontal and vertical components of TKE ( $k_H = ({1}/{2}) \overline {u_{H} u_{H}} = ({1}/{2}) ({1}/{2} \overline {uu} + ({1}/{2}) \overline {vv} )$ ; $k_w = ({1}/{2}) \overline {ww}$ ), where adding two times (3.6) with (3.7) results in (3.4) ( $k = 2k_H + k_w$ ). The sources and sinks in these equations are as follows: $P_H = ({A}/{2})( \overline {uu} + \overline {vv}$ ) and $P_w$ are the production rates of the horizontal and vertical components of TKE; $R_H = ({1}/{2\rho _{0}}) (\overline {ps_{11}} + \overline {ps_{22}})$ and $R_w$ are the pressure–strain correlations, which sum to zero (i.e. $2R_H = -R_w$ ) due to the incompressibility condition (3.1); and $\epsilon _{H} = ({\nu }/{2}) (\overline {\partial _{m} u_{1} \partial _{m} u_{1}} + \overline {\partial _{m} u_{2} \partial _{m} u_{2}})$ and $\epsilon _{w}$ are the dissipation rates of the horizontal and vertical components of TKE. Equation (3.8) is the volume-averaged budget for the vertical buoyancy flux, where $P_B$ is a source term due to forcing (equivalent to production by the mean-velocity gradient in stratified wall-bounded flows); $R_B$ is the pressure scrambling term; $2 N^2 k_w$ is a source due to $k_w$ (equivalent to production by the mean-density gradient in stratified wall-bounded flows); $-2 N^2 k_p$ is a sink due to $k_p$ (equivalent to buoyant destruction in stratified wall-bounded flows); and $\epsilon _B$ is a ‘dissipation’ term for vertical buoyancy flux (this term, however, is not sign definite). The transport terms are exactly zero because our flow is statistically homogeneous in all three spatial directions, which is why only the volume-averaged sources and sinks remain in (3.4)–(3.8).

3.3. Reynolds stress and buoyancy flux budgets and their relationship with mixing coefficient

We can visualise the energy exchange in our system under statistically stationary conditions following (3.5)–(3.7). In figure 3, ingoing arrows indicate sources and outgoing arrows indicate sinks of each of the three energy buckets. Forcing directly produces $k_H$ and $k_w$ (black arrows). There is also intercomponent exchange by pressure–strain redistribution, which under stable conditions, and with the three-dimensional forcing in our DNS, will convert $k_H$ into $k_w$ to reduce the large-scale anisotropy due to buoyancy effects (blue arrow). (Note that the conversion would have been from $k_w$ to $k_H$ for cases with $P_w/P_k=1$ as discussed in the previous section.) The buoyancy flux converts $k_w$ into $k_p$ (orange arrow) and, finally, there are dissipative losses (red arrows).

Figure 3. Energy exchange diagrams under the statistically stationary conditions of our DNS set-up, following (3.5)–(3.7). Ingoing arrows indicate sources and outgoing arrows indicate sinks of the three energy buckets (horizontal and vertical components of TKE and TPE). Black arrows indicate direct production by the forcing term, the blue arrow indicates the pressure–strain redistribution term, the orange arrow indicates the vertical buoyancy flux and red arrows indicate the dissipation terms.

Figure 4. Volume- and time-averaged values of ${\textit{Re}}_{L}$ and $\textit{Fr}_{k}$ from the DNS dataset of Yi & Koseff (Reference Yi and Koseff2022a ) (purple circles) and an additional higher Reynolds number simulation at strong stability (purple star).

Before presenting the key results, we first illustrate the parameter space spanned by the DNS dataset of Yi & Koseff (Reference Yi and Koseff2022a ) by plotting the volume- and time-averaged values of ${\textit{Re}}_{L}$ and $\textit{Fr}_{k}$ from the 42 simulations in this dataset in figure 4 (purple circles). Lines of constant ${\textit{Re}}_{b} = {\textit{Re}}_{L} \textit{Fr}_{k}^{2}$ are marked by the diagonal dashed lines (increasing in order of magnitude from $1$ to $10^{4}$ ). The dataset spans roughly three orders of magnitude in $\textit{Fr}_{k}$ , covering near neutral to very strongly stable conditions. The simulations with $\textit{Fr}_{k} \gt 1$ have ${\textit{Re}}_{L} \approx 100$ , whereas the simulations with $\textit{Fr}_{k} \leqslant 1$ have ${\textit{Re}}_{b}$ between $10$ and $100$ , where ${\textit{Re}}_{b}$ decreases for simulations with smaller values of $\textit{Fr}_{k}$ as a wider range of anisotropic scales are simulated at the expense of isotropic scales. We choose to plot the results as a function of $\textit{Fr}_{k}$ for simplicity and continuity with how they are originally presented by Yi & Koseff (Reference Yi and Koseff2022a ), but we should note that Reynolds number variations do affect the quantities of interest as we have discussed in § 2. To illustrate this multiparameter dependence, we also present the results from a higher Reynolds number simulation with strong stability ( $\textit{Fr}_{k}\sim \mathcal{O}(10^{-2})$ ) (purple star). Key simulation input and output parameters are summarised in table 1 with C-series simulations involving a fixed value of the forcing coefficient $A$ in time, whereas the K-, D- and V-series simulations involve time-varying values of $A$ to maintain specific values of $k$ , $\epsilon _{k}$ and $k_{w}$ , which are reported in the table. The time-varying forcing approach builds on that of Bassenne et al. (Reference Bassenne, Urzay, Park and Moin2016), and the details of how $A(t)$ is chosen can be found in Appendix B of Yi & Koseff (Reference Yi and Koseff2022a ). A key benefit of this time-varying forcing scheme is that the simulations more rapidly reach statistically stationary conditions relative to a constant forcing scheme. We illustrate this benefit further in relation to the growth of energy in large horizontal-scale motions termed ‘shear modes’ in strongly stable, forced turbulence simulations in Appendix A.

Table 1. Volume- and time-averaged input and output parameters for the simulations. The C-series simulations involve a fixed value of $A$ in time, whereas the K-, D- and V-series simulations involve time-varying values of $A$ to maintain specific values of $k$ , $\epsilon _{k}$ and $k_{w}$ , respectively, which are provided in the table. All simulations used $N = 64$ Fourier modes in each spatial direction, except for simulation V2_128, which used $N = 128$ modes. The C-, K- and D-series simulations all had values of $\kappa _{\textit{max}} \eta \approx 2$ , whereas the lowest value for the V-series simulations was $\kappa _{\textit{max}} \eta \approx 1.25$ .

Figure 5. Steady-state, volume- and time-averaged budgets of (a) $k_H$ , (b) $k_w$ and (d) $B$ as a function of the turbulent Froude number ( $\textit{Fr}_k$ ). The three panels correspond to (3.6)–(3.8). The mixing coefficient $\varGamma$ is plotted in panel (c) as a function of $\textit{Fr}_k$ with the dissipation anisotropy factor $c_{3} = \epsilon _{w}/\epsilon _{k}$ shown in colour. The results from the higher Reynolds number simulation (V2_128) are denoted by filled stars.

Next, we consider the volume-averaged budgets of the horizontal and vertical components of the TKE ( $k_H$ and $k_w$ ), and the vertical buoyancy flux ( $B$ ) as a function of the turbulent Froude number ( $\textit{Fr}_k$ ) in figure 5(a,b,d). We also connect these budgets to the relationship between the mixing coefficient ( $\varGamma$ ) and $\textit{Fr}_k$ in figure 5(c). In the following paragraphs, we denote weak stratification where we observe monotonic increase of $R_{w}/P_{k}$ (relative importance of pressure–strain redistribution of $k_{w}$ to total mechanical production), moderate stratification where $R_{w}/P_{k}$ plateaus along with $\varGamma$ , and strong stratification where pressure scrambling (source/sink in the vertical buoyancy flux budget) becomes positive and where $\varGamma$ begins to decrease. These choices are made given that the specific $\textit{Fr}_{k}$ thresholds associated with these dynamical changes in the volume-averaged budgets will be sensitive to how the turbulence is forced and the Reynolds number.

In figures 5(a) and 5(b), we plot the steady-state, volume- and time-averaged budgets of $k_H$ and $k_w$ normalised by the total production of TKE ( $P_k$ ). First, we consider the budget for $k_H$ (figure 5 a). When stratification is weak (large Froude numbers), the pressure–strain correlation is weak (blue exes), and there is a balance between production and dissipation (black and red triangles). As stratification increases (moving from large to small Froude numbers), production of $k_H$ increases in magnitude relative to total production (black triangles) since $k_H$ becomes larger than $k_w$ due to damping by buoyancy and consequent reduced production of the latter. With this change, the pressure–strain correlation becomes a sink of $k_H$ , and it reaches a plateau around $\textit{Fr}_{k} \approx 0.7$ before starting to decrease in magnitude for $\textit{Fr}_{k} \lt 0.1$ . Dissipation of $k_H$ initially becomes smaller until $\textit{Fr}_{k} \approx 0.7$ before increasing in magnitude for $\textit{Fr}_{k} \lt 0.7$ .

Next, we consider the budget for $k_w$ (figure 5 b). Like with $k_H$ , when stratification is weak (large Froude numbers), the pressure–strain correlations (blue exes) and buoyancy flux (orange stars) are negligible, and there is a balance between production and dissipation. As stratification increases (moving from large to small Froude numbers), vertical motions are damped by buoyancy ( $k_w \lt k_H$ ) and direct production of $k_w$ decreases. As large-scale anisotropy becomes important in the system, the pressure–strain correlation becomes a source of $k_w$ , increasing in magnitude until reaching a plateau at $\textit{Fr}_{k} \approx 0.7$ before decreasing for $\textit{Fr}_{k} \lt 0.1$ . Relatedly, the buoyancy flux increases in magnitude until $\textit{Fr}_{k} \approx 0.7$ , overtaking dissipation as the main sink of $k_w$ , before sharply decreasing in magnitude for $\textit{Fr}_{k} \lt 0.1$ . Dissipation of $k_w$ decreases as stratification is increased.

In figure 5(d), we plot the steady-state, volume- and time-averaged budget of the vertical buoyancy flux ( $B = ({g}/{\rho _{0}}) \overline {w\rho }$ ) as a function of the turbulent Froude number ( $\textit{Fr}_k$ ). The terms have been normalised by the sum of all source terms to keep their values in the range of $\pm 1$ . When stratification is weak (large Froude numbers), the buoyancy flux is generated primarily by the source due to $k_w$ (black triangles) and secondarily by the source due to forcing (black circles), and it is destroyed primarily by pressure scrambling (blue exes) and secondarily by dissipation (red triangles). As stratification increases (moving from large to small Froude numbers), the source due to forcing and the dissipation weaken, and a dominant balance between three terms arises: generation by source due to $k_w$ , destruction by pressure scrambling, and sink due to $k_p$ (orange stars). When $\textit{Fr}_{k} \approx 0.1$ , the pressure scrambling term becomes negligible, and there is a balance between source due to $k_w$ and sink due to $k_p$ . However, as stratification increases further ( $\textit{Fr}_{k} \lt 0.1$ ), the source due to $k_w$ weakens, and pressure scrambling switches signs and begins to generate buoyancy flux to ensure that the buoyancy flux remains positive, which is expected for stably stratified systems. The emerging dominant balance as one goes to even lower $\textit{Fr}_k$ than we simulated seems to be between the pressure scrambling term correlating vertical velocity and density to produce buoyancy flux and buoyancy destruction working to reduce the correlation and flux. This balance also emerged in DNS of stably stratified wall-bounded flow away from the wall and at high stabilities reported by Shah & Bou-Zeid (Reference Shah and Bou-Zeid2014) (see their figure 20d for $\textit{Ri}_b=0.968$ ).

Finally, in figure 5(c), we plot the mixing coefficient ( $\varGamma$ ) as a function of $\textit{Fr}_k$ with the dissipation anisotropy factor $c_{3} = \epsilon _{w}/\epsilon _{k}$ from § 2 shown in colour. Doing so illustrates how the range of isotropic scales vary simultaneously with $\textit{Fr}_{k}$ , highlighting the multiparameter variations that occur across the set of simulations (see figure 4). Because we keep the simulation resolution fixed while varying the stability in the original set of simulations (circles), an increasing range of anisotropic scales are simulated at the cost of decreasing range of isotropic scales, which explains why simulations with weak stability have values of $c_{3}$ close to $1/3$ (indicating small-scale isotropy), while simulations with stronger stability have systematically smaller values of $c_{3}$ (indicating departures from small-scale isotropy). The three rightmost scaling regimes ( $\varGamma \sim \textit{Fr}_{k}^{-2}$ , $\varGamma \sim \textit{Fr}_{k}^{-1}$ and $\varGamma \approx \textrm{constant}$ ) have been theoretically discussed by Garanaik & Venayagamoorthy (Reference Garanaik and Venayagamoorthy2019), but the leftmost scaling of $\varGamma \sim \textit{Fr}_{k}^{0.4}$ is an empirical fit to the results of simulations V1, V2, V3 and V4 (four leftmost circles). The deviation of the higher Reynolds number simulation (filled star) is discussed in the subsequent paragraph. As stratification increases, $\varGamma$ increases until it plateaus at a value of $\varGamma \approx 0.2$ for $\textit{Fr}_{k} \approx 0.7$ . This coincides with when the pressure–strain redistribution becomes equally important in generating $k_w$ and reaches a plateau (figure 5 b). As stratification increases further, $\varGamma \approx 0.2$ is maintained until $\textit{Fr}_k \lt 0.1$ when it begins to decrease. This is connected with a decrease in the pressure–strain redistribution, which is the main mechanism for generating $k_w$ for $\textit{Fr}_k \lt 0.1$ , and there is an accompanied decrease in the vertical buoyancy flux (figure 5 b). As redistribution from $k_H$ to $k_w$ weakens, leading to a diminished vertical buoyancy flux, we observe in the vertical buoyancy flux budget that the source due to $k_w$ term diminishes in magnitude relative to the sink due to $k_p$ term, and the pressure scrambling term changes sign and becomes a source of vertical buoyancy flux (figure 5 c).

We note that increasing the Reynolds number while keeping stability relatively fixed (e.g. compare simulations V2 and V2_128, corresponding to $\textit{Fr}_{k} \approx 0.031$ and $\textit{Fr}_{k} \approx 0.025$ , respectively; results from simulation V2_128 are denoted by filled stars) results in a more effective pressure–strain redistribution from $k_{H}$ to $k_{w}$ and weakens the loss of $k_{H}$ by dissipation (figure 5 a). This more effective generation of $k_{w}$ is accompanied by greater losses via enhanced conversion of $k_{w}$ into $k_{p}$ by the buoyancy flux and viscous dissipation (figure 5 b). Regarding the larger buoyancy fluxes associated with increasing Reynolds number, we observe that this change is associated with a relative strengthening of the source due to $k_{w}$ compared with pressure scrambling, which is due to the greater conversion of $k_{H}$ into $k_{w}$ by more effective pressure–strain redistribution (figure 5 d). These changes result in a greater value of $\varGamma$ as the Reynolds number is increased (figure 5 c), which can be explained by considering (2.6). We find that increasing the Reynolds number while keeping stability relatively fixed results in (i) negligible changes to the fraction of the TKE production that directly generates $k_{w}$ ( $P_{w}/P_{k}$ remains largely unchanged), (ii) enhanced pressure–strain redistribution into $k_{w}$ relative to TKE production ( $R_{w}/P_{k}$ increases) and (iii) more of the TKE loss occurring via dissipative losses of $k_{w}$ ( $c_{3}$ increases). Visually, we can see this by considering figure 2(a) (note that $P_{w}/P_{k} \sim \mathcal{O}(10^{-2})$ for simulations V2 and V2_128), where simultaneous increases in $R_{w}/P_{k}$ and $c_{3}$ result in larger values of $\varGamma = \textit{Ri}_{f}/(1 - \textit{Ri}_{f})$ (i.e. movement towards the top right corner of the realisable region that is left of the orange vertical line). Regarding the transition between $\varGamma \approx \textrm{constant}$ and $\varGamma \sim Fr^{0.4}$ , we note that the latter relationship is an empirical fit that needs to be revisited based on simulations that vary $\textit{Fr}_{k}$ while keeping ${\textit{Re}}_{b}$ fixed, which would reveal the relationship between $\varGamma$ and $\textit{Fr}_{k}$ at strong stability for different values of $c_{3}$ .

Bou-Zeid et al. (Reference Bou-Zeid, Gao, Ansorge and Katul2018) in studying the relationship between $\textit{Ri}_f$ and the gradient Richardson number ( $\textit{Ri}_g=N^2/S^2$ ) observed a similar relationship where $\textit{Ri}_f$ first increases with $\textit{Ri}_g$ (increasing stability), but then reaches a plateau, where the former regime was referred to as the ‘ $B$ -limited regime’, indicating a limitation in the magnitude of the vertical buoyancy flux due to the background density gradient and its ability to convert TKE to TPE, and the latter regime was referred to as the ‘ $R_w$ -limited regime’, indicating a limitation in the magnitude of the vertical buoyancy flux due to the pressure-redistribution term. These ideas are connected with the left and right flanks of the relationship between $\varGamma$ and $\textit{Ri}_g$ as described by Caulfield (Reference Caulfield2021). The results from the higher Reynolds number simulation at strong stability highlights the need to explore the effects of Reynolds number in both the ‘ $B$ -limited regime’ and the ‘ $R_{w}$ -limited regime’ for a more complete physical understanding of mixing associated with stably stratified turbulence.

In this section, we have summarised the key findings of Yi & Koseff (Reference Yi and Koseff2022a ), connecting the shape of $\varGamma$ as a function of some measure of dynamic stability (here we have used $\textit{Fr}_k$ ) to changes in the generation of $k_w$ by pressure–strain redistribution and the associated change in the generation of $B$ by the pressure scrambling term for very strong stability. We also briefly explored the role of the Reynolds number through an additional simulation, which (when combined with our analysis in § 2) suggests that greater values of $\varGamma$ will result when Reynolds number is increased while keeping $\textit{Fr}_{k}$ fixed for strong stability. Given that Yi & Koseff (Reference Yi and Koseff2023) also reported sign changes of the pressure scrambling terms at strong stability for three different types of shear-forced, stably stratified turbulence simulations, we attempt to explain this seemingly more general feature of stably stratified turbulence in the next section by further studying these two pressure-related quantities.

4.1. Pressure decomposition

To diagnose these pressure effects further, we consider the Poisson equation governing the pressure fluctuations by taking the divergence of (3.2),

(4.1) \begin{equation} \partial _{j} \partial _{j} p = -\rho _{0}\partial _{j} u_{m} \partial _{m} u_{j} - g\partial _{z}{\rho }. \end{equation}

We then decompose the pressure fluctuations into nonlinear and buoyancy components ( $p_{\textit{NL}}$ and $p_{{B}}$ ) that separately satisfy the two right-hand side source terms:

(4.2) \begin{align} \partial _{j} \partial _{j} p_{\textit{NL}} &= -\rho _{0}\partial _{j} u_{m} \partial _{m} u_{j}, \end{align}

(4.3) \begin{align} \partial _{j} \partial _{j} p_{{B}} &= - g\partial _{z}{\rho }. \end{align}

As mentioned in § 3, the shear-forced datasets of Yi & Koseff (Reference Yi and Koseff2023) involve an additional right-hand side term in (4.1) due to the non-divergence of the shear forcing terms. This would require solving for a third set of pressure fluctuations associated with the terms involving the mean shear. The current dataset, however, allows us to study just the competition between the nonlinear and buoyancy components of the pressure fluctuations. When running the DNS, we solved the full Poisson equation for the pressure (4.1). Here, we separately compute $p_{\textit{NL}}$ and $p_{{B}}$ offline as a postprocessing step to probe their responses to increasing stability.

Figure 6. Normalised (a) pressure variance, (b) pressure–strain correlation and (c) pressure scrambling term as a function of $\textit{Fr}_{k}$ . The results from the higher Reynolds number simulation (V2_128) are denoted by filled stars.

In figure 6(a), we plot the normalised pressure variance as a function of the turbulent Froude number ( $\textit{Fr}_k$ ). With $p = p_{\textit{NL}} + p_{{B}}$ , the volume-averaged pressure variance is $\overline {p^{2}} = \overline {p_{\textit{NL}}^{2}} + \overline {p_{{B}}^{2}} + 2 \overline {p_{\textit{NL}} p_{{B}}}$ . The total pressure variance (blue exes) has been plotted using the sum of the three right-hand side terms to verify the offline postprocessing of the nonlinear and buoyancy components of the pressure fluctuations. The total pressure variance is dominated by the nonlinear component (purple squares), but the buoyancy component (blue circles) increases in magnitude with stability. The nonlinear and buoyancy components are anticorrelated, and this anticorrelation strengthens with stability (orange triangles).

Next, we consider the normalised pressure–strain correlations ( $R_{w}/P_{k}$ ) in figure 6(b), which at its peak amounts to ${\sim} 25 \,\%$ of the TKE production. Recall that positive values indicate transfer of $k_{H}$ into $k_{w}$ and negative values indicate transfer of $k_{w}$ into $k_{H}$ . For weak stratification, the nonlinear component makes up most of the pressure–strain correlations, but the total contribution of the pressure–strain term to the dynamics is small since the turbulence is nearly isotropic in these simulations. The buoyancy component increases with stratification until $\textit{Fr}_{k} \approx 0.7$ , contributing to the redistribution of energy from $k_{H}$ into $k_w$ , which promotes large-scale return-to-isotropy. For $\textit{Fr}_{k} \lt 0.7$ , however, the buoyancy component begins to decrease and eventually changes sign for $\textit{Fr}_{k} \approx 0.3$ , becoming a sink of $k_w$ . This conversion from vertical to horizontal TKE goes counter to the return-to-isotropy expectation as it redistributes energy from the less energetic component $k_{w}$ (that in our simulations still has a source associated with the forcing) to the more energetic ones $k_{H}$ . We find that the decrease in the total pressure–strain correlation for $\textit{Fr}_{k} \lt 0.1$ is driven by both a weaker exchange of $k_H$ into $k_w$ by the nonlinear component and a stronger exchange of $k_w$ into $k_H$ by the buoyancy component. Interestingly, when considering the higher Reynolds number simulation with $\textit{Fr}_{k}\approx 0.025$ , we note that while the buoyancy pressure–strain redistribution maintains a similar magnitude as the lower Reynolds simulation with $\textit{Fr}_{k}\approx 0.031$ , we observe that the net pressure–strain redistribution increases by approximately $50\,\%$ due to a significant increase of the nonlinear component. We additionally note that this significant change in the correlation is associated with relatively negligible differences of the pressure variance components between the lower and higher Reynolds number simulations (figure 6 a).

Finally, we consider the normalised pressure scrambling term in figure 6(c). For weak stratification, the nonlinear component makes up most of the pressure scrambling term and represents a large negative fraction (approximately $-70\,\%)$ of the total production of covariance. Therefore, the total pressure scrambling term acts as a sink of vertical buoyancy flux. However, as stratification increases, the nonlinear component begins to weaken, and this is accompanied by a largely monotonic increase of the buoyancy component. Focusing on very strong stratification ( $\textit{Fr}_k\lt 0.1$ ), we note that the sign change of the total pressure scrambling term is due to the compensation between the nonlinear and buoyancy components, and as the buoyancy component of the pressure scrambling term starts to dominate, the total pressure scrambling term becomes an increasingly important means of generating vertical buoyancy flux (see also figure 5 d). Comparing the lower and higher Reynolds number simulations ( $\textit{Fr}_{k} \approx 0.031$ and $\textit{Fr}_{k} \approx 0.025$ ), we observe that the nonlinear component of the pressure scrambling increases slightly in magnitude while the buoyancy component of the pressure scrambling decreases in magnitude, leading to a 24 % reduction in the net pressure scrambling term.

4.2. Interpreting the nonlinear and buoyancy pressure fluctuations and their correlations

To develop an understanding of the nonlinear and buoyancy pressure fluctuations, we return to (4.2) and (4.3), which relate the Laplacian of the nonlinear and buoyancy pressure fluctuations to their respective right-hand side source terms ( $f_{\textit{NL}}(\boldsymbol{x},t) = -\rho _{0} \partial _{j} u_{m} \partial _{m} u_{j}$ and $f_{{B}}(\boldsymbol{x},t) = -g\partial _{z} \rho$ ). Because our problem set-up is triply homogeneous, the flow variables can be represented as a sum of Fourier modes (denoted by the circumflex symbol), the Fourier coefficients of $\partial _{j} \partial _{j} p$ are given by $-\kappa ^2 \hat {p}(\boldsymbol{\kappa })$ with $\boldsymbol{\kappa }$ being the wavenumber vector and $\kappa ^{2}$ being the squared wavenumber magnitude. Therefore, we can express the Fourier coefficients of the nonlinear and buoyancy components of the pressure fluctuations as $\hat {p}_{\textit{NL}}(\boldsymbol{\kappa },t) = -\hat {f}_{\!\textit{NL}}/\kappa ^2$ and $\hat {p}_{{B}}(\boldsymbol{\kappa },t) = -\hat {f}_{\!B}/\kappa ^2$ , where $\hat {f}_{\!\textit{NL}}(\boldsymbol{\kappa },t)$ and $\hat {f}_{\!B}(\boldsymbol{\kappa },t)$ are the Fourier coefficients of $f_{\textit{NL}}$ and $f_{{B}}$ . One additional constraint on the pressure fields is that their volume-average values are zero, which requires $\hat {p}_{\textit{NL}}(\boldsymbol{\kappa } = 0,t) = 0$ and $\hat {p}_{{B}}(\boldsymbol{\kappa } = 0,t) = 0$ .

Relatedly, we can also consider this relationship between the right-hand side terms of (4.2) and (4.3) and the pressure fluctuations in physical space. By using Green’s function for the 3-D Poisson equation, we can write the pressure fluctuations as

(4.4) \begin{align} p_{\textit{NL}}(\boldsymbol{x},t) &= \iiint - \frac {f_{\textit{NL}}(\boldsymbol{x}',t) - \overline {f}_{\!\textit{NL}}(t)}{4\pi \lvert \boldsymbol{x}' - \boldsymbol{x} \rvert }\, {\textrm{d}} \boldsymbol{x}', \end{align}

(4.5) \begin{align} p_{{B}}(\boldsymbol{x},t) &= \iiint - \frac {f_{{B}}(\boldsymbol{x}',t) - \overline {f}_{\!B}(t)}{4\pi \lvert \boldsymbol{x}' - \boldsymbol{x} \rvert }\, {\textrm{d}} \boldsymbol{x}', \end{align}

which involve convolution integrals of Green’s function $-1/(4 \pi r)$ ( $r$ is the distance between two points marked by $\boldsymbol{x}$ and $\boldsymbol{x}'$ ) and the two right-hand side source terms with their volume-averaged values subtracted, $f_{\textit{NL}} - \overline {f}_{\!\textit{NL}}$ and $f_{{B}} - \overline {f}_{\!B}$ , which is necessary to ensure a zero volume-average for the pressure fields. The pressure fields are defined using convolution integrals of the right-hand side terms with the Green’s function because they are governed by elliptic partial differential equations (PDEs), where the pressure at some point $\boldsymbol{x}$ depends on the values of the right-hand side terms over the whole domain. We note that by the convolution theorem, we can also arrive at $\hat {p}_{\textit{NL}}(\boldsymbol{\kappa },t)=-\hat {f}_{\!\textit{NL}}/\kappa ^{2}$ and $\hat {p}_{{B}}(\boldsymbol{\kappa },t)=-\hat {f}_{\!B}/\kappa ^{2}$ from (4.4) and (4.5), where the Fourier coefficients of Green’s function are $-1/\kappa ^{2}$ .

Next, we make one final manipulation to $f_{\textit{NL}}$ by re-expressing it using the rate-of-strain and rate-of-rotation tensors so that $\partial _{j} u_{m} = S_{mj} + R_{mj}$ and $\partial _{m} u_{j} = S_{jm} + R_{jm}$ . With these substitutions, we find

(4.6) \begin{equation} f_{\textit{NL}} = -\rho _{0}\partial _{j} u_{m} \partial _{m} u_{j} = -\rho _{0} (S_{jm} S_{jm} - R_{jm} R_{jm}) = -\rho _{0} \left ( S_{jm} S_{jm} - \frac {1}{2} \omega _{q} \omega _{q}\right )\!, \end{equation}

where we have simplified two terms by noting that a double contraction of symmetric and anti-symmetric tensor sums to zero, and rewriting the rate-of-rotation tensor in terms of the vorticity vector (see a discussion of this decomposition and its relation to the nonlinear pressure fluctuations by Vlaykov & Wilczek Reference Vlaykov and Wilczek2019). We also note that

(4.7) \begin{equation} \overline {f}_{\!\textit{NL}} = \rho _{0} \left ( \frac {1}{2} \overline {\omega _{q} \omega _{q}} - \overline {S_{jm} S_{jm}} \right ) = 0, \end{equation}

since $({1}/{2}) \overline {\omega _{q} \omega _{q}} = ({1}/{2}) (\overline {\partial _{j} u_{i} \partial _{j} u_{i}}) = \overline {S_{jm} S_{jm}}$ for homogeneous turbulence (see the discussion of pseudo and true dissipation in § 5.3 of Pope Reference Pope2000).

Here, we briefly summarise the motivation for taking the above-stated steps. We are seeking to understand the pressure–strain and pressure scrambling terms, which requires us to first understand what causes the pressure fluctuations in our problem. To do this, we considered the pressure Poisson equation given in (4.1). If the right-hand side source terms only involve a single Fourier mode, it is possible to state that the pressure fluctuations scale proportionally to the right-hand side terms following $-\partial _{j} \partial _{j} p \propto p \propto -f$ , but because the right-hand side source terms consist of contributions from a range of scales for turbulent flows, this statement cannot be applied generally in physical space. Rather, we have to resort to the exact relationship between the pressure fields and the right-hand side terms of (4.1) by considering the Green’s function solution for the pressure fluctuations ((4.2) and (4.3)), which led to (4.4) and (4.5) that involve convolution integrals of the right-hand side terms of (4.1) with a weighting factor that drops off as $1/r$ . Through the convolution theorem, these non-local products in physical space can be stated as local products in Fourier space ( $\hat {p}_{\textit{NL}}(\boldsymbol{\kappa },t) = -\hat {f}_{\!\textit{NL}}/\kappa ^2$ and $\hat {p}_{{B}}(\boldsymbol{\kappa },t) = -\hat {f}_{\!B}/\kappa ^2$ ), where the factor of $1/\kappa ^2$ indicates that larger scales of the right-hand side terms contribute more significantly to the pressure fluctuations than smaller scales. Next, we noted that the volume average of $f_{\textit{NL}}$ is zero (and so is the volume average of $\partial _{z} \rho$ due to statistical homogeneity of the density fluctuations in the constant mean density gradient simulations analysed here). Finally, given that the convolution integrals ((4.4) and (4.5)) filter the integrands such that the large scales are weighted more favourably (i.e. contribute more to the integral), we can interpret the positive and negative pressure fluctuations as follows:

1. (i) nonlinear pressure fluctuations, (a) $p_{NL} \gt 0$ , regions where large-scale weighted strain dominates, (b) $p_{NL} \lt 0$ , regions where large-scale weighted rotation dominates;

2. (ii) buoyancy pressure fluctuations, (a) $p_B \gt 0$ , regions where large-scale weighted unstable density gradients dominate, (b) $p_B \lt 0$ , regions where large-scale weighted stable density gradients dominate.

In the next subsection, these interpretations will prove useful in elucidating the physics of the quadrant analyses we perform on the vertical pressure–strain term in (3.7) and the vertical pressure scrambling term in (3.8).

4.3. Quadrant analyses of the buoyancy flux, pressure–strain redistribution and pressure scrambling dynamics

We now proceed to study the vertical density flux ( $w\rho$ ) and the decomposed pressure correlations that can help clarify the flow physics resulting in TKE redistribution in or out of the vertical TKE component ( $p_{\textit{NL}} \partial _{z} w$ , $p_{{B}} \partial _{z} w$ ) and the generation or destruction of buoyancy flux ( $p_{\textit{NL}} \partial _{z} \rho$ , $p_{{B}} \partial _{z} \rho$ ) by applying quadrant analyses (e.g. Wallace Reference Wallace2016). For the vertical density flux, quadrants 1 and 3 are associated with down-gradient fluxes, which convert $k_w$ to $k_p$ , and quadrants 2 and 4 are associated with counter-gradient fluxes, which convert $k_p$ to $k_w$ . For the two pressure–strain correlations (with $k_{H} \gt k_{w}$ for our simulations), quadrant 1 and 3 events convert $k_H$ to $k_w$ , promoting large-scale isotropy, whereas quadrant 2 and 4 events convert $k_w$ to $k_H$ , promoting large-scale anisotropy. For the two pressure scrambling terms, quadrant 1 and 3 events generate vertical buoyancy flux, strengthening the conversion of $k_w$ to $k_p$ , whereas quadrant 2 and 4 events destroy vertical buoyancy flux, weakening the conversion of $k_w$ to $k_p$ . We summarise these ideas and provide brief descriptions of these five sets of correlations in table 2.

In figure 7(a,b,c), we consider the joint probability distribution functions (p.d.f.s) of the vertical density flux ( $w\rho$ ) for three different values of $\textit{Fr}_k$ : weakest stability $\textit{Fr}_k \approx 5.13$ , where $B$ is negligible and pressure scrambling is negative; strong stability $\textit{Fr}_k \approx 0.25$ , where $B$ is still near its maximum but pressure scrambling is zero; and strongest stability $\textit{Fr}_k \approx 0.016$ , where $B$ has started to diminish and pressure scrambling is positive. The dark contours correspond to the five p.d.f. values ranging between $10^{-1}$ and $10^{-5}$ , with one decade spacing in between. The $x$ - and $y$ -axis variables are normalised by the root-mean-squared (r.m.s.) values from each simulation.

For stably stratified turbulence, the averaged vertical buoyancy flux $B = ({g}/{\rho _{0}}) \overline {w\rho }$ is expected to be positive, indicating a net conversion of $k_w$ to $k_p$ (the positive sign can be rationalised through an analogous argument typically used to explain why the Reynolds shear stress $\overline {uw} \lt 0$ for sheared flows; see also the orange stars in figure 5 b). In line with this, we observe that the two variables are positively correlated with down-gradient fluxes (quadrants 1 and 3) occurring more frequently than counter-gradient fluxes (quadrants 2 and 4), as indicated by the larger probability mass contained in the down-gradient flux quadrants than in the counter-gradient flux quadrants (see the text near the four corners). As stability increases (left to right), we observe that counter-gradient fluxes occur more often, and the overall correlation weakens (inner p.d.f. contours become more circular). Since counter-gradient buoyancy fluxes represent parcels returning towards their equilibrium elevation (where the mean density is equal to that parcel density), when there is weak mixing, enhanced counter-gradient fluxes result in weaker correlation between $w$ and $\rho$ and implies weaker vertical exchange of buoyancy.

Table 2. Descriptions of the contributions from the four quadrants to the five correlations of interest: (i) vertical density flux, (ii) nonlinear pressure–strain term, (iii) buoyancy pressure–strain term, (iv) nonlinear pressure scrambling term and (v) buoyancy pressure scrambling term.

Figure 7. (a–c) Joint probability distribution functions (p.d.f.s) and (d–f) covariance integrands of the normalised vertical density flux ( $w\rho /(w_{\textit{rms}} \rho _{\textit{rms}})$ ) for (a,d) very weak, (b,e) strong and (c, f) very strong stability. For panels (a–c), the text in each quadrant indicates the probability within each quadrant, whereas for panels (d–f), the text in each quadrant indicates that quadrant’s contribution towards the overall correlation.

Another way to explain this phenomenon is to consider the turbulent Froude number as a ratio of two time scales $\textit{Fr}_{k} = \epsilon /(Nk) \sim \tau _{B}/\tau _{L}$ , where $\tau _{B} \sim 1/N$ is the buoyancy time scale characterising the reversible exchange between $k_{w}$ and $k_{p}$ associated with density perturbations, and $\tau _{L} \sim k/\epsilon _{k}$ is the large-eddy time scale associated with the overturning of large eddies in the system. For $\textit{Fr}_{k} \gt 1$ , the large eddies overturn faster than the time needed for reversible exchange between $k_{w}$ and $k_{p}$ , allowing irreversible mixing associated with eddy turnover. For $\textit{Fr}_{k} \sim \mathcal{O}(1)$ , the large eddies overturn at a comparable rate needed for reversible exchange between $k_{w}$ and $k_{p}$ . Finally, for $\textit{Fr}_{k} \lt 1$ , the large eddies overturn more slowly than the reversible exchange between $k_{w}$ and $k_{p}$ , allowing parcels to return to their equilibrium position before they can fully mix.

If we recall that Rotta-style return-to-isotropy models take $\tau _{L}$ as the time scale associated with the pressure–strain redistribution, we can also interpret the turbulent Froude number as a measure of the relative importance of two reversible exchange processes: (i) reversible exchange between $k_{H}$ and $k_{w}$ by pressure–strain redistribution occurring over $\tau _{L}$ ; (ii) reversible exchange between $k_{w}$ and $k_{p}$ by vertical buoyancy flux occurring over $\tau _{B}$ . For $\textit{Fr}_{k} \gt 1$ , reversible exchange between $k_{H}$ and $k_{w}$ occurs more rapidly and thus more effectively than the reversible exchange between $k_{w}$ and $k_{p}$ . However, there is negligible overall mixing of density because the background density gradient is negligible, and the turbulence is weakly anisotropic, resulting in negligible pressure–strain redistribution. For $\textit{Fr}_{k} \sim \mathcal{O}(1)$ , reversible exchange between $k_{H}$ and $k_{w}$ by pressure–strain redistribution, and reversible exchange between $k_{w}$ and $k_{p}$ by vertical buoyancy flux occur at comparable rates. The background density gradient is significant and damps $k_{w}$ relative to $k_{H}$ , which results in large-scale anisotropy, but pressure–strain redistribution occurs rapidly enough to sustain $k_{w}$ , which in turn generates vertical buoyancy flux and sustains the reversible exchange between $k_{w}$ and $k_{p}$ (see figures 5 b and 5 d). For $\textit{Fr}_{k} \lt 1$ , the reversible exchange of $k_{H}$ and $k_{w}$ occurs more slowly than the reversible exchange between $k_{w}$ and $k_{p}$ by vertical buoyancy flux. This results in a reduced generation of $k_{w}$ by pressure–strain redistribution, which decreases $B$ because of a reduction in its primary generation mechanism (i.e. source due to $k_{w}$ , $2N^{2} k_{w}$ ; see figure 5 d). This then limits the conversion of $k_{w}$ to $k_{p}$ , preventing the complete decay of vertical perturbations.

Next, we quantify the contributions of the down-gradient and counter-gradient flux events to the overall flux by plotting the covariance integrand associated with these two joint p.d.f.s in figure 7(d,e, f). For each bin, the covariance integrand is computed by taking the values of the two variables and the associated value of the joint p.d.f., and multiplying them together. It thus represents the contribution of each bin to the overall flux. Down-gradient fluxes are shown in red and counter-gradient fluxes are shown in blue. As stability increases (left to right), we note that the down-gradient fluxes become less positively correlated, whereas the counter-gradient fluxes become more negatively correlated, resulting in a weaker positive correlation of the overall flux with increasing stability. Since both variables ( $w$ and $\rho$ ) are normalised by their root-mean-squared values, summing across the four quadrants results in the correlation coefficient rather than the overall flux, which is important to keep in mind to avoid incorrectly concluding that the buoyancy flux is stronger for the weakest stability simulation than for the strongest stability simulation. The correlation coefficient can be interpreted as a transport efficiency (Li & Bou-Zeid Reference Li and Bou-Zeid2011). Therefore, this result indicates that the vertical density flux is more effective (positively correlated) for weaker stability than for stronger stability, but we note from figure 5(b) that the buoyancy flux is negligible for the weakest stability case ( $\textit{Fr}_k \approx 5.13$ ), exhibits its peak value for the strong stability case ( $\textit{Fr}_k \approx 0.25$ ) and weakens but remains finite for the strongest stability case ( $\textit{Fr}_k \approx 0.016$ ).

Figure 8. (a) Probability mass of down-gradient and counter-gradient events (red and blue triangles, respectively), (b) correlation of down-gradient and counter-gradient events and the overall flux (orange triangles), (c) efficiency of the vertical density flux, and (d) normalised vertical buoyancy flux associated with the down-gradient, counter-gradient and overall flux as a function of $\textit{Fr}_k$ . The efficiency is defined as the net flux (sum over all quadrants) divided by the sum over just the down-gradient fluxes (sum over only quadrants 1 and 3). The orange stars are the net vertical buoyancy flux, and $P_k$ is the rate of production of TKE. The results from the higher Reynolds number simulation (V2_128) are denoted by filled stars.

To further understand how the contributions from the down-gradient and counter-gradient fluxes vary with increasing stability, we plot the following quantities as a function of $\textit{Fr}_k$ in figure 8: (a) the probability mass of down-gradient (red triangles) and counter-gradient events (blue triangles); (b) the correlation of down-gradient and counter-gradient events (blue triangles) and the correlation of the overall flux (orange stars); (c) the efficiency of the vertical buoyancy flux defined as the net vertical density flux normalised by the sum of only the down-gradient contributions; (d) the normalised buoyancy flux magnitudes associated with the down-gradient and counter-gradient events and the overall flux. We have chosen to combine the down-gradient and counter-gradient quadrants together to simplify our plots since the probability mass and correlations for the individual down-gradient and counter-gradient quadrants were quite similar to one another for all stability values, as expected given that our simulation set-up is homogeneous in the vertical direction, and invariant to vertical reflection of the density gradient and gravity force.

In figure 8(a), we see that for weak stratification, down-gradient fluxes occur approximately $75\,\%$ of the time and counter-gradient fluxes occur approximately $25\,\%$ of the time, and as stability is increased, we see that counter-gradient fluxes occur more frequently, such that the down-gradient and counter-gradient fluxes approach equal likelihoods. In figure 8(b), we plot the correlation coefficient associated with the down-gradient, counter-gradient and overall flux (red and blue triangles and orange stars) as a function of $\textit{Fr}_k$ . For weak stability, the down-gradient fluxes are strongly correlated with a value near $0.75$ , and the overall flux has a slightly weaker correlation due to the weak but finite negative correlations associated with the counter-gradient fluxes. With increasing stability, the down-gradient fluxes become more weakly correlated, whereas the counter-gradient fluxes become more strongly correlated, which drives the correlation of the overall flux towards zero and indicates less effective vertical density exchange. In figure 8(c), we plot another measure of the efficiency of the vertical density flux defined as the net flux normalised by the contributions from just the down-gradient fluxes. We note that for weak stability, the efficiency of the vertical density flux is approximately $0.9$ , but it monotonically weakens with increasing stability due to increasing contributions from counter-gradient fluxes. Lastly, in figure 8(d), we plot the normalised vertical buoyancy flux magnitudes associated with the down-gradient, counter-gradient and overall fluxes (red and blue triangles and orange stars). Because the turbulence is statistically stationary and triply homogeneous, we note that the orange symbols are equivalent to the mixing coefficient $\varGamma$ . When the stratification is weak relative to turbulence ( $\textit{Fr}_k\gt 1$ ), we note that increasing the background stratification results in sharp increases in the down-gradient fluxes, which are not accompanied by similar increases in the counter-gradient fluxes. As stratification and turbulence both become important ( $\textit{Fr}_k \sim \mathcal{O}(1)$ ), increasing the background stratification results in further enhancement of the down-gradient fluxes, but now this is accompanied by increases in the counter-gradient fluxes as well. This indicates simultaneous exchanges between $k_w$ and $k_p$ in both directions (i.e. conversion of $k_w$ into $k_p$ and conversion of $k_p$ into $k_w$ ) that are comparable in magnitude but with a net sum that still results in a net conversion of $k_{w}$ into $k_{p}$ . Lastly, when stratification becomes very strong ( $\textit{Fr}_k\lt 0.1$ ), increasing the background stratification results in a weakening of the down-gradient fluxes (red triangles), but the contributions from the counter-gradient fluxes remain largely unchanged (blue triangles), which indicates that the transfer from $k_w$ to $k_p$ becomes weaker, but the transfer from $k_p$ to $k_w$ remains largely unchanged.

Regarding the greater buoyancy flux magnitudes observed for the higher Reynolds number simulation (V2_128) ( $\textit{Fr}_{k}\approx 0.025$ ), we note that the correlations and the efficiency values reported in figures 8(b) and 8(c) are similar for this case to those for the lower Reynolds simulation ( $\textit{Fr}_{k}\approx 0.31$ ), which indicates that this change must be driven by the enhanced pressure–strain redistribution of $k_{H}$ to $k_{w}$ in the higher Reynolds number simulation (see figure 5 b). We also note that while approaching the limit of $\textit{Fr}_k \ll 1$ might lead to down-gradient and counter-gradient contributions becoming increasingly similar, resulting in a net zero flux, the turbulent flow would still involve exchanges among the horizontal, vertical and potential energy components, a net downscale cascade of energy, and dissipation at the viscous and diffusive scales. This is physically very different from the canonical set-up of adiabatically displacing fluid parcels in a stably stratified fluid that is initially at rest, which also involves a net zero buoyancy flux (no net exchange between $k_w$ and $k_p$ ) due to exact compensation between down-gradient and counter-gradient fluxes, but which lacks other physical aspects that are present for three-dimensional, stably stratified, turbulent flows (e.g. net downscale energy cascade and dissipative losses).

Figure 9. (a–c) Probability mass of the four quadrants of the pressure–strain correlations (Q1, red triangles; Q2, orange triangles; Q3, blue squares; Q4, purple circles); (d–f) correlation of the events that redistribute $k_{H}$ into $k_{w}$ (red triangles), events that redistribute $k_{w}$ into $k_{H}$ (blue triangles) and overall pressure-strain redistribution (black exes); (g–i) efficiency of the pressure-strain correlation; (j–l) normalised pressure–strain correlations associated with conversion of $k_{H}$ into $k_{w}$ , conversion of $k_{w}$ into $k_{H}$ and the overall redistribution (same symbols/colors convention as panels d–f). All variables are plotted as a function of $\textit{Fr}_{k}$ . The results from the higher Reynolds number simulation (V2_128) are denoted by filled stars.

Next, we seek to further understand how the contributions from the four quadrants of the pressure–strain correlations contribute to the overall redistribution of $k_{H}$ into $k_{w}$ by plotting the following quantities as a function of $\textit{Fr}_{k}$ in figure 9: (panels a,b,c) probability mass of the four quadrants of the pressure-strain correlations (Q1, red triangles; Q2, orange triangles; Q3, blue squares; Q4, purple circles); (panels d,e, f) correlations of the events that redistribute $k_{H}$ into $k_{w}$ (red triangles), events that redistribute $k_{w}$ into $k_{H}$ (blue triangles) and overall pressure–strain redistribution (black exes); (panels g,h,i) efficiency of the pressure–strain correlation; (panels j,k,l) normalised pressure–strain correlations associated with conversion of $k_{H}$ into $k_{w}$ , conversion of $k_{w}$ into $k_{H}$ , and the overall redistribution. The associated contour plots of the joint p.d.f.s and covariance integrands are shown in figures 13 and 14 in Appendix B.

We note that the peak pressure–strain correlation value is reached at $\textit{Fr}_{k} \approx 0.3$ (panel d), which coincides with the maximum probability mass values associated with Q3 redistribution events ( $k_{H}$ into $k_{w}$ ; blue squares in figure 9 a) and minimum values associated with Q2 and Q4 redistribution events ( $k_{w}$ into $k_{H}$ ; orange triangles and purple circles in figure 9 a). This also coincides with the maximum efficiency of the overall pressure–strain correlation, shown by the peak in $\eta _{p,\partial _{z} w}$ defined as the net pressure-strain correlation normalised by the sum over just the Q1 and Q3 redistribution events, which convert $k_{H}$ into $k_{w}$ (panel g). Focusing just on the nonlinear component of the pressure–strain correlations, we note that the peak correlation persists until stronger stability ( $\textit{Fr}_{k} \approx 0.1$ ) (panel e), and associatedly, the efficiency of the nonlinear pressure–strain correlation and normalised net redistribution ( $R_{w}/P_{k}$ ) both exhibit weaker declines for $\textit{Fr}_{k} \lt 0.1$ compared with the total pressure–strain correlation (compare panels h and k with panels g and j). These differences can be explained by considering the buoyancy component of the pressure–strain correlations (panels f,i,l). At weak stratification, Q1 and Q3 redistribution events (conversion of $k_{H}$ into $k_{w}$ ) dominate over Q2 and Q4 redistribution events (conversion of $k_{w}$ into $k_{H}$ ) (panel c), but the overall correlation (panel f) and normalised redistribution (panel l) are negligible due to the buoyancy effects being weak. As stratification increases until $\textit{Fr}_{k} \approx 0.7$ , we note that the probability masses of the four quadrants change insignificantly, but with stronger buoyancy effects, the correlation and normalised redistribution values now begin to contribute meaningfully to the overall redistribution, where, on average, the buoyancy pressure–strain correlation converts $k_{H}$ into $k_{w}$ . With further increasing stratification, the probability masses of the Q1 and Q3 redistribution events decrease, while the probability masses of the Q2 and Q4 redistribution events increase. This is related to a strengthening negative correlation associated with the events that convert $k_{w}$ into $k_{H}$ (blue triangles; panel f), which drives the overall correlation to change signs and become increasingly negative (black exes). Because the overall term changes sign, we compute the efficiency of the buoyancy pressure–strain correlation in two ways, where the positive values of the overall buoyancy pressure–strain component are normalised by the sum of the Q1 and Q3 events (red diamonds), and the negative values of the overall buoyancy pressure-strain component are normalised by the sum of the Q2 and Q4 events (blue diamonds). The peak efficiency of ${\sim} 0.4$ associated with redistribution from $k_{H}$ into $k_{w}$ occurs at $\textit{Fr}_{k} \approx 1$ , whereas the peak efficiency hovers at approximately $0.25$ for redistribution from $k_{w}$ into $k_{H}$ for $\textit{Fr}_{k} \lt 0.1$ . Before moving on, we reiterate that unlike the prevailing description of the nonlinear pressure–strain correlation as promoting large-scale isotropy, which is associated with a definite sign for the term, the buoyancy pressure–strain correlation exhibits more complex behaviour that promotes large-scale isotropy for weak to moderate stability, but then begins to promote large-scale anisotropy for stronger stability. Given that its overall effect is sensitive to the strength of the background stratification, this hints at similarities with and connections to the vertical buoyancy flux, which comprises both significant down-gradient fluxes, converting $k_{w}$ into $k_{p}$ , and counter-gradient fluxes, converting $k_{p}$ into $k_{w}$ . In § 4.4, we will mathematically demonstrate the link between the buoyancy pressure–strain redistribution and vertical buoyancy flux and explain this sign-change behaviour associated with increasing stability.

Figure 10. (a–c) Probability mass of the four quadrants of the pressure scrambling correlations (Q1, red triangles; Q2, orange triangles; Q3, blue squares; Q4, purple circles); (d–f) correlation of the events that generate buoyancy flux $B$ (red triangles), events that destroy buoyancy flux (blue triangles) and overall pressure scrambling term (black exes); (g–i) efficiency of the pressure scrambling term; (j–l) normalised pressure scrambling term associated with generation of $B$ , destruction of $B$ , net effect on $d_{t}B$ (same symbols/colours convention as panels d–f). All variables are plotted as a function of $\textit{Fr}_{k}$ . The results from the higher Reynolds number simulation (V2_128) are denoted by filled stars.

Now, we seek to understand how the contributions from the four quadrants of the pressure scrambling correlations contribute to the generation or destruction of the vertical buoyancy flux. In figure 10, we plot the following quantities as a function of $\textit{Fr}_{k}$ : (panels a,b,c) probability mass of the four quadrants of the pressure scrambling correlations (Q1, red triangles; Q2, orange triangles; Q3, blue squares; Q4, purple circles); (panels d,e, f) correlations of the events that generate $B$ (red triangles), events that destroy $B$ (blue triangles) and overall pressure scrambling correlation (black exes); (panels g,h,i) efficiency of the pressure scrambling correlation; (panels j,k,l) normalised pressure scrambling correlations associated with generation of $B$ (red), destruction of $B$ (blue) and the net effect (black). The associated contour plots of the joint p.d.f.s and covariance integrands are shown in figures 15 and 16 in Appendix B.

Recall that the net pressure scrambling term changes signs at $\textit{Fr}_{k} \approx 0.1$ as stability increases from near neutral conditions (cf. figure 10 j). This is associated with relatively small changes in the probability mass associated with Q1 and Q2 events but decreasing mass associated with Q4 events (destruction of $B$ ) and increasing mass associated with Q3 events (generation of $B$ ). These changes are broadly reflected also in the nonlinear and buoyancy components as well (see figure 10 b,c). Connectedly, the overall correlation of the pressure scrambling term (black exes in figure 10 d) also switches signs at $\textit{Fr}_{k} \approx 0.1$ , and this is due to both an increasing positive correlation of the sum of Q1 and Q3 events, which generate vertical buoyancy flux, and a reduced correlation of the sum of Q2 and Q4 events, which destroy vertical buoyancy flux. Looking at the correlation of the nonlinear component of pressure scrambling (panel e), we note that it is always negatively correlated, but this weakens due to an increasing contribution from positively correlated, nonlinear pressure scrambling events. The correlation of the buoyancy component of pressure scrambling is always positive (panel f), and the negatively correlated buoyancy pressure scrambling events remain negligible until $\textit{Fr}_{k} \lt 0.1$ .

Next, we consider the efficiencies associated with the three pressure scrambling correlations (panels g,h,i). For the total pressure scrambling term (figure 10 g), we define the efficiency by normalising the net pressure scrambling correlation by the sum of Q2 and Q4 events when the correlation is negative (blue diamonds), and by the sum of Q1 and Q3 events when the correlation is positive (red diamonds). The overall pressure scrambling is most efficient at neutral stability with an efficiency slightly above $60\,\%$ , but after the sign-change at $\textit{Fr}_{k} \approx 0.1$ , the efficiency generally increases with stability. For the nonlinear pressure scrambling term (figure 10 h), we note that the efficiency is largest for weak stratification, and monotonically decreases with increasing stability. For the buoyancy pressure scrambling term (figure 10 i), the efficiency increases with stability and reaches a maximum value of approximately $80\,\%$ before decreasing for $\textit{Fr}_{k} \lt 0.1$ . In closing, we stress that the overall buoyancy pressure scrambling term is sign definite with negligible contributions from the negatively correlated contributions until stability becomes strong ( $\textit{Fr}_{k} \lt 0.1$ ). We will seek to mathematically understand this behaviour in the next section.

4.4. Analytical expressions for the buoyancy component of pressure correlations and further physical interpretation

In this section, we seek to mathematically explain two particular observations about the buoyancy pressure correlations. First, as shown in figure 6(b), the buoyancy component of the pressure–strain correlation changes sign at strong stability and begins to redistribute $k_{w}$ into $k_{H}$ . This behaviour is in contrast to the physical picture underlying the Rotta-type, return-to-isotropy models that are often used to close the Reynolds stress equations (though non-Rotta models have been proposed for the buoyancy-linked redistribution; e.g. Hanjalić & Launder Reference Hanjalić and Launder2022). Second, the buoyancy component of the pressure scrambling seems sign definite (see figure 6 c), and it generally increases monotonically with stratification. This behaviour contrasts with the classic view of pressure scrambling terms as a sink of passive scalar fluxes, which decorrelate the velocity and scalar fields.

To explain these effects, we return to (4.3) and consider the Fourier coefficients of the buoyancy component of the pressure fluctuations,

(4.8) \begin{equation} \hat {p}_{{B}}(\boldsymbol{\kappa },t) = \frac {i g \kappa _{z} \hat {\rho }(\boldsymbol{\kappa },t)}{\kappa ^{2}}, \end{equation}

where $i = \sqrt {-1}$ is the imaginary unit and $\hat {\rho }(\boldsymbol{\kappa },t)$ are the Fourier coefficients associated with the density field. Next, we define $g_{1}(\boldsymbol{x},t) = ({1}/{\rho _{0}}) p_{{B}} \partial _{z} w$ and $g_{2}(\boldsymbol{x},t) =({g}/{\rho _{0}^{2}}) p_{{B}} \partial _{z} \rho$ , which are the unaveraged buoyancy components of the pressure–strain and pressure scrambling terms, respectively. Taking the forward Fourier transform of $g_{1}$ and $g_{2}$ in the three spatial dimensions and applying the Dirac delta function (cf. Appendix D of Pope Reference Pope2000) to collapse one of the integrals, we arrive at the following expressions:

(4.9) \begin{align} \hat {g}_{1}(\boldsymbol{\kappa },t) &= \frac {g}{\rho _{0}} \int {\textrm{d}}\boldsymbol{\kappa }' \left [-\frac {\kappa _{z}'(\kappa _{z} - \kappa _{z}')}{\kappa ^{\prime 2}} \right ] \hat {\rho }(\boldsymbol{\kappa }',t) \hat {w}(\boldsymbol{\kappa }-\boldsymbol{\kappa }',t), \end{align}

(4.10) \begin{align} \hat {g}_{2}(\boldsymbol{\kappa },t) &= \left ( \frac {g}{\rho _{0}} \right )^{2} \int {\textrm{d}}\boldsymbol{\kappa }' \left [-\frac {\kappa _{z}'(\kappa _{z} - \kappa _{z}')}{\kappa ^{\prime 2}} \right ] \hat {\rho }(\boldsymbol{\kappa }',t) \hat {\rho }(\boldsymbol{\kappa }-\boldsymbol{\kappa }',t), \end{align}

where $\hat {w}(\boldsymbol{\kappa },t)$ are the Fourier coefficients associated with the vertical velocity field. This derivation involves similar steps as in § 6.4.2 of Pope (Reference Pope2000), where the nonlinear advection term of the Navier–Stokes equations is expressed in Fourier space. Because we are interested in the volume-averaged effect of the buoyancy components of the pressure–strain and pressure scrambling terms, we set $\boldsymbol{\kappa } = 0$ in (4.9) and (4.10) and find

(4.11) \begin{align} R_{w,\textit{B}} &= \hat {g}_{1}(\boldsymbol{\kappa } = 0,t) = \frac {g}{\rho _{0}} \int {\textrm{d}}\boldsymbol{\kappa }' \left (\frac {\kappa _{z}'}{\kappa '} \right )^{2} \hat {\rho }(\boldsymbol{\kappa }',t) \hat {w}^{*}(\boldsymbol{\kappa }',t), \end{align}

(4.12) \begin{align} R_{B,\textit{B}} &= \hat {g}_{2}(\boldsymbol{\kappa } = 0,t) = \left ( \frac {g}{\rho _{0}} \right )^{2} \int {\textrm{d}}\boldsymbol{\kappa }' \left (\frac {\kappa _{z}'}{\kappa '} \right )^{2} \hat {\rho }(\boldsymbol{\kappa }',t) \hat {\rho }^{*}(\boldsymbol{\kappa }',t), \end{align}

where we have used the conjugate symmetry of Fourier coefficients of real functions to rewrite $\hat {w}(-\boldsymbol{\kappa }',t) = \hat {w}^{*}(\boldsymbol{\kappa }',t)$ and $\hat {\rho }(-\boldsymbol{\kappa }',t) = \hat {\rho }^{*}(\boldsymbol{\kappa }',t)$ .

Figure 11. Weighting factor $(\kappa _{z}/\kappa )^{2} = 1/[1 + (\kappa _{H}/\kappa _{z})^{2}]$ in (4.11) and (4.12) plotted as a function of (a) the horizontal and vertical wavenumbers ( $\kappa _{H}$ and $\kappa _{z}$ ) and (b) an aspect ratio based on the horizontal and vertical wavenumbers ( $\kappa _{H}/\kappa _{z}$ ).

There are several key aspects of (4.11) and (4.12). First, the volume-averaged buoyancy component of the pressure–strain correlation ( $R_{w,\textit{B}}$ ) (4.11) involves a weighted integral of the buoyancy flux cospectrum, where positive values of the integrand are associated with down-gradient vertical buoyancy fluxes, and negative values of the integrand are associated with counter-gradient vertical buoyancy fluxes. This reveals that down-gradient vertical buoyancy fluxes are connected to redistribution of $k_{H}$ into $k_{w}$ and that counter-gradient vertical buoyancy fluxes are connected to redistribution of $k_{w}$ into $k_{H}$ . Second, the volume-averaged buoyancy component of the pressure scrambling term ( $R_{B,\textit{B}}$ ) involves a weighted integral of the density variance spectrum (related to TPE). The integrand is quadratic, implying that the term is positive semidefinite (zero in the neutral limit). Note that without the wavenumber-dependent weighting factors, the integrals in (4.11) and (4.12) would be related to $\overline {w\rho }$ and $\overline {\rho \rho }$ . Third, the weighting factor $ (\kappa _{z}/\kappa )^{2}$ , which appears in both expressions, more strongly weights contributions from pancake structures ( $\kappa _{z} \gg \kappa _{H}$ , where $\kappa ^{2} = \kappa _{x}^{2} + \kappa _{y}^{2} + \kappa _{z}^{2} = \kappa _{H}^{2} + \kappa _{z}^{2}$ with $\kappa _{H} = \sqrt {\kappa _{x}^{2} + \kappa _{y}^{2}}$ ), which are prominent dynamical features of stratified turbulent flows (e.g. Lin & Pao Reference Lin and Pao1979; Hopfinger Reference Hopfinger1987; Riley & Lelong Reference Riley and Lelong2000; Billant & Chomaz Reference Billant and Chomaz2001; Lindborg Reference Lindborg2006; Shah & Bou-Zeid Reference Shah and Bou-Zeid2014; Caulfield Reference Caulfield2021; Kimura & Sullivan Reference Kimura and Sullivan2024). To aid this interpretation, we plot the weighting factor $(\kappa _{z}/\kappa )^{2}$ as a function of $\kappa _{H}$ and $\kappa _{z}$ in figure 11(a) and also as a function of an aspect ratio $\kappa _{H}/\kappa _{z}$ in figure 11(b) (or also $l_z/l_H$ , where $l_H$ and $l_z$ are the horizontal and vertical length scales associated with the horizontal and vertical wavenumbers, $\kappa _{H} = 2\pi /l_H$ and $\kappa _{z} = 2\pi /l_z$ ). This latter plot uses a re-expressed statement of the weighting factor where the numerator and denominator have been divided through by $\kappa _{z}^{2}$ . In both plots, the weighting factor approaches unity when $\kappa _{z} \gg \kappa _{H}$ (or equivalently, $l_z/l_H \ll 1$ ). The weighting factors in (4.11) and (4.12) reveal that the effect of stratification manifests in an anisotropic manner in Fourier space, which would be missed when only accounting for the wavenumber magnitude. Lastly, we note that (4.11) and (4.12) apply generally to all homogeneous stratified flows independent of forcing (e.g. decaying, forced, sheared, etc.) and across all stability conditions (neutral, stable, unstable). This fact helps explain why some relationships between $R_w$ , $R_B$ and $\varGamma$ have been observed across different types of stratified turbulent flows (e.g. Yi & Koseff Reference Yi and Koseff2022a , Reference Yi and Koseff2023).

With these points in mind, we return to figures 6(b) and 6(c) to interpret the buoyancy components of the pressure–strain and pressure scrambling terms. First, the buoyancy component of the pressure–strain correlation begins to decrease at $\textit{Fr}_{k} \approx 0.7$ and becomes increasingly negative for $\textit{Fr}_{k} \lt 0.3$ (see figure 6 b). Since the volume-averaged buoyancy flux remains positive for all values of $\textit{Fr}_k$ (i.e. $\overline {w\rho } \gt 0$ ), we note from (4.11) that $R_{w,\textit{B}} \lt 0$ must result from stronger counter-gradient buoyancy fluxes associated with layered structures that become more favourably weighted after the multiplication by $(\kappa _{z}/\kappa )^{2}$ . This means that, on average, counter-gradient buoyancy fluxes associated with layered motions contribute to redistribution of $k_{w}$ into $k_{H}$ , opposing large-scale return-to-isotropy. Relatedly, we observed more frequent and stronger contributions from counter-gradient fluxes in figure 8 with increasing stability. Because $R_{w,\textit{B}} \lt 0$ for $\textit{Fr}_{k} \lt 0.3$ is the result of volume and time averaging, this implies that there are persistent, large-scale, counter-gradient buoyancy fluxes, a result that was also observed by Almalkie & de Bruyn Kops (Reference Almalkie and de Bruyn Kops2012) for forced, homogeneous, stably stratified turbulence, and by Kaltenbach, Gerz & Schumann (Reference Kaltenbach, Gerz and Schumann1994) for stably stratified, homogeneous shear turbulence simulations for $\textit{Ri}_{g} = 0.5$ and $\textit{Pr} = 1$ . In the literature, these features have been treated distinctly from persistent, counter-gradient buoyancy fluxes at small scales (e.g. Holt, Koseff & Ferziger Reference Holt, Koseff and Ferziger1992; Kaltenbach et al. Reference Kaltenbach, Gerz and Schumann1994; Gerz & Schumann Reference Gerz and Schumann1996). Second, the buoyancy component of the pressure scrambling term is positive semidefinite and increases monotonically with stratification (see figure 6 c). From (4.12), we see that this is due to the integrand of $R_{B,{B}}$ being quadratic and vanishing in the limit of neutral stability.