Nonlocality of Mean Scalar Transport in Two-Dimensional Rayleigh-Taylor Instability Using the Macroscopic Forcing Method

The importance of nonlocality of mean scalar transport in 2D Rayleigh-Taylor Instability (RTI) is investigated. The Macroscopic Forcing Method (MFM) is utilized to measure spatio-temporal moments of the eddy diﬀusivity kernel representing passive scalar transport in the ensemble averaged ﬁelds. Presented in this work are several studies assessing the importance of the higher-order moments of the eddy diﬀusivity, which contain information about nonlocality, in models for RTI. First, it is demonstrated through a comparison of leading-order models that a purely local eddy diﬀusivity is insuﬃcient in capturing the mean ﬁeld evolution of the mass fraction in RTI. Therefore, higher-order moments of the eddy diﬀusivity operator are not negligible. Models are then constructed by utilizing the measured higher-order moments. It is demonstrated that an explicit model based on the Kramers-Moyal expansion of the eddy diﬀusivity kernel is insuﬃcient. An implicit model construction that matches the measured moments is shown to oﬀer improvements relative to the local model in a converging fashion.


Introduction
Rayleigh-Taylor Instability (RTI) is a phenomenon that occurs when a heavy fluid is accelerated into a light fluid.Specifically, RTI occurs when the following are present: 1) a density gradient, 2) an acceleration (associated with the body force) in the direction opposite that of the density gradient, and 3) a perturbation at the interface of the two fluids.RTI is present in many scientific and engineering applications such as supernovae (Gull 1975) and inertial confinement fusion (ICF) (Zhou 2017;Lindl 1995).In the case of ICF, RTI occurs when a perturbation forms between the outer heavy ablator and the inner light deuterium gas, which causes premature mixing in the target, thereby greatly reducing the efficiency of the process.Thus, RTI is of great interest to scientists and engineers, especially in the context of ICF.
During a typical ICF experiment design process, a Reynolds-Averaged Navier-Stokes (RANS) approach is often utilized to model the role of hydrodynamic instabilities such as RTI.This is despite the fact that RTI can be more accurately predicted using high-fidelity methods like direct numerical simulations (DNS) (Youngs 1994;Cook & Dimotakis 2001;Cook & Zhou 2002;Cabot & Cook 2006;Mueschke & Schilling 2009) and large eddy simulations (LES) (Darlington et al. 2002;Cook et al. 2004;Cabot 2006).Motivation for development of RANS models for various engineering applications like ICF can be understood by considering the computational cost of each method.DNS requires resolution of the smallest turbulent scales, and LES the energy-containing scales, which are still much smaller than the macroscopic physics (i.e., averaged fields) of engineering interest.On the other hand, by design, RANS must only resolve macroscopic scales, thereby requiring much lower computational cost.Thus, RANS models are commonly used in engineering practice, especially in design optimization, where hundreds of thousands of simulations are often performed.Such is especially the case in designing targets for ICF experiments (Casey et al. 2014;Khan et al. 2016).Due to the utility of RANS in such applications, the need for predictive RANS models remains salient.
Models of varying complexities have been applied to the RTI problem.Among the most commonly-used types are two-equation models.One such model is the ubiquitously-used - model (Launder & Spalding 1974).Particularly, Gauthier & Bonnet (1990) introduced algebraic relations for some closures to satisfy realizability constraints for the model to be valid under the strong gradients of RTI.Another popular two-equation model is the - model; a version was introduced by Dimonte & Tipton (2006) for RTI.One appeal of the - model is its inclusion of a transport equation for turbulence lengthscale  (in place of the transport equation for  in -) that can be related to the initial interface perturbation.The self-similarity of turbulent RTI is leveraged to set the model coefficients.
These two-equation models rely on the gradient diffusion approximation for the turbulent mass flux closure.The gradient diffusion approximation rests on the assumption that turbulence transports quantities in a manner similar to Fickian diffusion.Importantly, this approximation implies purely local dependence of the mean turbulent flux on the mean gradient, ignoring history effects and gradients at nearby points in space.However, this approximation may not be valid for mean scalar transport.Specifically, the turbulent mass flux contains features that the gradient diffusion approximation cannot capture (Morgan & Greenough 2015;Denissen et al. 2014), so a local coefficient may not be enough to scale the mean gradient to model turbulent mass flux.
Nonlocality in RTI has been studied in experiments and simulations.Clark et al. (1997) analyzed data from turbulent RTI experiments and compared the pressure-strain correlation and pressure production due to turbulent mass flux, suggesting spatial nonlocality of pressure effects.DNS studies by Ristorcelli & Clark (2004) and experiments by Mueschke et al. (2006) have also examined nonlocality of RTI in the context of two-point correlations.Thus, the nonlocal nature of RTI is well-known, and work has been done to capture this nonlocality in models.For example, two-point closures to account for nonlocality in RTI have been developed by several authors for RANS (Clark & Spitz 1995;Steinkamp et al. 1999b,a;Pal et al. 2018;Kurien & Pal 2022) and LES (Parish & Duraisamy 2017).While these works attempt to address the effects of nonlocality in RTI, they do so without directly studying the form of the nonlocal operator.
Several authors have studied ways to directly measure the nonlocal eddy diffusivity in other canonical flows.One such approach involves application of the Green's function.The Green's function approach starts from analytical derivations of relations between turbulent fluxes and mean gradients, which was done by Kraichnan (1987).Hamba (1995) then introduced a reformulation of these relations appropriate for numerical computation of nonlocal eddy diffusivities, which has been applied to study channel flow (Hamba 2004) and, most recently, homogeneous isotropic turbulence (HIT) (Hamba 2022).
A different approach to determining nonlocal eddy diffusivities is the Macroscopic Forcing Method (MFM) by Mani & Park (2021).In contrast to the Green's function approach, MFM is derived by considering arbitrary forcing added directly to the transport equations with its formulation rooted in linear algebra.Additionally, MFM offers extensions to the Green's function approach by utilization of forcing functions that are not of the form of a Dirac delta.
Harmonic forcing has been utilized to derive analytical fits to nonlocal operators in Fourier space (Shirian & Mani 2022).Additionally, forcing polynomial mean fields using inverse MFM offers a computationally economical path for determination of spatio-temporal moments of the eddy diffusivity operator in conjunction with the Kramers-Moyal expansion as opposed to computation of the moments from a full MFM analysis through post-processing (Mani & Park 2021).Previous works using MFM have revealed turbulence operators for a variety of flows.Shirian & Mani (2022) and Shirian (2022) measured nonlocal operators in space and time in HIT.Though the spatial nonlocal operator was measured in HIT, it was applied to a turbulent round jet and was shown to match experiments more closely than the purely local Prandtl mixing-length model.MFM has also been applied to turbulent wall-bounded flows, including channel flow (Park & Mani 2023b) and separated boundary layers (Park et al. 2022;Park & Mani 2023a), to measure the anisotropic but local eddy diffusivity.In those flows, incorporation of the MFM-measured anisotropic eddy diffusivity improved RANS model predictions significantly, and remaining model errors were attributed to missing nonlocal effects.
It is with a motivation towards RANS model improvement that the present work seeks to understand nonlocality of closure operators governing turbulent scalar flux transport in RTI using MFM.Note that it is not intended for MFM to supplant current RANS models.Instead, MFM is an analysis tool that can be used to assess models and discover the necessary characteristics for accurate models.Here, MFM allows for direct measurement of nonlocal closure operators, which has not yet been done in RTI.This new knowledge of nonlocality of the mean scalar transport closure operator in RTI will aid in the development of improved RANS models used for studying ICF.
It is important to note that this work presents MFM measurements for a simplified RTI problem: the flow is two-dimensional, incompressible, and low-Atwood number, and only passive scalar mixing is considered.Since the eddy diffusivity is not universal, the MFM measurements of its moments presented here cannot be directly extended to more complex RTI.However, valuable insight into trends in the eddy diffusivity for mean scalar transport in RTI can be gained in this work.This follows the common process for developing turbulence models, where models are first designed for simpler flows then tested on and adjusted for more complex flows.In this work, MFM is performed on a simplified RTI problem to give a preliminary look into the eddy diffusivity of Rayleigh-Taylor-type flows, but future work will involve extensions to more complex flow characteristics that are closer to the practical flow observed in ICF capsules.The intent of this work is to present MFM as a tool for determining characteristics of the eddy diffusivity of a flow (i.e., its nonlocality and the importance of its higher order moments) that a model should satisfy in order to accurately predict mean scalar transport.The current work will inform future studies with additional complexities, including three-dimensionality, finite Atwood number, compressibility, and coupling with momentum.This work is organized as follows.First, an overview of RTI is covered briefly in §2.Next, §3 gives an overview of the mathematical methods used in this work, including: 1) the generalized eddy diffusivity and its approximation via a Karmers-Moyal expansion; 2) MFM and its application for finding the eddy diffusivity moments; 3) self-similarity analysis.Simulation details, including the governing equations and the computational approach, are given in §4.Finally, results of several studies on the importance of higher-order eddy diffusivity moments as well as assessments of suggested operator forms incorporating nonlocality of the eddy diffusivity for mean scalar transport in RTI are presented in §5.The results show that nonlocality of the eddy diffusivity is important in mean scalar transport of the RTI problem studied here, and RANS models incorporating this nonlocality result in more accurate predictions than leading-order models.

Brief overview of RTI
RTI is characterized by spikes (heavy fluid moving into light fluid) and bubbles (light fluid into heavy fluid).The mixing widths of these spikes and bubbles are denoted as ℎ  and ℎ  , respectively, and the mixing half-width is defined as ℎ = 1 2 (ℎ  + ℎ  ).The behaviors of these quantities in RTI are dependent on the Atwood number, defined as (2.1) Here,   and   are the densities of the heavy and light fluids, respectively.In the limit of low-Atwood number and late time, the mixing layer width is expected to reach a self-similar state of growth that scales quadratically with time: where  is the mixing width growth rate.The mixing width growth rate can also be viewed as the net mass flux through the midplane (Cook et al. 2004).In this case,  can also be written as where ℎ is the time derivative of ℎ.In the limit of self-similarity, these two definitions of  are expected to converge to the same value.
In a simulation, ℎ can be measured as where   is the mass fraction of the heavy fluid (therefore,   = 1 −   is the mass fraction of the light fluid), and ⟨ * ⟩ denotes averaging over realizations and homogeneous direction .An alternative definition used in works such as Cabot & Cook (2006) and Morgan et al. (2017) is This definition is particularly useful, since it allows ℎ to be determined solely based on the RANS field.That is, there is no closure problem in determining ℎ with this definition.Thus, this is the ℎ reported in this work.From these two definitions, a mixedness parameter  can be defined, which can be interpreted as the ratio of mixed to entrained fluid (Youngs 1994;Morgan et al. 2017): (2.6) In the limit of self-similarity,  is expected to approach a steady-state value.
A metric for turbulent transition is the Taylor Reynolds number: where  = 1 2 ⟨ ′   ′  ⟩ is the turbulence kinetic energy, and  is the effective Taylor microscale, approximated by (2.8) Here, the turbulent lengthscale  can be approximated as 1 5 the mixing layer width (Morgan et al. 2017).The large-scale Reynolds number can also be examined (Cabot & Cook 2006): where ℎ 99 is the mixing width based on 1%-99% mass fraction.Dimotakis (2000) determined that the criterion for turbulent transition is when   > 100 or   > 10, 000.

Model problem
In this work, a two-dimensional (2D), nonreacting flow with two species-a heavy fluid over a light fluid-is considered, with gravity pointing in the negative -direction.It must be noted that the behavior of 2D RTI is significantly different from three-dimensional (3D) RTI, the latter of which is more relevant to problems of engineering interest.It is well known that while 2D RTI is unsteady and chaotic, it is not strictly turbulent, since turbulence is a characteristic of 3D flows.In addition, 2D RTI has a faster late-time growth rate, develops larger structures, and is ultimately less well-mixed.These differences have been studied in RTI by Cabot (2006) and Young et al. (2001) and in Richtmeyer-Meshkov instability by Olson & Greenough (2014).
For this study, 2D RTI is chosen as the model problem instead of 3D RTI, since it is a good simplified setting for understanding nonlocality in RTI through the lens of MFM.Specifically, 2D RTI simulations are much less computationally expensive than those of 3D RTI, and MFM requires many simulations to attain statistical convergence.Thus, 2D RTI remains the focus of this work, with the hope that the understanding of nonlocality in this flow could be extended to nonlocality in 3D RTI.
In this 2D problem,  is the homogenous direction.In addition, there is no surface tension, the Atwood and Mach () numbers are finite but small, and the Peclet () number is finite but large.

Generalized eddy diffusivity and higher-order moments
In this work, the effect of nonlocality on mean scalar transport is of interest, so analysis begins with the scalar transport equation under the assumption of incompressibility: where u is the velocity vector and   is the molecular diffusivity of the heavy fluid.
After Reynolds decomposition and averaging, this becomes In this work, large  (the ratio of advective transport rate to diffusive transport rate) and small  are assumed.The former assumption means molecular diffusion is negligible, and the latter yields ⟨  ⟩ = 0, allowing the advective term to drop.Equation 3.2 becomes The term ⟨ ′  ′  ⟩ is the turbulent scalar flux, and this is the unclosed term that needs to be modeled.As mentioned previously, one reason the gradient diffusion approximation used to model this term is inaccurate is that it assumes locality of the eddy diffusivity.This assumption can be removed by instead considering a generalized eddy diffusivity that is nonlocal in space and time, as demonstrated by Romanof (1985) and Kraichnan (1987).For 2D RTI, such a model reduces to Here,  is the spatial coordinate in averaged space and  is the time at which the turbulent scalar flux is measured,  ′ is all points in averaged space, and  ′ is all points in time.This definition is exact for passive scalar transport, including in the case studied in this work.This nonlocal eddy diffusivity can also be viewed as a two-point correlation.This was first described by Taylor (1922) in homogeneous turbulence.Through Lagrangian statistical analysis, Taylor derived the following relation between diffusivity and velocity correlations: (3.5) Work by Shende et al. (2023) has shown that MFM recovers this Lagrangian formulation for eddy diffusivity in homogeneous flows.It should be noted that the above definition is not valid for inhomogeneous RTI (again, the exact definition of eddy diffusivity for the studied flow is the one in equation 3.4), but the intent here is to provide another interpretation of MFM more aligned with the well-understood two-point correlations.
The eddy diffusivity kernel can be approximated by Taylor-series-expanding the scalar gradient locally about  and , which results in the following Kramers-Moyal-like expansion for the turbulent scalar flux as done by Kraichnan (1987), Hamba (1995), andHamba (2004): Here,   are the eddy diffusivity moments; the first index, , denotes order in space, while the second, , denotes order in time.This is the form presented in Mani & Park (2021) and Liu et al. (2023).
When the eddy diffusivity kernel is purely local,  (,  ′ , ,  ′ ) =  00 ( −  ′ )( −  ′ ). (3.11) In this case,  00 is the only surviving moment, while all higher-order moments in space and time are zero.Any non-zero higher-order moment therefore characterizes the nonlocality of the eddy diffusivity kernel.Thus, this expansion implies explicitly a model form for the turbulent scalar flux that incorporates nonlocality of the eddy diffusivity.Truncating the expansion provides an approximation of ⟨ ′  ′  ⟩ but with the caveat that the expansion may not converge.This will be discussed in more detail later in §5.3.1.
Each   provides more information about the eddy diffusivity kernel with increasing order.For example,  00 represents the volume of the kernel in space-time.The coefficient corresponding to one higher-order in space,  10 , provides information about the centroid of the kernel in space. 20 contains information about the moment of inertia of the kernel in space,  01 contains information about the centroid of the kernel in time, and so on.

The Macroscopic Forcing Method
MFM is a method for numerically determining closure operators in turbulent flows (Mani & Park 2021).Much like a rheometer measures the molecular viscosity of a fluid by imposing a shear force on the flow, MFM forces the transport equation in a turbulent flow and extracts the closure operator from its response.Unlike the molecular viscosity, which is a material property, the turbulent closure operator is a property of the flow, so MFM measurements of one flow cannot be generalized for all flows; the MFM-measured closure for one flow cannot be applied exactly as it is to a different flow.However, MFM measurements of one flow can reveal characteristics of the turbulent closure that are expected be true for a family of similar flows.
Specifically, MFM can be used to determine the RANS closure operator, as shown in the pipeline diagram in figure 1.In MFM, two simulations are run at once: the donor and the receiver simulation.In this work, the donor simulation numerically solves the multicomponent Navier-Stokes equations in equations 4.1 -4.4.The receiver simulation "receives"   from the donor simulation and uses it to solve the scalar transport equation with a forcing : (3.12) Ultimately, forcings on the receiver simulation effect a response from the flow, and measuring this response allows for determination of the eddy diffusivity.In particular, these forcings are macroscopic.Here, macroscopic quantities are defined as fields that are unchanged by Reynoldsaveraging.Mathematically, the macroscopic forcing is such that  = .This macroscopic nature is crucial to the method, since it does not disturb the underlying mixing process, which allows for measurement of the closure operator without changing it.For details, see Mani & Park (2021).
In actuality, the Inverse Macroscopic Forcing Method (IMFM) is used to determine eddy diffusivity moments.That is, instead of the forcings being chosen, certain mean mass fraction fields are chosen.Numerically, mean mass fractions are enforced in each realization, so the averages (denoted by * ) described here are in , the homogeneous direction in space.The forcing needed to maintain the chosen   is determined implicitly along the process and is not directly used in the analysis.
As an illustration, the measurement of  00 can be considered.According to equation 3.6, choosing   =  (for  between −1/2 and 1/2) results in    = 1, and all other higher-order derivatives are zero.Thus, choosing this   in each realization results in the realization-and spatially-averaged measurement −⟨ ′  ′  ⟩ =  00 .Measurement of higher-order moments involves similar choices of   but requires information from lower order moments.For example, measuring  10 involves choosing   =  2 , which results in −⟨ ′  ′  ⟩ =  00 +  10 .Here,  00 comes from the simulation using   = .Thus,  10 is computed by subtracting  00 from the ⟨ ′  ′  ⟩ measurement from the simulation using Specifically, the following desired mean mass fractions are used for each moment for  between −1/2 and 1/2: From these   , the needed forcing in each timestep is numerically determined: where the superscript  denotes the timestep number,   desired is the mean mass fraction desired as outlined in equations 3.13 -3.16, and Δ is the timestep size.This MFM forcing bears some resemblance to other forcings used in the literature, such as interaction by exchange with the mean (IEM) (Sawford 2004;Pope 2001).One main difference between forcings in such methods and MFM is that the purpose of the latter is to drive the flow to a specified mean gradient, which allows for measurement-not enforcement -of the eddy diffusivity moments.In other words, in MFM for scalar transport, the input is a mean scalar gradient, and the output is the eddy diffusivity moment; in IEM and similar methods, the input is a desired moment (e.g., in IEM, the input moment is ⟨ 2 ⟩) and the output is a mixing model.In addition, methods such as IEM use microscopic forcings, while MFM uses macroscopic forcings, which is a distinguishing characteristic of the latter method.
To determine  00 ,  10 ,  01 , and  20 , four separate simulations are needed.For each of these simulations, the moments can be calculated using measurements of the turbulent scalar flux as follows: where   denotes the −⟨ ′  ′  ⟩ measured from the receiver simulation using the forcing corresponding the moment   .

Self-similarity analysis
We perform our analysis in the self-similar regime.First, we define a self-similar coordinate: so that ⟨  ⟩ is only a function of .Note that  requires a definition of ℎ().From the previous discussion on the self-similarity of RTI, an appropriate definition is ℎ() =  2 .Through self-similar analysis of equation 3.6, the eddy diffusivity moments and turbulent scalar flux can be normalized.Details of this process can be found in the Appendix.

Algebraic fit to mixing width
Recall that ℎ() =  2 is used in the self-similarity analysis.This is valid only for late time, so the subsequent analyses in this work are all done in this self-similar timeframe.Usually,  can be determined from ℎ( ) 2 , where ℎ() is computed from the simulation via equation 2.5.However, due to the convergence and statistical errors as well as the existence of a virtual time origin,  2 is not a good representation of ℎ() measured in the DNS.Instead, a fitting coefficient  * and virtual time origin  * are determined to make a shifted quadratic fit to ℎ() from the simulation: (3.23) With this fit, the normalizations of the turbulent scalar flux and moments become (3.28)For exact self-similarity, plots of the measured   against  must be independent of time.This expectation sets a criterion to assess the extent to which ideal self-similarity is achieved.Plots and assessment of the self-similar collapse of the measurements presented in this work are in the Appendix.

Governing equations
The governing equations solved in this work are the compressible multicomponent Navier-Stokes equations, which involve equations for continuity, diffusion of mass fraction   of species  (characterized by its binary molecular diffusivity   ), momentum transport, and transport of specific internal energy : Here,   is the material derivative   +      ,  is density,  is velocity,  is pressure, and  is gravitational acceleration, active in the − direction.The viscous stress tensor    and heat flux vector   are respectively defined as Here,  is the dynamic viscosity,  is the thermal conductivity,  is temperature, and ℎ  is the specific enthalpy of species .
Component pressures and temperatures of each species are determined using ideal gas equations of state.Under the assumption of pressure and temperature equilibrium, an iterative process is performed to determine volume fractions   that allow for computation of partial densities and energies.More details on the hydrodynamics equations and computation of component quantities can be found in Morgan et al. (2018).
Finally, total pressure is determined as the weighted sum of component pressures: In general, in these compressible equations,   are not passive scalars.However, the component equations of state are scaled so that a consistent hydrostatic pressure gradient is maintained across the mixing layer.Thus, in this work,   are effectively passive.

Computational approach
Simulations for 2D RTI are run using the Ares code, a hydrodynamics solver developed at Lawrence Livermore National Laboratory (LLNL) (Morgan & Greenough 2015;Bender et al. 2021).Ares employs an arbitrary Lagrangian-Eulerian (ALE) method based on the one by Sharp & Barton (1981), in which the governing equations (equations 4.1 to 4.4) are solved in a Lagrangian frame and then remapped to an Eulerian mesh through a second-order scheme.The spatial discretization is a second-order non-dissipative finite element method, and time advancement is a second-order explicit predictor-corrector scheme.
The Reynolds number (more specifically, the kinematic viscosity ) is set through a numerical Grashof number, such that Here, Δ is the grid spacing; in the simulations, a uniform mesh is used, and Δ = Δ = Δ.To ensure that the unsteady structures are properly resolved and for the simulation to appropriately be considered a DNS,  should be kept small.A  that is too large results in a simulation with dissipation dominated by numerics rather than the physics.Morgan & Black (2020) found that past  ≈ 12 in the Ares code, numerical diffusivity dominates molecular diffusivity.For our simulations, we use  = 0.05 and  = 1, the latter of which is in line with the DNS by Cabot & Cook (2006).These choices give a  of 10 −9  2 /.The Schmidt number , defined as /  , is set to unity, so   = 10 −9  2 /.The Mach number,  =   , where  is the speed of sound, characterizes compressibility effects of the flow. is set by the specific heat ratio , which is 5/3 in the simulations im this work.The maximum  is measured at the last timestep to be approximately 0.03, which is ascertained to be small enough to assume incompressibility.The Peclet number  characterizes the advection versus diffusion rate and is defined as .Here, a   and a   are reported, which use a large-scale   and the Taylor Reynolds number   , respectively.In the presented simulations,  = 1.The two  are computed in post-processing:   is approximately 8, 000, and   is approximately 54.Both are below the criterion established by Dimotakis (2000), suggesting that the simulated flow is transitional or pre-transitional.
The number of cells in each simulation is 2049 × 2049.The width   of the domain is 1, and the height   is 1.The boundary conditions are periodic in  and no slip and no penetration in .
Initially, the velocity field is zero, temperature is 293.15K, and pressure is 1 atm.A tophat perturbation based on the ones used by Morgan & Greenough (2015) and Morgan (2022) is imposed on the density field at the interface of the heavy and light fluids: where  1, and  2, are phase shift vectors randomly taken from a uniform distribution, and   is the length in  of the domain.Here, the minimum and maximum wavenumbers are set to  min = 8 and  max = 256, respectively.The stop condition of the simulations is when ℎ is approximately half the domain size in .This corresponds to the nondimensional simulation time  of 30.84.  is defined as   0 , where and ℎ 0 is the dominant lengthscale determined by the peak of the initial perturbation spectrum.
Before the MFM analysis was conducted, the results of the donor simulations were examined.In figure 2a, mixedness is observed to reach a value of around 0.6 but appears to not have converged yet.Figure 2b shows the two definitions of  over time.The first definition,  = ℎ/ 2 , reaches a value of about 0.05 by the end of the simulation, but it does not appear to be converged.The second definition,  = ℎ 2 /4ℎ, is oscillatory, due to the sensitivity of the time derivative to noise, and it appears to fluctuate about a value of approximately 0.04.It is acknowledged that this behavior indicates that the RTI simulated here is only weakly turbulent.However, it is observed that the flow is still self-similar at late times.The contour plot of ⟨  ⟩ in figure 3 exhibits parallel contour lines after  ≈ 17, indicating self-similarity at those times.It is also shown in  the Appendix that the mean concentration and normalized turbulent scalar flux profiles exhibit self-similar collapse after  ≈ 17, so the presented self-similar analysis is valid.
Figure 4 shows a plot of the algebraic fit for ℎ, described in equation 3.23.For the simulations presented here,  * is 0.0046 and  * is −1.6 × 10 3 .The plot shows a strong quadratic dependence of ℎ on  at late time, as ℎ fit matches DNS at  ⪆ 17, validating the self-similar ansatz of ℎ ∼  2 .
To further ensure the simulations are working as desired, the flow fields of the donor and receiver simulations can be examined qualitatively.The   contours at the last timesteps of each simulation are shown in figure 5.The receiver simulation shown is the one used to compute  00 (where ⟨  ⟩ = ).Self-similar RTI turbulent mixing is observed at this timestep, where the characteristic heavy spikes are sinking into the lighter fluid and the light bubbles rise into the heavier fluid.Both simulations have the same velocity fields, since the receiver simulation "receives" the velocity field from the donor simulation.In contrast with the donor simulation, which has a stark black-and-white difference between the heavy and light fluids, there is a grey gradient of density in the receiver simulation due to the imposed mean scalar gradient.The fluctuations of   in each simualtion are also compared in figure 6.The  ′  contours are not identical but are qualitatively very similar.In both simulations,  ′  is constrained to the mixing layer.Based on these observations, it is concluded that the simulations are visually working as intended.

Eddy diffusivity moments
Figure 7 shows normalized MFM measurements of the eddy diffusivity moments  00 ,  10 ,  01 , and  20 averaged over 1, 000 realizations and the homogeneous direction .Some expected characteristics of the measured moments are observed: (i) The leading order moment is over two magnitudes larger than the molecular diffusivity (10 −9  2 /).The scaled higher order moments shown are all at least one magnitude larger than the molecular diffusivity.
(ii)  00 is symmetric and always positive.
(iii)  10 is antisymmetric.This antisymmetry can be understood by interpreting  10 / 00 as the centroid of the eddy diffusivity kernel.Physically, for  > 0, it is expected that the mean scalar gradient at the center of the mixing layer (at a negative distance away) has more influence on the turbulent scalar flux than the mean scalar gradient at the outer edges, since the mixing layer is growing outwards.This makes the eddy diffusivity kernel skewed more towards the center of the domain, so  10 < 0 for  > 0. A similar effect occurs for  < 0, which results in  10 > 0.
(iv)  01 is symmetric and always negative.The latter must be true for the flow to depend on its history (it does not violate causality).
(v)  20 is symmetric and always positive, as is characteristic of moment of inertia of a positive kernel.
Based on the magnitudes of the normalized moments, some initial observations on importance of each moment can be made. 00 has the highest magnitude of all the moments, which is expected since it is the leading-order moment.The magnitudes of  10 /ℎ and  01 / are on the order of 10% of the magnitude of  00 , which suggests that they are non-negligible.On the other hand, the magnitude  20 /ℎ 2 is on the order of 1% that of  00 , so  20 is likely not an important moment to include in modeling RTI.More robust studies will be presented in the following sections to determine importance of each of the eddy diffusivity moments.
It is also observed that there is statistical error in the measurements.Due the chaotic nature of RTI, the moment measurements contain statistical error, but this error can be reduced by averaging many realizations.To demonstrate statistical convergence of the measurements, plots of  00 averaged over different numbers of realizations are included in figure 8.As the number of realizations increases, the plots become smoother, and it is found that after 1, 000 realizations, the rate of reduction in statistical error slows down significantly.Averaging over this number of realizations results in a smooth  00 measurement and higher-order moment measurements with acceptably less statistical error.Additionally, the higher the order of the moment, the slower its rate of statistical convergence.Recall that determination of higher-order moments requires information from lower-order moments.For example, in determining  01 ,  00 is subtracted from  01 , the turbulent scalar flux measurement in the simulation associated with  01 .Naturally, there is statistical error associated with both  01 and  00 .However, the error in  00 is amplified by , so the overall statistical error of  01 increases with time.This statistical error "leakage" occurs for all higher-order moments.The higher the order of the moment, the worse the statistical error, since information from more lower-order moments is needed, and so more statistical error is accumulated and amplified.The relatively high statistical error of the higher-order moments makes it challenging to study their importance.Particularly, taking derivatives of quantities with high statistical error amplifies the error, so smoother measurements are desired.In this work, the moment measurements are smoothed using a Savitzky-Golay filter function in Matlab with a polynomial order of unity and window size of 191.These smoothed moments are shown in figure 9.While it is possible to design an alternative formulation of MFM that removes leakage of statistical error from low-order moments to higher-order moments (see Lavacot et al. 2022), for this 2D study and for the order of moments considered here, the statistical convergence is sufficient.
Using these measurements, nonlocal timescales and lengthscales (   and    , respectively) can be defined: (5.1) Note that this analysis can only be done for −1 ⩽  ⩽ 1, since the moments are analytically zero outside the mixing layer.Nondimensionally, the nonlocal timescale is    =    / 0 , and the nonlocal lengthscale is    =    /ℎ.Contour plots of the nondimensionalized nonlocal timescale and lengthscale are in figure 10.Note that    scales as , so profiles of    / against  are also plotted in figure 11 in the self-similar time regime ( > 17).The scaled profiles collapse and have a centerline value of approcimately 0.1.This means that the mean fluxes at some time  are affected by mean scalar gradients 0.1 earlier.Figure 12 shows the minimum nonlocal lengthscale is at the centerline, where    ≈ 0.09.The maximum lengthscales occur near the outer edges of the mixing layer: at around  = ±0.87,   ≈ 0.27.This indicates that mean fluxes at the mixing layer edges depend mostly on mean scalar gradients about a quarter of a mixing width away, while mean fluxes at the centerline depend on mean scalar gradients about one tenth of a mixing width away; nonlocality appears to be stronger at the mixing layer edges than at the centerline.These nonlocal properties of the eddy diffusivity for RTI could not be predicted without direct measurement of the eddy diffusivity moments, which has been made possible through MFM.

Comparison of terms in turbulent scalar flux expansion
To aid in the determination of which moments are important for a RANS model, a comparison of the terms in the expansion of the turbulent scalar flux (equation 3.6) is presented.These terms involve gradients of ⟨  ⟩.Instead of using ⟨  ⟩ directly from the DNS, a fit to ⟨  ⟩ is used, since the statistical error in the raw measurement gets amplified by derivatives in .That is, the quantities of interest are sufficiently converged for plotting but not for operations involving  derivatives.Thus, an analytical fit to ⟨  ⟩ is obtained as follows: where the integral is determined numerically, and  and  are fitting coefficients.The coefficients  2 = 1.2 and  = 0.36 are found to give good agreement to the mean concentration profile from DNS, as shown in figure 13.The terms on the right hand side of equation 3.6 are plotted against the DNS measurement of the turbulent scalar flux in figure 14.Clearly, the  00 term is not enough to capture the turbulent scalar flux.It is observed that the  01 term is significant in magnitude in the middle of the domain, and the  10 term carries importance at the outer edges of the mixing layer.The term associated with the highest-order moment that was measured,  20 also appears to be of similar magnitude as the other moments, indicating it may also carry important information about nonlocality of the eddy diffusivity.These preliminary findings indicate nonlocality is certainly important for accurate modeling of mean scalar transport in this RTI problem, since the higher-order terms in equation 3.6 appear non-negligible compared to the leading-order term.It may be tempting to ascribe physical reasons for the behavior of the terms plotted in figure 14, but this is not so straightforward, especially since the full eddy diffusivity kernel for this problem as not yet been measured.Further, it would be inappropriate to draw conclusions about importance of each eddy diffusivity moment in a RANS model, since the operator form must be scrutinized first.A faulty operator form could give misleading implications about certain eddy diffusivity moments.It turns out that a simple superposition of these terms, which would represent a truncation of equation 3.6, does not accurately represent the true eddy diffusivity kernel and actually leads to divergence of predictions, so such an operator form would not be appropriate; this will be covered more in depth later in §5.3.Nevertheless, the results shown here are strong evidence of nonlocality of the eddy diffusivity kernel for the RTI simulated here.

Comparison of leading-order model against a local model
To demonstrate the shortcomings of models using purely local coefficients, an MFM-based leading-order model and the - RANS model are compared.The intent of this study is not to immediately propose a "better" RANS model to replace -, nor is it to suggest the MFMbased leading-order model is more accurate than the - model.In fact, it is expected that the MFM-based leading-order model will perform poorly, since it does not include important higherorder moments of eddy diffusivity.Instead, this study emphasizes the necessity of higher-order moments and shows how MFM can reveal incorrect model forms.
In particular, a 1D - simulation is run, and the eddy diffusivity and mean concentration profiles are extracted from the results to be compared to those of the MFM-based model using the measured  00 that was presented in §5.1.The - simulation used in this section is implemented in Ares, and details of the implementation are in Morgan & Greenough (2015) and Morgan (2018).Note that the - simulation is used here for illustration purposes and should not be confused with the 2D DNS simulations used to obtain our MFM moments.
The MFM-measured  00 is used for the leading-order MFM-based model: (5.4) To solve this,  00 MFM is obtained from the smoothed MFM measurements and transformed to self-similar coordinates.The resulting  00 MFM is a function of  =  ℎ DNS , where ℎ DNS =  * DNS ( −  * DNS ) 2 is an algebraic fit to the mixing width from the DNS.The equation is then solved semi-analytically in conjunction with the mean mass fraction evolution equation in selfsimilar coordinates: (5.6) The - model uses the gradient diffusion approximation for the turbulent flux: where   =   ⟨⟩ √ 2.  is one of the model coefficients set by similarity constraints derived by Dimonte & Tipton (2006).Particularly, this work uses the coefficient calibration detailed in (Morgan & Greenough 2015), and the coefficients are chosen to achieve the same  as the DNS.Here,   is unity and   is 2.47.The - RANS model is solved in spatio-temporal coordinates, and the ⟨⟩, , and  obtained from the solution are used to compute   and, consequently,  00 - , which is purely local.For a meaningful comparison with the MFM-based model,  00 - is transformed to  00 - according to the self-similar coordinate  =  Figure 15a shows the mean concentration profiles computed using each of the two models.As expected, the MFM-based leading-order model performs poorly, not capturing the slope of the DNS profile, since that model only uses the leading-order eddy diffusivity moment and incorporates no information about nonlocality of the eddy diffusivity.The - model exhibits divergence from DNS at the outer edges of the mixing layer, since it is designed to predict a linear ⟨  ⟩ profile.However, it does capture the slope of the DNS profile, despite it also using a leading-order closure.In addition, it is observed in figure 15b that the MFM-measured  00 MFM is significantly lower in magnitude than  00 - .Here, MFM reveals that the - model is using an incorrect model form, since the  00 it is using does not match the MFM measurement.In fact, the - model is using this higher-magnitude coefficient in order to compensate for the error in model form and achieve a linear mean concentration profile with a slope that matches DNS.Despite this compensation, the - model still disagrees with DNS results at the outer edges of the mixing layer, which are important to capturing the average reaction rate in reacting flows like in ICF.A more accurate RANS model would more closely match the eddy diffusivity moments measured by MFM.As results will show shortly, the gap between the leading-order MFM-based model ⟨  ⟩ and the DNS measurement would be bridged by inclusion of higher-order moments, which would introduce information about the nonlocality of the eddy diffusivity.

Assessment of nonlocal operator forms
In this section, two RANS operator forms using information about the nonlocality of the eddy diffusivity are presented.These are the explicit and implicit operator forms; the former is a truncation of the turbulent scalar flux expansion 3.6, and the latter will be presented shortly.It must be stressed that the intention of the following studies is not to propose a new RANS model.Ultimately, a RANS model should not depend on direct MFM measurements that can only be retrieved from impractically many DNS.Instead, these studies are performed to further assess the importance of each of the eddy diffusivity moments, determine which combinations of moments best enhance the performance of a RANS model, and examine the differences between the explicit and implicit operator forms.The aim of these studies is to inform development of more predictive RANS models for RTI, not to suggest that these are the exact models that should be used.
In addition,  20 will not be included in the following studies.This is mainly due to the high statistical error in the measurement that makes it difficult to ascertain whether errors in the results are due to this statistical error or solely the addition of the moment to the model.From the comparison of terms in §5.2.1, it is expected that  20 is not as important as  10 and  01 to include in a RANS model.This should be tested in future work when a more statistically-converged measurement is achieved for  20 , ideally in a 3D analysis.

Explicit operator form
The explicit operator form is a truncation of the expansion of the turbulent scalar flux, as defined in equation 3.6.Hamba (1995) and Hamba (2004) have examined this form in the context of shear flows.Transformation of this expansion to self-similar coordinates and substitution into 3.3 results in which can be solved numerically for ⟨  ⟩.The   used in the numerical solve are the smoothed, normalized moments.To determine which eddy diffusivity moments are important in constructing RANS models for RTI, different combinations of   terms are kept in equation 5.8, and the results are compared to DNS.In the numerical solve, equation 5.8 is discretized on a staggered mesh, and derivatives are computed using central finite differences.A matrix-vector equation is assembled and solved for ⟨  ⟩ with Dirichlet boundary conditions.
Figure 16 shows the turbulent scalar fluxes computed using the explicit operator form, and figure 17 shows the corresponding mean concentration profiles.Again, it is apparent that the leading-order moment is not enough to capture the turbulent scalar flux.The combination using  00 ,  10 , and  01 -the moments deemed most important in §5.2.1-gives the best match to the DNS measurement.
It is particularly remarkable that a converged turbulent scalar flux can be obtained using  00 ,  10 , and  01 .As mentioned previously, it is known that equation 3.6 may not converge.That is, the expansion must be taken to infinite terms to remove error; truncating the expansion can result in significant error.This is analogous to a Kramers-Moyal expansion, which cannot be approximated adequately by more than two terms, after which it requires infinite terms for convergence (Pawula 1967;Mauri 1991).To understand how adding terms to equation 3.6 can result in greater error, one can consider the eddy diffusivity kernel associated with each term.The leading-order moment is associated with a delta function kernel, as it is purely local.However, when equation 3.4 is replaced by equation 3.6, an integral operator is replaced with a high-order differential operator.This means that the nonlocal effects are approximated by derivatives of delta functions; see Liu et al. (2023) for more details.It has been shown that, in general, the eddy diffusivity kernel is not a superposition of finite delta functions, as it is smooth (Mani & Park 2021;Liu et al. 2023).Therefore, truncation of the expansion does not match the shape of the eddy diffusivity kernel, leading to errors in prediction of the turbulent scalar flux.While the  00 ,  10 , and  01 combination did not diverge, adding  20 does lead to divergent results for these reasons, so this combination is not presented here.
Another issue with the explicit operator form is its numerical implementation.In spatio- temporal space, some terms associated with higher-order moments involve mixed derivatives (e.g., the term  01  2   ), which would undergo another spatial gradient when substituted into equation 3.3.Such terms are difficult to handle numerically.In this work, the model is implemented in the more convenient self-similar space, but, ultimately, a spatio-temporal model would be developed, as it is more practical.It is thus pertinent to work towards a better method to incorporate nonlocal information in a RANS model that does not encounter the Kramers-Moyal-like convergence issue and is easier to implement.

Implicit operator form and the Matched Moment Inverse (MMI)
In this section, an implicit operator form is introduced as a solution to both the increasing error when adding terms from the turbulent scalar flux expansion and implementation challenges associated with the explicit operator form.Recall that the explicit operator form fails to match the shape of the eddy diffusivity kernel without infinite terms of the turbulent scalar flux expansion.In this implicit operator form, the aim is to match the shape of the eddy diffusivity kernel, instead of using the truncated expansion for the turbulent scalar flux.Using the four moments that have been measured, this model form is where   (, ) are model coefficients fitted corresponding to each of the eddy diffusivity moments   measured using MFM.The bracketed operator on the left hand side is the Matched Moment Inverse (MMI) operator.The way this model form is designed to match the eddy diffusivity kernel shape is detailed in Liu et al. (2023).In addition, this form is significantly easier to implement numerically in spatio-temporal space, since it can be directly time-integrated using explicit methods.In this way, it is also easy to add more terms with higher-order moments, as it simply requires extension of the operator.
In self-similar coordinates, this becomes where it is found through self-similar analysis that (5.12) (5.14) The coefficients are determined through a process illustrated as follows in spatio-temporal coordinates for simplicity.If one wants to construct a model in the form of equation 5.9, four equations must be formulated to determine the four coefficients.This is done by using measurements from the four simulations used to determine the four moments  00 ,  10 ,  01 , and  20 .For example, the first equation results from substitution of  00 for −⟨ ′  ′  ⟩ and the associated desired ⟨  ⟩  ; the remaining three equations follow, using the other three moments: 1 (5.18) This system of equations is then rearranged into a matrix equation  MMI a = b, which is solved for the coefficients in vector a = ( 00 ,  10 ,  01 ,  20 )  .Note that this matrix equation is constructed over every point in space and time, so a = a(, ).In this work, analysis is done in self-similar coordinates, in which a = a().If one wishes to construct a model with different moments, the MMI operator and equations must be modified accordingly.For example, a model using only  00 and  10 would have an MMI operator of the form 1 +  10   and use only the first two equations (with the  01 and  20 terms removed).Thus, models using different combinations of moments would use different MMI coefficients   .
To summarize, for this implicit operator form, the following system of equations is solved in self-similar coordinates: where L MMI is the MMI operator constructed using some combination of moments, such as in equation 5.10.Numerically, the following system is solved: (5.21) (5.22)where  is the matrix representing the numerical MMI operator, and D  is the matrix representing the numerical derivative with respect to .This can be rewritten as a block matrix-vector multiplication: where b is a vector representing the right-hand side of equations 5.22, with the proper boundary conditions enforced.In this study, zero gradient boundary conditions are used for the turbulent scalar flux, and Dirichlet boundary conditions are used for the mean concentration.The system is solved using finite differences on a staggered mesh.
Presented in this work are the determinants of the MMI matrix and resulting   for two different combinations of moments.Figure 18 shows that with the combination of  00 ,  10 , and  01 , the determinant of the MMI matrix is positive for all , so   are all well-behaved.This is indicative of good model form.On the other hand, figure 19 shows that with the combination of  00 and  01 , the determinant of the MMI matrix crosses zero, so   contain singularities which effect poor RANS predictions (observable in plots presented later).Singular matrices arising in the MMI solve for a certain form of the implicit operator form may indicate that form is poor, in the sense that the combination of moments does not make a good RANS model.Since MMI appears to be sensitive to the information it takes in to determine the implicit operator coefficients, one must take special care and choose a model form that avoids this issue.It is found that MMI determinant zero-crossings do not occur for any of the moment combinations tested in this work other than the  00 and  01 combination, but it may happen with combinations of other higher-order moments not measured here.
Turbulent scalar fluxes computed using the implicit operator form are shown in figure 20.The implicit operator form's turbulent scalar flux prediction using just  00 is identical to that of the explicit operator form, by construction.It is apparent that adding either  10 or  01 alone is insufficient.As noted earlier, adding  01 leads to a particularly poor prediction due to singular MMI matrices at some .The best match to DNS is attained using the combination of  00 ,  10 , and  01 .In fact, it is evident that the implicit operator form using  00 ,  10 , and  01 predicts the turbulent scalar flux more accurately than the explicit operator form using the same moments.This is because the implicit operator form is designed to match the shape of the eddy diffusivity kernel, and the explicit operator form may not be accomplishing this.
These trends in the explicit and implicit operator forms can be observed again in the predictions of the mean concentration profile, shown in figure 21.Particularly, the implicit operator form using  00 ,  10 , and  01 gives a very good prediction of the mean concentration that nearly overlaps the DNS measurement.For a clearer comparison of the explicit and implicit operator forms using these moments, figure 22 shows the derivatives of the DNS-and model-computed ⟨  ⟩.The implicit operator form predicts a magnitude and shape closer to the DNS measurement than the explicit operator form does.In particular, the implicit operator form captures the shape of the tails much better than the explicit operator form.

Conclusion
In this assessment, it is determined that nonlocality must be considered in developing more predictive models for RTI.The studies presented in this work are facilitated using MFM, a numerical tool for precisely measuring closure operators.Four of the eddy diffusivity moments of RTI ( 00 ,  10 ,  01 , and  20 ) are measured, and it is demonstrated that the higher-order moments, which contain information about the nonlocality of the eddy diffusivity kernel, should not be neglected when constructing models for RTI.
Specifically, it is determined that  00 ,  10 , and  01 are the most important moments for constructing a model for RTI.Two methods for constructing RANS models using these moments are presented.First, an explicit operator form, based on a Kramers-Moyal-like expansion derived by taking the Taylor series expansion of the scalar gradient in the generalized eddy diffusivity, is described and tested.While incorporation of higher-order moments in the explicit operator form results in more accurate predictions than a leading-order model, there exist several issues.One problem is that the expansion used for the explicit operator form may not converge, so addition of higher-order moments leads to less accurate predictions.Another problem is that the explicit operator form is difficult to implement numerically.
Thus, an implicit operator form is presented to address these issues with the explicit operator form.Since an implicit operator form involves an invertible matrix operator, it is easier to implement than an explicit operator form.In addition, the proposed implicit operator form is  designed to match the shape of the eddy diffusivity kernel via the MMI operator, in contrast to the explicit operator form, which truncates a non-converging Kramer-Moyal expansion.It is shown that the implicit operator form exhibits a marked improvement in predictions over the explicit operator form.Incorporation of nonlocality into RANS models via these operator forms comes with several challenges.For one, development of any new model must consider scalar realizability.While this is not thoroughly explored in this work, since an actual model is not yet proposed, one approach to preserve realizability is suggested by Braun & Gore (2021), where the turbulent scalar flux is rewritten as an advection-like term and added to the original advection term in order to enforce physical mean component mass fractions; a conservative numerical scheme maintains realizability.Further, the new model must be tested on more complex RTI for it to be useful in practical settings such as ICF.This includes assessment of the model for 3D, finite Atwood, and compressible (Richtmeyer-Meshkov) flows.The model should be tested in the same validation cases as other models for RTI, such as the tilted rig (Denissen et al. 2014) and gravity reversal (Banerjee et al. 2010).Based on these evaluations, which may also involve new MFM measurements where the method is extended to more complex flow regimes, the new model can be amended, as is the usual process of turbulence model development.This is left for future work, when a new model is developed based on the findings presented here.
One obstacle encountered in these studies is the inherent statistical error in the DNS com-  putations.The higher-order moments particularly contain high statistical error due to buildup of error from the lower-order moments on which they depend.Because of this, it is admittedly difficult to draw definite conclusions about the effect of higher-order moments beyond the firstorder moments.That is, due to the statistical error, it is currently unclear if inclusion of moments beyond first-order in a RANS model would significantly improve its predictions, or if just the first-order temporal and spatial moments along with the leading-order moment are sufficient.This motivates development of a technique to accelerate the statistical convergence of these higher-order moments.Such a method could also be used to study the effect of other higher-order moments that were not measured in this present work, since they would have suffered from high statistical error with the current method.
It must be stressed that the results in this work are for 2D RTI and should not be directly applied to 3D RTI.As noted previously, the third spatial dimension significantly impacts the turbulent physics of RTI.In particular, 3D RTI has a lower growth rate than 2D, so lower magnitudes of the eddy diffusivity moments are expected in 3D.Despite the quantitative difference in physics between 2D and 3D, they are qualitatively similar in the RANS space, so trends in the shapes of the eddy diffusivity moments are expected to persist in 3D.In other words, the form of the turbulent scalar flux closure in 3D is expected to be the same as in 2D, but the coefficients would be different.These expected trends are yet to be confirmed, and future work should involve applying MFM to 3D RTI.
Through this work, an understanding of nonlocality in 2D RTI has been developed.It has been shown that incorporation of information about the nonlocality of the eddy diffusivity may greatly improve the accuracy of a RANS model.This work demonstrates this by testing operators using MFM measurements of the nonlocal eddy diffusivity.In practice, a RANS model for RTI would not have to rely on these MFM measurements directly; one would not have to perform many MFM simulations to construct a model.In other words, MFM should be seen as a diagnostic tool rather than the means for building the actual model.The ultimate goal is to develop an improved, more predictive model for RTI by incorporating nonlocal information, which the present work has demonstrated to be significant for accurate prediction of mean scalar transport in 2D RTI.
To nondimensionalize the eddy diffusivity moments, equation 3.6 is substituted into equation 3.3: This reveals nondimensionalizations for the eddy diffusivity moments.The prefactors to the derivatives of ⟨  ⟩ on the right hand side are denoted as the normalized eddy diffusivity moments   .The turbulent scalar flux scales with the leading-order term in equation 3.6.Substitution of the nondimensionalization for  00 (equation 3.25) into the leading-order term in equation 3.6 and transformation to self-similar coordinates gives the scaling for the turbulent scalar flux: Figure 23 shows the unscaled moments measured directly from the MFM simulations.The profiles are taken from the portion of the simulation where the flow is self-similar ( ⪆ 17).It is obvious that without normalizing the moments as described above there is no self-similar collapse.The moments are scaled and plotted against  in figure 24

Figure 2 :
Figure 2: Self-similarity parameters computed from a donor simulation.

Figure 5 :
Figure 5:   contours (black: light, white: heavy) and velocity vector fields (red arrows) of the (a) donor simulation and (b) receiver simulation with  enforcing ⟨  ⟩ = .These snapshots are taken at the last timestep.

Figure 7 :
Figure 7: Moments of eddy diffusivity kernel normalized by appropriate length and timescales.Data averaged over 1,000 realizations and homogeneous direction .

Figure 9 :
Figure 9: Smoothed moments (dashed red) over raw MFM measurements of moments (solid black).The moments are taken from the mean data at the last timestep of the simulations and are transformed to self-similar space.

Figure 12 :
Figure 12: Nondimensional nonlocal lengthscale profiles at different times for  > 17.Lighter lines correspond to later times; darker lines correspond to earlier times.

Figure 14 :Figure 15 :
Figure 14: Comparison of terms in the expansion of the turbulent scalar flux.Black: DNS measurement of turbulent scalar flux, solid blue:  00 term, dashed pink:  10 term, dash-dotted green:  01 term, dotted orange:  20 term.
Figure 16: Turbulent scalar flux predictions using the explicit operator form.Captions of each plot list moments used in the model.
Figure 17: Mean concentration profile predictions using the explicit operator form.Captions of each plot list moments used in the model.
Figure 18: (a) Determinant of  MMI over  and (b, c, d) MMI coefficients   over  for the implicit operator form using  00 ,  10 , and  01 .
Figure 19: (a) Determinant of  MMI over  and (b, c) MMI coefficients   over  for the implicit operator form using  00 and  01 . MMI is singular at the  at which its determinant crosses zero.

Figure 20 :
Figure 20: Turbulent scalar flux predictions using the implicit operator form.Captions of each plot list moments of eddy diffusivity used in the model.
Figure 21: Mean concentration profile predictions using the implicit operator form.Captions of each plot list moments of eddy diffusivity used in the model.
Figure23shows the unscaled moments measured directly from the MFM simulations.The profiles are taken from the portion of the simulation where the flow is self-similar ( ⪆ 17).It is obvious that without normalizing the moments as described above there is no self-similar collapse.The moments are scaled and plotted against  in figure24to demonstrate the selfsimilar collapse.The normalized turbulent scalar flux and mean concentration profiles are shown in figures 25 and 26, also showing self-similar collapse.

Figure 23 :Figure 24 :
Figure 23: Unscaled moments at different times over .Dimensions of the moments are as follows: 00 is  2 /,  10 is  3 /,  01 is  2 , and  20 is  4 /.Red profiles are from early times, and cyan profiles are from late times, when the flow is self-similar.Lighter lines correspond to later times; darker lines correspond to earlier times.

Figure 25 :
Figure 25: Self-similar collapse of turbulent scalar flux.− ⟨ ′  ′ ⟩ is dimensionless.Red profiles are from early times, and cyan profiles are from late times, when the flow is self-similar.Lighter lines correspond to later times; darker lines correspond to earlier times.

Figure 26 :
Figure 26: Mean concentration profiles at different times.Red profiles are from early times, and cyan profiles are from late times, when the flow is self-similar.Lighter lines correspond to later times; darker lines correspond to earlier times.