A Priori Direct Numerical Simulation Analysis of the Closure of Cross-Scalar Dissipation Rate of Reaction Progress Variable and Mixture Fraction in Turbulent Stratified Flames

The cross-scalar dissipation rate of reaction progress variable and mixture fraction εcξ~\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{{\varepsilon_{c\xi } }}$$\end{document} plays an important role in the modelling of stratified combustion. The evolution and statistical behaviour of εcξ~\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{{\varepsilon_{c\xi } }}$$\end{document} have been analysed using a direct numerical simulation (DNS) database of statistically planar turbulent stratified flames with a globally stochiometric mixture. A parametric analysis has been conducted by considering a number of DNS cases with a varying initial root-mean-square velocity fluctuation u′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u^{\prime }$$\end{document} and initial scalar integral length scale ℓϕ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell_{\phi }$$\end{document}. The explicitly Reynolds averaged DNS data suggests that the linear relaxation model for εcξ~\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{{\varepsilon_{c\xi } }}$$\end{document} is inadequate for most cases, but its performance appears to improve with increasing initial ℓϕ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell_{\phi }$$\end{document} and u′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u^{\prime }$$\end{document} values. An exact transport equation for εcξ~\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{{\varepsilon_{c\xi } }}$$\end{document} has been derived from the first principle, and the budget of the unclosed terms of the εcξ~\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{{\varepsilon_{c\xi } }}$$\end{document} transport equation has been analysed in detail. It has been found that the terms arising from the density variation, scalar-turbulence interaction, chemical reaction rate and molecular dissipation rate play leading order roles in the εcξ~\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{{\varepsilon_{c\xi } }}$$\end{document} transport. These observations have been justified by a scaling analysis, which has been utilised to identify the dominant components of the leading order terms to aid model development for the unclosed terms of the εcξ~\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{{\varepsilon_{c\xi } }}$$\end{document} transport equation. The performances of newly proposed models for the unclosed terms have been assessed with respect to the corresponding terms extracted from DNS data, and the newly proposed closures yield satisfactory predictions of the unclosed terms in the εcξ~\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{{\varepsilon_{c\xi } }}$$\end{document} transport equation.


Introduction
Stratified mixture combustion occurs when there is a limited mixing time between the unburned reactants such that some premixing takes place but not to the extent of homogeneity. This type of combustion is often encountered in industrial combustors since it offers benefits in terms of improved energy efficiency and reduced pollutant emissions (Lipatnikov, 2017). Stratified combustion may also occur unintentionally by not allowing a sufficient mixing time in premixed combustors. Therefore, it is necessary to have reliable models for numerical simulations which will be pivotal to the design of future generation combustors that exploit the advantages offered by turbulent stratified mixture combustion.
A complete thermo-chemical description of stratified combustion requires a passive scalar (e.g. mixture fraction ) to describe the local mixture composition and an active scalar (e.g. reaction progress variable c ) to determine the progress of the chemical reaction. The normalised mixture fraction assumes a value of zero in the pure air stream and unity in the pure fuel stream, and is defined as (Bilger, 1988): where Y F and Y O are the fuel and oxidiser mass fractions respectively, s = (Y O ∕Y F ) st is the mass stoichiometric ratio, and a subscript F∞ ( O∞ ) indicates the mass fraction value of fuel (oxygen) in a pure fuel (air) stream. Likewise, the reaction progress variable c quantifies the extent of completion of the chemical reaction by rising from zero in the unburned reactants to a value of unity in the fully burned products, and can be defined as (Hélie and Trouvé, 1998): where st = Y O∞ ∕(sY F∞ + Y O∞ ) is the stoichiometric mixture fraction.
Presumed probability density function (PPDF) (Libby and Williams, 2000;Ribert et al., 2005;Robin et al., 2006), flamelet based tabulated chemistry (Darbyshire et al., 2010) and Flamelet Generated Manifold (Nguyen et al., 2010;Fiorina et al., 2015;Marincola et al., 2013) based modelling approaches require solving the transport equations of the Favre averaged active and passive scalar variances, as well as their covariance, which in turn need closures of their respective dissipation rates. Several previous analyses (Mura et al., 2007;Chakraborty, 2010a, 2011a;Chakraborty et al., 2011) focussed on modelling the scalar dissipation rates of mixture fraction and reaction progress variable in turbulent stratified mixture combustion, but relatively limited effort has been directed to the closure of cross-scalar dissipation rate. Mura et al. (2007) and Nguyen et al. (2010) proposed algebraic closures of the cross-scalar dissipation rates of fuel mass fraction and mixture fraction fluctuations ̃ Y = D∇Y �� F ⋅ ∇ �� ∕ (where D is the molecular diffusivity, and q , q = q∕ and q �� = q −q are the Reynolds averaged, Favre-averaged and Favre fluctuation of a general quantity q respectively, with being the gas density) which plays a key role in the context of Libby-Williams (L-W) model (Libby and Williams, 2000). Following this, Malkeson and Chakraborty (2011b) derived the exact transport equation of ̃ Y and proposed closures for the unclosed terms based on a-priori DNS analysis. Flame Generated Manifold (Nguyen et al., 2010;Marincola et al., 2013;Fiorina et al., 2015) closures of stratified mixture combustion utilise reaction progress variable c rather than Y F as the active scalar, and therefore it is important to understand the behaviour of ̃ c in the flame brush of turbulent stratified flames. However, to date, the closure of the reaction progress variable and mixture fraction cross scalar dissipation rate ̃ c = D∇c �� ⋅ ∇ �� ∕ has not been addressed in detail in the existing literature (Malkeson and Chakraborty, 2010a, b) despite its importance in the closure of the transport equation of the Favre averaged reaction progress variable c , the covariance c′′ ′′ and flame displacement speed S d in turbulent stratified mixture combustion . Although Y F and c are both active scalars, their gradients and correlations with ∇ can vary significantly depending on the extent of mixture stratification. Thus, it is important to analyse ̃ c independently since closures of ̃ Y do not automatically apply to ̃ c . Thus, the main objectives of this paper are to: (1) analyse the statistical behaviour of ̃ c within the flame brush in turbulent stratified flames.
(2) evaluate the effectiveness of the linear relaxation algebraic model at predicting the variation of ̃ c within the flame brush.
(3) identify the leading order contributors in the transport equation of ̃ c for turbulent stratified mixture combustion for different turbulence intensities and mixture stratification conditions.
(4) provide adequate model expressions for the unclosed terms of the ̃ c transport equation.
To meet the objectives, a three-dimensional Direct Numerical Simulation (DNS) database of six statistically planar turbulent stratified flames has been considered (Brearley et al., 2020). This database consists of cases with a globally stochiometric mixture and an initially bimodal distribution of equivalence ratio in the unburned gas. The cases vary by initial values of root-mean-square velocity fluctuation normalised by the laminar burning velocity of the stoichiometric mixture u � ∕S b( =1) and the initial scalar integral length scale normalised by the velocity integral length scale ∕ .
The rest of the paper will be organised in the following manner. The mathematical background and numerical frameworks pertaining to this analysis are presented in the next two sections. This will be followed by the discussion of the results before the summary of the main finding, and conclusions are drawn.

Mathematical Background
The cross-scalar dissipation rate ̃ c can be modelled using an algebraic expression in the context of the linear relaxation (LR) methodology (Mura et al., 2007;Malkeson and Chakraborty, 2010a), where the cross-scalar dissipation rate is modelled as: where C is a model constant of the order of unity, ̃= u �� j ∕ x k u �� j ∕ x k ∕ is the Favre-averaged turbulent kinetic energy dissipation rate, k = u �� j u �� j ∕2 is the Favre averaged turbulent kinetic energy and is the kinematic viscosity. Alternatively, an additional modelled transport equation for ̃ c can be solved. The transport equation for ̃ c can be by considering the individual transport equations for the reaction progress variable  and the mixture fraction: where ̇c is the reaction rate of reaction progress variable c , which takes the following form when c is defined using Eq. 2 Malkeson and Chakraborty, 2010b): with ̇F being the reaction rate of the fuel. The term A is the cross-scalar dissipation contribution, which is defined as Malkeson and Chakraborty, 2010b): It is worth noting that the statistical behaviour of N c = D(∇c ⋅ ∇ ) determines the behaviour of A . The cross-scalar dissipation rate values are often smaller than the values of scalar dissipation rate of the reaction progress variable (Malkeson and Chakraborty, 2010a;Inanc et al., 2022). However, the cross-scalar dissipation rate appears explicitly in the expression of displacement speed S d = |∇c| −1 (Dc∕Dt) in stratified flames Malkeson and Chakraborty, 2010b) and this quantity accounts for the effects of relative orientations of ∇c and ∇ on flame propagation (i.e. back-supported or front-supported flame propagation) which may not be negligible for stratified flames (e.g. Marincola et al., 2013;Inanc et al., 2020). Moreover, scalar dissipation rates N = D∇ ⋅ ∇ , N c = D∇c ⋅ ∇c and cross-scalar dissipation rate N c = D∇c ⋅ ∇ affect the value of in stratified flames .
Differentiating Eqs. 4 and 5 with respect to x i one obtains transport equations for c∕ x i and ∕ x i , respectively. Multiplying the c∕ x i transport equation by ∕ x i , and multiplying the ∕ x i transport equation by c∕ x i , then summing the resulting equations leads to a transport equation for the instantaneous cross-scalar dissipation rate N c = D∇c ⋅ ∇ . The transport equation for the cross-scalar dissipation rate of the Favre averaged scalars can be derived in a similar fashion. The transport equation of the Favre-averaged reaction progress variable c is used to obtain the transport equation of c∕ x i , and the transport equation of the Favre-averaged mixture fraction ̃ is used to obtain the transport equation of ̃∕ x i . Multiplying the c∕ x i transport equation by ̃∕ x i , and multiplying the ̃∕ x i transport equation by c∕ x i , then summing the resulting transport equations together yields the transport equation for D ∇c ⋅ ∇̃ . The final step of the derivation is to subtract the two derived equations to yield � c = � N c −D∇c ⋅ ∇̃ . For further information, interested readers are referred to Malkeson and Chakraborty (2011b) where a full mathematical description of the derivation of the transport equation of ̃ Y is presented, which is identical to the derivation of the transport equation of ̃ c . The ̃ c transport equation is shown below after ignoring the terms arising due to diffusivity variations: The first term on the left-hand side represents the temporal change of ̃ c , while the second term represents the effects of mean advection. On the right-hand side, D 1 is the contribution of molecular diffusion, −D 2 is the contribution of molecular dissipation, T 1 is the contribution of turbulent transport, T 2 arises due to density variation, T 3 is the contribution of scalar gradient alignment with fluid dynamic strain rates, and T 4 is the combined contribution of chemical reaction rate and mixture inhomogeneity. The terms T 1 , T 2 , T 3 , T 4 and −D 2 are unclosed and require modelling. The mathematical definitions of these terms are given as: The statistical behaviours of T 1 , T 2 , T 3 , T 4 and −D 2 will be extracted from the explicitly averaged DNS data and will be discussed in Sect. 4 of this paper. (9(i)) 1 3

Numerical Implementation
The simulations of turbulent stratified flames were simulated using the DNS code SENGA + (Jenkins and Cant, 1999). The simulations deal with statistically planar turbulent stratified flames propagating into an inhomogeneous methane-air mixture. A modified single-step chemistry proposed by Tarrazo et al. (2006) representative of methane-air combustion has been considered for the current analysis where the activation temperature and heat of combustion are taken to be functions of equivalence ratio. It has been shown elsewhere Chakraborty, 2010a, b, 2011a) that the laminar burning velocity variation with equivalence ratio is accurately captured using this thermo-chemistry. This simplification of chemical representation allows for a detailed parametric analysis in terms of u � ∕S b( =1) and ∕ (which was reported in detail by Brearley et al. (2020)) with a reasonable amount of computational cost. The relative costs between simple and detailed chemistry DNS were discussed in detail by Keil et al. (2021a). It has been demonstrated recently by Keil et al. (2021a, b) that the flame propagation statistics obtained from simple chemistry DNS are found to be qualitatively similar to the results obtained from detailed chemistry DNS. Moreover, the quantitative differences between displacement speed statistics between the simple and detailed chemistry results are comparable to the uncertainty involved with the different possible definitions of reaction progress variable. As most of the heat release takes place in premixed mode in stratified mixture combustion, it is expected that the findings of Keil et al. (2021a, b) are likely to be also valid for stratified mixture combustion. It was indeed demonstrated by  that the statistics of scalar gradients obtained from the thermochemistry used in this paper remain in good qualitative agreement with detailed chemistry DNS. Finally, several previous simple chemistry based analyses including DNS studies (Hélie and Trouvé, 1998;Ribert et al., 2005;Robin et al., 2006;Mura et al., 2007;Chakraborty, 2010a, b, 2013) gave rise to significant advancements in the fundamental understanding and modelling of stratified mixture combustion and the same approach has been adopted in this paper.
In all cases, a grid size of 800 × 400 × 400 grid points uniformly distributed in a computational domain of size 70.2 st × 30.1 st 2 has been considered, where is the thermal flame thickness of the stoichiometric mixture with T ad( =1) , T 0 and T being the adiabatic flame temperature of the stoichiometric mixture, unburned gas temperature, and the instantaneous dimensional temperature, respectively. The subscript L refers to the values in the 1D unstretched laminar premixed stoichiometric methane-air flame. The aforementioned grid spacing satisfies the requirement of having at least ten grid points within st . The mean direction of flame propagation runs parallel to the x direction (long side) of the computational domain, with opposite-facing boundaries on the unburned and burned side of the computational domain being taken to be partially non-reflecting and specified in accordance with the Navier-Stokes Characteristic Boundary Conditions technique (Poinsot and Lele, 1992). The transverse boundaries are taken to be periodic. All the spatial differentiations for the internal grid points are carried out using a tenth-order central difference scheme for the internal grid points but the order of accuracy gradually reduces to a one-sided second-order scheme at the two nonperiodic boundaries. The time advancement has been conducted using an explicit low storage third-order Runge-Kutta scheme (Wray, 1990). The initial velocity fluctuations were specified using a pseudo-spectral method (Rogallo, 1981) for a prescribed integral length scale and rms of velocity fluctuations u ′ that follows the Batchelor-Townsend spectrum (Batchelor and Townsend, 1948). The mixture inhomogeneity in the unburned gas is initialised by a bi-modal distribution of equivalence ratio = 1 − st ∕ st (1 − ) following the methodology proposed by Eswaran and Pope (1988) for the prescribed values of mean equivalence ratio , root-mean-square (rms) equivalence ratio fluctuation ′ , and integral length scale of equivalence ratio fluctuations. All species are assumed to be perfect gases with a Lewis number of unity. The heat release parameter = T ad( =1) − T 0 ∕T 0 is taken to be 4.5 for all cases considered here, which represents unburned reactants preheated to 415 K for stoichiometric methane-air mixtures. Standard values are considered for Prandtl number Pr = 0.7 and the ratio of specific heats = 1.4 . The turbulence length scale to stoichiometric flame thickness ratio is taken to be ∕ st = 3.0 for all cases. The initial mean equivalence ratio is taken to be unity (i.e., = 1.0 ) and the initial rms equivalence ratio is taken to be � = 0.35. The stratified mixture combustion is expected to remain within the flammability limit and thus the range of equivalence ratio encountered in this analysis remains within the flammability range (i.e. 0.6 < < 1.5 ) for methane-air combustion for the choices of = 1.0 and � = 0.35.
The reacting scalar and temperature fields within the flame are initialized by 1D unstretched laminar premixed flame solution corresponding to = 1.0 and = 4.5 . Table 1 shows the initial values of the simulation parameters u k ∕ 0 (where 0 and 0 are the unburned gas density and viscosity, respectively and k is the turbulent kinetic energy evaluated over the whole domain) is calculated based on the integral length scale and ranges from 42.0 to 84.0 for initial u � ∕S b( =1.0) = 4.0 to 8.0 , respectively. All the cases considered here nominally belong to the thin reaction zones regime combustion (Peters, 2000). For these initial conditions, 2.41 and 1.44 grid points reside within the Kolmogorov length scale for initial u � ∕S L = 4.0 and 8.0 respectively, with this number increasing as the simulation progresses due to the decay of the turbulence. Under decaying turbulence, the simulations should be run for t sim ≥ max T turb , T chem , where the eddy turnover time T turb = ∕u � , and the chemical time scale t chem = st ∕S b( =1) . In all the cases, T chem remains greater than T turb , and all cases considered here have been run for t sim = 2.0T chem , which is about 2.67T turb and 5.33T turb for initial u � ∕S L = 4.0 and 8.0 respectively. This simulation duration remains either comparable to or greater than several previous analyses (Hélie and Trouvé, 1998;Haworth et al., 2000;Jimenez et al., 2002;Chakraborty, 2010a,b, 2011a,b). It is important to note that the cross-scalar dissipation statistics are extracted in this paper at the end of the simulation time when the quasi-stationary state has been obtained in terms of temporal evolutions of flame surface area and turbulent burning velocity, as previously shown by Brearley et al. (2020) for this database. The cross-scalar dissipation rate statistics and model performances reported in this paper remain qualitatively similar at different time instants since the quasi-stationary state is obtained (e.g. roughly halfway through the simulation).

3
To obtain the Reynolds/Favre-averaged values of a general quantity (i.e. Q and Q ), the quantities of interest are arithmetically averaged over the statistically homogeneous transverse planes (i.e. y − z directions) normal to the direction of mean flame propagation (i.e. the x direction) following previous analyses (Hélie and Trouvé, 1998;Swaminathan and Bray, 2005;Chakraborty, 2010a,b, 2011a,b;Chakraborty et al., 2011). For statistically planar flames, c is a unique function of the coordinate in the direction of the mean flame propagation (i.e. the x direction in this case) and thus all the results in Sect. 4 will be represented as a function of c . It has been assessed that halving the sample size in the transverse direction does not significantly alter the results.
Three additional simulation cases where initial u � ∕S b( =1) = 10 from the database reported by Brearley et al. (2020) were used to assist in the statistical analysis of ̃ c and the model development. These results have been omitted for conciseness and because they do not enhance the discussion due to their similarity with the results for initial u � ∕S b( =1) = 8.0.

Flame-Turbulence Interaction
The distributions of normalised mixture fraction ∕ st and reaction progress variable c in the central x − y plane when the statistics are extracted are shown in Fig. 1 for the cases considered here. All the cases considered here are representative of combustion in the thin reaction zones regime (Peters, 2000) and accordingly, local occurrences of flame thickening can be observed from the reaction progress variable contours due to penetration of energetic, turbulent eddies within the preheat zone. However, the reaction zone thickness remains smaller than the Kolmogorov length scale and thus this zone remains unperturbed by the small-scale turbulence but gets wrinkled by the large-scale turbulent motion. It can further be seen from Fig. 1 that the mixture inhomogeneity is predominantly obtained on the unburned gas side of the flame and the effects of stratification decay towards the burned gas side of the flame due to the increase in mass diffusivity with temperature. This behaviour has implications on the distribution of co-variance c′′ ′′ and cross-scalar dissipation rate ̃ c , which will be discussed in the following sub-section.

Statistical Behaviours of Co-Variance and Cross-Scalar Dissipation Rate
The variations of the variances c′′2 , Ỹ′′2 F and ̃ ′′2 and co-variances of c′′ ′′ and Ỹ′′ F ′′ with c across the flame brush are shown in Fig. 2 for all cases. It can be seen from Fig. 2 that the magnitudes of ̃ ′′2 , c′′ ′′ and Ỹ′′ F ′′ increase with increasing initial ∕ value. The mean scalar dissipation rate of mixture fraction in the volume of inhomogeneous mixture can be taken to scale as ⟨N ⟩ ∼ D⟨ �2 ⟩∕ 2 (Hélie and Trouvé, 1998;Malkeson and Chakraborty, 2011a;Brearley et al., 2020) with being the Taylor micro-scale of mixture fraction distribution. For moderate values of turbulent Reynolds number, remains of the same order as , and thus ⟨N ⟩ assumes high values for small values of . This leads to a more rapid mixing rate for smaller values of and that is why the mixture fraction fluctuations and ̃ ′′2 remain small for the cases with small initial values of . This is reflected in the small magnitudes of ̃ ′′2 for small initial values of and accordingly small magnitudes of c′′ ′′ and Ỹ′′ F ′′ are obtained for small values of . This behaviour is particularly strong for large values of u ′ (e.g. initial u � ∕S b( =1) = 8.0 case) because of enhanced mixing, which leads to decreases in magnitudes of ̃ ′′2 , c′′ ′′ and Ỹ′′ Figure 2 shows that the variations of c′′ ′′ and Ỹ′′ F ′′ within the flame brush can be significantly different from each other. In the unburned gas Y �� F = Y F∞ �� and thus at the leading edge Ỹ′′ F ′′ and Ỹ′′2 F take a value of ̃ ′′2 for Y F∞ = 1.0 , whereas c′′ ′′ vanishes on the unburned gas side of the flame because c �� = 0 in the unburned gas. Similarly, c �� = 0 in the burned gas and therefore c′′ ′′ vanishes on the burned gas side of the flame brush but Ỹ′′ F ′′ can assume non-zero values in the burned gas. The maximum value of Ỹ′′2 F is obtained at the middle of the flame brush for all cases. However, it is evident from Fig. 2 that the distributions of c′′2 and Ỹ′′2 F are markedly different. The reaction progress variable variance c′′2 attains the maximum value close to the middle of the flame brush (i.e., c ≈ 0.5 ) for all cases but its maximum value remains smaller than 0.25 (i.e., max 0.02 × � c ��2 < 5 × 10 −3 ). In the limit of infinitely fast chemistry (i.e., Da ≫ 1 ), the PDF of c can be considered to be bimodal with 1 3 impulses at c = 0.0 and c = 1.0 (Bray et al., 1985), which yields � c ��2 =c(1 −c) and max c��2 = 0.25 (i.e. max 0.02 ×c ��2 = 5 × 10 −3 ). This further suggests that � c ��2 =c(1 −c) cannot be used for the low Damköhler number (i.e., Da < 1 ) stratified mixture combustion considered here where � c ��2 <c(1 −c) and the PDF of c cannot be approximated by a bimodal distribution with impulses at c = 0.0 and c = 1.0 . Although the distributions of c′′ ′′ and Ỹ′′  ′′ for turbulent stratified mixtures was presented elsewhere  and the transport equation of c′′ ′′ can be obtained by replacing Y F with c in the transport equation of Ỹ′′ F ′′ , which yields: It can be seen from Eq. 10 that the closure of ̃ c is needed to solve this equation. The variations of ̃ c and ̃ Y with c for all cases considered here are shown in Fig. 3. The predictions of the linear relaxation model are also shown for C c = 3.0 and C Y = 3.0 which provides a reasonable order of magnitude agreement with DNS data for most cases. The results do not reveal any significant relation between C and the parameters tested in this study ( u � ∕S b( =1) , and ∕ ). The model performs relatively better for cases with high initial u � ∕S b( =1) . The linear relaxation model only considers the turbulent time scale and ignores the contribution of the chemical timescale. Thus, this model performs relatively better for low Damköhler number cases where the effects of chemistry have a reduced influence. The linear relaxation model also performs better for large values of ∕ , indicating that the existence of considerable mixture inhomogeneity is a requirement for the model to perform well. Moreover, this model also fails to capture the qualitative behaviour and the correct sign of ̃ c for small and moderate values of (see cases A-D in Fig. 3). The findings from Fig. 3 suggest that the linear relaxation model cannot be a sufficiently general model for the closure of ̃ c in turbulent stratified mixture combustion. Figure 3 also shows the variation ̃ Y and its equivalent linear relaxation model for the sake of comparison. Due to the fact that Y F and are closely related (i.e. Y F = Y F∞ ) in the unburned gas, ̃ Y assumes non-zero value in the unburned gas and decays within the flame brush for all cases. Furthermore, ∇Y �� F and ∇ �� remain significantly correlated inside the flame which is evident from the mostly positive values of ̃ Y . The quantities ∇c �� and ∇ �� also remain somewhat correlated inside the flame, which yields non-zero values of ̃ c within the flame brush. However, the variations of ̃ c and ̃ Y are qualitatively different which indicates that the correlation between ∇c �� and ∇ �� is fundamentally different to that between ∇Y �� F and ∇ �� . The behaviour of ̃ Y is found to be consistent with previous findings by Malkeson & Chakraborty (2010a, 2011b. The findings of Fig. 3 indicate that it might be necessary to solve a modelled transport equation for the closures of ̃ c and ̃ Y . The closures of the unclosed terms of the transport equation of ̃ Y have been discussed elsewhere (Malkeson and Chakraborty, 2011b) and this analysis will henceforth focus on the closures of the unclosed terms of the ̃ c transport equation. (10)

Statistical Behaviours of the Unclosed Terms of the ̃ c Transport Equation.
The variations of the terms appearing on the right-hand side of the ̃ c transport equation (see Eq. 8), across the flame brush for all cases are shown in Fig. 4. It is evident from Fig. 4 that the magnitude of molecular diffusion D 1 remains negligible in all cases. The contribution of turbulent transport T 1 plays a significant role only for the small u � ∕S b( =1) cases but its importance diminishes as the turbulence intensity increases. The remaining terms (i.e. T 2 , T 3 , T 4 , −D 2 ) play leading order roles in the evolution of ̃ c and have similar orders of magnitude for all cases considered here. The contribution of density variation T 2 assumes mostly positive values for a major part of the flame brush. The term representing the scalar turbulence interaction T 3 takes a predominantly negative value throughout the flame brush in these cases with a similar magnitude to that of T 2 . The combined contribution of chemical reaction rate and mixture inhomogeneity T 4 assumes predominantly positive values towards the unburned gas side of the flame brush before becoming weakly negative towards the burned gas side of the flame brush. The molecular dissipation term −D 2 exhibits predominantly negative values but local positive values can be discerned on the burned gas side of the flame brush for case B and on the unburned gas side of the flame brush for case E. In cases (A-D) the terms T 3 , T 4 − D 2 play leading order roles but in these cases T 4 predominantly assumes negative values for a major part of the flame brush.
Replacing c by Y F ( ̇c by ̇F ) and setting A = 0 in Eqs. 8 and 9 yields the transport equation of ̃ Y (see Malkeson and Chakraborty (2011b) for the full equation). The variations of the terms appearing on the right-hand side of the ̃ Y transport equation across the flame brush for all cases are shown in Fig. 5. A comparison between Figs. 4 and 5 reveals that the evolution of ̃ Y is fundamentally different to that of ̃ c . Most notably, T 3 acts as a source term and −D 2 acts as a sink term in the unburned gas for the ̃ Y transport with approximately similar magnitudes. However, the magnitude of the sink term exceeds that of the source term in all cases due to the nature of decaying turbulence. The behaviours of the terms of the ̃ Y transport equation were discussed in detail in Malkeson and Chakraborty (2011b) and thus will be not be elaborated further in this analysis. The closures of the unclosed terms of the Ỹ′′  The observations made from Fig. 4 can further be substantiated using a scaling analysis where the density is scaled by the unburned gas density 0 , the mean velocity is scaled using a reference mean velocity scale U ref , the fluctuating velocity components are scaled by S b , spatial derivatives of the mean quantities are scaled by , whereas spatial derivatives of  Swaminathan and Bray (2005). The assumptions yield the following scaling estimates for the terms shown in Fig. 4: Alternatively, following Tennekes and Lumley (1972) and Mantel and Borghi (1994), if the turbulent velocity fluctuations are scaled using u ′ , the mean gradients are scaled using the integral length scale , the length scale associated with the gradient of fluctuating velocity and mass fractions are scaled using the Taylor micro-scale T , and the fluctuations of density gradient and reaction rate gradient are scaled with respect to st , it is possible to obtain: Malkeson and Chakraborty (2011b) argued that the second derivatives of the mass fraction and mixture fraction fluctuations are scaled with respect to the dissipation cut-off scale D , which can be taken to be equal to the Obukhov-Corrsin scale D 3 ∕̃ 3∕4 in the thin reaction zones regime (Peters, 2000). As the Schmidt number Sc remains of the order of unity (i.e. Sc = 0.7 ), it is possible to scale D 21 and D 22 as: In order to demonstrate the applicability of the scaling relations given by Eqs. 11i-iii, the variations of different components of T 1 , T 2 , T 3 , T 4 , (−D 2 ) with c are exemplarily shown in Fig. 6 for cases E and F. The same qualitative behaviour has been observed for the other cases, so they are not explicitly shown here for the sake of conciseness. It can be seen from Fig. 6  . Moreover, Eqs. 11i and ii reveal that the magnitudes of D 1 and T 1 become increasingly insignificant in comparison to T 2 , T 3 , T 4 and −D 2 with increasing Re t , and the leading order contributions arise from the terms T 2 , T 3 , T 4 and (−D 2 ). This conclusion is consistent with the observations made from Fig. 4. Figure 6 further shows that T 2 ≈ T 21 , and this is supported by Eqs. 11i and 11ii indicating the weakening of T 22 with increasing Re and/or Da . Equation 11 indicates that T 32 acts as a leading order term, and the contributions T 31 and T 33 weaken with increasing Re and/or Da . It can be seen from Fig. 6 that that T 32 is the dominant contributor to T 3 but in this case T 31 and T 33 play non-negligible roles due to small values of Da (i.e. Da < 1 ). Equation 5 suggests that T 41 is expected to dominate over T 42 , which is consistent with the results shown in Fig. 6. Finally, Eqs. 11i and 11iii show that D 21 is the dominant term in D 2 , with the contribution of D 22 weakening with increasing Re and/or Da , which is consistent with the variations shown in Fig. 6. The scaling estimates based on Eq. 11 and their validation based on simulation results will play a key role for the purpose of model development for T 1 , T 2 , T 3 , T 4 and (−D 2 ) for the ̃ c transport equation.
Modelling of the Turbulent Transport Term T 1 . It can be seen from Eq. 9i that the modelling of T 1 depends on the closure of T 11 because T 12 and T 13 are closed in the context of second-moment closure as the scalar fluxes u ′′ j c ′′ and u ′′ j ′′ are needed for the purpose of solving the transport equations of c and ̃ , respectively. Furthermore, it can be seen from Eq. 11 that the magnitudes of T 12 and T 13 are expected to be smaller than that of T 11 for high values of Re t . Thus, the modelling of T 11 is the key to the closure of T 1 and therefore, the modelling of T 11 will be discussed in this sub-section. Equation 9i further indicates that T 11 closure depends on the modelling of u ′′ j c . The simplest closure that can be used for u ′′ j c is the usual gradient hypothesis which can be expressed as: where t = Ck 2 ∕̃ is the eddy viscosity and c is the appropriate turbulent Schmidt number. The predictions of Eq. 12 with c = 1.0 are compared to u ′′ 1 c (i.e. the only non-zero component of u ′′ j c ) extracted from DNS data in Fig. 7 for the cases considered here. It can be seen from Fig. 7 that Eq. 12 does not adequately capture both qualitative and quantitative behaviours for a major part of the flame brush especially for small values of u � ∕S b( =1) towards the burned gas side of the flame brush. The different sign of the prediction of Eq. 12 in comparison to u ′′ 1 c extracted from DNS data is indicative of the counter-gradient transport. The effects of counter-gradient transport become strong within the  (Veynante et al., 1997) and therefore the discrepancy between u ′′ j c obtained from DNS data and the prediction of Eq. 12 are particularly prominent from the middle to the burned gas side of the flame brush. Therefore, Eq. 12 cannot be considered as a general model for u ′′ j c , and an ideal model for u ′′ j c should be able to capture both gradient and counter-gradient transport. In order to meet this objective, a new for u ′′ 1 c has been proposed in the following manner: where The mean burned gas density and the mean laminar burning velocity are evaluated as (Domingo et al., 2002): Here, S b ( ) and b ( ) are the laminar burning velocity and burned gas density expressed as a function of mixture fraction , respectively, P( ) is the probability density function (PDF) of , and the lower and upper limits of mixture fraction are given by l and u , respectively. The PDF of can be modelled using a -function, which was discussed elsewhere (Poinsot and Veynante, 2001;) and thus will not be repeated here. The term ensures that u ′′ 1 c assumes the sign depending on the nature of the transport (i.e. whether counter-gradient or gradient) behaviour of ( u �� 1 c �� ) because u ′′ 1 ′′ predominantly exhibits gradient type transport (i.e. − u �� 1 �� = ( t ∕ )̃∕ x 1 with being the turbulent Schmidt number), as is a passive scalar. The term � is representative of the Bray number (Veynante et al., 1997) for turbulent stratified combustion. It was discussed elsewhere (Veynante et al., 1997;Chakraborty and Cant, 2009a, b) that counter-gradient (gradient) type transport is obtained when the Bray number (i.e. ∝ � ) is greater (smaller) than unity (Veynante et al., 1997).
√k ≫ 1 suggests a situation where the flame normal acceleration due to thermal expansion overcomes the effects of turbulent fluctuations, which lead to counter-gradient transport (Veynante et al., 1997). By contrast, turbulent velocity fluctuations overwhelm the effects of thermal expansion effects for � Veynante et al., 1997). The involvement of g in Eq. 13i ensures that the first term on the right-hand side of Eq. 13i signifying gradient type transport does not play a major role for � √k ≫ 1 and the behaviour of u ′′ 1 c is governed by the second term on the right-hand side of Eq. 13i. By contrast, the first term on the right-hand side plays the dominant role for small values of � and ensures gradient type transport of u ′′ 1 c . The predictions of Eq. 13 are also shown in Fig. 7, which shows that Eq. 13 mostly captures the correct sign of u ′′ 1 c extracted from DNS data more successfully than the gradient hypothesis model (i.e. Equation 12). This also suggests that Eq. 13 allows for the prediction of counter-gradient transport, which Eq. 12, by design, is not able to capture. Although there is a scope for improvement of the performance of Eq. 13 for cases with small initial values of ∕ (e.g. ∕ = 1.0 cases A and B) irrespective of turbulence intensities, this model performs reasonably well for high initial values of ∕ (e.g. ∕ = 2.0 and 3.0 cases). It is important to note that the mixing rate increases with decreasing ∕ (Hélie and Trouvé, 1998;Brearley et al., 2020) and thus one obtains an almost homogeneous mixture as a result of mixing for small values of ∕ (e.g. cases A and B). Therefore, the inaccurate predictions of Eq. 13 for cases A and B may not have a major implication on the modelling of turbulent stratified mixture combustion. Moreover, Eq. 11 suggests that T 1 is not the leading order contributor in the ̃c transport equation and this can be substantiated from Fig. 4. Thus, the local discrepancy between Eq. 13 predictions of u ′′ 1 c may not have a major implication on the transport equationbasedclosure of ̃c . It is important to note that the success of the model given by Eq. 13 depends on the closure of ( u �� 1 c �� ) . The modelling of ( u �� 1 c �� ) was discussed elsewhere (Veynante et al., 1997;Cant, 2009a,b, 2015;Malkeson and Chakraborty, 2012) in detail and thus is not discussed in this paper.

Modelling of the Density Variation Term T 2 .
The variations of T 2 with c for cases A-F are shown in Fig. 8a-f, respectively. The density variation term T 2 predominantly shows positive values for all cases considered here. The order of magnitude of T 2 can be estimated using T 2 ∼̃c̃c 0 ∕ b − 1 upon scaling ̇c as ̇c ∼ c and −(D∕ ) p∕ x i ∕ x i can be scaled as: It is important to note that the qualitative behaviour of T 2 was found to be captured by ̃c̃c 0 ∕ b − 1 as dictated by the scaling analysis and a proportionality parameter K is needed for the quantitative prediction. This yields the following model for T 2 : It has been found that the proportionality parameter K is not sensitive to the ∕ value, but it is a function of the turbulence intensity (or local turbulent Reynolds number). This is empirically expressed as: where Re L = 0k 2 ∕ 0̃ is a measure of the turbulent local Reynolds number. The expression given by Eq. 14ii is one of the possibilities out of several possible ones but this functional form (i.e., K = a 1 + a 2 ∕ 1 + exp − a 3 Re L − a 4 ) was chosen because it ensures that K increases from 0.1 to an asymptotic value (i.e., K = 0.25 ) for Re L → ∞ . The exact values of a 1 , a 2 , a 3 and a 4 are chosen based on a regression analysis based on T 2 and ̃c̃c 0 ∕ b − 1 extracted from DNS data. The practical applicability of the model is not expected to be sensitive to a 3 and a 4 because Re L ≫ 15.0 is obtained for most practical applications. The predictions of Eq. 14 are compared with T 2 extracted from DNS data in Fig. 8, which reveals that the density variation term T 2 can be reasonably predicted by the expression given by Eq. 14.
Modelling of the Scalar-Turbulence Interaction Term T 3 . In order to model the scalar-turbulence interaction term T 3 , all the components of this term (i.e. T 31 , T 32 and T 33 ) need to be closed. Therefore, modelling of each of the 2 st with c along with the predictions of Eq. 14 for cases A-F (a-f) components of T 3 for the ̃c transport equation will be addressed in turn. The term T 31 is made up of two sub-terms T (i) 31 and T (ii) 31 as shown in Eq. 9iii. The statistical behaviours of T (i) 31 and T (ii) 31 are governed by the alignment of local strain rate eigendirections with reaction progress variable gradient and mixture fraction gradient respectively. It has been observed that the mixture fraction gradient aligns preferentially with the most compressive principal strain rate eigendirection. By contrast, reaction progress variable alignment with local strain rate eigendirection is dependent on the relative strengths of turbulent straining and the strain rate induced by heat release (Chakraborty and Swaminathan, 2007a;Chakraborty et al., 2009;Malkeson and Chakraborty, 2011c). The reaction progress variable gradient is shown to align with the most extensive principal strain rate when the strain rate induced by heat release overwhelms the effects of turbulent straining and vice versa (Chakraborty and Swaminathan, 2007a;Chakraborty et al., 2009;Malkeson and Chakraborty, 2011c). The terms can be taken to scale with u �� j �� ̃∕k and u �� j c �� ̃∕k based on previous modelling strategies adopted by Mantel and Borghi (1994) and Chakraborty et al. (2008), which yield the following model expression for T 31 : where C 1 and C 2 are the model parameters. According to Swaminathan and Bray (2005) the order of magnitudes of T (i) 31 and T (ii) 31 are given by: It is evident from Eqs. 16i and ii that the order of magnitudes of the model expressions are consistent with the order of magnitude estimates of T (i) 31 and T (ii) 31 presented in Eq. 11i. Moreover, the model expression given in Eq. 15 is also consistent with the order of magnitude estimates presented in Eq. 11ii when the scaling arguments of Tennekes and Lumley (1972) and Mantel and Borghi (1994) are adopted: The predictions of the model given by Eq. 15 are compared to T 31 obtained from DNS data in Figs. 9a-f for cases A-F, respectively, and it can be seen that Eq. 15 satisfactorily models T 31 when C 1 = 1.0 and C 2 = 0.15 are considered.
The term T 32 is the dominant component of the scalar-turbulence interaction term T 3 (see Fig. 6 and scaling estimates given by Eq. 11ii) and thus its modelling is crucial to the modelling of the ̃c transport equation. The variations of T 32 with c across the flame brush for cases A-F are shown in Figs. 10a-f, respectively, which shows that T 32 assumes both positive and negative values within the flame brush. The relative alignments of ∇c and ∇ with local principal strain rates effectively determine the behaviour and sign of T 32 . It is well-known (Batchelor, 1959;Gibson, 1968;Kerr, 1985;Ashurst et al., 1987;Ruetsch and Maxey, 1991;Leonard and Hill, 1991;Nomura and Elgobashi, 1992) that the passive scalar gradient (e.g. ∇ ) preferentially aligns with the most compressive principal strain rate eigendirection. However, the reactive scalar gradient (e.g. ∇c ) exhibits preferential collinear alignment with the most compressive principal strain rate eigendirection when turbulent straining dominates over the strain rate arising from the flame normal acceleration, whereas a preferential alignment between the reactive scalar gradient and the most extensive principal strain rate eigendirection is obtained when the strain rate originating from flame normal acceleration overcomes turbulent straining (Chakraborty and Swaminathan, 2007a;Chakraborty et al., 2009;Malkeson and Chakraborty, 2011c). It was demonstrated earlier by Malkeson and Chakraborty 2 st with c along with the predictions of Eq. 17 with C 3 = 0.028 and C 4 = 0.01∕ 1 + Ka L 0.5 for cases A-F (a-f) (2011c) that ∇ continues to align with the most compressive principal strain rate eigendirection even when ∇c aligns with the most extensive principal strain rate eigendirection. Thus, the values and signs of T 32 are determined by the relative competition of the strain rates due to turbulence and flame normal acceleration. This also needs to be reflected in the model expression of T 32 and accordingly the following model expression is proposed here: where C 3 and C 4 are the model parameters and Da L =kS b ∕̃b is the local Damköhler number with b = 2D 0 ∕S b being a measure of the thermal flame thickness. The first term on the right-hand side (i.e. C 3 �̃∕k�̃c √̃c ∕ √ �̃c � ) accounts for the alignment of scalar gradients with the most compressive principal strain rate eigendirection under the action of turbulent straining (~ ̃∕k ∼ u � ∕ ). The second term on the right-hand side of Eq. 17 (i.e.
represents the effects of preferential collinear alignment of ∇c with the most extensive principal strain rate eigendirection under the action of strain rate induced by flame normal acceleration, which scales with 0 ∕ b − 1 (S b ∕ b ) (Chakraborty and Swaminathan, 2007a;Chakraborty et al., 2009) Swaminathan and Bray (2005). The effects of the strain rate are expected to decrease with increasing Karlovitz number (Chakraborty and Swaminathan, 2007b;Chakraborty et al., 2008Chakraborty et al., , 2011, and thus the model parameter C 4 may have some dependence on the local Karlovitz number . The predictions of Eq. 17 are shown in Figs. 10a-f where C 3 = 0.028 and C 4 = 0.01∕ 1 + Ka L 0.5 are used. For these choices of C 3 and C 4 , the model given by Eq. 17 reasonably captures the behaviour of T 32 extracted from DNS data for most cases considered here, but there are local discrepancies between the model predictions and DNS data. However, it is worthwhile to recognise that Eq. 17 is the first attempt to model the term T 32 in the ̃c transport equation and there is a scope for further improvement of the modelling of T 32 . It is also important to note that C 4 decreases with increasing Ka L and thus the second term on the right-hand side of Eq. 17 tends to weaken as the combustion regime approaches the broken reaction zone with the increase in the value of Karlovitz number. The variations of T 33 with c across the flame brush for cases A-F are shown in Fig. 11a-f respectively. Figure 11 indicates that T 33 assumes predominantly negative values for a major part of the flame brush, except for a small region on the burned and unburned gas sides of the flame brush in cases B and C, respectively (see Fig. 11b and c). For the statistically planar flames, T 33 takes the following form: Thus, the behaviour of T 33 is expected to be affected by ̃c ∼ D( c �� ∕ x 1 )( �� ∕ x 1 ) and ũ 1 ∕ x 1 , and therefore T 33 is modelled as: where C T 3 is a model parameter. The order of magnitude of Eq. 19 takes the following form according to Swaminathan and Bray (2005): st with c along with the predictions of Eq. 18 with C T 3 = 0.03 for cases A-F (a-f) A comparison between Eqs. 11i and 20 reveals that the model expression given by Eq. 19 is consistent with the order of magnitude of T 33 . Moreover, the modelled expression given by Eq. 19 is also consistent with the order of magnitude estimate predicted by Eq. 11ii, as shown below: The predictions of Eq. 19 are compared to T 33 extracted from DNS data in Figs. 11a-f for cases A-F, respectively when C T 3 = 0.03 . It can be seen from Fig. 11a-f that C T 3 = 0.03 yields satisfactory agreement with DNS data for all the cases considered here although Eq. 19 locally overpredicts the magnitude of the negative value at the middle of the flame brush in case A. It is worthwhile to note that Eq. 19 provides a simple expression as a first attempt towards the modelling of the ̃c transport equation, and therefore there will be further scope for improvement in the future.

Modelling of the Combined Reaction and Dissipation Contribution
The reaction rate and the molecular dissipation contributions of the ̃c transport are often modelled by their combined contribution (Mantel and Borghi, 1994;Chakraborty et al., 2008Chakraborty et al., , 2011. A similar approach has been adopted here following Malkeson and Chakraborty (2011b) for the ̃Y transport equation because T 4 and −D 2 may have larger magnitudes than that of T 4 − D 2 . This is particularly important because an imbalance created by the modelling inaccuracies of individual models for T 4 and −D 2 may not capture the qualitative behaviours of T 4 − D 2 . The variation of T 4 − D 2 with c has been found to be best qualitatively captured by the variation of the quantity ̃c throughout the flame brush in all cases. Moreover, T 4 − D 2 is dimensionally inconsistent with ̃c and a relevant time scale must be considered, which is estimated to be ̃∕ � ��2 , which yields the following model of T 4 − D 2 : Here, C D1 and C D2 are the model parameters.
The parameter The order of magnitude of the model expression given by Eq. 22 can be estimated in the following manner according to the scaling arguments proposed by Swaminathan and Bray (2005): It can be seen from Eqs. 22 and 11i the order of magnitude of the modelled expression is consistent with the order of magnitudes of T 4 and −D 2 . On scaling ̃ as: ̃∼ 1∕t D , where t D can be expressed as: t D ∼ 2 D ∕D , the right-hand side of Eq. 23 yields an order of magnitude estimate of u � 2 ∕l 2 Re 1∕2 t , which is consistent with the leading order scaling estimate of the combined contribution of (T 4 − D 2 ) as predicted by Eqs. 11ii.
The performance of the model given by Eq. 22 for all cases considered here is shown in Fig. 12a-f for cases A-F, respectively. It demonstrates that the model satisfactorily captures the behaviour and order of magnitude of (T 4 − D 2 ) for all cases when C D1 = 0.35 and C D2 = 0.4 are taken although this model (i.e. Equation 22) does not adequately capture (0.2m −c) and − 0.40 ̃c̃∕ � ��2 1 + (c∕{1 +c}) 2 m can be considered to be the modelled expressions of T 4 and −D 2 , respectively, which is not explicitly shown here for the sake of brevity.

Conclusions
A Direct Numerical Simulations (DNS) database of freely propagating statistically planar turbulent flames propagating into stratified mixture has been utilised in order to analyse the statistical behaviour of the cross-scalar dissipation rate ̃c transport. The statistical behaviours of the different terms of the ̃c transport equation have been explained based on physical arguments and the relative contributions of these terms to the overall ̃c transport have been explained based on detailed scaling arguments. The modelling of the unclosed terms of the ̃c transport equation has been considered in the context of Reynolds Averaged Navier-Stokes (RANS) simulations. It has been found that the density variation term T 2 , scalar-turbulence interaction term T 3 , reaction rate contribution T 4 and the molecular dissipation term −D 2 remain the leading order contributors for all cases considered here. It has been observed in all cases that the magnitude of the turbulent transport term T 1 remains small in comparison to the leading order contributors. Suitable model expressions have been identified for the contributions of the unclosed terms (i.e. T 1 , T 2 , T 3 , T 4 and −D 2 ) of the ̃c transport equation based on a priori analysis of DNS data. These model expressions are proposed in such a manner that the underlying physics and the order of magnitude of the unclosed term in question are adequately addressed. The predictions of the proposed models have been assessed with respect to the corresponding quantities extracted from DNS data. The new models are found to predict the unclosed terms of the ̃c transport equation satisfactorily for all the cases considered here. It is worthwhile to consider that this analysis is one of the first attempts to model the unclosed terms of the ̃c transport equation and thus there is some scope for improvement in the future for some of the model expressions. Admittedly, the closure of the unclosed terms of the ̃c transport equation involves a number of model parameters but the number of parameters used for the closures of the unclosed terms of the ̃c transport equation is comparable to those used for the closure of the ̃c (Mantel and Borghi, 1994;Chakraborty et al., 2008Chakraborty et al., , 2011 and ̃Y ) transport equations. Moreover, the effects of the differential diffusion rate of heat and mass and detailed chemical kinetics are not addressed in the current analysis. In order to attain a more comprehensive understanding of the cross-scalar dissipation rate transport, the effects of detailed chemistry and differential diffusion will need to be studied in detail in the future for higher values of turbulent Reynolds numbers. Undoubtedly, the models proposed here need to be implemented in RANS simulations. However, the performance of the ̃c modelling cannot be assessed in isolation because the modelling inaccuracies involved with k ,̃, ̃c and ̃ could potentially affect the predictions of RANS simulations in a complex manner and it is possible that some modelling inaccuracies will either augment each other or cancel with other. In a-priori analysis, such as the one carried out in this study, the true capability of the models is demonstrated but the interaction between the modelling error and numerical error (one of which originates due to discretisation error using a coarser RANS grid) warrants a-posteriori assessment of the models. It is important to note that only the model expressions for T 1 , T 31 and T 33 (i.e. Equations 12, 13, 15, 16) involve gradients of Favremean quantities which need to be evaluated on RANS grids. It is not possible to estimate the numerical error aspect due to discretisation error involving a coarser RANS grid in this type of a-priori analysis because the flame brush thickness scales with the integral length scale of turbulence and therefore the RANS grid size becomes comparable to the size of the simulation domain. Furthermore, the modelling and numerical errors interact not in a straightforward manner in RANS simulations. All of these aspects need to be accounted for a posteriori assessment of the cross-scalar dissipation rate closure based on RANS simulations for configurations in which experimental data is available for comparison with simulation results. This will also form the basis of future analyses on turbulent stratified mixture combustion.