In this study, the two-level simulation (TLS) modeling strategy developed earlier for momentum mixing is assessed in turbulent flows with passive scalar mixing and transport. Compared to the conventional algebraic and transport equation based subgrid-scale closures for large-eddy simulation (LES), where the effects of the unresolved small scales (SS) of motion on the large scales (LS) are modeled, in the TLS model, both large and small scales are explicitly evolved in a coupled two-scale decomposition method called TLS. Earlier studies demonstrated the ability of the SS model in the TLS to simulate high Reynolds number wall-bounded flows with a reasonable LS grid resolution. In the present study, this two-scale strategy is extended to model passive scalar mixing in turbulent channel flows. A comprehensive *a priori* assessment of the scalar SS model is performed using direct numerical simulation datasets for passive scalar transport in a fully developed turbulent channel flow at frictional Reynolds number $ R e \tau = 395$ and at Schmidt numbers $ S c =$ 1 and 4. This assessment is carried out to study the evolution of the large and the small scales of the scalar field in physical and spectral space and also to assess the SS prediction of quantities relevant for modeling of scalars such as scalar dissipation rate, counter-gradient transport, and backscatter. A TLS model is, then, used in *a posteriori* near-wall study using a hybrid TLS–LES approach to assess the two-scale scalar model.

## I. INTRODUCTION

Turbulent scalar transport is observed in several flows of interest, such as in stratified flow in water, heat and mass transfer in internal flows, in high-speed aerospace systems, and in pollutant dispersion. We focus on simulating passive scalar transport and mixing both due to the large-scale (LS) macromixing and the localized small-scale (SS) micromixing. In particular, micromixing effects of straining and shearing turbulent motions can lead to a large magnitude of the scalar gradient, enhanced molecular diffusion of the scalar field, and increased dissipation of scalar fluctuations. An accurate description of the micromixing is, therefore, critical to study the dynamics of the scalar field. However, both experimental and numerical investigations of scalar micromixing, for example, in liquids, suffer from the challenge of resolution requirement, as in such flows, the Schmidt number $ S c = \nu / \kappa \u223c O ( 100 )$, or higher, which makes the smallest length scale for the scalar field, the Batchelor scalar $ \eta B = \eta K / S c 1 / 2$ much smaller than the smallest flow length scale, and the Kolmogorov scale $ \eta K$.^{1} Here, *ν* is the kinematic viscosity and *κ* is the scalar diffusivity. The excessive resolution needed to capture scalar micromixing renders the direct numerical simulation (DNS) strategy, computationally intractable, thus requiring alternate strategies. This study assesses this two-scale modeling approach as a cost-effective simulation strategy.

The two-level simulation (TLS) approach is not the only modeled methodology to simulate unsteady scalar mixing and transport. Large-eddy simulation (LES) strategies with subgrid scalar closure models based on eddy diffusivity have predicted many features of scalar transport but have also shown some disagreements. For example, eddy diffusivity models^{2} disregard several key aspects of the subgrid-scale closure (SGS) physics such as the presence of counter-gradient transport, anisotropy, and backscatter. Alternate modeling approaches that alleviate some of these shortcomings have also been proposed, such as models based on the scale-similar approach,^{3} mixed formulation,^{4} the nonlinear model,^{5} models based on interscale energy transfer,^{6,7} models based on the hybrid Reynolds-averaged Navier-Stokes (RANS)-LES strategy^{8} such as the dynamic delayed detached eddy simulation,^{9} variants of the gradient model,^{10–12} models employing the generalized gradient diffusion hypothesis for defining the tensor diffusivity,^{13–15} the linear eddy model,^{16} and the multiscale modeling strategy.^{17}

Compared to the conventional LES strategy, which uses implicit grid filtering, the multiscale strategy^{18–21} employs scale decomposition to explicitly represent quantities in terms of its LS and SS components. While the LS is solved on the resolvable three-dimensional (3D) domain, the SS component is explicitly obtained by solving a modeled set of equations (i.e., using a subgrid simulation) instead of using the SGS eddy diffusivity type closure models. These scale decomposition methods primarily differ in their strategy to obtain the SS field and, here we focus on one such a strategy, the TLS model^{21} that has shown promise as an effective simulation approach for high Reynolds (*Re*) number isotropic^{21} and wall-bounded^{22} turbulent flows using relatively coarse grid resolutions. The earlier studies showed that many SS physics such as vorticity, counter-gradient transport, and backscatter of momentum can be captured.^{23} A hybrid TLS–LES approach has also been developed^{24} to allow for efficient computation of high *Re* wall-bounded turbulent flows where the TLS model is used in the near-wall region and the conventional LES is used in outer regions. The approach has been extended further to allow for a dynamic blending of the TLS and LES approaches in complex flows so that TLS remains active in regions of flow exhibiting features such as anisotropy, nonequilibrium, and separation/reattachment, and LES is used elsewhere.^{25}

In the present study, the model is examined specifically for passive scalar mixing and transport as a first step toward future application to more complex flows involving reactive scalars. We first summarize the formulation and the modeling assumptions of the TLS model for scalar and momentum mixing and then assess its ability using *a priori* analysis of DNS datasets ( $ R e \tau = 395$ and Schmidt numbers *Sc* = 1 and 4). For some application demonstration, we show *a posteriori* use of the TLS model as a near-wall model within the hybrid TLS–LES strategy for the same conditions but using a relatively coarser LS grid. We have considered only two values of *Sc* in this study to allow for both *a priori* and *a posteriori* assessment of the TLS model by comparing with the corresponding DNS. As stated before, DNS at higher values of *Sc* tends to be computationally expensive; therefore, in this study, we limit ourselves to only a moderate increase in the value of *Sc*. However, as discussed later, the effects of *Sc* are apparent on the statistical features of the passive scalar transport. Therefore, the chosen values of *Sc* can be considered adequate for the demonstration of the capabilities of the TLS model in capturing the effects of a moderate increase in *Sc*. Furthermore, a detailed comparison of modeling strategies relying on the different representation of the subgrid quantities is out of the scope of the present study. However, the effectiveness of the TLS model is demonstrated by comparing it with the widely used dynamic eddy diffusivity model^{26} and a no-model (NM) approach for the SGS scalar flux, while employing the same grid and the same numerical approach for the resolved LS fields. More detailed uses of hybrid TLS–LES for a range of $ R e \tau $ and *Sc* values and a detailed comparison with other models will be reported in the near future.

This article is arranged as follows. The governing equations and modeling approaches are described in Sec. II. The details of the computational setup and the numerical methodology are presented in Sec. III. A comprehensive assessment of the TLS model is discussed in Sec. IV. Finally, a summary of the key findings of this study, and future outlook are presented in Sec. V.

## II. MATHEMATICAL FORMULATION

The TLS approach and the resulting governing equations for both LS and SS fields have been described in detail in the past.^{21} Here, we repeat them focusing on passive scalar transport. A brief description of the hybrid TLS–LES strategy^{24} for near-wall modeling is also included for completeness.

### A. Governing equations

The incompressible Navier–Stokes equations along with passive scalar transport are given by

*u*is the velocity component,

_{i}*p*is the physical pressure divided by constant density

*ρ*, $\varphi $ is the passive scalar, and

*i*= 1, 2, and 3 are the Cartesian coordinate directions

*x*,

*y*, and

*z*, respectively. These equations are complete after specifications of proper initial and boundary conditions.

In DNS, the governing equations given by Eq. (2.1) are solved wherein all the relevant length- and time-scales associated with both the flow and the scalar fields are well resolved on the computational grid. However, for scenarios that occur at high *Re* and for flows where *Sc* > 1, DNS tends to be computationally intractable, and alternative cost-effective techniques, such as LES or, as studied here, TLS, are needed.

### B. Two-scale decomposition

The TLS model uses scale separation of any variable $ \psi ( x , t )$ into its LS and SS components. Here, the LS function is denoted by $ L \Delta $, and although it can be defined in many ways,^{19,21} in the TLS model, it is defined on the underlying LS grid $ G \Delta $. Thus,

where superscripts “L” and “S” denote large- and small-scale components of $ \psi ( x , t )$, respectively. The LS field is defined as^{21}

where the LS function $ L \Delta $ is the application of a local averaging operator $ S \Delta $ followed by an interpolation operator $ I \Delta $ with $ x k$ representing the nodes corresponding to the LS grid. The local averaging operator $ S \Delta $ can be defined in different ways. It depends upon the LS grid $ G \Delta $ and the approach to determine the discrete LS value of a field variable. In the present study, it is specified to be a sampling operator, where the LS value of a field variable is defined as the value of the variable at the nodes of $ G \Delta $, i.e., $ \psi L ( x k , t ) = \psi ( x k , t )$. Alternatively, in complex flows, a local volume-based averaging operator can be used to specify $ S \Delta $. The LS component of a field variable is uniquely defined once the operators $ S \Delta $ and $ I \Delta $ are specified.

Applying $ S \Delta $ on $ \psi ( x , t )$ yields a discrete representation of the LS field denoted by $ \psi L ( x k , t )$, which is defined at $ x k$. This is, then, interpolated into the SS (finer) grid in order to obtain a continuous representation of the LS field denoted by $ \psi L ( x , t )$. Once the continuous representation of the LS field is obtained, the corresponding SS field (denoted by *ψ ^{S}*) is obtained using Eq. (2.2). The SS grid is much finer than the LS grid and is chosen to ensure resolution of the relevant length scales in a manner similar to DNS.

Although the decomposition presented in Eq. (2.2) appears analogous to the classical Reynolds decomposition, there are some important differences.^{18,19,21} For example, the large-scale operation on the SS field is not zero [i.e., $ ( \psi S ) L \u2260 0$]; the SS part of any LS quantity is not zero: $ ( \psi L ) S \u2260 0$; the products always have nonzero LS and SS parts: $ ( \psi L \psi L ) S \u2260 0$ and $ ( \psi S \psi S ) L \u2260 0$; and the LS operator is not idempotent: $ ( \psi L ) L \u2260 \psi L$.

An example of scale separation is illustrated in Fig. 1, where a representative field $\varphi $ obtained from the DNS of a turbulent channel flow is used to obtain its LS and SS components. The discrete LS field is first obtained on 16 uniformly spaced points from the DNS data given on 256 points in a streamwise oriented line. Afterward, the discrete LS field is used to define the continuous LS field using a cubic B-spline interpolation, which is, then, used along with the DNS data in Eq. (2.2) to obtain the SS field. It can be observed that the high wavenumber (small scale) features are not in the LS field but present in the corresponding SS field. The goal of the TLS approach is to provide a modeled representation of the SS field in the 3D LS evolution to account for many of the key features in the DNS field.

### C. TLS equations

Applying the two-scale decomposition [Eq. (2.2)] to Eq. (2.1) leads to the following LS and SS equations:^{21}

Here, Eqs. (2.6) and (2.9) represent the scalar extension in the TLS approach.^{21} The terms $ F i L$ and $ G L$ couple the LS velocity and scalar fields to the SS velocity and scalar fields, respectively. These terms are given by

The combined LS and SS equations when expressed in 3D recover the full DNS equations given by Eq. (2.1), thus offering no computational advantage. In the TLS approach, however, a SS closure model^{21} is used to reduce the computational cost of *simulating* the SS, and this SS is included in the coarse-grained LS equation [Eqs. (2.4) and (2.6)]. The TLS closure for the velocity and the scalar fields here are identical to the earlier approach,^{21} and so, only some details are repeated here.

### D. Modeling of small scales

In the TLS model, a grid-within-grid strategy is employed (illustrated in Fig. 2), wherein the LS equations are solved on a 3D grid using a conventional solver, while the SS equations are solved on a collection of three orthogonal one-dimensional (1D) lines embedded within the LS grid. The SS grid is chosen to resolve all the smallest scales (e.g., Batchelor or Kolmogorov scale depending on *Sc* and *Re*) of interest along the respective direction, and a modeled form of the SS equations (described below) is solved on these 1D lines. In the original approach,^{21} the SS lines spanned the entire domain, but in the hybrid TLS–LES approach,^{24} the SS domain is localized to the regions of interest, for example, in a wall-bounded flow, the TLS model is only used in the near-wall region [Fig. 2(b)]. Note that the orientation of the three 1D orthogonal lines along Cartesian directions has been chosen to have an easier numerical implementation, but the governing equations are general and SS lines can be oriented such that the implementation suits a particular numerical method. For example, when an unstructured LS grid is used, the 1D lines for small scales can be chosen to be lines connecting adjacent cell centers.

Although Fig. 2(a) suggests a Cartesian strategy,^{21} the two-scale decomposition is independent of the grid structure of the LS or the need for SS orthogonal lines as long as the scale separation procedure [Eq. (2.2)] is followed. The hybrid TLS–LES model^{24} is an even more computational simplification of the full TLS model.^{21}

In TLS, the SS equations are solved along the three 1D orthogonal lines, which leads to multiple realizations of the SS field. For example, there is one equation for each SS velocity component, $ u i S$, along each orthogonal direction. Thus, in a given LS grid, there are three realizations of each of the velocity component, and similarly, three realizations of $ \varphi S$ along each orthogonal direction are obtained. These realizations are ensemble-averaged to recover the SS fields that are added to the LS terms in the LS governing equations given by Eqs. (2.4) and (2.6).

To solve the SS equations along each of the 1D lines in Fig. 2 requires assumptions that primarily impact how the derivatives in the cross-direction (to the 1D direction) are modeled. The rationale and the justification used for the velocity field were discussed earlier and are not repeated here, for brevity.^{21} In this study, we focus primarily on the assumptions used to obtain the modeled 1D SS scalar equations. These assumptions are as follows:

- the SS second-order derivative along the 1D line
*l*for the SS scalar field $ \varphi S$ is modeled as the averaged sum of the SS second-order derivatives along all three orthogonal directions through_{k}$ \u2202 2 \varphi S \u2202 x k \u2202 x k = 1 3 \u2211 j = 1 3 \u2202 2 \varphi S \u2202 x j 2 .$Here, “

*k*” is a free index and refers to the line*l*, which is parallel to the corresponding coordinate $ x k \u2009 ( k = 1 , \u2009 2 \u2009 and \u2009 3 )$. For example, along line_{k}*l*_{1}, SS second order derivatives along*x*_{2}and*x*_{3}directions are not available and therefore are lumped together with the SS second-order derivative along the*x*_{1}direction. A similar approach is used for the equations along lines*l*_{2}and*l*_{3}. - the SS contribution to the advection term from the SS scalar field is neglected in the transverse direction $ ( j \u2260 k )$ to the line
*l*,_{k}$ \u2202 \u2202 x j [ ( u j S + u j L ) ( \varphi S + \varphi L ) ] S = \u2202 \u2202 x j [ ( u j S ( l k ) + u j L ) ( \varphi S ( l k ) + \varphi L ) ] S .$

Although this seems a drastic assumption, it is noted earlier^{21} that the advection terms neglected on any line *l _{k}* are available (albeit, in an approximate sense) along the other two orthogonal lines.

The above approach makes some assumptions regarding the SS scalar field that will be revisited later when an extension to multiscalar mixing and reacting flows are considered.

With these assumptions, the 1D SS scalar equation along line *l _{k}* is given by

with the corresponding forcing term expressed as

### E. Hybrid TLS–LES method

To assess the TLS model for scalar mixing in an actual simulation, the hybrid TLS–LES strategy is employed for wall-bounded flow, as shown in Fig. 2(b). Here, the TLS model is applied only in the near-wall region and away from the wall, and a conventional LES with a SGS dynamic eddy diffusivity model^{26} is used. The transition between the TLS and the LES region is achieved by a smooth blending function $ K \u2261 K ( x , t ) \u2208 [ 0 , \u2009 1 ]$, which can be either statically^{24} or dynamically adjusted.^{25} Here, *K* = 0 or 1 corresponds to regions where only LES or TLS models are used, whereas for $ K \u2208 ( 0 , \u2009 1 )$, contributions from both LES and TLS models are used. The TLS–LES method is primarily designed for computationally efficient simulations since the use of the TLS model in the entire computational domain can be expensive. It has been validated for a wide range of flows, which include flow in a periodic channel, separating/reattaching flow, wakes, and mixing layer. The full description of TLS–LES implementation is presented in the cited references, and here, the same approach is employed for scalars and offers (albeit) a brief *a posteriori* assessment of the TLS model.

## III. COMPUTATIONAL SETUP AND APPROACH

The computational setup, the simulation parameters for DNS, and the numerical methodology employed to solve the equations are discussed in this section.

### A. Flow configuration

Passive scalar transport in a fully developed turbulent channel flow at a fixed friction Reynolds number $ R e \tau \u2261 u \tau h / \nu = 395$ is considered in this study (The other higher values of $ R e \tau $ will be considered later). Here, *h* is the channel half-height and $ u \tau $ is the wall-shear or friction velocity. The extent of the computational domain is $ 2 \pi h$, 2*h*, and *πh* in the streamwise (*x*), wall-normal (*y*) and spanwise (*z*) directions, respectively. Periodic boundary conditions in the streamwise and spanwise directions are used along with no-slip bottom and top walls with Dirichlet boundary conditions for the scalar field. To examine the effects of a moderate increase in the molecular diffusivity, two Schmidt numbers, *Sc* = 1 and 4, are studied here. As described in Sec. I, DNS tends to be computationally expensive for higher values of *Sc* at a fixed value of $ R e \tau $. Therefore, to allow for both *a priori* and *a posteriori* assessment using the DNS results, we have considered only a moderate increase in the value of *Sc*. Hereafter, $ x \u2261 ( x 1 , x 2 , x 3 ) \u2261 ( x , y , z )$ is used interchangeably.

For the *a priori* assessment, two DNS cases, denoted by A1-D and A2-D, are simulated with *Sc* = 1 and *Sc* = 4, respectively. The computational grid is uniformly spaced in the streamwise and spanwise directions, and in the wall-normal direction, a nominal stretching of the grid is performed to cluster more points in the near-wall regions. The grid resolution is chosen to ensure that the scalar field is resolved as per the value of *Sc*. Therefore, the total number of grid points used for case A2-D ( $ 512 \xd7 384 \xd7 512$) is 8 times the number of points used by case A1-D ( $ 256 \xd7 192 \xd7 256$). In terms of wall units, the grid resolution is $ \Delta x 1 + \u2248 5$ and $ \Delta x 3 + \u2248 2.5$ along the streamwise and spanwise directions for case A2-D. In the wall-normal direction, the minimum and maximum resolutions are about 0.1 and 2.7 wall units for this case with more than 25 points within $ x 2 + = 10$.

For the TLS–LES studies (labeled later as A1-T and A2-T for *Sc* = 1 and *Sc* = 4, respectively), LS grid resolution ( $ 64 \xd7 46 \xd7 64$) is used as in an earlier study.^{25} In terms of the wall units, the LS grid resolution is about 40 and 20 wall units in the homogeneous streamwise and spanwise directions, respectively. In the wall-normal direction, the minimum and maximum LS grid resolutions are about 5 and 30 units, which is coarser in comparison to a wall-resolved LES (typically 2 units).^{27} The SS grid resolution is set to 10 SS cells per LES cell along each direction for both the cases and is chosen to ensure that the scalar field is spatially resolved in a manner comparable to a DNS. The blending of the TLS and LES models is performed by a wall-normal coordinate based smooth and static blending function $ K ( x 2 )$,^{25} where the TLS model remains active up to $ x 2 + \u2248 80$ measured from the walls. The blending function is given by

Here, $ x 2 d = | x 2 \u2212 x 2 , w |$ with $ x 2 , w$ denoting the location of the wall in the vertical direction and $ d = 0.1 h , \u2009 c 1 = 2$, and $ c 2 = 0.2$. The same blending function is utilized for both the flow and the scalar fields. At each LS time step, the SS flow and scalar fields are evolved using a sub-cycling strategy, starting with a zero value for the SS fields. Further details of this strategy for the flow field are presented elsewhere,^{21,23} and the same strategy is followed for the evolution of SS scalar equations. To further assess the performance of the TLS model, two other modeling approaches are considered using the same LS grid. These include a no-model approach and a conventional eddy viscosity/diffusivity-based model.^{26,28} The no-model approach ignores the SGS contributions to the SGS stress tensor and the scalar flux terms. In the conventional approach, referred to as the locally dynamic kinetic energy model (LDKM), the eddy viscosity is obtained by solving an additional modeled transport equation for the SGS turbulent kinetic energy^{28} and the eddy diffusivity is determined using a dynamic approach.^{26} These additional cases are labeled as A1-L (A1-N) and A2-L (A2-N) for cases with *Sc* = 1 and *Sc* = 4, respectively, using the LDKM (no-model) approach.

The flow in the channel is driven by specifying a constant pressure gradient so that in the statistically stationary state, a pre-specified value of $ R e \tau $ is approximately attained. The flow is initialized by imposing three-dimensional perturbations on the mean velocity field and specifying a mean profile for the scalar field. The approach of the flow and the scalar fields to a statistically stationary state is assessed in terms of the time evolution of global quantities such as the coefficient of friction ( $ C f$), and Nusselt number (*Nu*), which are defined as

Here, $ \tau w$ and $ Q w$ denote the wall shear stress and wall scalar flux, respectively, $ U b$ is the bulk velocity, and $ \Delta \varphi $ denotes the difference in the value of the scalar field between the bottom and top walls. After the flow evolves to a statistically stationary state, the turbulence statistics are obtained by averaging along the homogeneous *x*_{1} and *x*_{3} directions and then obtaining a running average over time.

### B. Numerical methodology

For DNS, the governing equations for the flow and the scalar fields given by Eq. (2.1) are solved using a parallel 3D incompressible flow solver using the artificial compressibility method.^{28} The flow solver has been recently extended to allow for DNS and LES of passive and active scalar transport.^{29} The spatial discretization utilizes a fourth-order-accurate finite-difference method on a generalized curvilinear nonstaggered structured grid. The semidiscrete equations obtained after the spatial discretization are advanced in time using a second-order accurate backward-difference method. The artificial compressibility method utilizes a dual time-stepping algorithm to perform pseudo-time-integration of the transport equations for the pressure, momentum, and scalar fields. This is achieved by using a five-stage Runge–Kutta time-stepping scheme. The numerical method is energy conserving and has been extensively used in the past to study a wide range of turbulent flows, such as isotropic turbulence, channel flow, separating/reattaching flows, wakes, mixing layers, and stratified flows. The parallelization of the solver utilizes the standard domain decomposition procedure for which the message passing interface (MPI) library is used.

In the TLS and TLS–LES approaches, the SS equations for the flow and the scalar fields are evolved during each LS time step. Since the SS fields can be intermittent and have sharp gradients, the small-scale equations are integrated along 1D lines in time using an explicit, two-step component-wise total-variation-diminishing scheme.^{30} Furthermore, since the SS equations can be integrated independently between two LS time-steps at all locations, a hybrid programing paradigm that leverages features of the distributed (MPI) and shared (OpenMP) programing paradigms is used for cost-effective simulations.^{25}

Since the TLS requires the solution of the 1D SS equations, the computational cost of TLS is higher than that of a conventional LES using the same LS grid. For example, in the present case, the cost of the TLS–LES case is about 10 times the cost of the LDKM-based LES. However, the hybrid parallelization algorithm allows the use of OpenMP threads for an efficient solution of the independent 1D SS equations. In particular, with the use of 6 threads, the computational cost reduces from 10 times to about 2 times the cost of the LDKM-based LES case. Note that a corresponding wall-resolved LES for $ R e \tau = 395$ and *Sc* = 1 requires a grid of $ 128 \xd7 97 \xd7 128$. By using 2 OpenMP threads, the cost of the present TLS–LES case reduces to about 80% of the cost of such a wall-resolved LES.

## IV. RESULTS AND DISCUSSION

*A priori* assessment of the TLS model for the scalar field is first carried out using the DNS field in the physical and spectral space. Then, the implications of the modeling assumptions are discussed with a focus on the SS scalar physics predicted by the TLS model. Finally, *a posteriori* results of the near-wall TLS–LES strategy are discussed for the same conditions used in the *a priori* studies.

The *a priori* analysis along 1D lines *l*_{1} and *l*_{2} along the streamwise (*x*_{1}) and wall-normal (*x*_{2}) directions is only one discussed since the results along the spanwise (*x*_{3}) line *l*_{3} are similar to those along *l*_{1} and therefore not included for brevity. For analysis, 10 instantaneous snapshots segregated uniformly in time by 0.25 of the convective flow through time are used. While examining the results along *l*_{1}, four different wall-normal locations, namely, $ x 2 + =$ 10, 30, 100, and 395, are considered. These four locations approximately represent the regions of the flow corresponding to the “peak of the turbulence intensity,” the “buffer layer,” the “log layer,” and the “outer layer.” These four different regions are well known to exhibit different types of turbulence dynamics.^{31} For the *a priori* analysis, as stated before, the DNS data for a field are used to get the discrete LS field from which the continuous LS field is obtained using a cubic B-spline interpolation, followed by the SS field by using Eq. (2.2).

### A. Two-scale decomposition

The two-scale decomposition of the scalar field is illustrated in Fig. 3 for case A1-D where the LS and SS fields are obtained by using two different LS grid resolutions along the streamwise direction at $ x 2 + = 10$. With the LS^{4Δ} grid, the LS field adequately represents the spatial variation of the scalar field in comparison to the DNS field, as the magnitude of the SS field is significantly lower. However, the SS contributions become significant as the LS grid is coarsened. For example, the intensity of the SS fluctuations is about 3 times with the LS^{16Δ} grid compared to the LS^{4Δ} grid. While using the LS^{16Δ} grid, apart from capturing the small-scale variations, the SS field also shows some part of the coherent LS variations seen on the much finer LS grid. This ability of the SS field to capture (albeit, approximately) the content not resolved on the LS grid is characteristic of the two-scale approach (and other multiscale methods) and differs significantly from the classical eddy diffusivity type SGS model. The 1D SS fields are dictated by the LS dynamics through the term $ G L$ in Eq. (2.14), and as will be shown later, arbitrarily coarse LS grids cannot be employed.

The two-scale decomposition in the spectral space is shown in Fig. 4 where the 1D streamwise scalar spectra corresponding to the DNS, the LS, and the SS fields are shown at $ x 2 + = 10$ for the two cases using LS^{4Δ} and using two different LS grids for case A2-D. Similar to the behavior of the SS field in Fig. 3, the spectral contribution is increased [see Fig. 4(b)] as the LS grid is coarsened. In case A2-D (*Sc* = 4), the scalar spectrum shows a broader range of scales compared to the lower $ S c = 1$ case (A1-D), and across the entire range of wave numbers (length scales), a consistently higher value of the spectrum is observed, which is due to the enhanced level of scalar fluctuations.

The scale decomposition shows that the full spectral content of the DNS scalar field cannot be represented on the coarse-grained LS grid, but some of the missing content is included in the SS field. However, although the TLS model is less sensitive to the LS grid, an arbitrarily coarse LS grid cannot be used since the SS field is still obtained from a modeled equation.

### B. Assessment of TLS modeling assumptions

The SS modeling assumptions discussed in Sec. II D are now assessed along lines *l*_{1} and *l*_{2} as the results along line *l*_{3} exhibited similar behavior as along line *l*_{1} and are, therefore, not included for brevity.

#### 1. Modeling assumption for SS scalar diffusion

The first assumption (see Sec. II D) is the reduction of the 3D SS scalar diffusion term along a particular 1D line such that the difference between the SS second-order derivative of $\varphi $ in a particular 1D line *l _{k}* and the averaged sum of the SS second-order derivatives in all three directions is equal to zero, i.e.,

*V*= 0, where

_{k}Here, *k* is not a repeated index. A deviation from *V _{k}* = 0 implies the error associated with this assumption.

Figure 5 shows the probability density functions (PDFs) of *V _{k}* along

*l*

_{1}and

*l*

_{2}for both the cases using different LS grid resolutions. Along line

*l*

_{1}, the results are shown only in the near-wall region, i.e., $ x 2 + = 10$. At other locations, a similar shape of the PDF is observed. The PDFs from both the cases and at both the LS grid resolutions show qualitatively a similar behavior with a narrow peak centered around a value of zero with a wider tail, indicating that overall this assumption can be considered reasonable. However, the modeling assumption also implies that the dissipative influence of events with $ V k \u2260 0$ is excluded but as seen, such events have low probability. An increase in

*Sc*on the PDF of

*V*is noticeable by the increase in the decay rate of the PDF and the narrower peak, indicating an increase in the probability of events around

_{k}*V*= 0. These results suggest that this assumption may be better satisfied for higher

_{k}*Sc*, but further studies are warranted and will be considered in the future.

As discussed earlier for the momentum equation by Kemenov and Menon,^{21} the events corresponding to *V _{k}* = 0 can be segregated into two different groups for the scalar field as well. The first group corresponds to regions where the SS second-order derivatives of the scalar field are relatively small, which in turn leads to smaller values of

*V*. Such regions exhibit weaker fluctuations of the scalar field and therefore can be adequately represented by the resolved LS field. On the other hand, the second group corresponds to the flow regions where the SS second-order derivatives are nonnegligible and are approximately equal to each other characterized by intense turbulent fluctuations of the scalar field, high scalar dissipation, and a significant effect of the SS dynamics on the LS dynamics. The deviation from these regions would, then, suggest inaccuracy associated with this assumption.

_{k}To further examine SS second-order derivative terms, we consider the joint PDF of the modeled SS second-order derivative $ M v , k = 1 3 \u2211 j = 1 3 \u2202 2 \varphi S / \u2202 x j 2$ and the exact (directly obtained from the DNS) SS second-order derivative $ E v , k = \u2202 2 \varphi S / \u2202 x k 2$ along line *l _{k}*. These PDFs are shown in Fig. 6 along

*l*

_{1}at $ x 2 + = 10$ and along

*l*

_{2}for both the cases. Along both the lines

*l*

_{1}and

*l*

_{2}, the contours of the PDFs attain a peak near the origin, showing the overall validity of the assumption. Since the bisector of the first and the third quadrant corresponds to the model assumption

*V*= 0, the events away from the bisector indicate $ V k \u2260 0$, which are excluded by this assumption. However, along the bisector and away from the origin indicates events where the modeled SS second order derivatives can be well represented.

_{k}Along line *l*_{1}, the events with large magnitude $ E v , k$ correlate very well with the event corresponding to a large magnitude of $ M v , k$ with the same sign. The correlation coefficient between these terms is given in Table I where at all the wall-normal locations, the correlation coefficient is greater than 70%. The effect of increasing *Sc* is evident in Fig. 6 by the narrowing of the contours toward the bisector and an increase in the correlation coefficient (see Table I). The shape of the joint PDF along line *l*_{2} shows much higher probability along the bisector and a much higher correlation coefficient shown in Table I. This can be attributed to a higher magnitude of $ \u2202 2 \varphi L / \u2202 x 2 2$ compared to $ \u2202 2 \varphi L / \u2202 x 1 2$ and $ \u2202 2 \varphi L / \u2202 x 3 2$, which leads to a corresponding higher magnitude of $ \u2202 2 \varphi S / \u2202 x 2 2$. Since *l*_{2} explicitly resolves the derivatives along the *x*_{2} direction, the correlation of modeled and exact terms is higher along this line.

. | . | l _{1}
. | . | . | . | . |
---|---|---|---|---|---|---|

Case . | LS grid . | $ x 2 + = 10 $ . | $ x 2 + = 30 $ . | $ x 2 + = 100 $ . | $ x 2 + = 395 $ . | l _{2}
. |

A1-D | LS^{4Δ} | 0.72 | 0.79 | 0.85 | 0.9 | 0.98 |

LS^{8Δ} | 0.74 | 0.76 | 0.81 | 0.85 | 0.98 | |

LS^{16Δ} | 0.72 | 0.72 | 0.77 | 0.8 | 0.98 | |

A2-D | LS^{4Δ} | 0.81 | 0.81 | 0.86 | 0.9 | 0.98 |

LS^{8Δ} | 0.79 | 0.77 | 0.81 | 0.85 | 0.99 | |

LS^{16Δ} | 0.74 | 0.74 | 0.76 | 0.8 | 0.98 |

. | . | l _{1}
. | . | . | . | . |
---|---|---|---|---|---|---|

Case . | LS grid . | $ x 2 + = 10 $ . | $ x 2 + = 30 $ . | $ x 2 + = 100 $ . | $ x 2 + = 395 $ . | l _{2}
. |

A1-D | LS^{4Δ} | 0.72 | 0.79 | 0.85 | 0.9 | 0.98 |

LS^{8Δ} | 0.74 | 0.76 | 0.81 | 0.85 | 0.98 | |

LS^{16Δ} | 0.72 | 0.72 | 0.77 | 0.8 | 0.98 | |

A2-D | LS^{4Δ} | 0.81 | 0.81 | 0.86 | 0.9 | 0.98 |

LS^{8Δ} | 0.79 | 0.77 | 0.81 | 0.85 | 0.99 | |

LS^{16Δ} | 0.74 | 0.74 | 0.76 | 0.8 | 0.98 |

Overall, the *a priori* analysis suggests that this assumption is reasonable for the scalar field as well. Further assessment at other values of $ R e \tau $ and *Sc* and *a posteriori* simulations are planned.

#### 2. Modeling assumption for SS scalar advection

The SS model [see Eq. (2.13)] along each of the 1D lines neglects the advection contribution from the other two orthogonal lines and can be reexpressed in a nonconservative form using the continuity equation as

In this assumption, only the SS advection term along line *l _{k}* is accounted as it is available along this direction.

Figure 7 shows the contours of the joint PDF of the exact and modeled SS advective terms along line *l*_{1} and *l*_{2}. Along line *l*_{1}, the results are shown at the near-wall location $ x 2 + = 10$. The modeled term corresponds to the second term on the right-hand-side of Eq. (4.2), which is only a part of the scalar advection term. We can observe that along both lines *l*_{1} and *l*_{2}, and in both the cases, higher probabilities are observed for events with small values of the SS advective terms. As the magnitude of the SS advection term increases, the probability decreases. These PDFs also suggest that the modeled SS advection term admits values higher or lower in magnitude compared to the most probable values of the exact SS advection term.

The effect of higher *Sc* is evident in Fig. 7, where now the deviation away from the bisector is higher for the higher *Sc* case and suggests that the range of fluctuations of the scalar field is larger in all directions. This is confirmed by the correlation coefficients shown in Table II. The correlation coefficient reduces on line *l*_{2} compared to line *l*_{1}, and an increase in *Sc* leads to a reduction in the correlation coefficient.

. | . | l _{1}
. | . | . | . | . |
---|---|---|---|---|---|---|

Case . | LS grid . | $ x 2 + = 10 $ . | $ x 2 + = 30 $ . | $ x 2 + = 100 $ . | $ x 2 + = 395 $ . | l _{2}
. |

A1-D | LS^{4Δ} | 0.61 | 0.45 | 0.43 | 0.51 | 0.12 |

LS^{8Δ} | 0.61 | 0.41 | 0.39 | 0.44 | 0.05 | |

LS^{16Δ} | 0.57 | 0.4 | 0.41 | 0.47 | −0.03 | |

A2-D | LS^{4Δ} | 0.32 | 0.24 | 0.23 | 0.2 | 0.07 |

LS^{8Δ} | 0.37 | 0.27 | 0.24 | 0.2 | 0.03 | |

LS^{16Δ} | 0.41 | 0.29 | 0.28 | 0.23 | 0.0 |

. | . | l _{1}
. | . | . | . | . |
---|---|---|---|---|---|---|

Case . | LS grid . | $ x 2 + = 10 $ . | $ x 2 + = 30 $ . | $ x 2 + = 100 $ . | $ x 2 + = 395 $ . | l _{2}
. |

A1-D | LS^{4Δ} | 0.61 | 0.45 | 0.43 | 0.51 | 0.12 |

LS^{8Δ} | 0.61 | 0.41 | 0.39 | 0.44 | 0.05 | |

LS^{16Δ} | 0.57 | 0.4 | 0.41 | 0.47 | −0.03 | |

A2-D | LS^{4Δ} | 0.32 | 0.24 | 0.23 | 0.2 | 0.07 |

LS^{8Δ} | 0.37 | 0.27 | 0.24 | 0.2 | 0.03 | |

LS^{16Δ} | 0.41 | 0.29 | 0.28 | 0.23 | 0.0 |

These observations suggest that in general, this assumption may not be completely adequate for the SS scalar model. However, since only a part of the total advection term in Eq. (4.2) is modeled along each 1D line, contributions from the other lines may compensate for the low correlation observed. For this paper, however, we continue with the 1D model as formulated focusing more on the predictions of well-known scalar quantities that are important for LES closures.

### C. Scalar dissipation rate

The scalar dissipation rate (*χ*) is a measure of the rate of the scalar mixing process,^{31} which smoothens the scalar gradient. In high *Re* turbulent flows, the velocity field fluctuates at multiple scales, leading to a complex gradient of the scalar field, which in turn enhances the mixing process. Apart from *χ*, another quantity of interest is the unresolved variance of the scalar field:^{32} [ $ ( \varphi \u2033 2 ) L = ( \varphi 2 ) L \u2212 ( \varphi L ) 2$]. In the TLS model, $ ( \varphi \u2033 2 ) L$ can be explicitly computed, but since *χ* involves scalar gradients, only part of the SS component of *χ* can be captured.

The scalar dissipation rate is given by

Figure 8(a) shows the PDF of *χ* at the four wall-normal locations. The PDF shape attains a narrow peak with a wide tail. The extent of the tail is reduced with the distance from the wall. At all the locations, the PDF of *χ* approaches a lognormal distribution, $ f = 1 \chi \sigma 2 \pi exp \u2009 [ \u2212 ( ln \u2009 \chi \u2212 \mu ) 2 2 \sigma 2 ]$. A reference lognormal distribution is included in Fig. 8(a), where the parameters *μ* and *σ* are determined using the average values of the mean and standard deviation of *χ* at the four locations. The approach to a lognormal distribution is well known and has been reported in past studies.^{33–35}

It is apparent that *χ* is significantly influenced by the level of turbulence intensity and exhibits a wide range of variations. Using the TLS decomposition, *χ* can be expressed as $ \chi = \chi L + \chi S$, and the 1D streamwise spectra are shown in Fig. 8(b). It is apparent that at both the LS grid resolutions, the LS and SS spectra together reasonably capture the behavior of the full DNS spectrum. In particular, with a decrease in the LS grid resolution to $ LS 16 \Delta $, the contribution of the SS spectrum *χ ^{S}* increases significantly compared to the contribution of the finer grid $ LS 4 \Delta $. However, although

*χ*can be captured, it is still somewhat approximate since the SS equations are solved on a 1D line with the underlying TLS modeling assumptions discussed in Sec. IV B. We can use the correlation of the exact

^{S}*χ*denoted by $ \chi S , DNS$ and the TLS modeled

^{S}*χ*denoted by $ \chi S , TLS$ to examine this feature. The DNS and TLS approximation of these quantities are given by

^{S}*k*is not a repeated index and corresponds to line

*l*.

_{k}It is apparent from Eq. (4.4) that $ \chi S , TLS$ can only represent *χ ^{S}* partly along the 1D line in the TLS model. However, as the TLS model employs three orthogonal 1D lines, it can still capture parts of

*χ*in the other two orthogonal directions. This is evident from the contours of joint PDF of $ \chi S , TLS$ and $ \chi S , DNS$ shown in Fig. 9 along lines

^{S}*l*

_{1}and

*l*

_{2}. Along line

*l*

_{1}at the near-wall location ( $ x 2 + = 10$), the contours of the joint PDF are spread around the diagonal. Qualitatively, the shape of the PDF is similar at both values of

*Sc*; however, the higher

*Sc*case shows a reduced level of spread around the diagonal. This spread occurs because in the present case, there is a large wall-normal gradient, i.e., $ \u2202 \varphi / \u2202 x 2 \u226b \u2202 \varphi / \u2202 x 1$; therefore, $ \chi S , TLS$ has significant error along

*l*

_{1}. However, the contours of the joint PDF are along the bisector in Fig. 9(b) for line

*l*

_{2}, indicating that along line

*l*

_{2}, significant parts of

*χ*are captured by the TLS model.

^{S}These results are further evident in the correlation coefficients shown in Table III. In particular, the correlation coefficient along line *l*_{1} is the lowest at $ x 2 + = 100$ in both the cases and decreases as the LS grid size is increased, which indicates higher levels of error in the modeling of *χ ^{S}*. On the other hand, the correlation coefficient is nearly unity along line

*l*

_{2}. These results further demonstrate that the use of three orthogonal 1D lines provides some in-built corrective features in the 1D SS model and may be important especially when large SS anisotropy is prevalent.

. | . | l _{1}
. | . | . | . | . |
---|---|---|---|---|---|---|

Case . | LS grid . | $ x 2 + = 10 $ . | $ x 2 + = 30 $ . | $ x 2 + = 100 $ . | $ x 2 + = 395 $ . | l _{2}
. |

A1-D | LS^{4Δ} | 0.66 | 0.65 | 0.59 | 0.72 | 0.99 |

LS^{8Δ} | 0.6 | 0.62 | 0.55 | 0.65 | 0.99 | |

LS^{16Δ} | 0.54 | 0.61 | 0.53 | 0.61 | 0.99 | |

A2-D | LS^{4Δ} | 0.78 | 0.75 | 0.62 | 0.76 | 0.99 |

LS^{8Δ} | 0.73 | 0.7 | 0.58 | 0.69 | 0.99 | |

LS^{16Δ} | 0.67 | 0.66 | 0.55 | 0.65 | 0.99 |

. | . | l _{1}
. | . | . | . | . |
---|---|---|---|---|---|---|

Case . | LS grid . | $ x 2 + = 10 $ . | $ x 2 + = 30 $ . | $ x 2 + = 100 $ . | $ x 2 + = 395 $ . | l _{2}
. |

A1-D | LS^{4Δ} | 0.66 | 0.65 | 0.59 | 0.72 | 0.99 |

LS^{8Δ} | 0.6 | 0.62 | 0.55 | 0.65 | 0.99 | |

LS^{16Δ} | 0.54 | 0.61 | 0.53 | 0.61 | 0.99 | |

A2-D | LS^{4Δ} | 0.78 | 0.75 | 0.62 | 0.76 | 0.99 |

LS^{8Δ} | 0.73 | 0.7 | 0.58 | 0.69 | 0.99 | |

LS^{16Δ} | 0.67 | 0.66 | 0.55 | 0.65 | 0.99 |

The ability of the TLS model to capture *χ ^{S}* is further examined in terms of its correlation with the estimated square of the LS gradient magnitude

*ζ*, which is estimated as $ \zeta L = [ ( \varphi \varphi ) L \u2212 \varphi L \varphi L ] / [ n \Delta ] 2$. Figure 10 shows the joint PDF of

^{L}*ζ*and

^{L}*χ*for the two cases along line

^{S}*l*

_{1}using two different LS grid resolutions. It is evident that even with a much finer LS grid ( $ LS 4 \Delta $), the correlation is much weaker with a significant spread away from the diagonal bisector. The correlation gets further weak when the LS grid is coarsened to LS

^{16Δ}, and the shape of the joint PDF is also altered, indicating that the scale-resolved

*ζ*cannot determine the SS dynamics. This is also evident quantitatively from the correlation coefficient (not shown here, for brevity). It appears that a better representation of the SS dissipation rate is required to capture the SS scalar dynamics in the LS field. Such a representation is available within the TLS model.

^{L}### D. Turbulent scalar transport and backscatter

The SS turbulent transport and backscatter are key from a modeling perspective as they determine the effects of the SS dynamics and the interscale coupling. Therefore, an assessment of these properties in the TLS model is carried out in this section.

#### 1. SS turbulent transport

In conventional LES, a widely popular approach of modeling the SGS scalar flux is to relate it to the resolved gradient of the scalar field by using an eddy diffusivity^{2} through $ \lambda j sgs = \u2212 \kappa t \u2202 \varphi \xaf \u2202 x j$, where *κ _{t}* is the eddy diffusivity. Here, the overbar denotes the spatial filtering operation used in an LES. The eddy diffusivity

*κ*can be modeled algebraically through

_{t}^{36}$ \kappa t = C \varphi \Delta \xaf 2 | S |$, where $ C \varphi $ is a model coefficient, $ \Delta \xaf$ is the LES filter width, and $ | S \xaf |$ is the magnitude of the resolved strain-rate tensor $ S \xaf i j = 1 2 ( \u2202 u \xaf i \u2202 x j + \u2202 u \xaf j \u2202 x i )$.

The model coefficient $ C \varphi $ can be specified to be a constant^{37,38} as in the Smagorinsky^{39} type scalar model or can be computed dynamically.^{36,40} The dynamic evaluation can lead to $ C \varphi < 0$, implying a support for counter-gradient SGS scalar transport. However, in actual simulations when $ C \varphi < 0$ occurs locally, it is typically set to 0 to avoid unphysical effects.^{36}

Since TLS model does not employ the notion of eddy diffusivity and since past studies have shown that it can account for both co- and counter-gradient SS turbulent momentum transport processes,^{23} we reexamine the SS turbulent transport of the scalar field $\varphi $ using the correlation of $ \lambda j L$ with respect to the LS gradient of the scalar field, i.e., $ \u2202 \varphi L / \u2202 x j$. The LS equation for the scalar field given by Eq. (2.6) can be reexpressed as

where $ \lambda j L$ is the corresponding SS scalar flux, a term analogous to $ \lambda j sgs$ used in the filtered LES equations. It is given by

Figure 11 shows the joint PDF of $ \lambda 2 L$ and $ \u2202 \varphi L / \u2202 x 2$ for both the cases. The results are shown along line *l*_{1} at two different wall-normal locations, $ x 2 + = 10$ and $ x 2 + = 395$, using a LS grid of size $ LS 16 \Delta $. It is evident that the SS scalar flux can be both co- and counter-gradients in nature (a similar behavior is also observed for the other LS grid sizes and at the other wall-normal locations and therefore not shown). The effect of *Sc* on the shape of the joint PDF is particularly noticeable in the near-wall region, as the resolved scalar gradient is larger for higher *Sc*, leading to an enhanced SS scalar transport.

The SS scalar transport is further assessed using two notional eddy diffusivities, which is defined as

Here, $ \kappa t F$ is based on the conventional SGS closure, whereas $ \kappa t L$ uses the model relation between SS scalar flux and the LS gradient of the scalar field. The model coefficient $ C \varphi $ is specified to be 0.0196 following a past study.^{38} It is important to note that the notion of the SS scalar flux in the TLS model is different from the SGS scalar flux in the eddy diffusivity type models.

Figure 12 shows PDFs of the two eddy diffusivities for the two cases at $ x 2 + = 10$ and $ x 2 + = 395$. As expected and according to Eq. (4.7), $ \kappa t F \u2265 0$, whereas $ \kappa t L$ can have positive or negative values. The shape of the PDF of $ \kappa t F$ approaches a lognormal distribution in both the cases and at both the locations, indicating the presence of a large variation of $ \kappa t F$. The effect of *Sc* is seen as a higher mean value of $ \kappa t F / \kappa $ in the higher *Sc* case at all the wall-normal locations (see Table IV), which is due to differences in the grids between the two cases.

Case . | LS grid . | $ x 2 + = 10 $ . | $ x 2 + = 30 $ . | $ x 2 + = 100 $ . | $ x 2 + = 395 $ . |
---|---|---|---|---|---|

A1-D | LS^{4Δ} | 18.9 | 7.5 | 3.6 | 1.2 |

LS^{8Δ} | 75.7 | 29.9 | 14.3 | 4.9 | |

LS^{16Δ} | 303.2 | 120.4 | 57.2 | 19.7 | |

A2-D | LS^{4Δ} | 19.6 | 7.9 | 3.7 | 1.4 |

LS^{8Δ} | 78.4 | 31.5 | 15.0 | 5.5 | |

LS^{16Δ} | 313.8 | 125.6 | 60.0 | 22.0 |

Case . | LS grid . | $ x 2 + = 10 $ . | $ x 2 + = 30 $ . | $ x 2 + = 100 $ . | $ x 2 + = 395 $ . |
---|---|---|---|---|---|

A1-D | LS^{4Δ} | 18.9 | 7.5 | 3.6 | 1.2 |

LS^{8Δ} | 75.7 | 29.9 | 14.3 | 4.9 | |

LS^{16Δ} | 303.2 | 120.4 | 57.2 | 19.7 | |

A2-D | LS^{4Δ} | 19.6 | 7.9 | 3.7 | 1.4 |

LS^{8Δ} | 78.4 | 31.5 | 15.0 | 5.5 | |

LS^{16Δ} | 313.8 | 125.6 | 60.0 | 22.0 |

The PDF of $ \kappa t L$ exhibits wider tails and narrower peaks in both the cases and at both the wall-normal locations. Unlike $ \kappa t F$, which shows a lognormal behavior, the shape of the PDF of $ \kappa t L$ can be represented with either the Tsallis distribution,^{41} or the stretched exponential distribution,^{42} which are typically used to represent turbulence statistics,^{43,44} or the behavior of the velocity gradient.^{42,45} In particular, the peak is much narrower in the near-wall region. The effects of *Sc* are evident in terms of a marginally higher value of the probability of events corresponding to a larger magnitude of $ \kappa t L$. The results show that the TLS model supports both co- and counter-gradient types of SS turbulent scalar transport (which are observed in applications such as stratified flows,^{46,47} turbulent heat transfer,^{48} and chemically reacting flows^{49–51}). Even though the PDF shows the possibility of $ \kappa t L$ having positive or negative values locally, the mean value of $ \kappa t L / \kappa $ summarized in Table V tends to be mostly positive. In comparison to the mean value of $ \kappa t F$, the magnitude of the mean value of $ \kappa t L$ is always smaller, and the effect of *Sc* is much more apparent. This is because unlike $ \kappa t L , \u2009 \kappa t F$ is independent of the scalar gradient [see Eq. (4.7)].

Case . | LS grid . | $ x 2 + = 10 $ . | $ x 2 + = 30 $ . | $ x 2 + = 100 $ . | $ x 2 + = 395 $ . |
---|---|---|---|---|---|

A1-D | LS^{4Δ} | 1.1 | 17.1 | 25.4 | 17.2 |

LS^{8Δ} | 1.3 | 4.4 | 16.3 | 19.7 | |

LS^{16Δ} | 2.0 | −9.7 | 15.1 | −0.9 | |

A2-D | LS^{4Δ} | 5.0 | 23.0 | 24.7 | 13.5 |

LS^{8Δ} | 11.3 | 11.5 | 15.3 | 24.0 | |

LS^{16Δ} | 7.2 | −23.1 | 4.8 | 30.7 |

Case . | LS grid . | $ x 2 + = 10 $ . | $ x 2 + = 30 $ . | $ x 2 + = 100 $ . | $ x 2 + = 395 $ . |
---|---|---|---|---|---|

A1-D | LS^{4Δ} | 1.1 | 17.1 | 25.4 | 17.2 |

LS^{8Δ} | 1.3 | 4.4 | 16.3 | 19.7 | |

LS^{16Δ} | 2.0 | −9.7 | 15.1 | −0.9 | |

A2-D | LS^{4Δ} | 5.0 | 23.0 | 24.7 | 13.5 |

LS^{8Δ} | 11.3 | 11.5 | 15.3 | 24.0 | |

LS^{16Δ} | 7.2 | −23.1 | 4.8 | 30.7 |

#### 2. SS backscatter

The SS scalar dissipation rate ( $ \u03f5 \varphi L = \u2212 \lambda j L \u2202 \varphi L \u2202 x j$) is an important quantity for the scalar transport as it affects the turbulence-scalar interaction and the scalar dissipation processes.^{52} It appears in the transport equation of the variance of the SS scalar field $ ( \varphi \u2033 2 ) L$ given by

If the SS scalar flux is modeled as in a conventional LES using the eddy diffusivity based approach, then $ \u03f5 \varphi L = \kappa t \u2202 \varphi L \u2202 x j \u2202 \varphi L \u2202 x j \u2265 0$, as the eddy diffusivity $ \kappa t \u2265 0$. A positive value of $ \u03f5 \varphi L$ physically means that this term acts as a sink in Eq. (4.8), which dissipates the scalar fluctuations, and acts as a forward scatter. However, as shown in Sec. IV D 1, both co- and counter-gradient SS turbulent transport processes of the scalar field are possible. To further understand the characteristics of the interscale interactions, we examine the behavior of $ \u03f5 \varphi L$.

Figure 13 shows the PDF of $ \u03f5 \varphi L$ for the two cases obtained using two different LS grid resolutions along line *l*_{1} at the near-wall location ( $ x 2 + = 10$). The shape of the PDF in both the cases and using both the LS grid resolutions is positively skewed with a narrow peak near $ \u03f5 \varphi L = 0$ and a wider tail corresponding to positive and negative values indicating the presence (with low probability) of large magnitudes of positive and negative values. The increase in *Sc* leads to a narrower PDF with a decrease in the probability of backscatter (negative) events compared to the forward (positive) events. Furthermore, even though there is a considerable presence of backscatter in both the cases, the mean value of $ \u03f5 \varphi L > 0$, implying a dominance of the forward cascade process in the statistical sense. As the LS grid is coarsened, the tails of the PDF get wider, implying an increase in the probability of the events corresponding to large-magnitude positive (forward) and negative values (backscatter) of the SS scalar dissipation.

The aspects of the interscale transfer mechanism for $ \u03f5 \varphi L$ are examined further in the spectral space. Figure 14 shows the streamwise spectra of $ \u03f5 \varphi L$ and $ \u03f5 \u2212 L$. Here, $ \u03f5 \u2212 L = ( \u03f5 \varphi L \u2212 | \u03f5 \varphi L | ) / 2$ denotes the backscatter part of the scalar dissipation. It can be observed that $ \u03f5 \u2212 L$ has contribution across all the length scales. A higher level of dissipation is evident as expected across the entire range of scales in the higher *Sc* case. The spikes in the spectra correspond to the locations on the LS grid, where $ \u03f5 \varphi L$ is identically equal to zero because of the employed LS function. Due to this, $ | \u03f5 \varphi L |$ has a discontinuity in its derivative, causing the observed spikes.

To further understand the characteristics of $ \u03f5 \u2212 L$ from the modeling perspective, its correlation with the magnitude of the resolved LS gradient of the scalar field ( $ | \u2207 \varphi L |$) is examined. Since the eddy diffusivity based closure for the unresolved scalar flux utilizes $ \u2207 \varphi L$, the aforementioned correlation can indicate if the existing LES closures can be improved to allow for the inclusion of backscatter.

Figure 15 shows contours of the joint PDF of $ \u03f5 \u2212 L$ and $ | \u2207 \varphi L |$ for both the cases at the inner ( $ x 2 + = 10$) and outer ( $ x 2 + = 395$) locations using LS^{16Δ} as the LS grid size. (The qualitative behavior of the joint PDF for other LS grid resolutions and at the two other wall-normal locations ( $ x 2 + = 30$ and 100) is the same.) The shape of the joint PDF and the value of the correlation coefficient are nearly independent of *Sc* and the LS grid resolution. In particular, the correlation is enhanced in the outer region, which is also evident from the change in the shape of the PDF at $ x 2 + = 395$ compared to $ x 2 + = 10$. Overall, the correlation coefficient saturates to about –0.4 in both cases (not shown, for brevity) and at all the wall-normal locations except the near-wall region where the magnitude of the correlation is reduced to around –0.25. Such a behavior has also been observed for the variation of the backscatter part of the turbulent kinetic energy with respect to the resolved strain-rate magnitude in the past studies.^{23,53}

The *a priori* analysis presented here indicates that the TLS model extended to the scalar field has the requisite properties to capture some of the key SS features considered important for scalar mixing. Of course, the true assessment is when this model is used in an actual (*a posteriori*) simulation, and in this paper, we provide a demonstration of the TLS model in a TLS–LES study of these DNS cases. More detailed studies are under way and will be reported in the future.

### E. *A posteriori* demonstration using the TLS–LES approach

The TLS model for the scalar field is used as a near-wall model within the hybrid TLS–LES strategy to simulate the same two DNS cases but at a much lower LS resolution ( $ 64 \xd7 46 \xd7 64$). The SS grid resolution is set to 10 SS cells per LES cell along each direction for both the cases, which is adequate for spatial resolution of the scalar field in a manner comparable to a DNS. These cases are labeled as A1-T and A2-T, and the details of the numerical setup are described earlier in Sec. III A. Furthermore, to demonstrate the effectiveness of the TLS model as a near-wall model for the flow and the scalar fields, two other modeling strategies for LES, namely, the LDKM^{28} and the no-model (NM), are used to perform additional simulations. These cases are labeled as A1-L (A1-N) and A2-L (A2-N) for *Sc* = 1 and 4, respectively, for the LDKM (NM) based LES. Note that the same LS grid is used in these LES cases with the near-wall grid resolution of $ x 2 + \u2248 5$, which is coarser in comparison to a wall-resolved LES, where typically $ x 2 + < 2$ is required to capture the near-wall dynamics.^{27}

The ability of the TLS model is first demonstrated in terms of the prediction of the wall relevant quantities such as the friction coefficient $ C f$ and the Nusselt number *Nu* [see Eq. (3.2)], which are summarized in Table VI. Overall, a good agreement of the TLS–LES results with the corresponding DNS results is observed. It is further observed that the TLS–LES results are better compared to the results from the LDKM and NM approaches. This is expected because the near-wall grid resolution is coarser in these additional LES cases.

Case . | Sc
. | $ R e \tau $ . | $ C f \xd7 10 3 $ . | Nu
. |
---|---|---|---|---|

A1-D | 1 | 398 | 6.3 | 15.6 |

A2-D | 4 | 403 | 6.4 | 34.5 |

A1-T | 1 | 395 | 6.4 | 16.7 |

A2-T | 4 | 393 | 6.2 | 36.4 |

A1-L | 1 | 433 | 7.7 | 18.5 |

A2-L | 4 | 437 | 7.7 | 43.5 |

sA1-N | 1 | 403 | 8.0 | 17.5 |

A2-N | 4 | 400 | 8.0 | 38.8 |

Case . | Sc
. | $ R e \tau $ . | $ C f \xd7 10 3 $ . | Nu
. |
---|---|---|---|---|

A1-D | 1 | 398 | 6.3 | 15.6 |

A2-D | 4 | 403 | 6.4 | 34.5 |

A1-T | 1 | 395 | 6.4 | 16.7 |

A2-T | 4 | 393 | 6.2 | 36.4 |

A1-L | 1 | 433 | 7.7 | 18.5 |

A2-L | 4 | 437 | 7.7 | 43.5 |

sA1-N | 1 | 403 | 8.0 | 17.5 |

A2-N | 4 | 400 | 8.0 | 38.8 |

The effectiveness of the TLS–LES modeling strategy is further evident from the vertical profile of the normalized mean streamwise velocity ( $ u 1 +$) and the scalar ( $ \varphi +$) fields shown in Fig. 16, where the results from all the models for *Sc* = 1 are compared with the reference results. In particular, the NM case underpredicts both $ u 1 +$ and $ \varphi +$, whereas the LDKM case underpredicts $ u 1 +$ compared to the DNS case. This reflects in the prediction of $ C f$ and *Nu* shown in Table VI. As stated in Sec. I, a detailed comparison of different modeling strategies for the SGS scalar flux is out of the scope of the present study. Therefore, we limit the analysis of the results to assess the capabilities of the TLS model in comparison to the reference DNS and proposed correlations in the past studies,^{55–58} which is discussed next.

The vertical profiles of $ u 1 +$ and $ \varphi +$ obtained from the TLS–LES cases for both the values of *Sc* are shown in Fig. 17. The results are compared with the corresponding DNS cases and the reference numerical^{54} and analytical expression.^{57} As expected, due to the passive nature of the scalar, the increase in *Sc* does not affect the variation of $ u 1 +$, and an excellent agreement is observed with respect to the reference DNS results. The profiles of $ \varphi +$ on the other hand clearly show the effects of an increase in *Sc*, which is again captured very well by the TLS–LES approach in good agreement with the reference results. In particular, the increase in *Sc* leads to a sharp gradient in the scalar field in the near-wall region, which is also evident from a significant increase in the value of *Nu* shown in Table VI. The effect of *Sc* on the profile of $ \varphi +$ is similar to the past studies of passive scalar transport in wall-bounded turbulent flows.^{59,60}

The vertical profiles of the second-order turbulence statistics corresponding to the flow and the scalar fields are compared in Fig. 18. These profiles correspond to the streamwise turbulence intensity ( $ u 1 \u2032 +$), the vertical turbulent momentum flux [ $ \u2212 ( u 1 \u2032 u 2 \u2032 ) +$], the intensity of the scalar fluctuations ( $ \varphi \u2032 +$), and the vertical turbulent scalar flux [ $ \u2212 ( u 2 \u2032 \varphi \u2032 ) +$] fields. The variation of $ u 1 \u2032 +$ obtained from the TLS–LES cases shows a very good agreement with the reference results, whereas the values of $ \u2212 ( u 1 \u2032 u 2 \u2032 ) +$ are underpredicted. However, the peak locations of both these quantities and the trends in vertical variation are overall in good agreement with the reference results. As expected, there is no effect of *Sc* on these quantities due to the passive nature of the scalars.

The effects of an increase in *Sc* are apparent in the variation of $ \varphi \u2032 +$ and $ \u2212 ( u 2 \u2032 \varphi \u2032 ) +$. In particular, the peak value of $ \varphi \u2032 +$ increases with an increase in *Sc*, and the peak location moves closer to the wall. The vertical variation of $ \varphi \u2032 +$ from the TLS–LES agrees very well with the DNS results for *Sc* = 1. However, at *Sc* = 4, while the peak location shows a very good agreement, an overprediction is evident in the peak value. Furthermore, away from the wall where conventional LES is used, the values of $ \varphi \u2032 +$ are underpredicted. The effects of an increase in *Sc* are also noticeable in the variation of $ \u2212 ( u 2 \u2032 \varphi \u2032 ) +$, where the values are consistently higher in the *Sc* = 4 case compared to the *Sc* = 1 case for $ x 2 + > 3$. In addition, $ \u2212 ( u 2 \u2032 \varphi \u2032 ) +$ approaches unity in regions away from the wall, as it is normalized by the wall scalar flux, which is the sum of turbulent and molecular fluxes. Note that the contribution of the molecular fluxes to the total flux tends to zero in the core of the flow.^{60} The layer in which there is a significant diffusive contribution to the total scalar flux, i.e., the diffusive sublayer, gets thinner with an increase in *Sc*. These trends are captured by both TLS–LES cases. The observed differences in the second-order statistics of the scalar field at higher *Sc* can be attributed to a very coarse LS grid used in the present study.

The ability of the TLS–LES strategy to capture the effects of *Sc* can also be examined in terms of the nondimensional transfer coefficient *K*^{+} of the scalar field defined as

Here, $ \u27e8 \xb7 \u27e9$ denotes averaging along the homogeneous *x*_{1} and *x*_{3} directions and in time. Figure 19 compares the variation of *K*^{+} obtained from DNS and TLS–LES cases with respect to the correlations reported in the past studies. These correlations are given by

The first two correlations^{55,56} are obtained from experimental data, and the last correlation is from a past numerical study.^{58} We can observe that the effects of *Sc* on *K*^{+} are captured very well by the TLS–LES approach.

The results presented in this section show that the TLS–LES strategy where the TLS is used as a near-wall model is able to capture the key statistics of the passive scalar transport in the turbulent channel flow.

## V. CONCLUSIONS

In this study, a two-scale simulation modeling approach is assessed for its ability to model passive scalar transport and mixing in a wall-bounded fully developed turbulent channel flow (at $ R e \tau = 395$) at two different *Sc* numbers (1 and 4). Several features of the TLS model are addressed in an *a priori* study by using DNS datasets. The capabilities of the TLS model are examined in terms of the physical and spectral space characteristics, the sensitivity of the predictions to the TLS modeling assumptions, and the ability to capture critical parameters of modeling interest such as the scalar dissipation rate, co-/counter-gradient scalar transport, and interscale interactions.

The scale-decomposition strategy employed by the TLS model shows that the large- and small-scale variations of the scalar field are appropriately represented in both the physical and spectral space, thus highlighting a key feature of the TLS model to approximately capture the variation of the scalar across the entire range of scales even when the LS field is (relative to DNS) coarse-grained. The TLS model also shows an appropriate response to a moderate increase in the value of the Schmidt number.

Further analysis showed that since the flow has a high level of anisotropy in the near wall-region and has a large gradient of the scalar field in the wall-normal direction, the SS scalar dissipation rate along the streamwise oriented 1D line is only partly captured, but its contribution on the line in the wall-normal direction is captured very well. Both *Sc* number cases show a substantial local presence of counter-gradient transport and backscatter, which are usually ignored by eddy diffusivity based LES closures.

The *a priori* assessment is revisited in a *a posteriori* study of the TLS model as a near-wall model within the hybrid TLS–LES strategy, and the overall results are still consistent given the rather coarse-grained LS simulation that is utilized. Still, further studies are warranted, including more *a posteriori* assessment of the scalar mixing modeling in much higher *Re* flows (as done earlier for a range of high *Re* problems focusing only on momentum mixing^{23,25}) that are not easily accessible to DNS. Finally, the extension of this approach to fuel-air mixing and combustion is also being assessed.

## ACKNOWLEDGMENTS

This work was supported by the Office of Naval Research (ONR) under Grant No. N00014-16-1-2577. The computation time provided by the DOD HPC center at the Navy DSRC and at the SimCenter of the University of Tennessee at Chattanooga is highly appreciated.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## References

*International Symposium on Turbulence and Shear Flow Phenomena (*