Numerical experiments on turbulent entrainment and mixing of scalars

Abstract Numerical experiments on the turbulent entrainment and mixing of scalars in a incompressible flow have been performed. These simulations are based on a scale decomposition of the velocity field, thus allowing the establishment from a dynamic point of view of the evolution of scalar fields under the separate action of large-scale coherent motions and small-scale fluctuations. The turbulent spectrum can be split into active and inactive flow structures. The large-scale engulfment phenomena actively prescribe the mixing velocity by amplifying inertial fluxes and by setting the area and the fluctuating geometry of the scalar interface. On the contrary, small-scale isotropic nibbling phenomena are essentially inactive in the mixing process. It is found that the inertial mechanisms initiate the process of entrainment at large scales to be finally processed by scalar diffusion at the molecular level. This last stage does not prescribe the amount of mixing but adapts itself to the conditions imposed by the coherent anisotropic motion at large scales. The present results may have strong repercussions for the theoretical approach to scalar mixing, as anticipated here by simple heuristic arguments which are shown able to reveal the rich dynamics of the process. Interesting repercussions are also envisaged for turbulence closures, in particular for large-eddy simulation approaches where only the large scales of the velocity field are resolved.

parallel-to-the-interface scales. These multi-scale spatially evolving cascade mechanisms have been finally rigorously explained in Cimarelli et al. (2021) by means of the theoretical framework provided by the generalized Kolmogorov equation.
One of the main limitations in the understanding of turbulence entrainment and mixing is the lack of dynamical information of the process, i.e. the cause and the effect of the elementary processes composing the flow system. Most of the studies mentioned above deal with the analysis of flow realization snapshots, thus losing possible information about the dynamic causality of the phenomena that naturally occur over time. In this context, it would be important to assess the previously shown multi-scale mechanisms of turbulent entrainment and mixing from a dynamic point of view. In the present work, we address this issue by analysing the evolution of the process under the separate action of large and small velocity scales. To this aim, numerical experiments are performed where the velocity coupling of the passive scalar equation with the Navier-Stokes equations is modified to alternatively suppress the role of large and small velocity scales in the process of scalar mixing. Hence, this approach allows us to clearly assess their contribution from a dynamical point of view. To the authors' knowledge this is the first attempt to address the dynamic causality of large and small velocity scales in turbulent mixing and entrainment of a scalar.
The paper is organized as follows. The numerical settings and the main flow features of the temporal planar jet are shown in § 2. Some heuristic arguments in support of the relevance of the different scales of turbulence in the scalar mixing process are reported in § 3. The numerical experiments are described in § 4 and their results are presented in § 5. In § 6, some insights into the essential processes at the basis of turbulent entrainment and mixing are envisaged. Finally, the work is closed by concluding remarks in § 7.

Numerical settings
The selected base flow for the numerical experiments on turbulent entrainment and mixing is a turbulent planar temporal jet (da Silva & Pereira 2008;van Reeuwijk & Holzner 2014;Cimarelli et al. 2021). This type of flow represents a paradigm for free shear flows by reproducing essential flow features such as an inhomogeneous mean shear, a turbulent/non-turbulent interface region and a simultaneous presence of two interacting well-separated ranges of scales: the large-scale motions reminiscent of the Kelvin-Helmholtz instability and the small-scale fluctuations of the developed turbulent motion. The main advantage of this type of flow with respect to other paradigms, e.g. spatially developing jets, mixing layers and wakes (Stanley, Sarkar & Mellado González 2002), is the two-dimensional statistical spatial homogeneity. Indeed, the flow develops in time rather than in the streamwise direction thus gaining a statistical spatial homogeneity and losing the statistical homogeneity in time.
The flow is governed by the Navier-Stokes equations coupled with the evolution of a passive scalar: Equations (2.1) are written in a dimensionless form by using the initial jet width H for the length scales, the initial jet velocity U 0 for the velocity scales and the initial jet scalar concentration C 0 for the scalar scales. When not specifically stated, all the results of the present work are reported dimensionless by using H, U 0 and C 0 as reference scales.
Here Re = U 0 H/ν is the Reynolds number, where ν is the kinematic viscosity, while Sc = ν/α is the Schmidt number, where α is the diffusivity of the scalar. In (2.1), index i = 1, 2, 3 corresponds to the streamwise, spanwise and cross-flow (u, v, w) velocities and (x, y, z) directions. Equations (2.1) are solved using the open-source code Incompact3d (Laizet & Lamballais 2009). The spatial discretization is based on a high-order compact finite-difference scheme, sixth-order accurate, and the time integration is performed using a third-order Runge-Kutta scheme. Periodic boundary conditions are applied in all directions.
The initial condition is a fluid layer that is quiescent and characterized by a null concentration of the scalar except for a thin region −H/2 < z < H/2 where the streamwise velocity and the scalar are non-zero and homogeneously distributed in the streamwise and spanwise directions: where σ 0 = H/35 is the initial momentum thickness. In order to facilitate a rapid transition to turbulence, a perturbation is superimposed to the initial condition of all three components of the velocity field consisting of a uniform random noise whose amplitude is 0.01U 0 . The simulations have been performed for Re = 3000 and Sc = 1. The numerical domain is a cuboid of size 36H × 36H × 24H discretized with 900 × 900 × 600 equally spaced points, thus leading to the same resolution in all three spatial directions. The worst condition for the resolution is reached at the centreline during the initial turbulent transition for t = 20 where x/η cl ≈ 4.5 with η cl the centreline value of the Kolmogorov scale. As shown in the following, in the self-similar decay of the flow, the Kolmogorov scale increases from this minimum both in time η cl ∝ √ t and with cross-flow position, thus highlighting that the spatial resolution employed is appropriate (see also the Appendix where a validation of the numerical solution is reported). Finally, a constant time step t = 0.01 has been used leading to a Courant-Friedrichs-Lewy (CFL) number that on average is CFL ≈ 0.3.
Averages, hereafter denoted as · , have been computed by taking advantage of the statistical properties of the flow, i.e. spatial homogeneity in the streamwise and spanwise directions and symmetry in the cross-flow direction, and by using an ensemble of two independent flow realizations obtained using different random perturbations in the initial conditions. The customary Reynolds decomposition of the flow in a mean and fluctuating field is adopted: where U = u and Θ = θ are the average fields and υ and ϑ are the fluctuating ones.

Flow features and self-similarity
The evolution of the flow is shown in figure 1 by means of iso-contours of the scalar field at different time instants. From the initial condition, the flow field exhibits first Kelvin-Helmholtz-like instabilities and then transitional mechanisms giving rise to turbulent fluctuations. As shown in figure 2(a), these production mechanisms of turbulent fluctuations exceed the dissipative ones giving rise to an increase of the turbulent kinetic energy and of the scalar variance content up to t ≈ 15 and t ≈ 20 where the maxima of ϑϑ cl (t) = ϑϑ (z = 0, t) and k cl (t) = υ i υ i (z = 0, t)/2 are reached, respectively. Hence, for t > 20, turbulence starts to decay. As shown by the instantaneous behaviour of the scalar field in figure 1, the flow, while decaying, performs entrainment and mixing producing increasingly larger scales of motion and increasing its width. Note that also at the latest stages of the decay, the large-scale structures of the main flow instability are still visible.
The decay of the flow is characterized by a self-similar behaviour. The only relevant parameters in the self-similar regime are the volume flux Q 0 ≡ U dz, the cross-sectional scalar content C 0 ≡ Θ dz and the time t. The former two are invariant for this flow type. The conserved quantities, Q 0 = const. and C 0 = const., can be easily derived from the cross-stream integration of the mean streamwise momentum equation and of the mean scalar equation, respectively. These three relevant parameters immediately imply that the self-similar decay is characterized by an increase of the length scales as √ Q 0 t and by a decrease of the velocity and scalar scales as √ Q 0 /t and C 0 / √ Q 0 t, respectively. In figure 2(b), we show the behaviour of the mean centreline velocity and scalar concentration, U cl (t) = U(z = 0, t) and Θ cl (t) = Θ(z = 0, t), and of the jet half-widths h Ω and h θ , defined as the centreline distance where the mean enstrophy and scalar profiles reduce to 2 % their centreline value. From the plot it is clear that the flow has reached a dynamic equilibrium for t > 30 when all the observables start to follow a self-similar behaviour. Note that the scalar and enstrophy interfaces measured by means of h Ω and h θ almost collapse. In accordance with these scaling of velocities and lengths, the self-similar regime is characterized by a constant Taylor Reynolds number, Re λ = √ 2k cl /3λ cl /ν = 60, where λ cl = √ 10νk cl / cl is the centreline value of the Taylor microscale and cl the centreline turbulent dissipation (see the inset of figure 2a).

A heuristic view of entrainment and mixing
The evolution of scalar fields is governed by inertial and diffusive transport phenomena. The mixing and entrainment properties of the scalar are, hence, the result of complex dynamical interactions between these two processes. Indeed, scalar entrainment is directly accomplished by diffusive mechanisms only, but, from a dynamical point of view, their strength significantly depends on the properties of the inertial scalar fluxes especially in turbulent flows. It is then reasonable to use information from both the inertial and diffusive phenomena to argue about the rate of scalar entrainment and mixing. Based on heuristic arguments, we try here give an estimate of these two processes by considering general scale-dependent properties of the flow motion. Indeed, in the case of Schmidt number close to unity, both mixing due to turbulence convection and scalar diffusion occur at similar length and time scales. These predictions are then used to produce alternative scalings for the rate of scalar entrainment and mixing.
The diffusive flux across iso-concentration surfaces can be generally estimated as where C is a characteristic scalar concentration and S is the area of the scalar interface indirectly induced by a turbulent motion with cutoff scale (Catrakis, Aguirre & Ruiz-Plancarte 2002). On the other hand, the turbulent scalar flux across iso-concentration surfaces is directly determined by the convective action of the turbulent field and can be estimated as where u is the characteristic velocity difference of a turbulent motion with cutoff scale . Let us recall that this reasoning applies for flows of unity Schmidt number where the velocity and scalar cutoff length scales are similar. Equations (3.1) and (3.2) highlight that the essential role of turbulent motion is to influence the convoluted scalar surfaces by prescribing the interface area S and its fluctuation u , thus influencing both the diffusive and inertial scalar fluxes. Based on the assumption that iso-surfaces in turbulence are fractal-like, it is possible to derive a useful relation for estimating the area of complex turbulence surfaces: where S 0 and L are some characteristic area and length and D is the fractal dimension of the surface. In this context, a fundamental result has been obtained by Sreenivasan, Ramshankar & Meneveau (1989). Indeed, by assuming that all fluxes across iso-surfaces are independent of Reynolds number, Sreenivasan et al. (1989) found that an estimate of the fractal dimension of turbulent interfaces is D = 7/3 in good agreement with experimental results. By inserting the estimate of the surface area (3.3) in the models for the diffusive and inertial scalar fluxes (3.1) and (3.2), we obtain the following relations: (3.5) It is useful now to introduce the relative speed of a scalar iso-concentration surface: where V b is the boundary entrainment rate and V e is the entrainment velocity (da Silva et al. 2014). By considering this entrainment parametrization, the scalar concentration flux can be globally estimated as It is now relevant to determine a scaling of V θ by using information from the inertial and diffusive fluxes. By equating this expression with the estimates (3.4) and (3.5), we can define two alternative scalings of the entrainment velocity: (3.9) Interestingly, these two relations do not depend on the characteristic scalar concentration, thus suggesting the relevance of studying the effects of the structure of the velocity field on the mixing properties of scalar fields as is done in the following sections.

Numerical experiments
As shown in the previous analysis, all the scales of the turbulent motion contribute to the wrinkling of the scalar interface; however, of general relevance is the understanding of the separated role played by the large (engulfment) and small (nibbling) scales. Indeed, it is reasonable to argue that the large-scale field is responsible for large-scale events of scalar transport and, hence, of deformation and stretching of the scalar interface. These large-scale phenomena are in general a result of motions related to the flow instability, such as the Kelvin-Helmholtz instability for the case of temporal jets. On the other hand, the small-scale motion is expected to be more universal at least for high Reynolds numbers. It is again reasonable to argue that small turbulent scales produce a less convoluted scalar interface but composed of smaller and smaller ripples. It is then relevant to establish how the diffusive and convective laws of the scalar interface (3.4) and (3.5) change under the separate deformation action of large-and small-scale turbulent fluctuations. In the present work, we aim at dynamically establishing the above mentioned scale dependency of turbulent mixing and entrainment by performing numerical experiments where the scalar field is separately driven by large-and small-scale motions. First, a filtering operator is defined and applied to the velocity field: and G is the kernel of a generic filter in space. The filtering formalism allows us to decompose the turbulent motion in the space of scales as This decomposition of the velocity field is used to simulate the different dynamics of two independent scalar fields,θ and θ , separately driven by the large and small scales of the turbulence motion, respectively. The system of equations solved is thus highlighting that the evolution of the scalar fieldθ is driven by a flow motion composed of large-scale fluctuations,θ =θ[ū]. On the contrary, the scalar field θ evolves under the action of a flow motion composed of small-scale fluctuations, θ = θ [U + u ].
As can be seen, for the scalar field θ the action of the mean flow has been explicitly added since u = 0, while, for the scalar fieldθ, the mean is implicitly included because ū = U. It is important to notice that the filtering operation is applied for the sole decomposition of the velocity field in large and small scales that separately drive the evolution of the scalar fieldsθ and θ , respectively. Hence, the two scalar fieldsθ and θ are the result of a different dynamical evolution imposed by the separate action of large- and small-scale motions and should not be confused with the filtering decomposition of the scalar field itself, i.e. θ / =θ + θ . The initial conditions are taken from the unfiltered solution at t = 40 thus allowing us to study the different evolution ofθ and θ in the fully turbulent self-similar conditions of the flow. A two-dimensional Gaussian filter is employed in the two homogeneous directions, with characteristic length Δ = 2.2λ cl . Hence, the filter length increases during the flow decay as √ t by following the behaviour of the Taylor microscale λ cl . By analysing other values of filter lengths, we found that the selected value Δ = 2.2λ cl allows us to nicely split the large-scale motions typical of the main instability of the flow from the more homogeneous and isotropic small-scale turbulent field (da Silva, Lopes & Raman 2015).
As shown in figure 3, the large-scale spanwise rolls reminiscent of the Kelvin-Helmholtz instability are carried by the large-scale fieldū. Superimposed to this large-scale motion, very anisotropic turbulent fluctuations are also observed in the field u. On the contrary, no signature of large-scale spanwise rolls is found in the small-scale field u . This scenario is confirmed by the streamwise velocity spectra E k (k x ) evaluated at the jet interface z = h θ and shown in figure 4(a). Indeed, the large-scale fieldū is found to capture the large-scale spectral behaviour including the peak at small wavenumbers, while the small-scale field u reproduces the small-scale spectral behaviour. In particular, the crossover scale between the two fields is found to take place well within the inertial range where the turbulent spectrum follows the classic k −5/3 power law. Hence, the more homogeneous and isotropic small-scale fluctuations of the flow find their support in the field u .
To better characterize the large-and small-scale fields, the probability density function of the velocity decomposition evaluated at the jet interface is shown in figure 4(b). The positively skewed distribution of the cross-flow velocity fluctuations at the interface is found to be almost entirely carried by the large-scale fieldū. On the contrary, the small-scale field u exhibits an almost symmetric distribution with a highly intermittent character. All these features are at the basis of the significantly different entrainment and mixing properties of the scalar field experiments that are discussed in the following section.
In closing this section, let us point out that other types of filter kernels, including also the filtering in the cross-flow direction, have been considered. The detailed study of the effects of the filter type, of the filter lengths and of the initial conditions on the evolution of the scalar fields described by the system of (4.4) is reported in Boga (2020). For the purpose of the present work, no physically relevant differences have been observed. In the same work, different filter lengths for the scalar fieldsθ and θ have also been considered to enforce a separation of scales that due to the relatively low Reynolds number considered is otherwise limited. Again, no physically relevant phenomena are observed.

Results
We start by considering the scalar field evolution under the separate action of the largeand small-scale velocity fields. As shown in figure 5, the evolution of the scalar fieldθ roughly resembles that of the unfiltered case θ in terms of both large-scale pattern and spreading of the jet. Indeed, the main difference between the two fields is given by the reduced amount of mixing, i.e. large-scale engulfment events are present but not combined with small-scale nibbling. As a result, at the final stages of the evolution for t = 140, the region occupied by the scalar fieldθ is almost the same as that for θ but with the difference that unmixed flow regions of low-concentration levels are present within it. A completely different scenario appears for the scalar field θ . In this case the large-scale pattern is completely different with respect to that of θ, with large-scale engulfment events being completely suppressed. The low-concentration regions initially present within the jet for t = 50 are absorbed during the evolution by the small-scale nibbling and diffusion, thus leading to a significantly more mixed turbulent jet at the final stage of the evolution for t = 140. The completely different evolution of the two scalar fieldsθ and θ is also found to produce very distinct scalar-interface topologies for t = 140. A less convoluted interface with small-scale wrinkling is produced by the scalar field θ in contrast to a highly convoluted surface with almost absent small-scale wrinkling for the scalar fieldθ .
We analyse now the evolution of the three scalar fields θ,θ and θ from a statistical point of view. As shown in figure 6, the small-scale-driven scalar field θ shows a slower decay of the mean centreline scalar concentration together with a more rapid decay of its variance. On the contrary, the large-scale-driven scalar fieldθ shows a behaviour that more closely matches that of the unfiltered field in terms of mean scalar concentration decay, i.e. only a slightly slower decay is observed (see figure 6a). This matching of the unfiltered behaviour is, however, not observed when considering the scalar concentration variance. As shown in figure 6(b), a significantly slower decay of the scalar fluctuations is indeed observed.   The origin of these behaviours can be studied by considering the scalar variance budgets for θ,θ and θ : where it is possible to recognize the turbulence production and dissipation terms and the turbulent and diffusive fluxes As shown in figure 7(a), the turbulence production of scalar fluctuations is largely reduced for the scalar field ϑ driven by the small-scale velocity field. This behaviour can be easily explained by considering that turbulence production is a large-scale anisotropic phenomenon and, hence, is largely inhibited when coupling a scalar field with an almost isotropic small-scale velocity field. Accordingly, the scalar fieldθ being driven by a very anisotropic large-scale velocity field exhibits a turbulence production that is almost unaltered with respect to that of ϑ. The behaviour of turbulence dissipation is shown in the inset of figure 7(a). At the initial stages of the decay, the dissipations of bothθ and ϑ show a reduction with respect to that of the scalar field ϑ. On the other hand, for t ≥ 70 only the dissipation of the scalar field ϑ shows a reduction while the dissipation ofθ is found to almost recover the dissipation values of ϑ. These behaviours of turbulence production and dissipation can be used to explain the previously observed decays of the mean and variance of the scalar concentration (see again figure 6). By noting that the turbulence production of scalar fluctuations is a phenomenon that subtracts intensity from the mean field and releases it to the fluctuations, it is immediately clear that for the scalar field θ , the large reduction in the production of scalar fluctuations leads to a significantly slower decay of the mean scalar concentration Θ (t). This reduction in the production of scalar fluctuations is not balanced by an equivalently large decrease of dissipation thus leading to a reduction of the net source of scalar fluctuations, i.e. π − θ < π − θ . As a consequence, a more rapid decay of the intensity of the scalar fluctuations ϑ ϑ (t) is observed. The opposite scenario occurs for the scalarθ driven by the anisotropic large-scale velocity field. In this case the production of scalar fluctuations is only marginally reduced, thus justifying the slightly slower decay of the mean scalar concentrationΘ(t). This slight reduction of production is combined with a larger reduction of dissipation especially at the first stages of the evolution, thus leading to an increase of the net source of scalar fluctuations, i.e.π −˜ θ > π − θ . This behaviour justifies the observed slower decay of the intensity of the scalar fluctuations θθ (t). The formalism of the scalar variance budget allows us also to address the inertial and diffusive phenomena of spatial spreading of the scalar fluctuations from the production to the dissipative regions of the flow, i.e. the turbulent and diffusive fluxes ϕ u and ϕ α . As shown in figure 7(b), the scalar fluctuations are transported from the production region towards the inner and outer regions of the jet, z/h θ/2 < 1 and z/h θ/2 > 1, respectively. This transport is essentially performed by inertial mechanisms for both ϑ andθ, i.e. the scalar fluctuations driven by the unfiltered field and by the large-scale velocity field, respectively. Indeed, the turbulent flux is always far greater than the diffusive one: ϕ u ϕ α and ϕ u φ α . The scenario is, however, different for the scalar fluctuations ϑ driven by the small-scale velocity field. Also in this case the fluctuations are mainly transported by inertial mechanisms, ϕ u > ϕ α , but their intensity is markedly reduced: ϕ u ϕ u . On the contrary, the diffusion mechanisms are found to surprisingly remain almost unaltered, i.e. ϕ α ≈ ϕ α . Accordingly, we argue that the reduction of the inertial transport mechanisms is not induced by a decrease of the scalar fluctuation intensity given by the contraction of the turbulence production but simply by a less efficient coupling of the small-scale velocity field with the scalar fluctuations, i.e. | ϑ ϑ υ 3 | | ϑϑ υ 3 |. As a result, the percentage contribution of diffusion with respect to advection in the spatial spreading of the scalar fluctuations ϑ is largely increased. As expected, the described different evolution of the mean and fluctuating scalar field for the three numerical experiments also has an impact on the global mixing and entrainment properties of the scalar jet. In figure 8(a), the entrained volume V, defined as the volume of fluid where the scalar concentration is θ > θ th with θ th = 0.02Θ cl , is shown as a function of time for the three numerical experiments. During the flow decay, the entrained volume is significantly reduced for the scalar driven by the small-scale velocity field, V < V. On the other hand, a smaller reduction of the entrained volume is observed for the scalar driven by the large-scale velocity field,Ṽ < V. This difference between the two scalar fieldsθ and θ is enhanced when considering the spreading of the scalar jet. As shown in the inset of 8(a), the position of the scalar interface is only slightly reduced for the scalar driven by the large-scale velocity field,h θ ≈ h θ , in contrast to the scalar driven by the small-scale velocity field where the position of the scalar interface advances very slowly, h θ h θ . The different entrainment and propagation of the scalarsθ and θ can be explained by recalling the behaviours of the turbulent and diffusive fluxes. Contrary to the scalarθ, the scalar driven by the small-scale velocity field θ is characterized by a significantly reduced  inertial flux, |ϕ u | |ϕ u |, thus explaining the very slow spreading of the jet, h θ h θ . We argue that, associated with this reduction of inertial fluxes, a reduced amount/intensity of entrapment events is experienced by the scalar θ . As a result, the presence and the size of null-concentration blobs entrapped within the scalar jet z < h θ are markedly reduced. At the same time the diffusion processes for the scalar θ remain almost unaltered, ϕ α ≈ ϕ α . Hence, we argue that the time scale of the diffusion processes is small enough to effectively complete the entrainment of the reduced amount of entrapped null-concentration blobs. This picture is confirmed by figure 5, where the instantaneous flow patterns of the scalar field θ show a significantly more mixed flow region within the scalar jet z < h θ . These behaviours are not observed for the scalar driven by the large-scale velocity fieldθ since in this case engulfment events are largely reproducedφ u ϕ u while retaining the same intensity of the diffusion processesφ α ≈ ϕ α . Indeed, the pattern of the scalarθ shown in figure 5 exhibits a number of unmixed flow regions within the scalar jet z <h θ . This different behaviour of entrainment for the two scalars explains why the reduction of the entrained volume for the scalar θ , V < V, is not as large as the reduction of the spreading of the jet, h θ h θ . The probability density function of the three scalars ϑ,θ and ϑ actually confirms the above mentioned entrainment and mixing properties. As shown in the inset of figure 8(b), the probability density function at the jet centreline highlights a symmetric distribution of the scalar fluctuations for ϑ and a very wide and asymmetric distribution for ϑ and ϑ. In particular, this asymmetry is essentially given by an increased probability of large negative scalar fluctuations that can be clearly associated with the presence of unmixed flow regions of low scalar concentration even at the jet centreline due to very-large-scale engulfment events. In this respect, it can be noticed that the probability density functions of ϑ andθ have a negative tail that extends up to ϑ/Θ = −1 corresponding to θ = 0 and, hence, to completely unmixed flow regions. This behaviour is not observed for ϑ , thus confirming that this field lacks engulfment events and, at the same time, exhibits a mixed flow in the scalar jet core.
As shown in figure 8(b), the scenario is completely different at the scalar interface, z = h θ , z =h θ and z = h θ . In this case, all three scalar fields ϑ,θ and ϑ show a peak of probability for ϑ/Θ = −1 corresponding again to θ = 0. Indeed, the presence of entrainment at the scalar interface is revealed by a wide positive tail of the probability distribution. Interestingly, for the scalars ϑ andθ this positive tail is very wide and flat, thus revealing that at the scalar interface a wide range of scalar fluctuation intensities is almost equally probable. We argue that this plateau is again due to large-scale engulfment mechanisms that intersect the mean scalar interface and contain a wide range of scalar values. On the contrary, the scalar ϑ driven by the small-scale velocity field exhibits a positive tail characterized by a decrease of the probability with the intensity of the scalar fluctuations. Indeed, in this case engulfment is absent and the small scales of nibbling can contain only a fraction of the entire range of scalar values.

Theoretical framework
The different mixing and entrainment properties of the three numerical experiments highlighted in the previous section can now be analysed in terms of heuristic arguments in order to reveal the essential features of the process. In § 3, we have shown how the complex dynamics of turbulent fluctuations has an impact on the entrainment and mixing rates of a scalar field by establishing the geometrical and fluctuating properties of the scalar surface through which a scalar concentration flux occurs following the well-established Fick's law of diffusion. We apply now these arguments to the three scalar concentration experiments here developed.
In the case of a temporal planar jet and considering the scalar iso-concentration surface θ = 0.02Θ cl , a proper choice of the characteristic scales for the diffusive and inertial scalar fluxes (3.4) and (3.5) is where h θ is the jet half-width defined by the mean scalar interface itself, while η θ int and u η int are the Batchelor scale and the Kolmogorov velocity evaluated at the scalar interface z = h θ . Without loss of generality we have chosen the scalar iso-concentration surface θ = 0.02Θ cl , but let us highlight that the same theoretical framework applies for all scalar surfaces. With these characteristic scales, the diffusive and inertial fluxes become and it is easy to show that the two fluxes equal each other: An alternative scaling is given by considering an outer velocity scale in the assumptions (6.1a-d), i.e. u = υ int , (6.4) where υ int is the standard deviation of the velocity fluctuations at z = h θ . In this case, the inertial flux differs from the diffusive one and can be written as (6.5) In the temporal jet, the entrainment hypothesis (3.7) can be written as where S 0 is the (x-y)-plane surface of the jet. Hence, the different estimates of the scalar surface velocity become thus leading to the following scaling: For the scalarθ driven by the large-scale velocity field, the above assumptions apply by simply replacing the characteristic scales, i.e.
thus leading to the following velocity scales: On the contrary, for the scalar θ driven by the small-scale velocity field, the choice of the characteristic scales needs more attention. Indeed, due to the absence of large-scale engulfment events, the largest characteristic length for the scalar surface is given by the filter length Δ (see the instantaneous scalar topology shown in figure 5). Accordingly, the characteristic scales are and the estimates of the scalar surface velocity become

Comparison with the numerical experiments
To verify these scalings, it is important to note that the entrainment velocity V e is null in the temporal jet, since d dt U dz = 0, (6.13) and, hence, the relative speed of scalar interface defined in (3.6) is entirely determined by the boundary entrainment rate V θ = V b . The boundary entrainment velocity can be computed as the rate of variation of the volume of fluid entrained V: (6.14) In figure 9(a), the behaviour of the boundary entrainment velocity for the three numerical experiments is shown. In accordance with the flow analysis reported in § 5, the scalar field θ exhibits only a slight reduction of the entrainment rate, while, for the scalar field θ , this reduction is larger. As shown in the inset of figure 9(a), after an initial transient that is more marked for the scalar field θ , a self-similar scaling V b ∝ t −1/2 applies for t > 80. The heuristic approach defined in the previous section is based on the idea that iso-surfaces are determined by the structure of turbulence and are fractal-like (Sreenivasan et al. 1989). The estimate of the corresponding iso-surface area for the three numerical experiments is reported in figure 9(b). Let us recall that for the assumptions made in the previous section, the fractal-like surfaces read respectively, for the scalar fields θ,θ and θ . As shown in figure 9(b), both the scalar fieldsθ and θ experience a reduction of the scalar interface area. In theθ scalar field, this reduction is given by the lack of small-scale wrinkles at the scalar interface due to the absence of small-scale nibbling. Conversely, in the θ scalar field, this reduction is caused by the lack of large-scale convolutions of the scalar interface due to the absence of large-scale engulfment. Since this reduction is larger for the scalar field θ , S <S < S , we argue that large-scale engulfment more efficiently increases the interface area than  Figure 10. (a,b) Predictions of the boundary entrainment rate based on (6.7a,b), (6.10a,b) and (6.12), respectively, for the three scalar fields θ,θ and θ . (c,d) Refined predictions of the boundary entrainment rate based on (6.20) and (6.25).
small-scale nibbling. Let us notice that we have verified these scalings by measuring the actual iso-surface area of θ = θ th and a good agreement has been found. The behaviour of the scalar interface areas supports the observed reduction of entrainment rate for the scalar fieldsθ and θ (see again figure 9a). But, actually, the physics at the basis of the rate of entrainment is more complex and eludes the simple indication given by the interface area. Indeed, the heuristic models for the entrainment rate (6.7a,b), (6.10a,b) and (6.12), despite being based on the fractal surface area, actually show a completely different behaviour. Let us try to grasp the essential reasons.
In figure 10(a,b), the predictions of the entrainment velocity given by (6.7a,b), (6.10a,b) and (6.12) are shown. For the scalars θ andθ, respectively driven by the unfiltered and large-scale velocity fields, the predictions are reasonable in terms of both behaviour and reduction of the entrainment rate for the scalarθ. Compare these behaviours with the measurements reported in figure 9(a). On the other hand, the prediction for the scalar θ lacks in capturing the reduction of the entrainment rate that, on the contrary, is estimated to be the largest one. Despite being characterized by the largest predicted reduction of the scalar iso-surface area, see figure 9(b), equations (6.12) predict a large entrainment rate because of a smaller value of the Batchelor scale and of a larger value of the velocity fluctuation intensity, i.e. These inequalities are simply due to the lower propagation rate of the scalar field θ with respect to θ andθ. As a result, the interface of the scalar field θ experiences a velocity field associated with a much more inner location with respect to that experienced by the other two scalar interfaces since h θ h θ ≈ h θ . At inner positions, the variance and dissipation of velocity fluctuations are larger than those at the outer regions of the jet (see (Cimarelli et al. 2021)) and the operation of filtering out the large-scale fluctuations is not able to balance this flow inhomogeneity thus justifying inequalities (6.16a,b). Overall, we have that the scalar θ driven by the small-scale velocity field experiences a reduction of the scalar interface area S <S < S that is overcome by a larger increase of the diffusion velocity α/η θ int α/η θ int ≈ α/η θ int and of the velocity fluctuations υ int υ int ≈ υ int . For this reason, (6.7a,b), (6.10a,b) and (6.12) actually predict a larger entrainment rate for the scalar θ .

Heuristic refinement of the entrainment and mixing assumptions
Evidently, the heuristic assumptions at the basis of (6.7a,b), (6.10a,b) and (6.12) miss some relevant physical processes. We try here to grasp the essential aspects.
Let us start with the scalar diffusion. The hypothesis (6.2a,b), (6.17) certainly works for the scalar fields θ andθ but not for θ . Indeed, for the latter, the decrease in the scalar surface area, S <S < S , is overcome by a larger increase of the characteristic scalar gradient, (6.18) due both to a slower decay of the mean scalar concentration Θ cl Θ cl ≈ Θ cl (see figure 6a) and to an inner location of the scalar interface that in turn leads to η θ int η θ int ≈ η θ int . It appears that the mean scalar concentration is a valid indicator of the scalar fluctuation intensity for the fields θ andθ but not for θ . Indeed, the increase of mean scalar concentrationΘ cl is associated with an increase also of the scalar fluctuation intensity θθ cl (see figure 6(a,b). On the contrary, the increase of mean scalar concentration Θ cl is associated with a decrease of the scalar fluctuation intensity ϑ ϑ cl (see again figure 6a,b). In practical applications, this different behaviour could be ascribed for example to the presence of external agents or chemical reactions. Accordingly, since the scalar fluctuation intensity could not be directly related to the behaviour of the mean scalar concentration, it is important to directly take it into account in the scalar diffusion expression: (6.19) This expression leads to the following velocities: (6.20) respectively, for the scalar fields θ ,θ and θ . As shown in figure 10(c), these are actually able to predict the reduction of entrainment for both the scalar fieldsθ and θ . We consider now the inertial scalar flux. The hypothesis (6.5), (6.21) certainly works for the scalar fields θ andθ but not for θ . Indeed, for the latter, the decrease in the scalar surface area, S <S < S , is overcome by a larger increase of the velocity fluctuations, (6.22) due to the inner location of the scalar interface for θ . Hence, contrary to measurements, a larger velocity of propagation of the scalar θ is predicted. This behaviour actually suggests that the global intensity of the velocity fluctuations is not a complete indicator for the transport of scalar fluctuations. Evidently, the scalar transport is not accomplished efficiently by all the scales of the turbulent spectrum. The main difference between the intensity of the velocity fluctuationsυ int and υ int is given by the velocity scales involved. The intensityυ int is produced by an anisotropic field of motion characterized by large-scale engulfment events while the intensity υ int is produced by essentially isotropic small-scale fluctuations. It is then possible to argue that only the anisotropic part of the velocity field significantly contributes to the net propagation of scalars in agreement with a number of recent works (Cimarelli et al. 2015(Cimarelli et al. , 2021Buxton, Breda & Dhall 2019). Accordingly, we introduce the isotropy indicator This parameter tends to 1 when the intensity of the fluctuations is increasingly contained in the cross-flow component and diminishes when the anisotropic state departs from the one-component condition by reaching 1/3 corresponding to isotropic conditions. In velocity fields characterized by large-scale engulfment events, γ > 1/3, and, hence, it can be used to identify the portion of the intensity of the velocity fluctuations that is active in the scalar propagation. Accordingly, the inertial scalar flux can be rewritten as (6.24) and the entrainment velocities read respectively, for the scalar fields θ ,θ and θ . As shown in figure 10(d), these predictions of the scalar entrainment velocity are actually able to predict the reduction of entrainment for both the scalar fieldsθ and θ . Let us finally notice that these predictions are based on scaling arguments and, hence, are intended to give a good qualitative rather than a quantitative description of the boundary entrainment rate of figure 9(a).

Conclusions
The entrainment and mixing processes of scalar fields are investigated for the symmetries of an incompressible temporal planar jet. For Schmidt numbers close to unity, both mixing due to turbulence and scalar diffusion occur with similar time scales. In particular, scalar diffusion processes create a flux through scalar iso-concentration surfaces whose geometrical and fluctuating properties are prescribed by the complex dynamics of turbulence, thus forming coupled mechanisms of entrainment which elude a clear understanding. The definition of two scalar fields separately driven by large-and small-scale fluctuations is a sound starting point to address this issue. To this aim, numerical experiments have been conducted in order to establish from a dynamical point of view the scale dependency of turbulence mixing and entrainment of a passive scalar. Three independent scalar fields have been solved, θ,θ and θ , respectively driven by the unfiltered motion and by the large-and small-scale velocity fields defined by the filter length Δ = 2.2λ cl . The large-scale field is found to capture the large-scale anisotropic engulfment events while the small-scale field retains only the small-scale motions commonly associated with nibbling. We found that the evolution of the scalar fieldθ is almost unaltered with respect to that of θ . The propagation of the scalar field is only slightly reducedh θ ≈ h θ and the main difference is given by an increased number of unmixed flow regions within the jet core that leads to a reduced amount of entrained volume,Ṽ < V. Accordingly, it is found that both the inertial and diffusive scalar fluxes are almost unaltered,φ u ≈ ϕ u andφ α ≈ ϕ α , thus showing that anisotropic large-scale fluctuations are responsible for almost the entire mechanism of scalar mixing. On the contrary, the small-scale velocity field driving the scalar θ does not capture engulfment events. In this case, the inertial fluxes are markedly reduced, ϕ u ϕ u , thus showing that small-scale isotropic fluctuations do not contribute significantly to the process. As a result, the spreading of the jet is markedly reduced, h θ h θ . On the other hand, the diffusive mechanisms are almost unaltered, ϕ α ≈ ϕ α , and, together with the absence of large-scale engulfment, allows for a more complete mixing of the scalar jet core by diffusion thus leading to a reduction of the entrained volume V < V that is less pronounced than the reduction of the spreading of the jet h θ h θ . Overall, we found that inertial fluxes guided by large-scale anisotropic phenomena of engulfment initiate the processes of scalar entrainment and mixing and essentially commit their strength to be further processed through a cascade process by nibbling mechanisms at small scales. Hence, this last stage does not prescribe the amount of entrainment and mixing, but essentially adapts itself to the conditions imposed by the large-scale features of the velocity field. This is also demonstrated by the fact that the evolution of the scalar fieldsθ and θ diverges from the very beginning, t = 40, where only the velocity field is significantly changed. We also argue that the velocity stirring of the scalar iso-concentration surfaces is more effectively accomplished by large-scale engulfment events than by small-scale nibbling. Indeed, we estimate that large-scale convolutions of the scalar iso-surface are responsible for a larger increase of the iso-surface area than small-scale wrinkles.
From a heuristic point of view, classical arguments lead to an entrainment velocity that is proportional to the intensity of the velocity fluctuations measured at the scalar surface times its area. These arguments actually fail in predicting the behaviour of the scalar field θ . Indeed, due to the lower rate of propagation, the scalar interface for θ experiences very large velocity fluctuations, υ int υ int ≈ υ int , but the velocity of entrainment is actually reduced. This phenomenon actually suggests that only the large-scale portion of  Table 1. Numerical settings of the present simulation compared to those of the reference simulation by Cimarelli et al. (2021). For the simulation by Cimarelli et al. (2021), the time step is variable and prescribed by the condition CFL < 0.3; hence, the value reported is the average time step used. It is important to note that the values of the Reynolds number and of the domain extension here reported for the simulation of Cimarelli et al. (2021) differ from those reported in their work because the reference length H used for non-dimensionalization was defined as half the initial jet width contrary to the entire initial width considered here (see § 2.1).
the turbulent spectrum is actually active for scalar propagation. Accordingly, we found that the small-scale isotropic structure of the velocity field driving the scalar θ is essentially inactive for the scalar entrainment despite its intensity being large. On the contrary, the large-scale anisotropic motion driving the scalarθ is almost entirely active, thus leading to a larger velocity of propagation even if the intensity of its fluctuations is smaller. The adoption of a simple isotropy indicator to model the active fraction of the velocity fluctuations confirms this conjecture. The large-scale coherent structure of engulfment imposes the rate of entrainment more efficiently with respect to small-scale eddies by activating intense inertial fluxes and by setting the geometry and the area of the scalar surface.
It is worth noting that the present results are relevant also for reduced-order approaches to scalar mixing, in particular for large-eddy simulation where only the large scales of the velocity field are resolved. The present numerical experiments show that the scalar mixing processes are almost unaltered when driven solely by the large scales of the velocity field. Only a lack of scalar mixing in the jet core is observed that, however, can be recovered by simply adding a turbulent diffusivity to the scalar equation using classical gradient-diffusion hypothesis.
In closing this work, let us mention that the considered Reynolds number is rather low, Re λ = 60. By increasing the Reynolds number, the ratio between inviscid and viscous phenomena may change as long as a larger separation of scales is achieved. This aspect should be taken into account for future investigations. work. The two simulations differ in the numerical schemes, the domain extension and the spatial and temporal resolution adopted. The numerical code used in Cimarelli et al. (2021) is less accurate being based on a fourth-order-accurate spatial discretization and a third-order Adams-Bashforth scheme for time integration (see Verstappen & Veldman 2003). On the other hand, the domain lengths in the streamwise and spanwise directions and the spatial resolution are improved with respect to the present settings (see table 1). Despite these differences, the agreement between the two numerical datasets is very good. As shown in figure 11, the temporal behaviours of the mean centreline velocity U cl and of the mean centreline scalar concentration Θ cl predicted by the two simulations are essentially the same. A very good agreement is also observed in the cross-flow distribution of the numerical solutions as shown in figure 12 where the mean and variance profiles of velocity and scalar concentration are reported in self-similar variables, i.e. by re-scaling velocity with U cl , scalar concentration with Θ cl and cross-flow position with the jet half-width z 1/2 defined as the centreline distance at which the mean velocity profile reduces to half its centreline value.