# Lagrangian model for passive scalar gradients in turbulence

Xiaolong Zhang1  Maurizio Carbone2,3  Andrew D. Bragg1<EMAIL_ADDRESS>1Department of Civil and Environmental Engineering, Duke University, Durham,
NC 27708, USA 2Max Planck Institute for Dynamics and Self-Organization, Am
Faßberg 17, 37077 Göttingen, Germany 3Theoretical Physics I, University of
Bayreuth, Universitätsstr. 30, 95447 Bayreuth, Germany

###### Abstract

The equation for the fluid velocity gradient along a Lagrangian trajectory
immediately follows from the Navier-Stokes equation. However, such an equation
involves two terms that cannot be determined from the velocity gradient along
the chosen Lagrangian path: the pressure Hessian and the viscous Laplacian. A
recent model handles these unclosed terms using a multi-level version of the
recent deformation of Gaussian fields (RDGF) closure (Johnson & Meneveau,
Phys. Rev. Fluids, 2017). This model is in remarkable agreement with DNS data
and works for arbitrary Taylor Reynolds numbers $\Rey_{\lambda}$. Inspired by
this, we develop a Lagrangian model for passive scalar gradients in isotropic
turbulence. The equation for passive scalar gradients also involves an
unclosed term in the Lagrangian frame, namely the scalar gradient diffusion
term, which we model using the RDGF approach. However, comparisons of the
statistics obtained from this model with direct numerical simulation (DNS)
data reveal substantial errors due to erroneously large fluctuations generated
by the model. We address this defect by incorporating into the closure
approximation information regarding the scalar gradient production along the
local trajectory history of the particle. This modified model makes
predictions for the scalar gradients, their production rates, and alignments
with the strain-rate eigenvectors that are in very good agreement with DNS
data. However, while the model yields valid predictions up to around
$\Rey_{\lambda}\approx 500$, beyond this, the model breaks down.

###### keywords:

##  1 Introduction

Scalar transport in turbulence plays significant roles in various practical
applications, ranging from geophysical and environmental problems such as the
advection and dispersion of pollutants in the atmosphere (Nironi et al., 2015;
Mazzitelli & Lanotte, 2012), mixing of nutrients in the ocean (Bhamidipati et
al., 2020; Chor, 2020), and chemical reactions in industrial flows (Dimotakis,
2005; Hill, 1976). Scalar transport is also important from the perspective of
fundamental turbulence research, with previous studies showing that the scalar
field is a sensitive detector of the structures in turbulent flows, such that
their study has yielded insights into the physics of turbulent flows
themselves (Tong & Warhaft, 1994). The scalar field may feedback on the
velocity field under certain conditions, such as for stratified turbulence or
thermally driven turbulence (Zhang et al., 2022; Lohse & Xia, 2010). The focus
of this paper is, however, on the case of passive scalars.

The transport of scalars in turbulent flows is challenging to understand both
because of the complexity of the underlying turbulent flow that advects,
stretches and compresses the scalar field, and also because of the molecular
diffusion to which it is subject that leads to non-trivial differences
compared to the transport of fluid particles (Ottino & Ottino, 1989; Warhaft,
2000). Indeed, while it has often been assumed that the statistical properties
of passive scalars in turbulent flows should reflect the analogous properties
of the underlying velocity field, e.g. a similarity in the statistical
distribution of the turbulent kinetic energy and scalar dissipation rates,
this is not in general the case. For example, it has been found that
Kolmogorov’s hypothesis of local isotropy of the small scales of a turbulent
flow is strongly violated when applied to passive scalars, with both
experiments and direct numerical simulations (DNS) finding that the skewness
of the scalar derivative remains of the order of unity when a large-scale mean
scalar gradient is imposed, while it would be zero for a locally isotropic
flow (Sreenivasan & Antonia, 1977; Sreenivasan, 1991; Pumir, 1994; Mestayer et
al., 1976). This strong violation of small-scale isotropy of the scalar field
is often attributed to the presence of ramp-cliff structures in the scalar
field through which large and small scales of the scalar field are directly
connected (Buaria et al., 2021a; Shraiman & Siggia, 2000). It has also been
found that intermittency in the scalar field is even stronger than that for
the velocity field (Watanabe & Gotoh, 2004). Indeed, even for a stochastic
model where the velocity field has Gaussian statistics, the scalar field has
non-Gaussian statistics (Tong & Warhaft, 1994; Kraichnan, 1994; Falkovich et
al., 2001). Another profound difference is that while the velocity field
exhibits a dissipation anomaly (i.e. the averaged turbulent kinetic energy
dissipation rate is independent of viscosity for high Reynolds numbers), the
scalar field does not, with the scalar dissipation rate decreasing as $\sim
1/\log(Sc)$ as the Schmidt number $Sc$ is increased (Buaria et al., 2021b).
Interestingly, however, recent DNS results have shown that while the expected
correspondence between the analogous velocity and scalar statistics is not
observed for $Sc\leq 1$, it is recovered for sufficiently large $Sc$ (even for
$Sc=7$) where a viscous-convective sub-range emerges in the scalar field
(Shete et al., 2022).

In general, the properties of a passive scalar field depend upon both the
Reynolds number $\Rey$ and Schmidt number $Sc$, and in many practical
problems, both $\Rey$ and $Sc$ are large. For $Sc>1$, the smallest scale (in a
mean-field sense) in the scalar field is thought to be the Batchelor scale
(Batchelor, 1959) $\eta_{B}=Sc^{-1/2}\eta$, where $\eta$ is the Kolmogorov
length scale, and therefore resolving flows with high $\Rey$ and $Sc$ is very
challenging using DNS (as well as experiments) due to the spatial and temporal
resolution constraints. For problems where the small-scale properties of the
scalar field are important, large eddy simulations are not helpful. It is
therefore highly desirable to develop models for the small scales of the
scalar field that are capable of handling large ranges of $\Rey$ and $Sc$, as
well as being computationally efficient. One possibility is to develop
Lagrangian models for the scalar gradients in turbulent flows, inspired by the
corresponding models for the velocity gradient that have been highly
successful both in terms of making predictions and leading to new insights
into the small-scale dynamics of turbulent flows (Meneveau, 2011). Analogous
models for scalars could be used to explore and understand the small-scale
dynamics of scalar fields at $\Rey$ and $Sc$ that are currently far out of
reach using either DNS or experiments.

Lagrangian models for the velocity gradients are derived from the Navier-
Stokes equations, but they require modelling/approximations for the pressure
Hessian and viscous terms which are unclosed in the reference frame of a
single fluid particle trajectory. Various models have been proposed, including
the restricted Euler model (Vieillefosse, 1982), the tetrad model (Chertkov et
al., 1999), the recent fluid deformation model (Chevillard & Meneveau, 2006),
as well as closures based on random Gaussian fields (Wilczek & Meneveau, 2014;
Johnson & Meneveau, 2016) and, more generally, on tensor function
representation of the unclosed terms (Leppin & Wilczek, 2020). With the
exception of the restricted Euler model, these models predict the steady-
state, non-trivial properties of the velocity gradients in turbulent flows,
including the preferential alignment of the vorticity with the intermediate
strain-rate eigenvector, intermittency, and the multifractal scaling of the
moments of the velocity gradients. Unfortunately, most of the models are only
capable of predicting flows with low to moderate $\Rey$. Recently, the multi-
level recent deformation of Gaussian fields (ML-RDGF) model has been
developed, which was shown to predict DNS data accurately for arbitrary
Reynolds numbers (Johnson & Meneveau, 2017).

In the equation for a passive scalar gradient along a fluid particle
trajectory, the scalar gradient diffusion term is unclosed. A closed model for
scalar gradients in turbulent flows was previously derived based on a simple
linear relaxation model for the scalar gradient diffusion (Martín et al.,
2005), with the velocity gradient in the equation specified using a restricted
Euler model that was modified to include a linear damping term to model
viscous effects (Martın et al., 1998). Although the model showed general
qualitative agreement with DNS data, there were significant quantitative
inaccuracies for a number of key quantities, including significant errors in
the predictions for the alignments of the scalar gradients with the strain-
rate eigenvectors, and significant underprediction of large fluctuations of
the scalar gradients. More advanced models also have been proposed that use
the recent fluid deformation approximation (Chevillard & Meneveau, 2006) to
model the scalar gradient diffusion term (Hater et al., 2011; Gonzalez, 2009).
These models also yielded qualitatively reasonable predictions, but suffered
from quantitative inaccuracies and can only make predictions for relatively
low Reynolds numbers due to their use of the recent fluid deformation
approximation (Chevillard & Meneveau, 2006). Our work significantly advances
these models in two ways. First, the velocity gradients will be specified
using the much more sophisticated ML-RDGF model that also allows for
predictions to be made at arbitrarily large $\Rey$. Second, the recent
deformation of Gaussian fields closure (Johnson & Meneveau, 2016) will be used
to provide a more sophisticated closure for the scalar gradient diffusion
term. This closure leads to nonlinear terms that can regulate the growth of
the scalar gradients in regions of intense stretching, which can be vital to
preventing blow-ups from occurring in the model. The closure also incorporates
the effect of $Sc$ on the scalar gradient dynamics. In this way, the model is
capable of predicting the effect of both $\Rey$ and $Sc$ on the scalar
gradient dynamics.

##  2 Model for scalar gradients

###  2.1 Governing equations for instantaneous and characteristic variables

For an incompressible, Newtonian fluid, the equations governing the evolution
of the fluid velocity gradient $\bm{A}\equiv\bm{\nabla u}$ and scalar gradient
$\bm{B}\equiv\bm{\nabla}\phi$ are

$\displaystyle D_{t}\bm{A}$ $\displaystyle=-\bm{A\cdot
A}-\bm{H}+\nu\nabla^{2}\bm{A}+\bm{F_{A}},$ (1) $\displaystyle D_{t}\bm{B}$
$\displaystyle=-\bm{A}^{\top}\bm{\cdot B}+\kappa\nabla^{2}\bm{B}+\bm{F_{B}},$
(2)

where $D_{t}\equiv\partial_{t}+\bm{u\cdot\nabla}$ is the Lagrangian
derivative, $\bm{H}\equiv\bm{\nabla\nabla}p$ is the pressure Hessian, $p$ is
the fluid pressure (normalized by the fluid density), $\nu$ is the fluid
kinematic viscosity, $\kappa$ is the scalar diffusivity, and
$\bm{F_{A}},\bm{F_{B}}$ are forcing terms (assumed to be known/prescribed). In
the Lagrangian frame, the spatial derivative terms in (1) and (2) are unknown,
and therefore these terms must be modeled in terms of functionals of the
associated quantity, e.g $\nu\nabla^{2}\bm{A}$ must be modelled as some
functional of $\bm{A}$.

Concerning (1), various closure approaches have been developed, including the
Recent Fluid Deformation Approximation (RFDA) (Chertkov et al., 1999), Random
Gaussian Fields Closure (RGFC) (Wilczek & Meneveau, 2014), and the Recent
Deformation of Gaussian Fields (RDGF) closure (Johnson & Meneveau, 2016). More
recently, a multi-level version of the RDGF model has emerged (referred to as
ML-RDGF) that provides a model for $D_{t}\bm{A}$ that is valid for arbitrary
Reynolds numbers. These closures for (1) all lead to a model for $\bm{A}$ that
can generate steady-state statistics. This is in contrast to the Restricted
Euler (RE) model (Vieillefosse, 1982; Meneveau, 2011) that ignores the
anisotropic contribution to $\bm{H}$ and sets $\nu=0$, leading to a model for
$\bm{A}$ that exhibits a finite-time singularity. Subsequent studies found
that the addition of viscous effects to the RE model is not sufficient to
prevent a finite-time singularity; the anisotropic pressure Hessian must also
be accounted for. This could point to a potential challenge in closing (2); if
the closure approximation for $\kappa\nabla^{2}\bm{B}$ is not sufficiently
accurate then the predictions from the resulting model could also generate
finite-time singular solutions for $\bm{B}$ since the equation for $\bm{B}$
contains no pressure Hessian to regulate the amplification term $\bm{B\cdot
A}$ that causes $\|\bm{B}\|$ to grow when $\bm{B\cdot}(\bm{B\cdot A})<0$. In
this sense, developing a suitable closure for (2) may be more challenging than
that for (1).

The RGFC model (and also the RDGF and ML-RDGF models, since they are
extensions of the RGFC model) does not directly approximate the unclosed terms
in (1) (unlike the RFDA model) but instead closes the equation for
$d_{t}\bm{\mathcal{A}}$, which is defined as the solution to (1) averaged over
the subset of fluid particles that experience the same value of $\bm{A}$. The
statistics of $\bm{\mathcal{A}}$ correspond to the statistics of $\bm{A}$,
however, the advantage is that the terms requiring closure in the equation for
$d_{t}\bm{\mathcal{A}}$ are conditionally averaged quantities, and therefore
statistical approaches based on the properties of random Gaussian fields may
be employed to close the terms in the equation. Such a procedure cannot be
done when directly closing the terms in (1) since these involve instantaneous
quantities.

Mathematically, this procedure that replaces $D_{t}\bm{A}$ with
$d_{t}\bm{\mathcal{A}}$ may be formalized as follows. Suppose that $\bm{A}$
evolves in a phase-space with time-independent coordinates
$\bm{a}\in\mathbb{R}^{3\times 3}$. The PDF of $\bm{A}$ is then defined as
$\mathcal{P}(\bm{a},t)\equiv\langle\delta(\bm{A}-\bm{a})\rangle$ which solves
the Liouville equation

$\displaystyle\partial_{t}\mathcal{P}=-\bm{\nabla_{a}\cdot}\Big{(}\mathcal{P}\Big{\langle}D_{t}\bm{A}\Big{\rangle}_{\bm{A}=\bm{a}}\Big{)},$
(3)

where $\langle\cdot\rangle_{\bm{A}=\bm{a}}$ denotes an ensemble average
conditioned on $\bm{A}=\bm{a}$. We now introduce the time-dependent
characteristic variable $\bm{\mathcal{A}}$ defined via

$\displaystyle d_{t}\bm{\mathcal{A}}$
$\displaystyle\equiv\bm{G}(\bm{\mathcal{A}},t),$ (4)
$\displaystyle\bm{G}(\bm{a},t)$
$\displaystyle\equiv\Big{\langle}D_{t}\bm{A}\Big{\rangle}_{\bm{A}=\bm{a}}.$
(5)

Whereas $\bm{A}$ evolves according to the instantaneous Navier-Stokes
equation, $\bm{\mathcal{A}}$ evolves according to the conditionally averaged
Navier-Stokes equation. Nevertheless, since $\langle
D_{t}\bm{A}\rangle_{\bm{A}=\bm{a}}$ only involves a partial average over the
flow, the field $\bm{G}$ will in general exhibit nonlinear dependence on
$\bm{a}$ and $t$, and hence the trajectories $\bm{\mathcal{A}}$ generated by
(4) will still vary chaotically in time for a turbulent flow.

We may then define the PDF
$\varrho(\bm{a},t)\equiv\langle\delta(\bm{\mathcal{A}}-\bm{a})\rangle$ which
solves

$\displaystyle\partial_{t}\varrho$
$\displaystyle=-\bm{\nabla_{a}\cdot}\Big{(}\varrho\Big{\langle}d_{t}\bm{\mathcal{A}}\Big{\rangle}_{\bm{\mathcal{A}}=\bm{a}}\Big{)}.$
(6)

From the definition of the characteristic variable we have

$\displaystyle\Big{\langle}d_{t}\bm{\mathcal{A}}\Big{\rangle}_{\bm{\mathcal{A}}=\bm{a}}=\Big{\langle}\bm{G}(\bm{\mathcal{A}},t)\Big{\rangle}_{\bm{\mathcal{A}}=\bm{a}}=\bm{G}(\bm{a},t)=\Big{\langle}D_{t}\bm{A}\Big{\rangle}_{\bm{A}=\bm{a}},$
(7)

where the second equality follows since $\bm{G}(\bm{a},t)$ is not random. In
view of this, if $\varrho(\bm{a},0)=\mathcal{P}(\bm{a},0)$ then it follows
that $\varrho(\bm{a},t)=\mathcal{P}(\bm{a},t)\,\forall t$, and hence the
statistics of $\bm{A}$ may be obtained via solutions to (4).

Based on the definitions above we have

$\displaystyle d_{t}\bm{\mathcal{A}}$
$\displaystyle=-\bm{\mathcal{A}\cdot\mathcal{A}}-\langle\bm{H}\rangle_{\bm{\mathcal{A}}}+\nu\langle\nabla^{2}\bm{A}\rangle_{\bm{\mathcal{A}}}+\langle\bm{F_{A}}\rangle_{\bm{\mathcal{A}}},$
(8)

where we have used the short-hand notation

$\displaystyle\langle\cdot\rangle_{\bm{\mathcal{A}}}=\langle\cdot\rangle_{\bm{A}=\bm{a}}\Big{|}_{\bm{a}=\bm{\mathcal{A}}}.$
(9)

While closing the equation for $D_{t}\bm{A}$ requires approximating
instantaneous quantities, closing $d_{t}\bm{\mathcal{A}}$ requires
approximating conditionally averaged quantities. This is a major advantage
since it means that powerful statistical approaches can be used to develop
systematic closures for $d_{t}\bm{\mathcal{A}}$. This then is the approach
adopted in the RGFC and RDGF models, where
$\langle\bm{H}\rangle_{\bm{\mathcal{A}}}$ and
$\langle\nabla^{2}\bm{A}\rangle_{\bm{\mathcal{A}}}$ are closed under the
approximation that the velocity field $\bm{u}$ has Gaussian statistics.

A similar approach to this may be adopted for (2), but now the averaging must
be conditional on both $\bm{A}=\bm{a}$ and $\bm{B}=\bm{b}$, where
$(\bm{a},\bm{b})\in(\mathbb{R}^{3\times 3},\mathbb{R}^{3})$ are the time-
independent coordinates of the phase-space in which $\bm{A}$ and $\bm{B}$
evolve. We define

$\displaystyle d_{t}\bm{\mathcal{B}}$
$\displaystyle\equiv\bm{J}(\bm{\mathcal{A}},\bm{\mathcal{B}},t),$ (10)
$\displaystyle\bm{J}(\bm{a},\bm{b},t)$
$\displaystyle\equiv\Big{\langle}D_{t}\bm{B}\Big{\rangle}_{\bm{A}=\bm{a},\bm{B}=\bm{b}},$
(11)

and then based on (2) this leads to

$\displaystyle d_{t}\bm{\mathcal{B}}$
$\displaystyle=-\bm{\mathcal{A}}^{\top}\bm{\cdot\mathcal{B}}+\kappa\langle\nabla^{2}\bm{B}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}+\langle\bm{F_{B}}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}},$
(12) $\displaystyle\langle\cdot\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}$
$\displaystyle=\langle\cdot\rangle_{\bm{A}=\bm{a},\bm{B}=\bm{b}}\Big{|}_{\bm{a}=\bm{\mathcal{A}},\bm{b}=\bm{\mathcal{B}}}.$
(13)

An approach based on the RGFC or RDGF models may then be used to close the
term $\langle\nabla^{2}\bm{B}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}$
under the assumption that the scalar field $\phi$ has Gaussian statistics.
Just as the solutions to (4) can be used to construct the statistics of
$\bm{A}$, so also can the solutions to (10) be used to construct the
statistics of $\bm{B}$.

###  2.2 Model for the velocity gradients

In our model for $\bm{\mathcal{B}}$, the ML-RDGF model will be used to
prescribe $\bm{\mathcal{A}}$, and since our model for $\bm{\mathcal{B}}$ will
be based on the RDGF approach, we summarize the key ideas in this modeling
approach before applying them to deriving a closed equation for
$d_{t}\bm{\mathcal{B}}$

####  2.2.1 RDGF closure

The deformation tensor $\bm{D}(t,s)$ for a fluid particle with position
$\bm{x}^{f}(t|\bm{X},s)$ that satisfies $\bm{x}^{f}(s|\bm{X},s)=\bm{X}$ is
defined as

$\displaystyle\bm{D}(t,s)\equiv\frac{\partial}{\partial\bm{X}}\bm{x}^{f}(t|\bm{X},s),\quad
s\in[0,t],$ (14)

and evolves according to

$\displaystyle\partial_{t}\bm{D}=\bm{A\cdot D},\quad\bm{D}(0,0)=\mathbf{I},$
(15)

where $\mathbf{I}$ is the identity matrix. While the solution to this is given
by a time-ordered exponential, the RFDA approximates the solution by assuming
that $\bm{A}$ is constant over a recent-deformation timescale $\tau$, with the
deformation at times $s<t-\tau$ ignored. In this case, the approximate
solution is

$\displaystyle\bm{D}(t,s)$
$\displaystyle\approx\exp\Big{(}(t-s)\bm{A}\Big{)},\quad t-s\in[0,\tau].$ (16)

The idea then is as follows: the quantities $\bm{H}\equiv\bm{\nabla\nabla}p$
and $\nabla^{2}\bm{A}$ are evaluated along the fluid particle trajectory
$\bm{x}^{f}(t|\bm{X},s)$, and they may be related to their corresponding
values at the reference configuration $\bm{x}^{f}(s|\bm{X},s)=\bm{X}$ using
$\bm{D}$. Using the approximate solution for $\bm{D}$ in (16) and setting
$s=t-\tau$ leads to

$\displaystyle\langle\bm{H}\rangle_{\bm{\mathcal{A}}}$
$\displaystyle\approx\Big{\langle}\bm{D}^{-\top}\bm{\cdot}\Big{(}\frac{\partial^{2}}{\partial\bm{X}\partial\bm{X}}{p}\Big{)}\bm{\cdot}\bm{D}^{-1}\Big{\rangle}_{\bm{\mathcal{A}}}=\bm{\mathcal{D}}^{-\top}\bm{\cdot}\Big{\langle}\frac{\partial^{2}}{\partial\bm{X}\partial\bm{X}}{p}\Big{\rangle}_{\bm{\mathcal{A}}}\bm{\cdot}\bm{\mathcal{D}}^{-1},$
(17) $\displaystyle\langle\nabla^{2}\bm{A}\rangle_{\bm{\mathcal{A}}}$
$\displaystyle\approx\Big{\langle}\bm{D}^{-\top}\bm{\cdot}\Big{(}\frac{\partial^{2}}{\partial\bm{X}\partial\bm{X}}{\bm{A}}\Big{)}\bm{\cdot}\bm{D}^{-1}\Big{\rangle}_{\bm{\mathcal{A}}}=\bm{\mathcal{D}}^{-\top}\bm{\cdot}\Big{\langle}\frac{\partial^{2}}{\partial\bm{X}\partial\bm{X}}{\bm{A}}\Big{\rangle}_{\bm{\mathcal{A}}}\bm{\cdot}\bm{\mathcal{D}}^{-1},$
(18)

where $\bm{\mathcal{D}}=\exp(\tau\bm{\mathcal{A}})$, and $(\cdot)^{-\top}$
denotes the transpose of the inverse of a tensor. The RDGF model then uses the
RGFC approach (Wilczek & Meneveau, 2014) to approximate the conditional
averages in (17) and (18).

The basic motivation for the RDGF model is that the closure approximations for
$\langle\bm{H}\rangle_{\bm{\mathcal{A}}}$ and
$\langle\nabla^{2}\bm{A}\rangle_{\bm{\mathcal{A}}}$ can be improved by
applying the RGFC at $\bm{X},t-\tau$ rather than $\bm{x},t$. This is because
when the RGFC is applied at $\bm{X},t-\tau$ it is then transformed under the
flow map into something more realistic. For example, while the conditional
averages in (17) and (18) are approximated assuming that $\bm{u}$ has Gaussian
statistics, due to the Lagrangian transformation described by
$\bm{\mathcal{D}}$, the resulting approximations for
$\langle\bm{H}\rangle_{\bm{\mathcal{A}}}$ and
$\langle\nabla^{2}\bm{A}\rangle_{\bm{\mathcal{A}}}$ obtained through (17) and
(18) will, in general, correspond to those for a non-Gaussian field $\bm{u}$.

The final closed equation obtained using the RDGF closure has the form
(Johnson & Meneveau, 2016)

$\displaystyle
d\bm{\mathcal{A}}=\bm{\mathcal{N}}_{\mathcal{A}}\\{\bm{\mathcal{A}},\tau,\tau_{\eta}\\}dt+\bm{\Sigma\cdot}d\bm{\mathcal{W}},$
(19)

where the forcing term has been chosen to be
$\langle\bm{F_{A}}\rangle_{\bm{\mathcal{A}}}dt=\bm{\Sigma\cdot}d\bm{\mathcal{W}}$,
with $d\bm{\mathcal{W}}$ denoting a tensor-valued Wiener process, and
$\bm{\Sigma}$ denoting a diffusion tensor that depends on the coefficients
$D_{s}$ and $D_{a}$ which determine the growth rate of the mean-square values
of the strain-rate
$\bm{\mathcal{S}}\equiv(\bm{\mathcal{A}}+\bm{\mathcal{A}}^{\top})/2$ and
rotation-rate
$\bm{\mathcal{R}}\equiv(\bm{\mathcal{A}}-\bm{\mathcal{A}}^{\top})/2$ tensors.
For brevity, we do not include the details of the nonlinear operator
$\bm{\mathcal{N}}_{\mathcal{A}}\\{\cdot\\}$ which may be found in previous
works (Johnson & Meneveau, 2016).

The three unknown parameters $\tau,D_{a},D_{s}$ are obtained by an
optimization procedure that seeks those values for which the model satisfies
known constraints for isotropic turbulence, namely
$2\langle\|\bm{\mathcal{S}}\|^{2}\rangle=1/\tau_{\eta}^{2}$, and the two
homogeneity relations due to Betchov (1956),
$\langle\bm{\mathcal{A}}\bm{:\mathcal{A}}\rangle=0$,
$\langle(\bm{\mathcal{A}}\bm{\cdot\mathcal{A}})\bm{:\mathcal{A}}\rangle=0$.
The values obtained by this procedure are $\tau=0.1302\tau_{\eta}$,
$D_{s}=0.1014/\tau_{\eta}^{3}$, and $D_{a}=0.0505/\tau_{\eta}^{3}$.

####  2.2.2 ML-RDGF closure

The RDGF model described by (19) does not contain any dependence on
$\Rey_{\lambda}$. To address this, a multi-level version of the RDGF closure
(called ML-RDGF) was developed (Johnson & Meneveau, 2017). The key idea behind
this model is that in a turbulent flow, there exist velocity gradients at
different scales in the flow, and the velocity gradient dynamics at different
scales are coupled because of the energy cascade. Moreover, this coupling will
be influenced by the fact that the energy flux through the cascade is not
constant, but fluctuates in time and space. The ML-RDGF model extends the RDGF
model to take this into account by replacing (19) with

$\displaystyle
d\bm{\mathcal{A}}^{[n]}=\bm{\mathcal{N}}_{\mathcal{A}}\\{\bm{\mathcal{A}}^{[n]},\tau,\tau_{n}\\}dt-
d_{t}(\ln\tau_{n})\bm{\mathcal{A}}^{[n]}+\bm{\Sigma}^{[n]}\bm{\cdot}d\bm{\mathcal{W}}^{[n]},\quad
n=1,2,...N,$ (20)

in which the $\tau_{\eta}$ appearing in (19) has been replaced by the time-
dependent timescale $\tau_{n}(t)$, so that now $\tau=0.1302\tau_{n}$,
$D_{s}=0.1014/\tau_{n}^{3}$, and $D_{a}=0.0505/\tau_{n}^{3}$.

Equation (20) is actually a system of $N$ coupled equations, and
$\bm{\mathcal{A}}^{[n]}$ represents the velocity gradient at the $n^{th}$
level, corresponding to the velocity gradient filtered on some scale. For the
first level $n=1$, the timescale is fixed $\tau_{1}=\beta^{N-1}\tau_{\eta}$,
where $\beta=10$ is chosen in the model. For $n\geq 2$,
$\tau_{n}(t)\equiv\beta^{-1}\|\bm{\mathcal{S}}^{[n-1]}\|^{-1}$ (where
$\bm{\mathcal{S}}^{[n-1]}$ is the strain-rate associated with
$\bm{\mathcal{A}}^{[n-1]}$), such that the time evolution of
$\bm{\mathcal{A}}^{[n]}$ is coupled to the evolution at the larger scale where
the velocity gradient is $\bm{\mathcal{A}}^{[n-1]}$. The solution at level
$n=N$ then corresponds to the full (unfiltered) velocity gradient, i.e.
$\bm{\mathcal{A}}^{[N]}=\bm{\mathcal{A}}$. This multi-level model is capable
of predicting $\bm{\mathcal{A}}$ for the discrete Reynolds numbers
$\Rey_{\lambda}=\Rey_{\lambda}^{[n=1]}\beta^{N-1}$, where
$\Rey_{\lambda}^{[n=1]}$ is the Taylor Reynolds number corresponding to the
first level $n=1$, which was chosen in previous studies (Johnson & Meneveau,
2017) to be $\Rey_{\lambda}^{[n=1]}=60$. Then a modified timescale for the
second level, $\tau_{2}$, enables the model to predict flows at arbitrary
$\Rey_{\lambda}$.

The predictions for $\bm{\mathcal{A}}^{[N]}=\bm{\mathcal{A}}$ were shown to be
in excellent agreement with DNS and experimental data (Johnson & Meneveau,
2017), and revealed that the model makes robust predictions for the
intermittency of $\bm{\mathcal{A}}$ up to the highest Reynolds number
considered, $\Rey_{\lambda}=\textit{O}(10^{6})$.

###  2.3 Closure for scalar gradient equation based on RDGF

Based on its excellent performance, the ML-RDGF model will be used to specify
$\bm{\mathcal{A}}$ in the equation for the scalar gradient $\bm{\mathcal{B}}$.
The RDGF closure scheme will be used to close the scalar gradient diffusion
term $\langle\nabla^{2}\bm{B}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}$. A
multi-level version is not required since the effect of $\Rey_{\lambda}$ on
the model for $\bm{\mathcal{B}}$ will already be accounted for through the use
of the ML-RDGF to specify $\bm{\mathcal{A}}$ in the equation for
$\bm{\mathcal{B}}$.

Analogous to (17) and (18), under the recent fluid deformation approximation,
the scalar gradient diffusion term may be expressed as

$\displaystyle\langle\nabla^{2}\bm{B}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}\approx\bm{\mathcal{D}}^{-\top}\bm{\cdot}\Big{\langle}\frac{\partial^{2}}{\partial\bm{X}\partial\bm{X}}{\bm{B}}\Big{\rangle}_{\bm{\mathcal{A}},\bm{\mathcal{B}}}\bm{\cdot}\bm{\mathcal{D}}^{-1}.$
(21)

The Random Gaussian Fields Closure (RGFC) can be used to derive a closed
expression for the conditional average appearing in this expression, leading
through (21) to an RDGF closure for
$\langle\nabla^{2}\bm{B}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}$.

In the present context of the scalar field, the RGFC begins with the
assumption that $\phi$ has Gaussian statistics defined in terms of the
characteristic functional

$\displaystyle\Sigma^{\phi}[\lambda(\bm{x})]=\exp\Bigg{[}-\frac{1}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\lambda(\bm{x})R^{\phi}(\bm{x},\bm{x}^{\prime})\lambda(\bm{x}^{\prime})\,d\bm{x}\,d\bm{x}^{\prime}\Bigg{]},$
(22)

(time label is suppressed here and in what follows for simplicity) where
$\lambda(\bm{x})$ is the Fourier variable conjugate to $\phi(\bm{x})$, and
$R^{\phi}(\bm{x},\bm{x}^{\prime})\equiv\langle\phi(\bm{x})\phi(\bm{x}^{\prime})\rangle$,
which for isotropic turbulence has the form

$\displaystyle
R^{\phi}(\bm{x},\bm{x}^{\prime})=R^{\phi}(r)=\langle\phi^{2}\rangle
f_{\phi}(r),\quad\bm{r}\equiv\bm{x}-\bm{x}^{\prime},\quad r\equiv\|\bm{r}\|,$
(23)

where $f_{\phi}(r)$ is the scalar spatial correlation function. Following the
work of Wilczek & Meneveau (2014), the characteristic function for the scalar
gradient $\bm{B}\equiv\bm{\nabla}\phi$ is obtained from (22) as

$\displaystyle\Sigma^{\bm{B}}[\bm{\mu}(\bm{x})]=\exp\Bigg{[}-\frac{1}{2}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\bm{\mu}(\bm{x})\bm{\cdot}\bm{R}^{\bm{B}}(\bm{x},\bm{x}^{\prime})\bm{\cdot}\bm{\mu}(\bm{x}^{\prime})\,d\bm{x}\,d\bm{x}^{\prime}\Bigg{]},$
(24)

where $\bm{\mu}$ is the Fourier variable conjugate to $\bm{B}$, and
$\bm{R}^{\bm{B}}(\bm{x},\bm{x}^{\prime})\equiv\langle\bm{B}(\bm{x})\bm{B}^{\prime}(\bm{x}^{\prime})\rangle$,
which is related to $R^{\phi}(\bm{x},\bm{x}^{\prime})$ through

$\displaystyle\bm{R}^{\bm{B}}(\bm{x},\bm{x}^{\prime})=\langle\bm{\nabla}\phi(\bm{x})\bm{\nabla}^{\prime}\phi(\bm{x}^{\prime})\rangle=\bm{\nabla}\bm{\nabla}^{\prime}R^{\phi}(\bm{x},\bm{x}^{\prime}).$
(25)

For an isotropic scalar field where (23) applies, we then have

$\displaystyle\bm{R}^{\bm{B}}(\bm{x},\bm{x}^{\prime})=\bm{R}^{\bm{B}}(\bm{r})=\langle\phi^{2}\rangle\Bigg{[}\Bigg{(}\frac{f_{\phi}^{\prime}(r)}{r}-f_{\phi}^{{}^{\prime\prime}}(r)\Bigg{)}\frac{\bm{rr}}{r^{2}}-\frac{f_{\phi}^{\prime}(r)}{r}\mathbf{I}\Bigg{]},$
(26)

where prime superscripts denote differentiation with respect to $r$.

According to (24), the statistics of $\bm{B}$ are described by a Gaussian
characteristic functional when $\phi$ is assumed to be Gaussian. Just as in
previous works (Wilczek & Meneveau, 2014; Johnson & Meneveau, 2016) when
deriving a closure model for $D_{t}\bm{\mathcal{A}}$, it is acknowledged that
the statistics of $\bm{B}$ are not Gaussian in a real turbulent flow (nor are
they Gaussian even when $\bm{u}$ is Gaussian (Falkovich et al., 2001)). The
motivation for this choice is simply that it is the only option when deriving
a statistical field closure. It must be appreciated, however, that only with
respect to the closure of
$\langle\nabla^{2}\bm{B}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}$, are the
statistics of $\bm{B}$ approximated as being Gaussian; the model that follows
from this choice nevertheless generates statistics for $\bm{\mathcal{B}}$ that
are highly non-Gaussian (as will be shown later).

The specific term to be closed using RGFC is

$\displaystyle\Big{\langle}\frac{\partial^{2}}{\partial\bm{X}\partial\bm{X}}{\bm{B}}\Big{\rangle}_{\bm{a},\bm{b}}=\Big{\langle}\frac{\partial^{2}}{\partial\bm{X}\partial\bm{X}}{\bm{B}}(\bm{x}^{f}(t|\bm{X},s),t)\Big{\rangle}_{{\bm{A}}(\bm{x}^{f}(t|\bm{X},s),t)=\bm{a},{\bm{B}}(\bm{x}^{f}(t|\bm{X},s),t)=\bm{b}},$
(27)

which is evaluated at $\bm{a}=\bm{\mathcal{A}},\bm{b}=\bm{\mathcal{B}}$ in
(21). Just as the RFDA assumes
${\bm{A}}(\bm{x}^{f}(t|\bm{X},s),t)\approx{\bm{A}}(\bm{X},s)$ for
$t-s\in[0,\tau]$, we also assume
${\bm{B}}(\bm{x}^{f}(t|\bm{X},s),t)\approx{\bm{B}}(\bm{X},s)$ for
$t-s\in[0,\tau]$, and inserting this yields

$\displaystyle\Big{\langle}\frac{\partial^{2}}{\partial\bm{X}\partial\bm{X}}{\bm{B}}\Big{\rangle}_{\bm{a},\bm{b}}\approx\Big{\langle}\frac{\partial^{2}}{\partial\bm{X}\partial\bm{X}}{\bm{B}}(\bm{X},s)\Big{\rangle}_{{\bm{A}}(\bm{X},s)=\bm{a},{\bm{B}}(\bm{X},s)=\bm{b}},$
(28)

which puts the conditional average into a form to which the RGFC procedure
(Wilczek & Meneveau, 2014) can be applied.

Before proceeding, we note that although the argument of the conditional
average

$\Big{\langle}\frac{\partial^{2}}{\partial\bm{X}\partial\bm{X}}{\bm{B}}(\bm{X},s)\Big{\rangle}_{{\bm{A}}(\bm{X},s)=\bm{a},{\bm{B}}(\bm{X},s)=\bm{b}}$

does not contain $\bm{A}$ but only $\bm{B}$, the conditionality on
$\bm{A}(\bm{X},s)=\bm{a}$ cannot formally be removed. This is because since
the evolution of $\bm{B}$ depends upon $\bm{A}$ (see (2)), then the value of

$\frac{\partial^{2}}{\partial\bm{X}\partial\bm{X}}{\bm{B}}$

will not be uncorrelated from $\bm{A}$, in general. While a closure for the
full conditional average can be obtained using the RDGF approach, in order to
simplify the closure analysis the following approximation is made

$\displaystyle\Big{\langle}\frac{\partial^{2}}{\partial\bm{X}\partial\bm{X}}{\bm{B}}(\bm{X},s)\Big{\rangle}_{{\bm{A}}(\bm{X},s)=\bm{a},{\bm{B}}(\bm{X},s)=\bm{b}}\approx\Big{\langle}\frac{\partial^{2}}{\partial\bm{X}\partial\bm{X}}{\bm{B}}(\bm{X},s)\Big{\rangle}_{{\bm{B}}(\bm{X},s)=\bm{b}}.$
(29)

It is crucial to note, however, that the overall closure for
$\langle\nabla^{2}\bm{B}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}$ does
partially include the effect of the conditioning upon $\bm{\mathcal{A}}$,
since the right hand side of (21) contains $\bm{\mathcal{D}}$ and not $\bm{D}$
precisely because of the conditionality in the averaging operator on the left
hand side of (21). Therefore, the closure for
$\langle\nabla^{2}\bm{B}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}$ captures
some of the dependency of the scalar diffusion on the local velocity gradients
in the flow.

With the statistics of $\bm{B}$ prescribed using (24) and the covariance
tensor for $\bm{B}$ prescribed using (26), a closure for (29), and hence
$\langle\nabla^{2}\bm{B}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}$, may be
obtained following the same steps as in the work of Wilczek & Meneveau (2014).
The basic steps are as follows. First, using the approach described in
Appendix A of the work of Wilczek & Meneveau (2014), the unclosed quantity is
re-written in terms of a two-point quantity

$\displaystyle\Big{\langle}\frac{\partial^{2}}{\partial\bm{X}\partial\bm{X}}{\bm{B}}(\bm{X},s)\Big{\rangle}_{{\bm{B}}(\bm{X},s)=\bm{b}}=\lim_{r\to
0}\frac{\partial^{2}}{\partial\bm{r}\partial\bm{r}}\Big{\langle}\bm{B}(\bm{X}+\bm{r},s)\Big{\rangle}_{{\bm{B}}(\bm{X},s)=\bm{b}}.$
(30)

Applying the steps outlined in Appendix B of Wilczek & Meneveau (2014) to the
scalar field we obtain

$\displaystyle\Big{\langle}\bm{B}(\bm{X}+\bm{r},s)\Big{\rangle}_{{\bm{B}}(\bm{X},s)=\bm{b}}=\bm{R}^{\bm{B}}(\bm{r})\bm{\cdot}[\bm{R}^{\bm{B}}(\bm{0})]^{-1}\bm{\cdot
b},$ (31)

where

$\displaystyle[\bm{R}^{\bm{B}}(\bm{0})]^{-1}=\frac{\mathbf{I}}{\langle\phi^{2}\rangle
f_{\phi}^{{}^{\prime\prime}}(0)}.$ (32)

The resulting expression

$\displaystyle\Big{\langle}\frac{\partial^{2}}{\partial\bm{X}\partial\bm{X}}{\bm{B}}(\bm{X},s)\Big{\rangle}_{{\bm{B}}(\bm{X},s)=\bm{b}}=\lim_{r\to
0}\Bigg{(}\frac{\partial^{2}}{\partial\bm{r}\partial\bm{r}}\bm{R}^{\bm{B}}(\bm{r})\Bigg{)}\bm{\cdot}[\bm{R}^{\bm{B}}(\bm{0})]^{-1}\bm{\cdot
b},$ (33)

may then be computed using (26), and when the result is evaluated at
$\bm{b}=\bm{\mathcal{B}}$ and substituted into (21), the closed expression
obtained is

$\displaystyle\kappa\langle\nabla^{2}\bm{B}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}$
$\displaystyle\approx\delta_{\mathcal{B}}\Big{(}\bm{\mathcal{C}}_{R}^{-1}\bm{\cdot\mathcal{B}}+\bm{\mathcal{C}}_{R}^{-\top}\bm{\cdot\mathcal{B}}+\mathrm{tr}(\bm{\mathcal{C}}_{R}^{-1})\bm{\mathcal{B}}\Big{)},$
(34) $\displaystyle\delta_{\mathcal{B}}$
$\displaystyle\equiv\frac{\kappa}{3}\frac{f_{\phi}^{{}^{\prime\prime\prime\prime}}(0)}{f_{\phi}^{{}^{\prime\prime}}(0)},$
(35)

where
$\bm{\mathcal{C}}_{R}^{-1}\equiv\bm{\mathcal{D}}^{-1}\bm{\cdot}\bm{\mathcal{D}}^{-\top}$
is the inverse of the right Cauchy-Green tensor.

In the work of Johnson & Meneveau (2016), the coefficient analogous to
$\delta_{\mathcal{B}}$ in the closure for
$\langle\nabla^{2}\bm{A}\rangle_{\bm{\mathcal{A}}}$, namely
$\delta_{\mathcal{A}}$, was estimated based on the enstrophy production-
dissipation balance at steady-state. The same procedure can be applied to
approximate $\delta_{\mathcal{B}}$ based on the steady-state production-
dissipation balance $\langle\bm{S:BB}\rangle=-\kappa\langle\|\bm{\nabla
B}\|^{2}\rangle$, leading to

$\displaystyle\delta_{\mathcal{B}}$
$\displaystyle\approx\frac{1}{9\tau_{\eta}}\mathrm{tr}(\bm{\mathcal{C}}_{L})\gamma_{\mathcal{B}},$
(36) $\displaystyle\gamma_{\mathcal{B}}$
$\displaystyle\equiv\tau_{\eta}\frac{\langle\bm{S:BB}\rangle}{\langle\|\bm{B}\|^{2}\rangle},$
(37)

where
$\bm{\mathcal{C}}_{L}\equiv\bm{\mathcal{D}}\bm{\cdot}\bm{\mathcal{D}}^{\top}$
is the left Cauchy-Green tensor.

Whereas the RDGF model for $\langle\nabla^{2}\bm{A}\rangle_{\bm{\mathcal{A}}}$
is nonlinear in $\bm{\mathcal{A}}$ due to the contributions from
$\bm{\mathcal{C}}_{L}$ and $\bm{\mathcal{C}}_{R}$, the closure for
$\langle\nabla^{2}\bm{B}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}$ in (34)
is linear in $\bm{\mathcal{B}}$. This could potentially suggest an issue with
(34) since it is known that models for
$\langle\nabla^{2}\bm{A}\rangle_{\bm{\mathcal{A}}}$ that are linear in
$\bm{\mathcal{A}}$ can suffer from finite-time singularities, depending on the
initial conditions (Martın et al., 1998). However, given that the production
term in the equation for $\bm{\mathcal{B}}$ is
$-\bm{\mathcal{A}}^{\top}\bm{\cdot\mathcal{B}}$, then the dependence of the
closure in (34) on $\bm{\mathcal{C}}_{L}$ and $\bm{\mathcal{C}}_{R}$ may
indirectly prevent singular growth of $\|\bm{\mathcal{B}}\|$, at least in
regions where this growth is associated with large values of
$\|\bm{\mathcal{A}}\|$, because in those regions the growth of
$\bm{\mathcal{B}}$ could be modulated through $\bm{\mathcal{C}}_{L}$ and
$\bm{\mathcal{C}}_{R}$. Moreover, we performed tests where in the closure for
$\langle\nabla^{2}\bm{B}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}$, we set
$\bm{\mathcal{D}}=\mathbf{I}$, i.e. so that the recent deformation mapping was
removed. Simulations of the model using this blew up, showing that the
contributions in (34) involving $\bm{\mathcal{D}}$ do indeed play a key
indirect role in preventing singular growth of $\bm{\mathcal{B}}$.

In view of (36), any dependence of the closed expression for
$\kappa\langle\nabla^{2}\bm{B}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}$ on
$\kappa$ is contained within $\gamma_{\mathcal{B}}$, and this term must be
specified using data. However, in order for the model to depend upon $\kappa$
and hence $Sc$, then the data must specify $\gamma_{\mathcal{B}}$ as a
function of $Sc$, which is not desirable. An alternative approach is as
follows: The closure approximation in (34) is linear in $\bm{\mathcal{B}}$,
and $1/\delta_{\mathcal{B}}$ may be regarded as a linear relaxation timescale
for $\bm{\mathcal{B}}$, describing how the diffusion term
$\kappa\langle\nabla^{2}\bm{B}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}$
causes $\bm{\mathcal{B}}$ to relax to its equilibrium value. In view of this,
$1/\delta_{\mathcal{B}}$ should scale with the small-scale scalar timescale
$\tau_{\phi}$, and for $Sc\geq 1$ this timescale can be estimated using the
Batchelor scale (Batchelor, 1959; Donzis et al., 2005) as
$\tau_{\phi}\sim\tau_{\eta}Sc^{-1/3}$, while for $Sc<1$ the Corrsin scale
(Corrsin, 1951) may be used to obtain $\tau_{\phi}\sim\tau_{\eta}Sc^{-1/2}$ .
In view of this, the $Sc$ dependence may be explicitly accounted for in
$\delta_{\mathcal{B}}$ by using

$\displaystyle\delta_{\mathcal{B}}$
$\displaystyle\approx\frac{Sc^{-\xi}}{9\tau_{\eta}}\mathrm{tr}(\bm{\mathcal{C}}_{L})\alpha_{\mathcal{B}},$
(38) $\displaystyle\alpha_{\mathcal{B}}$
$\displaystyle\equiv\tau_{\eta}\frac{\langle\bm{S:BB}\rangle}{\langle\|\bm{B}\|^{2}\rangle}\Bigg{|}_{Sc=1},$
(39)

where $\xi=1/2$ for $Sc<1$, and $\xi=1/3$ for $Sc\geq 1$. From our DNS (see
§3) we obtain the value $\alpha_{\mathcal{B}}\approx-0.32$. However, in
anticipation of results to be shown later, we note that using this fixed value
for $\alpha_{B}$ yields a model whose predictions are not accurate (see figure
4). Therefore, a modification will be introduced wherein
$\alpha_{\mathcal{B}}$ is specified based on the local scalar gradient
production along the trajectory history of the particle, as in equation (46),
which dramatically improves the predictions from the model (cf. figure 6).

Analogous to the equation for $\bm{\mathcal{A}}$, the forcing term in the
scalar gradient equation is chosen to be
$\langle\bm{F_{B}}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}dt=\sigma
d\bm{\mathcal{W}}$, where now $d\bm{\mathcal{W}}$ is a vector-valued Wiener
process with increments defined by

$\displaystyle\langle d\bm{\mathcal{W}}\rangle$ $\displaystyle=\bm{0},$ (40)
$\displaystyle\langle d\bm{\mathcal{W}}d\bm{\mathcal{W}}\rangle$
$\displaystyle=\mathbf{I}dt.$ (41)

Since the equation for $\bm{\mathcal{B}}$ is linear, the statistics scale with
the forcing amplitude $\sigma$, and therefore $\sigma$ may be chosen
arbitrarily if the results generated by the model are suitably normalized.
Using this forcing term, the final form of the model equation is written as a
stochastic differential equation

$\displaystyle d\bm{\mathcal{B}}$
$\displaystyle\approx-\bm{\mathcal{A}}^{\top}\bm{\cdot\mathcal{B}}dt+\delta_{\mathcal{B}}\Big{(}\bm{\mathcal{C}}_{R}^{-1}\bm{\cdot\mathcal{B}}+\bm{\mathcal{C}}_{R}^{-\top}\bm{\cdot\mathcal{B}}+\mathrm{tr}(\bm{\mathcal{C}}_{R}^{-1})\bm{\mathcal{B}}\Big{)}dt+\sigma
d\bm{\mathcal{W}}.$ (42)

For applications where there is a mean scalar gradient, the forcing term could
be replaced by the term describing the production of fluctuating scalar
gradients due to the imposed mean scalar gradient.

Finally, just as the ML-RDGF model replaces the constant timescale for
$\bm{\mathcal{A}}^{[N]}$ with the fluctuating timescale $\tau_{N}(t)$, the
Kolmogorov timescale $\tau_{\eta}$ that appears in (38) should also for
consistency be replaced by $\tau_{N}(t)$ (since the ML-RDGF model is being
used to specify $\bm{\mathcal{A}}$ in (42)). In other words (38) is to be
replaced by

$\displaystyle\delta_{\mathcal{B}}$
$\displaystyle\approx\frac{Sc^{-\xi}}{9\tau_{N}(t)}\mathrm{tr}(\bm{\mathcal{C}}_{L})\alpha_{\mathcal{B}}.$
(43)

This ensures that for $Sc=1$, the local timescale on which both
$\bm{\mathcal{B}}$ and $\bm{\mathcal{A}}$ fluctuate is $O(\tau_{N}(t))$.

##  3 Numerical Simulations

We test the model predictions against data from Direct Numerical Simulation
(DNS) of a passive scalar field advected by an incompressible, three-
dimensional, statistically steady and isotropic turbulent velocity field. The
DNS code uses a standard Fourier pseudo-spectral method (Canuto, 1988) to
solve the discretized Navier-Stokes and passive scalar equations on a triply-
periodic cubic domain. The required Fourier transforms are executed in
parallel using the P3DFFT library (Pekurovsky, 2012) and the aliasing error is
removed via the $3/2$ rule (Canuto, 1988). The code is described in further
detail in Carbone et al. (2019).

The non-dimensional governing equation for the passive scalar in Fourier space
reads

$\displaystyle\partial_{t}\hat{\phi}+\textrm{i}\bm{k}\bcdot\widehat{\bm{u}\phi}=-\kappa\|\bm{k}\|^{2}\hat{\phi}+\hat{F},$
(44)

where the hat indicates a Fourier transform, “i” is the imaginary unit,
$\bm{u}(\bm{x},t)$ is the turbulent velocity field, $\phi(\bm{x},t)$ is the
passive scalar field and $\hat{\phi}(\bm{k},t)$ denotes its spatial Fourier
transform. For consistency with the stochastic model, the passive scalar field
is driven by a stochastic forcing, such that the complex forcing $\hat{F}$ is
Gaussian and white in time, with correlation

$\displaystyle\left\langle\hat{F}(\bm{k},t)\hat{F}(\bm{k}^{\prime},t^{\prime})\right\rangle=2\sigma_{0}^{2}\|\bm{k}\|^{2}\delta(\bm{k}+\bm{k}^{\prime})\delta(t-t^{\prime}),\;0<\|\bm{k}\|<k_{f}.$
(45)

The forcing is confined to the wavevectors $\bm{k}$ within a sphere of radius
$k_{f}$, and we choose $k_{f}=\sqrt{7}$ since it yields large-scale statistics
that are close to being isotropic. Finally, the spectrum of the forcing
$|\hat{F}|^{2}(k)$ scales as $\|\bm{k}\|^{2}$, compatible with energy
equipartition among the smallest Fourier modes. The constant parameter
$\sigma_{0}$ in (45) regulates the dissipation rate and it is adjusted to
yield an appropriate Kolmogorov scale $\eta$ prescribed by the spatial
resolution requirements. We simulate a flow with $\Rey_{\lambda}=100$ and
$Sc=1$. The spatial resolution is $k_{\max}\eta\simeq 3$, with $k_{\max}=N/2$
being the maximum resolved wavenumber, for all three simulations. The time
integration of equation (44) is performed by means of a second-order Runge-
Kutta scheme designed for stochastic differential equations (Honeycutt, 1992)
and the CFL number stays below 0.3.

Regarding the numerical simulations of the model, we solve the scalar gradient
equation (42) through a second-order predictor-corrector method (Kloeden &
Platen, 2018) with a time step $dt=0.05\tau_{\eta}$. Each level of the ML-RDGF
model equation was solved using its own appropriate time step, namely, level
$n$ was solved using $dt=0.05\langle\tau_{n}(t)\rangle$. Tests were performed
using smaller time steps and these tests indicated that the aforementioned
time steps were small enough to achieve convergence of the results.

We numerically solve the model equation (42) for $\Rey_{\lambda}=100,300,500$
and the results from these will be referred to as M1, M2, and M3 in the
results section. While M1 is designed to match the DNS, M2 and M3 will be used
to explore the model’s ability to capture the effect of $\Rey_{\lambda}$ on
the scalar gradient statistics.

##  4 Results and Discussion

###  4.1  Velocity Gradients

In our scalar gradient model, the velocity gradients are specified using the
ML-RDGF model(Johnson & Meneveau, 2017). We therefore begin by comparing the
predictions of this model against the DNS data in order to assess its accuracy
in predicting the statistics $\bm{\mathcal{A}}$, since any inaccuracies in
this model will in turn lead to inaccuracies in our model for the scalar
gradients.

\begin{overpic}[trim=14.22636pt 170.71652pt 0.0pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/A11_PDF_all_Re_lambda_withDNS.pdf}
\put(70.0,0.0){$a_{11}/a_{11,rms}$} \put(-13.0,60.0){\rotatebox{90.0}{PDF}}
\put(122.0,115.0){DNS} \put(122.0,107.0){M1} \put(122.0,98.0){M2}
\put(122.0,89.0){M3} \end{overpic} (a)  \begin{overpic}[trim=14.22636pt
170.71652pt 0.0pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/A12_PDF_all_Re_lambda_withDNS.pdf}
\put(70.0,0.0){$a_{12}/a_{12,rms}$} \put(-13.0,60.0){\rotatebox{90.0}{PDF}}
\put(122.0,115.0){DNS} \put(122.0,107.0){M1} \put(122.0,98.0){M2}
\put(122.0,89.0){M3} \end{overpic} (b)

Figure 1: PDFs of (a) longitudinal $a_{11}/a_{11,rms}$ and (b) transverse
$a_{12}/a_{12,rms}$ velocity gradients.

In Figure 1 we compare the PDFs of the longitudinal and transverse components
of the velocity gradient predicted by the ML-RDGF model with the DNS results.
(Recall that $\bm{a}$ and $\bm{b}$ are the phase-space variables conjugate to
$\bm{\mathcal{A}}$ and $\bm{\mathcal{B}}$, respectively, and we use the
subscripts “$rms$” and “$av$” to denote the root-mean-square and averaged
values of the variable under consideration). One can see that the PDFs for the
longitudinal gradients (Figure 1(a)) are negatively skewed while the PDFs for
the transverse gradients are symmetric (Figure 1(b)). The former is associated
with the self-amplification of the velocity gradients (Tsinober, 2001), while
the latter is a constraint due to isotropy. The predictions from M1 are in
excellent agreement with the DNS data, and the model makes realistic
predictions at the higher $\Rey_{\lambda}$ to which M2 and M3 correspond. The
ability of the ML-RDGF model to accurately predict the components of the
velocity gradient, and realistic intermittency trends with increasing
$\Rey_{\lambda}$ (and over a much larger range of $\Rey_{\lambda}$ than
considered here) were previously demonstrated in detail in the original ML-
RDGF paper of Johnson & Meneveau (2017).

Further insight into the ability of the ML-RDGF model to predict the velocity
gradient dynamics can be obtained by considering its predictions for the
velocity gradient invariants $Q\equiv-\bm{a}\bm{:a}/2$ and
$R\equiv-(\bm{a\cdot}\bm{a})\bm{:a}/3$. The invariant $Q$ measures the
relative strength of the local strain rate and vorticity in the flow, with
$Q>0$ denoting vorticity-dominated regions of the flow, while $R$ measures the
relative importance of the local strain self-amplification and enstrophy
production, with $R<0$ denoting regions dominated by enstrophy production
(Tsinober, 2001). The joint PDFs of $Q$ and $R$ from the ML-RDGF model for M1,
M2, M3 are presented in Figure 2, along with the DNS results for comparison.
As observed previously (Johnson & Meneveau, 2017), the ML-RDGF captures the
main features of the $Q,R$ joint-PDF well, including the signature sheared-
drop shape of the joint-PDF, which is preserved by the model as
$\Rey_{\lambda}$ is increased. The joint-PDF contours extend along the right
Viellefosse tail, and the PDF is concentrated in the quadrants $Q>0,R<0$ and
$Q<0,R>0$, consistent with previous studies (Johnson & Meneveau, 2016, 2017).
For M1, the joint-PDF is more compact compared to the DNS, indicating that the
model underpredicts the probability of large values of $Q,R$. The agreement
between the ML-RDGF model and DNS data was shown to be better in Johnson &
Meneveau (2017), however, their DNS data was for $\Rey_{\lambda}=430$,
indicating that the quantitive accuracy of the model is better at higher
$\Rey_{\lambda}$. Considering the results for M2 and M3 shows that as
$\Rey_{\lambda}$ increases, the contours of the joint-PDF spread out in the
$Q,R$ phase-space, showing that the model captures the effect of
$\Rey_{\lambda}$ on the intermittency in the flow.

\begin{overpic}[trim=0.0pt 170.71652pt 28.45274pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_1_2_JPDF_DNS.pdf}
\put(75.0,-1.0){$R/\tilde{Q}_{av}^{3/2}$}
\put(8.0,55.0){\rotatebox{90.0}{$Q/\tilde{Q}_{av}$}} \end{overpic} (a)
\begin{overpic}[trim=0.0pt 170.71652pt 28.45274pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_1_2_JPDF_Re_lambda_100.pdf}
\put(75.0,-1.0){$R/\tilde{Q}_{av}^{3/2}$}
\put(8.0,55.0){\rotatebox{90.0}{$Q/\tilde{Q}_{av}$}} \end{overpic} (b)

\begin{overpic}[trim=0.0pt 170.71652pt 28.45274pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_1_2_JPDF_Re_lambda_300.pdf}
\put(75.0,-1.0){$R/\tilde{Q}_{av}^{3/2}$}
\put(8.0,55.0){\rotatebox{90.0}{$Q/\tilde{Q}_{av}$}} \end{overpic} (c)
\begin{overpic}[trim=0.0pt 170.71652pt 28.45274pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_1_2_JPDF_Re_lambda_500.pdf}
\put(75.0,-1.0){$R/\tilde{Q}_{av}^{3/2}$}
\put(8.0,55.0){\rotatebox{90.0}{$Q/\tilde{Q}_{av}$}} \end{overpic} (d)

Figure 2: Logarithm of joint-PDFs of $Q/\tilde{Q}_{av}$ and
$R/\tilde{Q}_{av}^{3/2}$ (where
$\tilde{Q}_{av}\equiv\langle\|\mathcal{A}\|^{2}\rangle$) from (a) DNS, (b) M1,
(c) M2, (d) M3. Colours indicate the values of the logarithm of the PDF.

A more careful and quantitative comparison between the model predictions and
DNS data for $Q,R$ can be made by comparing the PDFs of $Q$ and $R$
separately. In Figure 3(a), one can see that in the DNS, the PDF of $Q$ is
strongly positively skewed, which is associated with the fact that the
vorticity is more intermittent than the strain-rate in turbulent flows (Yeung
et al., 2018). By contrast, the ML-RDGF model predicts PDFs for $Q$ that are
much more symmetric, with the model significantly underpredicting regions of
intense vorticity, and slightly underpredicting regions of intense straining.
These discrepancies are, however, mainly in the tails of the PDF, with the
model predictions in good agreement with the DNS for values of the PDF that
are $\geq O(10^{-2})$.

In Figure 3(b) it is seen that the model also predicts PDFs of $R$ that are
relatively symmetric compared with the DNS data for which the PDF is
negatively skewed. The model prediction for M2, which corresponds to
$\Rey_{\lambda}=300$, is much closer to the DNS data than those of M1 which
has the same $\Rey_{\lambda}=100$ as the DNS.

In view of these results, it is seen that while the ML-RDGF model predicts the
components of the velocity gradients very accurately (as well as the
alignments between the vorticity and strain-rate eigenvectors(Johnson &
Meneveau, 2017)), it does not predict the invariants of the velocity gradients
accurately when compared with the DNS, at least not for the $\Rey_{\lambda}$
considered here. This could in turn lead to inaccuracies in the scalar
gradient model, above and beyond any arising from the closure approximations
for the scalar gradient diffusion term. These results point to the need for
further refinements in the ML-RDGF model.

\begin{overpic}[trim=14.22636pt 170.71652pt 0.0pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_1_PDF_all_Re_lambda.pdf}
\put(80.0,-2.0){$Q/\tilde{Q}_{av}$} \put(-13.0,60.0){\rotatebox{90.0}{PDF}}
\put(122.0,115.0){DNS} \put(122.0,107.0){M1} \put(122.0,98.0){M2}
\put(122.0,89.0){M3} \end{overpic} (a)  \begin{overpic}[trim=14.22636pt
170.71652pt 0.0pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_2_PDF_all_Re_lambda.pdf}
\put(80.0,-2.0){$R/\tilde{Q}_{av}^{3/2}$}
\put(-13.0,60.0){\rotatebox{90.0}{PDF}} \put(122.0,115.0){DNS}
\put(122.0,107.0){M1} \put(122.0,98.0){M2} \put(122.0,89.0){M3} \end{overpic}
(b)

Figure 3: PDFs of (a) $Q/\tilde{Q}_{av}$ and (b) $R/\tilde{Q}_{av}^{3/2}$,
where $\tilde{Q}_{av}\equiv\langle\|\bm{\mathcal{A}}\|^{2}\rangle$.

###  4.2  Scalar Gradients

We now turn to consider the predictions from our new model for the scalar
gradients. In Figure 4, we plot the PDFs of $Q_{b}/Q_{b,av}$ (where
$Q_{b}\equiv\|\bm{b}\|^{2}$) which is proportional to the scalar dissipation
rate $\epsilon_{\phi}\equiv\kappa\|\bm{b}\|^{2}$, as well as the PDF of
$b_{1}$, the scalar gradient component in one of the (arbitrary, due to
isotropy) directions. Concerning the PDF of $Q_{b}/Q_{b,av}$, the results show
that while the model is in good qualitative agreement with the DNS data,
capturing the slowly decaying tail of the PDF, it significantly underpredicts
the values of the PDF. The model predictions for the PDF of $b_{1}$ are also
in significant error, underpredicting small to intermediate values of
$b_{1}/b_{1,rms}$, and significantly overpredicting large values of
$b_{1}/b_{1,rms}$, such that the overall shape of the PDF is not well captured
by the model.

An investigation into the cause of these significant underpredictions revealed
that the problem is due to the model generating extremely large values of
$\|\bm{\mathcal{B}}(t)\|^{2}$. An example of the time series of
$\|\bm{\mathcal{B}}(t)\|^{2}/\langle\|\bm{\mathcal{B}}(t)\|^{2}\rangle$
generated by the model at $\Rey_{\lambda}=100$ is shown in Figure 5, together
with the time series of
$\|\bm{\mathcal{S}}(t)\|^{2}/\langle\|\bm{\mathcal{S}}(t)\|^{2}\rangle$ for
comparison. Although the signal $\|\bm{\mathcal{S}}(t)\|^{2}$ exhibits
significant fluctuations about the mean, $\|\bm{\mathcal{B}}(t)\|^{2}$
exhibits infrequent but enormous fluctuations about the mean, which only get
stronger as $\Rey_{\lambda}$ is increased. Although
$\|\bm{\mathcal{B}}(t)\|^{2}$ would be expected to be more intermittent than
$\|\bm{\mathcal{S}}(t)\|^{2}$, one would not anticipate intermittent
fluctuations in $\|\bm{\mathcal{B}}(t)\|^{2}$ as large as these, nor are they
manifested in the DNS data, and therefore they seem to indicate an issue with
the model. Since the integral of the PDF of $Q_{b}$ over its sample-space is
one, then because the model vastly overpredicts the probability of extremely
large values of $\|\bm{\mathcal{B}}(t)\|^{2}$, it underpredicts the
probability of values in the sample-space range shown in figure 4.

\begin{overpic}[trim=14.22636pt 170.71652pt 0.0pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_3_PDF_all_Re_lambda_Fixed_alpha.pdf}
\put(80.0,-2.0){$Q_{b}/Q_{b,av}$} \put(-13.0,60.0){\rotatebox{90.0}{PDF}}
\put(122.0,115.0){DNS} \put(122.0,107.0){M1} \put(122.0,98.0){M2}
\put(122.0,89.0){M3} \end{overpic} (a)  \begin{overpic}[trim=14.22636pt
170.71652pt 0.0pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/b1_PDF_all_Re_lambda_withDNS_FIXED_alpha.pdf}
\put(80.0,-2.0){$b_{1}/b_{1,rms}$} \put(-13.0,60.0){\rotatebox{90.0}{PDF}}
\put(122.0,115.0){DNS} \put(122.0,107.0){M1} \put(122.0,98.0){M2}
\put(122.0,89.0){M3} \end{overpic} (b)

Figure 4: PDFs of (a) $Q_{b}/Q_{b,av}$, (b) $b_{1}/b_{1,rms}$. The results
from the model are obtained with the (uncorrected) model coefficient
$\alpha_{\mathcal{B}}$ specified by equation (39).

\begin{overpic}[trim=14.22636pt 170.71652pt 0.0pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/time_series_B2.pdf}
\put(80.0,2.0){$t/\tau_{\eta}$}
\put(-15.0,40.0){\rotatebox{90.0}{$\|\bm{\mathcal{B}}(t)\|^{2}/\langle\|\bm{\mathcal{B}}(t)\|^{2}\rangle$}}
\end{overpic} (a)  \begin{overpic}[trim=14.22636pt 170.71652pt 0.0pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/time_series_S2.pdf}
\put(80.0,2.0){$t/\tau_{\eta}$}
\put(-5.0,40.0){\rotatebox{90.0}{$\|\bm{\mathcal{S}}(t)\|^{2}/\langle\|\bm{\mathcal{S}}(t)\|^{2}\rangle$}}
\end{overpic} (b)

Figure 5: Time series of $\|\bm{\mathcal{B}}(t)\|^{2}$ and
$\|\bm{\mathcal{S}}(t)\|^{2}$, normalized by their mean values, generated from
the model with $\Rey_{\lambda}=100$ and the (uncorrected) coefficient
$\alpha_{\mathcal{B}}$ specified by equation (39).

That the model vastly overpredicts the probability of extremely large values
of $\|\bm{\mathcal{B}}(t)\|^{2}/\langle\|\bm{\mathcal{B}}(t)\|^{2}\rangle$
must be due to deficiencies in the closure for
$\kappa\langle\nabla^{2}\bm{B}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}$. In
particular, the values of
$\kappa\langle\nabla^{2}\bm{B}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}$
predicted by the closure approximation when $\|\bm{\mathcal{B}}\|$ is large
are too small to sufficiently counteract the scalar production term (although
apparently, they are large enough to prevent singularities in the model since
simulations of the model do not blow up). A relatively simple modification to
the model to address this deficiency is to modify its specification of the
coefficient $\alpha_{\mathcal{B}}$ in (39) so that it includes information on
the locally averaged scalar production, rather than simply the global mean
value. This is achieved by replacing (39) with

$\displaystyle\alpha_{\mathcal{B}}(t)$
$\displaystyle=\tau_{\eta}\int^{t}_{t-\tau_{1}}\bm{\mathcal{S}}(t^{\prime})\bm{:}\bm{\mathcal{B}}(t^{\prime})\bm{\mathcal{B}}(t^{\prime})\,dt^{\prime}\Bigg{/}\int^{t}_{t-\tau_{1}}\|\bm{\mathcal{B}}(t^{\prime})\|^{2}\,dt^{\prime},$
(46)

which can be computed when solving the model since it only depends on
$\bm{\mathcal{S}}$ and $\bm{\mathcal{B}}$ at previous times, which are known.
With this, the global average involved in (39) is replaced with a local time
average over the trajectory history of the particle. The time integral is
chosen to span $[t-\tau_{1},t]$ in view of the fact that
$\bm{\mathcal{S}}\bm{:}\bm{\mathcal{B}}\bm{\mathcal{B}}$ and
$\|\bm{\mathcal{B}}\|^{2}$ have timescales on the order of the integral
timescale, which in the ML-RDGF is specified by $\tau_{1}$.

The advantage of using (46) is that the coefficient $\delta_{\mathcal{B}}$ in
(42) will then depend upon the local scalar gradient dynamics, and in regions
where the production of $\|\bm{\mathcal{B}}\|^{2}$ is large,
$\delta_{\mathcal{B}}$ will also be large relative to its value in regions
where the production of $\|\bm{\mathcal{B}}\|^{2}$ is small. In other words,
using (46) introduces nonlinearity into the closure for
$\kappa\langle\nabla^{2}\bm{B}\rangle_{\bm{\mathcal{A}},\bm{\mathcal{B}}}$
with respect to its dependence on $\bm{\mathcal{B}}$, and this may help oppose
the extremely large fluctuations predicted by the original form of the model.
In practice, since (46) requires time-history information, the model is solved
for $t\leq\tau_{1}$ using (39), and then for $t>\tau_{1}$, (46) is used.

\begin{overpic}[trim=14.22636pt 170.71652pt 0.0pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_3_PDF_all_Re_lambda.pdf}
\put(80.0,-2.0){$Q_{b}/Q_{b,av}$} \put(-13.0,60.0){\rotatebox{90.0}{PDF}}
\put(122.0,115.0){DNS} \put(122.0,107.0){M1} \put(122.0,98.0){M2}
\put(122.0,89.0){M3} \end{overpic} (a)  \begin{overpic}[trim=14.22636pt
170.71652pt 0.0pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/b1_PDF_all_Re_lambda_withDNS.pdf}
\put(80.0,-2.0){$b_{1}/b_{1,rms}$} \put(-13.0,60.0){\rotatebox{90.0}{PDF}}
\put(122.0,115.0){DNS} \put(122.0,107.0){M1} \put(122.0,98.0){M2}
\put(122.0,89.0){M3} \end{overpic} (b)

Figure 6: PDFs of (a) $Q_{b}/Q_{b,av}$, (b) $b_{1}/b_{1,rms}$, based on using
(46) instead of (39) to specify $\alpha_{\mathcal{B}}$ in the model.

\begin{overpic}[trim=14.22636pt 170.71652pt 0.0pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_4_PDF_all_Re_lambda.pdf}
\put(73.0,-2.0){$R_{b}/\tilde{Q}_{av}^{1/2}Q_{b,av}$}
\put(-10.0,60.0){\rotatebox{90.0}{PDF}} \put(122.0,115.0){DNS}
\put(122.0,107.0){M1} \put(122.0,98.0){M2} \put(122.0,89.0){M3} \end{overpic}
(a)  \begin{overpic}[trim=14.22636pt 170.71652pt 0.0pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/CosineAngles_PDF_Re_lambda_100_withDNS.pdf}
\put(80.0,-2.0){$\bm{e}_{b}\bm{\cdot}\bm{e}_{i}$}
\put(-10.0,60.0){\rotatebox{90.0}{PDF}} \put(40.0,116.0){DNS, $i=1$}
\put(40.0,107.0){DNS, $i=2$} \put(40.0,98.0){DNS, $i=3$} \put(40.0,90.0){M1,
$i=1$} \put(40.0,81.0){M1, $i=2$} \put(40.0,72.0){M1, $i=3$} \end{overpic} (b)

Figure 7: PDFs of (a) $R_{b}/\tilde{Q}_{av}^{1/2}Q_{b,av}$, (b) the inner
product between the unit vector $\bm{e}_{b}\equiv\bm{b}/\|\bm{b}\|$ and the
eigenvectors $\bm{e}_{i}$ of $\bm{\mathcal{S}}$.

Figure 6 once again shows the PDFs of $Q_{b}$ and $b_{1}$, but this time using
(46) instead of (39) to specify $\alpha_{\mathcal{B}}$ in the model. Comparing
the results to those in figure 4, it can be seen that the new specification of
$\alpha_{\mathcal{B}}$ dramatically improves the predictions from the model,
being now in excellent agreement with the DNS data. The results also show the
impact of $\Rey_{\lambda}$ on $Q_{b}$ as predicted by the model, with the
probability of intermediate values of $Q_{b}/Q_{b,av}$ and $b_{1}/b_{1,rms}$
predicted to decrease as $\Rey_{\lambda}$ increases, while the probability of
large $Q_{b}/Q_{b,av}$ and $b_{1}/b_{1,rms}$ increases as $\Rey_{\lambda}$
increases. However, further tests of the model revealed that when
$\Rey_{\lambda}$ is increased much beyond $\Rey_{\lambda}=500$ the predictions
of the model become unrealistic, with extremely large values of the scalar
gradient occurring in the model, and the model can even blow up. Therefore,
although the use of (46) dramatically improves the performance of the model
over the range of $\Rey_{\lambda}$ considered, it is not sufficient to
guarantee that the model makes reasonable predictions for arbitrarily large
$\Rey_{\lambda}$. An investigation into the causes of the failure of the model
at high $\Rey_{\lambda}$ and possible remedies for this are left to future
work.

\begin{overpic}[trim=0.0pt 170.71652pt 28.45274pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_1_3_JPDF_DNS.pdf}
\put(75.0,-1.0){$Q_{b}/Q_{b,av}$}
\put(8.0,55.0){\rotatebox{90.0}{$Q/\tilde{Q}_{av}$}} \end{overpic} (a)
\begin{overpic}[trim=0.0pt 170.71652pt 28.45274pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_1_3_JPDF_Re_lambda_100.pdf}
\put(75.0,-1.0){$Q_{b}/Q_{b,av}$}
\put(8.0,55.0){\rotatebox{90.0}{$Q/\tilde{Q}_{av}$}} \end{overpic} (b)

\begin{overpic}[trim=0.0pt 170.71652pt 28.45274pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_1_3_JPDF_Re_lambda_300.pdf}
\put(75.0,-1.0){$Q_{b}/Q_{b,av}$}
\put(8.0,55.0){\rotatebox{90.0}{$Q/\tilde{Q}_{av}$}} \end{overpic} (c)
\begin{overpic}[trim=0.0pt 170.71652pt 28.45274pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_1_3_JPDF_Re_lambda_500.pdf}
\put(75.0,-1.0){$Q_{b}/Q_{b,av}$}
\put(8.0,55.0){\rotatebox{90.0}{$Q/\tilde{Q}_{av}$}} \end{overpic} (d)

Figure 8: Logarithm of joint-PDFs of $Q/\tilde{Q}_{av}$ and $Q_{b}/Q_{b,av}$
from (a) DNS, (b) M1, (c) M2, (d) M3. Colors indicate the values of the
logarithm of the PDF.

In figure 7(a) we show the PDF of the scalar gradient production
$R_{b}\equiv\bm{a:bb}$. The results show that the model predictions are in
good agreement with the DNS data, with some underpredictions for the largest
fluctuations. The model captures the strong negative skewness of the PDF that
is associated with the predominance of scalar gradient production over
destruction. The model also predicts the largest fluctuations in $R_{b}$
become more probable as $\Rey_{\lambda}$ is increased due to intermittency in
the flow. In figure 7(b) we show the PDF of the inner product between the unit
vector $\bm{e}_{b}\equiv\bm{b}/\|\bm{b}\|$ and the eigenvectors $\bm{e}_{i}$
(corresponding to the ordered eigenvalues) of $\bm{\mathcal{S}}$. The model
predicts these non-trivial alignments very well, capturing the strong
preferential alignment with the compressional eigendirection $\bm{e}_{3}$, and
misalignment with the intermediate eigendirection $\bm{e}_{2}$ and extensional
eigendirection $\bm{e}_{1}$. However, the model predicts a misalignment with
$\bm{e}_{1}$ that is a little too strong, and a misalignment with $\bm{e}_{2}$
that is a little too weak. Only the results for M1 are shown, as the results
from M2 and M3 are almost identical. The current model predictions for these
alignment PDFs are in much better agreement with the DNS data than those of
the model of Martín et al. (2005) which uses a much more simplistic closure
approximation for the scalar diffusion term.

###  4.3 Joint-PDFs of Velocity and Scalar Gradients

\begin{overpic}[trim=0.0pt 170.71652pt 28.45274pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_1_4_JPDF_DNS.pdf}
\put(70.0,-1.0){$R_{b}/\tilde{Q}_{av}^{1/2}Q_{b,av}$}
\put(8.0,55.0){\rotatebox{90.0}{$Q/\tilde{Q}_{av}$}} \end{overpic} (a)
\begin{overpic}[trim=0.0pt 170.71652pt 28.45274pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_1_4_JPDF_Re_lambda_100.pdf}
\put(70.0,-1.0){$R_{b}/\tilde{Q}_{av}^{1/2}Q_{b,av}$}
\put(8.0,55.0){\rotatebox{90.0}{$Q/\tilde{Q}_{av}$}} \end{overpic} (b)

\begin{overpic}[trim=0.0pt 170.71652pt 28.45274pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_1_4_JPDF_Re_lambda_300.pdf}
\put(70.0,-1.0){$R_{b}/\tilde{Q}_{av}^{1/2}Q_{b,av}$}
\put(8.0,55.0){\rotatebox{90.0}{$Q/\tilde{Q}_{av}$}} \end{overpic} (c)
\begin{overpic}[trim=0.0pt 170.71652pt 28.45274pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_1_4_JPDF_Re_lambda_500.pdf}
\put(70.0,-1.0){$R_{b}/\tilde{Q}_{av}^{1/2}Q_{b,av}$}
\put(8.0,55.0){\rotatebox{90.0}{$Q/\tilde{Q}_{av}$}} \end{overpic} (d)

Figure 9: Logarithm of joint-PDFs of $Q/\tilde{Q}_{av}$ and
$R_{b}/\tilde{Q}_{av}^{1/2}Q_{b,av}$ from (a) DNS, (b) M1, (c) M2, (d) M3.
Colors indicate the values of the logarithm of the PDF.

We now turn to consider the relationship between the velocity and scalar
gradients predicted by the model. In Figure 8 we show the joint-PDFs of the
velocity gradient invariant $Q/\tilde{Q}_{av}$ and scalar invariant
$Q_{b}/Q_{b,av}$ from the DNS and model. For M1, there is an excellent
qualitative agreement with the DNS data, with the model capturing the
elongation of the PDF (note however that the appearance of the elongation is
somewhat exaggerated due to the different axis ranges) along the horizontal
axis toward regions of large $Q_{b}/Q_{b,av}$, indicating that large
fluctuations in the scalar gradients are much more probable than they are for
the velocity gradients. The model also captures the exponential-like behaviour
of the isocontours of the PDF, whose shape indicates that large values of
$Q_{b}/Q_{b,av}$ tend to occur in regions where $Q/\tilde{Q}_{av}$ is small,
and vice-versa. The quantitative errors in M1 are mainly associated with the
variation of the joint-PDF along the $Q/\tilde{Q}_{av}$ axis, which can be
understood in terms of the ML-RGDF’s underprediction of large fluctuations of
$Q/\tilde{Q}_{av}$ at $\Rey_{\lambda}=100$, as already observed when
considering the PDF of $Q/\tilde{Q}_{av}$ in figure 3. Comparing the results
from M1, M2, and M3 shows that the model predicts that as $\Rey_{\lambda}$ is
increased the shape of the joint-PDF is preserved, but becomes stretched along
the two axes. Correspondingly the probability of regions with comparable
values of $Q/\tilde{Q}_{av}$ and $Q_{b}/Q_{b,av}$ is predicted to decrease as
$\Rey_{\lambda}$ increases.

\begin{overpic}[trim=0.0pt 170.71652pt 28.45274pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_3_4_JPDF_DNS.pdf}
\put(70.0,-1.0){$R_{b}/\tilde{Q}_{av}^{1/2}Q_{b,av}$}
\put(8.0,50.0){\rotatebox{90.0}{$Q_{b}/Q_{b,av}$}} \end{overpic} (a)
\begin{overpic}[trim=0.0pt 170.71652pt 28.45274pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_3_4_JPDF_Re_lambda_100.pdf}
\put(70.0,-1.0){$R_{b}/\tilde{Q}_{av}^{1/2}Q_{b,av}$}
\put(8.0,50.0){\rotatebox{90.0}{$Q_{b}/Q_{b,av}$}} \end{overpic} (b)

\begin{overpic}[trim=0.0pt 170.71652pt 28.45274pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_3_4_JPDF_Re_lambda_300.pdf}
\put(70.0,-1.0){$R_{b}/\tilde{Q}_{av}^{1/2}Q_{b,av}$}
\put(8.0,50.0){\rotatebox{90.0}{$Q_{b}/Q_{b,av}$}} \end{overpic} (c)
\begin{overpic}[trim=0.0pt 170.71652pt 28.45274pt
170.71652pt,scale={0.3},clip,tics=20]{Figures/x_3_4_JPDF_Re_lambda_500.pdf}
\put(70.0,-1.0){$R_{b}/\tilde{Q}_{av}^{1/2}Q_{b,av}$}
\put(8.0,50.0){\rotatebox{90.0}{$Q_{b}/Q_{b,av}$}} \end{overpic} (d)

Figure 10: Logarithm of joint-PDFs of $Q_{b}/Q_{b,av}$ and
$R_{b}/\tilde{Q}_{av}^{1/2}Q_{b,av}$ from (a) DNS, (b) M1, (c) M2, (d) M3.
Colors indicate the values of the logarithm of the PDF.

In Figure 9 we show the joint-PDFs of the velocity gradient invariant
$Q/\tilde{Q}_{av}$ and scalar production invariant
$R_{b}/\tilde{Q}_{av}^{1/2}Q_{b,av}$ from the DNS and model. This joint PDF
provides insights into how straining and vortical regions of the flow might
contribute differently to the scalar gradient production. The model accurately
reproduces the qualitative behaviour of the PDF seen in the DNS data,
including the elongation of the PDF along the
$R_{b}/\tilde{Q}_{av}^{1/2}Q_{b,av}<0$ direction, associated with the
predominance of scalar gradient production over destruction. The results also
show that events with strong scalar gradient production are more probable in
regions where $Q/\tilde{Q}_{av}<0$, despite the fact that $Q/\tilde{Q}_{av}$
itself has a positively skewed PDF. This is because the antisymmetric part of
the velocity gradient, and therefore the vorticity, does not directly
contribute to the invariant $R_{b}$, but only contributes indirectly through
its impact on the local alignment of $\bm{\mathcal{B}}(t)$ with
$\bm{\mathcal{S}}(t)$, whereas the strain-rate directly affects $R_{b}$. The
model underpredicts the probability of the largest fluctuations, which again
is principally due to the ML-RDGF underpredicting the intermittency of
$Q/\tilde{Q}_{av}$. Comparing the results from M1, M2, and M3 shows that the
model predicts that as $\Rey_{\lambda}$ is increased the shape of the joint-
PDF is largely preserved, except for being stretched along the axes due to the
increased intermittency of the flow. The model does seem to predict, however,
that the overall probability of the system being in the quadrant $Q<0,R_{b}>0$
reduces as $\Rey_{\lambda}$ increases. Future comparisons with DNS at higher
$\Rey_{\lambda}$ will be needed to assess the accuracy of this prediction. It
is possible that this is a defect in the model that is in some way related to
the failure of the model at higher $\Rey_{\lambda}$.

Finally, in Figure 10 we show the joint-PDFs of the scalar gradient invariant
$Q_{b}/Q_{b,av}$ and scalar production invariant
$R_{b}/\tilde{Q}_{av}^{1/2}Q_{b,av}$. The prediction of the model for M1 is in
very good agreement with the DNS, both qualitatively and quantitatively, with
only small deviations. As with the other joint-PDFs, comparing the results
from M1, M2, and M3 shows that the model predicts that as $\Rey_{\lambda}$ is
increased the shape of the joint-PDF of $Q_{b}/Q_{b,av}$ and
$R_{b}/\tilde{Q}_{av}^{1/2}Q_{b,av}$ is preserved, except for being stretched
along the axes due to the increased intermittency of the flow.

##  5 Conclusions

A Lagrangian model for passive scalar gradients in isotropic turbulence has
been developed, with the scalar gradient diffusion term closed using the RDGF
approach, which has recently been very successfully applied to close the
equation for the fluid velocity gradients along fluid particle trajectories.
This closure yields a diffusion term that is nonlinear in the velocity
gradients, but linear in the scalar gradients, and comparisons of the
statistics generated by the closed model with DNS data revealed large errors.
An investigation revealed that these large errors were due to the scalar
gradient model generating erroneously large fluctuations, possibly due to the
diffusion term being linear in the scalar gradient under the RDGF closure.
This defect was addressed by incorporating into the closure approximation
information regarding the scalar gradient production along the local
trajectory history of the particle. With this modification, the closed form of
the diffusion term is now a nonlinear functional of the scalar gradients, and
the resulting model is in very good agreement with the DNS data.

Since the ML-RDGF model of Johnson & Meneveau (2017) is used to specify the
velocity gradients in the scalar gradient equation, we begin by comparing its
predictions with DNS data. In agreement with the results of Johnson & Meneveau
(2017), the model very accurately predicts the longitudinal and transverse
components of the velocity gradients, and reproduces the key features of the
joint-PDF of $Q$ and $R$, the second and third invariants of the velocity
gradient tensor. However, a more quantitative test of the model predictions
for the PDFs of $Q$ and $R$ individually against DNS data revealed that at
least for $\Rey_{\lambda}=100$, the ML-RDGF significantly underpredicts the
probability of large positive values of $Q$ and large negative values of $R$,
which are primarily associated with regions of intense enstrophy and enstrophy
production, respectively. These inaccuracies could then impact the accuracy of
the scalar gradient model, and highlight the need for further improvements in
the ML-RDGF model (although the model may predict the statistics of $Q$ and
$R$ more accurately at higher $\Rey_{\lambda}$).

Comparisons between the scalar gradient model and DNS data for
$\Rey_{\lambda}=100$ showed very good agreement. In particular, the model
accurately predicts the squared magnitude of the scalar gradients (which are
proportional to the scalar dissipation rates), as well as the individual
components of the scalar gradients. The model also captures well the PDF of
the scalar production, including its strong negative skewess that is
associated with the predominance of scalar gradient production of destruction,
but slightly underpredicts the most extreme fluctuations of the scalar
gradient production and destruction. Next, the PDFs of the inner product
between the scalar gradient direction and the strain-rate eigendirections were
considered, which provide insights into the nontrivial statistical geometry of
the passive scalar and velocity gradient dynamics. The model is in excellent
agreement with the DNS data regarding the strong preferential alignment
between the scalar gradient and the compressional eigendirection. However, the
model predicts a misalignment with the extensional eigendirection that is a
little too strong, and a misalignment with the intermediate eigendirection
that is a little too weak.

The ability of the model to capture the statistical relationship between the
velocity and scalar dynamics was considered next, by considering various
joint-PDFs of the velocity and scalar gradient tensors. The results showed
excellent qualitative agreement with the DNS data, with some quantitative
errors that seem to be rooted in the ML-RDGF model underpredicting the
probability of extreme fluctuations in the $Q$ and $R$ invariants. The joint-
PDF of the squared magnitude of the scalar gradient with the scalar gradient
production term predicted by the model was in excellent qualitative as well as
quantitative agreement with the DNS.

The predictions of the model at $\Rey_{\lambda}$ greater than that of the DNS
were also considered, and the predictions are reasonable up to around
$\Rey_{\lambda}\approx 500$. However, beyond this, the model breaks down and
leads to extremely large scalar gradients that can even cause the numerical
simulations of the model to blow up. Therefore, while the modification to the
scalar gradient diffusion term that incorporates the scalar gradient
production along the local trajectory history of the particle leads to
excellent predictions from the model at lower $\Rey_{\lambda}$, it is not
sufficient to prevent the model from generating extremely large fluctuations
at high $\Rey_{\lambda}$ where intermittency in the velocity gradients can
lead to very large local scalar gradient production events. Therefore, as
anticipated earlier, developing an accurate model for scalar gradients in
turbulence is in some ways more complicated than that for velocity gradients,
because scalar gradient dynamics lack a mechanism similar to the pressure
Hessian that controls the growth of the velocity gradients. For scalar
gradients, the closure for the diffusion term is a delicate matter since this
term alone is dynamically responsible for preventing finite-time singularities
of the scalar gradients. A crucial point for future work is therefore to
understand in more detail how the diffusion term regulates the growth of the
scalar gradients, and developing a closure model that is sufficiently
sophisticated to capture this.

Another aspect to be explored is the influence of the Schmidt number $Sc$ on
the scalar gradients. While the model does capture a $Sc$ dependence in the
scalar gradient diffusion term, given that the model breaks down for $Sc=1$
when $\Rey_{\lambda}>500$, it is likely that the model will also break down
for $\Rey_{\lambda}<500$ when $Sc$ becomes sufficiently large, since both
regimes promote the intensification of the scalar gradients.

[Acknowledgements]This work used the Extreme Science and Engineering Discovery
Environment (XSEDE), which is supported by NSF grant number ACI-1548562 (Towns
et al., 2014). Specifically, the Expanse cluster was used under allocation
CTS170009.

[Funding]A.D. Bragg and X. Zhang acknowledge support through a National
Science Foundation (NSF) CAREER award # 2042346

[Declaration of interests]The authors report no conflict of interest.

## References

  * Batchelor (1959) Batchelor, George K 1959 Small-scale variation of convected quantities like temperature in turbulent fluid part 1. general discussion and the case of small conductivity. Journal of fluid mechanics 5 (1), 113–133.
  * Betchov (1956) Betchov, R 1956 An inequality concerning the production of vorticity in isotropic turbulence. Journal of Fluid Mechanics 1 (5), 497–504.
  * Bhamidipati et al. (2020) Bhamidipati, Neeraja, Souza, Andre N & Flierl, Glenn R 2020 Turbulent mixing of a passive scalar in the ocean mixed layer. Ocean Modelling 149, 101615.
  * Buaria et al. (2021a) Buaria, Dhawal, Clay, Matthew P, Sreenivasan, Katepalli R & Yeung, PK 2021a Small-scale isotropy and ramp-cliff structures in scalar turbulence. Physical Review Letters 126 (3), 034504.
  * Buaria et al. (2021b) Buaria, Dhawal, Clay, Matthew P, Sreenivasan, Katepalli R & Yeung, PK 2021b Turbulence is an ineffective mixer when schmidt numbers are large. Physical Review Letters 126 (7), 074501.
  * Canuto (1988) Canuto, Claudio 1988 Spectral methods and a maximum principle. Mathematics of Computation 51 (184), 615–629.
  * Carbone et al. (2019) Carbone, M., Bragg, A. D. & Iovieno, M. 2019 Multiscale fluid–particle thermal interaction in isotropic turbulence. Journal of Fluid Mechanics 881, 679–721.
  * Chertkov et al. (1999) Chertkov, Michael, Pumir, Alain & Shraiman, Boris I 1999 Lagrangian tetrad dynamics and the phenomenology of turbulence. Physics of fluids 11 (8), 2394–2410.
  * Chevillard & Meneveau (2006) Chevillard, Laurent & Meneveau, Charles 2006 Lagrangian dynamics and statistical geometric structure of turbulence. Physical review letters 97 (17), 174501.
  * Chor (2020) Chor, Tomas Luis Guimaraes 2020 Mixing and transport of materials in the Ocean Surface Boundary Layer. University of California, Los Angeles.
  * Corrsin (1951) Corrsin, Stanley 1951 On the spectrum of isotropic temperature fluctuations in an isotropic turbulence. Journal of Applied Physics 22 (4), 469–473.
  * Dimotakis (2005) Dimotakis, Paul E 2005 Turbulent mixing. Annual Review of Fluid Mechanics 37, 329–356.
  * Donzis et al. (2005) Donzis, DA, Sreenivasan, KR & Yeung, P Kc 2005 Scalar dissipation rate and dissipative anomaly in isotropic turbulence. Journal of Fluid Mechanics 532, 199–216.
  * Falkovich et al. (2001) Falkovich, G., Gawedzki, K. & Vergassola, M. 2001 Particles and fields in fluid turbulence. Rev. Mod. Phys. 73, 913–975.
  * Gonzalez (2009) Gonzalez, M 2009 Kinematic properties of passive scalar gradient predicted by a stochastic lagrangian model. Physics of Fluids 21 (5), 055104.
  * Hater et al. (2011) Hater, T, Homann, H & Grauer, R 2011 Lagrangian model for the evolution of turbulent magnetic and passive scalar fields. Physical Review E 83 (1), 017302.
  * Hill (1976) Hill, James C 1976 Homogeneous turbulent mixing with chemical reaction. Annual review of fluid Mechanics 8 (1), 135–161.
  * Honeycutt (1992) Honeycutt, Rebecca L. 1992 Stochastic Runge-Kutta algorithms. I. White noise. Phys. Rev. A 45, 600–603.
  * Johnson & Meneveau (2016) Johnson, Perry L & Meneveau, Charles 2016 A closure for lagrangian velocity gradient evolution in turbulence using recent-deformation mapping of initially gaussian fields. Journal of Fluid Mechanics 804, 387–419.
  * Johnson & Meneveau (2017) Johnson, Perry L & Meneveau, Charles 2017 Turbulence intermittency in a multiple-time-scale navier-stokes-based reduced model. Physical Review Fluids 2 (7), 072601.
  * Kloeden & Platen (2018) Kloeden, Peter E & Platen, Eckhard 2018 Numerical methods for stochastic differential equations. In Nonlinear dynamics and stochastic mechanics, pp. 437–461. CRC Press.
  * Kraichnan (1994) Kraichnan, Robert H 1994 Anomalous scaling of a randomly advected passive scalar. Physical Review Letters 72 (7), 1016\.
  * Leppin & Wilczek (2020) Leppin, Leonhard A. & Wilczek, Michael 2020 Capturing velocity gradients and particle rotation rates in turbulence. Phys. Rev. Lett. 125, 224501.
  * Lohse & Xia (2010) Lohse, Detlef & Xia, Ke-Qing 2010 Small-scale properties of turbulent rayleigh-bénard convection. Annual Review of Fluid Mechanics 42, 335–364.
  * Martın et al. (1998) Martın, Jesús, Dopazo, César & Valiño, Luis 1998 Dynamics of velocity gradient invariants in turbulence: restricted euler and linear diffusion models. Physics of Fluids 10 (8), 2012–2025.
  * Martín et al. (2005) Martín, Jesús, Dopazo, César & Valiño, Luis 2005 Joint statistics of the scalar gradient and the velocity gradient in turbulence using linear diffusion models. Physics of Fluids 17 (2), 028101.
  * Mazzitelli & Lanotte (2012) Mazzitelli, Irene & Lanotte, Alessandra S 2012 Active and passive scalar intermittent statistics in turbulent atmospheric convection. Physica D: Nonlinear Phenomena 241 (3), 251–259.
  * Meneveau (2011) Meneveau, Charles 2011 Lagrangian dynamics and models of the velocity gradient tensor in turbulent flows. Annual Review of Fluid Mechanics 43 (1), 219–245.
  * Mestayer et al. (1976) Mestayer, Patrice G, Gibson, Carl H, Coantic, Michel F & Patel, Anandkumar S 1976 Local anisotropy in heated and cooled turbulent boundary layers. The Physics of Fluids 19 (9), 1279–1287.
  * Nironi et al. (2015) Nironi, Chiara, Salizzoni, Pietro, Marro, Massimo, Mejean, Patrick, Grosjean, Nathalie & Soulhac, Lionel 2015 Dispersion of a passive scalar fluctuating plume in a turbulent boundary layer. part i: Velocity and concentration measurements. Boundary-layer meteorology 156 (3), 415–446.
  * Ottino & Ottino (1989) Ottino, Julio M & Ottino, JM 1989 The kinematics of mixing: stretching, chaos, and transport, , vol. 3. Cambridge university press.
  * Pekurovsky (2012) Pekurovsky, Dmitry 2012 P3dfft: A framework for parallel computations of fourier transforms in three dimensions. SIAM Journal on Scientific Computing 34 (4), C192–C209.
  * Pumir (1994) Pumir, Alain 1994 A numerical study of the mixing of a passive scalar in three dimensions in the presence of a mean gradient. Physics of Fluids 6 (6), 2118–2132.
  * Shete et al. (2022) Shete, Kedar Prashant, Boucher, David J., Riley, James J. & de Bruyn Kops, Stephen M. 2022 Effect of viscous-convective subrange on passive scalar statistics at high reynolds number. Phys. Rev. Fluids 7, 024601.
  * Shraiman & Siggia (2000) Shraiman, Boris I & Siggia, Eric D 2000 Scalar turbulence. Nature 405 (6787), 639–646.
  * Sreenivasan & Antonia (1977) Sreenivasan, KR & Antonia, RA 1977 Skewness of temperature derivatives in turbulent shear flows. The Physics of Fluids 20 (12), 1986–1988.
  * Sreenivasan (1991) Sreenivasan, Katepalli R 1991 On local isotropy of passive scalars in turbulent shear flows. Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 434 (1890), 165–182.
  * Tong & Warhaft (1994) Tong, Chenning & Warhaft, Z 1994 On passive scalar derivative statistics in grid turbulence. Physics of Fluids 6 (6), 2165–2176.
  * Towns et al. (2014) Towns, J., Cockerill, T., Dahan, M., Foster, I., Gaither, K., Grimshaw, A., Hazlewood, V., Lathrop, S., Lifka, D., Peterson, G. D., Roskies, R., Scott, J. R. & Wilkins-Diehr, N. 2014 Xsede: Accelerating scientific discovery. Computing in Science & Engineering 16 (5), 62–74.
  * Tsinober (2001) Tsinober, Arkady 2001 An informal introduction to turbulence. Kluwer Academic Publishers.
  * Vieillefosse (1982) Vieillefosse, P 1982 Local interaction between vorticity and shear in a perfect incompressible fluid. Journal de Physique 43 (6), 837–842.
  * Warhaft (2000) Warhaft, Zellman 2000 Passive scalars in turbulent flows. Annual Review of Fluid Mechanics 32 (1), 203–240.
  * Watanabe & Gotoh (2004) Watanabe, Takeshi & Gotoh, Toshiyuki 2004 Statistics of a passive scalar in homogeneous turbulence. New Journal of Physics 6 (1), 40.
  * Wilczek & Meneveau (2014) Wilczek, Michael & Meneveau, Charles 2014 Pressure hessian and viscous contributions to velocity gradient statistics based on gaussian random fields. Journal of Fluid Mechanics 756, 191–225.
  * Yeung et al. (2018) Yeung, P. K., Sreenivasan, K. R. & Pope, S. B. 2018 Effects of finite spatial and temporal resolution in direct numerical simulations of incompressible isotropic turbulence. Phys. Rev. Fluids 3, 064603.
  * Zhang et al. (2022) Zhang, Xiaolong, Dhariwal, Rohit, Portwood, Gavin, de Bruyn Kops, Stephen M & Bragg, Andrew D 2022 Analysis of scale-dependent kinetic and potential energy in sheared, stably stratified turbulence. Journal of Fluid Mechanics 946.