Statistics of Scalar Flux Transport of Major Species in Different Premixed Turbulent Combustion Regimes for H2-air Flames

The statistical behaviour of turbulent scalar flux and modelling of its transport have been analysed for both major reactants and products in the context of Reynolds Averaged Navier Stokes simulations using a detailed chemistry Direct Numerical Simulation (DNS) database of freely-propagating H2 −air flames (with an equivalence ratio of 0.7) spanning the corrugated flamelets, thin reaction zones and broken reaction zones regimes of premixed turbulent combustion. The turbulent scalar flux in the cases representing the corrugated flamelets and thin reaction zones regimes of combustion exhibit predominantly counter-gradient transport, whilst a gradient transport has been observed for the broken reaction zones regime flame considered here. It has been found that the qualitative behaviour of the various terms of the turbulent scalar flux transport equation for the major species such as H2, O2 and H2O in the cases representing the corrugated flamelets and thin reaction zones regimes of combustion are mostly similar, whilst the behaviour is markedly different for the case representing the broken reaction zone regime. However, the terms for the scalar flux transport equation for H2 and O2 show same signs whereas the corresponding terms for H2O show signs opposite to those for H2 and O2. The performances of the well-established existing models for the unclosed terms of the turbulent scalar flux transport equation have been found to be similar for H2, O2 and H2O Some of the existing models for turbulent flux, pressure gradient, molecular diffusion and reaction contributions have been found to yield reasonable performance for the cases representing the corrugated flamelets and thin reaction zones regimes but the existing closures for these terms have been found to be mostly inadequate for the broken reaction zones regime flames.


Introduction
One of the major challenges of turbulent scalar transport modelling in the context of Reynolds Averaged Navier Stokes (RANS) simulations involves the closure of turbulent scalar flux components. The most widely used closure of turbulent scalar flux assumes a gradient hypothesis. According to the gradient hypothesis, the turbulent scalar flux components ρu i Y α for compressible turbulent flows are modelled as [1]: where ρ is the fluid density, u i is the i th component of fluid velocity, Y α is the scalar in question (here Y α is taken to be the mass fraction of species α), μ t is the eddy viscosity, Sc t = μ t /ρD t is the turbulent Schmidt number with D t being the eddy diffusivity and q = ρq/ρ and q = q −q are the Favre mean and fluctuations of a general variable q with the overline referring to the Reynolds averaging operation. The closure of ρu i Y α using the gradient hypothesis has well-known limitations [2][3][4][5][6][7][8]. For example, ρu i Y α and (−∂Ỹ α /∂x i ) are not collinearly aligned in turbulent channel flows. The difficulty in the modelling of ρu i Y α is further exacerbated by the fact that the turbulent scalar flux shows a counter-gradient behaviour in turbulent premixed flames when the flame normal acceleration dominates over turbulent fluctuation . Due to the limitation of algebraic closures of ρu i Y α , it is often necessary to solve a modelled transport equation of turbulent scalar flux components ρu i Y α in the context of second-moment closure and probability density function based modelling methodologies [31][32][33]. Several analyses concentrated on the closures of turbulent fluxes of ρu i Y α [34] and pressure gradient contributions [35][36][37] to the turbulent scalar flux transport for passive scalar mixing in non-reacting flows. However, the physics of turbulent scalar flux transport is considerably different in turbulent premixed flames due to the heat release arising from chemical reaction and self-induced pressure gradients. The modelling of the unclosed terms of the transport equation of ρu i Y α for turbulent premixed flames was discussed by Bray et al. [9] for the flamelet regime of combustion, and the modelled transport equation of turbulent scalar flux was used by Lindstedt and Vaos [31] and Tian and Lindstedt [32] in RANS simulations of laboratory scale burners. The models for the unclosed terms have been assessed based on a-priori Direct Numerical Simulation (DNS) analysis by Nishiki et al. [20] for flames representing the corrugated flamelets regime [38]. Chakraborty and Cant [23][24][25] demonstrated that the characteristic Lewis number has a significant influence on the statistical behaviours of the unclosed terms of the turbulent scalar flux transport equation based on a-priori DNS analyses. Furthermore, the analysis of Chakraborty and Cant [24] proposed modifications to the existing model expressions for the unclosed terms of the turbulent scalar flux transport equation for non-unity Lewis number and also for the thin reaction zones [37] combustion regime. In a subsequent analysis [27] the same authors addressed the effects of turbulent Reynolds number on the closures of the various terms of the turbulent scalar flux transport equation. Recently, Lai et al. [39] also analysed the nearwall behaviours of the unclosed terms of the turbulent scalar flux transport equation based on DNS data of head-on quenching of statistically planar turbulent premixed flames with different turbulence intensities and characteristic Lewis numbers. Based on this analysis Lai et al. [39] proposed near-wall modifications to the models of the unclosed terms of the turbulent scalar flux transport equation. In this regard, it is worthwhile to mention that the usage of turbulent scalar flux transport equation is rare in LES, whereas the turbulent scalar flux transport equation is used in the probability density function (PDF) method coupled with second-moment closure [31,32]. Thus, an analysis of the turbulent scalar flux transport in the context of Large Eddy Simulations (LES) will be of limited relevance. As a result, this analysis focuses primarily on RANS modelling.
To date, all the analyses 39] on the turbulent scalar flux transport in premixed combustion modelling have been carried out in the context of simple chemistry and transport. With the advancement of high-performance computing, it is now possible to incorporate detailed chemical mechanisms and solve transport equations for several species in engineering simulations [40,41]. It is not possible to translate the results for the turbulent scalar flux for one major species to any other major species by a linear transformation. The spatial distributions of each major species are different and thus the correlations between velocity and scalar fluctuations are different for every major species. For this reason, it is necessary to analyse the turbulent scalar fluxes for different major species in the context of multi-species RANS simulations. Moreover, it is necessary to ascertain if the same models for the unclosed terms in the transport equation of ρu i Y α work equally well for all the major species α in different regimes of combustion. The present analysis addresses this gap in the existing literature by the analysis of turbulent scalar fluxes ρu i Y α of H 2 , O 2 and H 2 O, and the well-established sub-models of the transport equation of ρu i Y α in the context of RANS using a three-dimensional Direct Numerical Simulations (DNS) database of H 2 -air flames with an equivalence ratio of 0.7 (which ensures that the flames remain globally thermo-diffusively neutral with respect to flame speed response to flame stretch [42]). The simulation parameters for this DNS database have been chosen in such a manner that the cases considered here represent typical combustion situations within the corrugated flamelets, thin reaction zones and broken reaction zones regimes of premixed turbulent combustion. The turbulent scalar fluxes ρu i Y α of H 2 , O 2 and H 2 O, and the transport equation of these scalar fluxes have been analysed using the aforementioned DNS database. Moreover, the modelling of the unclosed terms of the turbulent scalar flux transport equations will be discussed for the mass fractions of H 2 , O 2 and H 2 O. 1 In this respect, the main objectives of this paper are: The modelling attempts for the turbulent scalar flux transport in the context of turbulent reacting flows are relatively scarce [9,20,24,27,39] and all the existing analyses on turbulent premixed flames have been carried out for the flames in the flamelets regime using simple chemical mechanism [9,20]. Until now no model has been proposed for the low-Damköhler number conditions in the broken reaction zones regime. In this paper, the turbulent scalar flux transport has been addressed in the context of detailed chemical mechanism and non-unity Lewis number for the very first time. It is not to be expected that this analsyis will provide solutions to the modelling challenges that are open in the literature for decades and therefore the current paper should be treated as a first step towards achieving the goal of having a unified model of turbulent scalar flux transport for all regimes of premixed combustion for a multi-species system.
The rest of the paper will be organised as follows. The mathematical background and numerical implementation pertaining to the current analysis are provided in the next two sections. Following that, results are presented and discussed. The main findings are summarised and conclusions are drawn in the final section of this paper.

Mathematical Background
The transport equation of the Favre-averaged mass fractionỸ α = ρY α /ρ of species α is given by: whereω Y and D α are the chemical reaction rate of species α and diffusivity of species α, respectively. The last term on the right-hand side of Eq. 2 indicates turbulent transport of species α and one requires a model for the turbulent scalar flux ρu j Y α . Using the transport equations of Y α and u i , it is possible to derive a transport equation for ρu i Y α , which takes the following form [9,20,[23][24][25]39]: where τ ij = μ[(∂u i /∂x j ) + (∂u j /∂x i )] − (2μ/3)δ ij (∂u k /∂x k ) is the viscous stress tensor, P is the pressure and μ is the dynamic viscosity. The term T 1 is associated with turbulent transport of ρu i Y α , whereas T 2 and T 3 represent generation of turbulent scalar flux ρu i Y α by mean species and velocity gradients respectively. The terms T 4 and T 5 denote the contributions of mean and fluctuating pressure gradients respectively. The combined contributions of T 6 and T 7 is referred to as the molecular diffusion term. The term T 8 is the chemical reaction rate contribution to the turbulent scalar flux ρu i Y α transport. The statistical behaviours of T 1 − T 8 and their modelling will be discussed for turbulent scalar flux transports of H 2 , O 2 and H 2 O in Section 4 of this paper.

Numerical Implementation
In order to analyse the statistical behaviours of ρu i Y α , and T 1 − T 8 a three-dimensional detailed chemistry (9 species and 19 chemical reactions) DNS database of statistically planar H 2 -air flames with an equivalence ratio of 0.7 (i.e. φ = 0.7) was used [43]. The simulations have been conducted using a well-known DNS code [43], which employs high-order finite differences (8 th order at the internal grid points and gradually reducing to a one-sided 4 th order scheme at the non-periodic boundary) and a high order Runge-Kutta (4 th order) scheme for spatial discretisation and explicit time advancement respectively. The reacting scalar field has been initialised using the steady laminar H 2 -air flame solution corresponding to φ = 0.7. The thermo-physical properties are taken to be temperature dependent and are expressed according to CHEMKIN polynomials. The unburned gas temperature T 0 is taken to be 300K, which leads to an unstrained laminar burning velocity S L = 135.6cm/s and heat release parameter τ = (T ad − T 0 )/T 0 = 5.71 (where T ad is the adiabatic flame temperature) under atmospheric pressure. The numerical implementation for this database has been discussed elsewhere [43] in detail and thus a brief description is provided here. Turbulent inflow and outflow boundaries are taken in the direction of mean flame propagation, and transverse boundaries are considered to be periodic. The inflow and outflow boundaries are specified using an improved Navier Stokes characteristic boundary conditions (NSCBC) technique [44]. The inflow turbulent velocity fluctuations are specified by scanning a plane through a frozen field of turbulent homogeneous isotropic incompressible velocity generated using a pseudo-spectral method [45] following the Passot-Pouquet spectrum [46]. The inflow values of normalised root-mean-square turbulent velocity fluctuation u /S L , turbulent length scale to flame thickness ratio l T /δ th , flame thickness to the Kolmogorov length scale ratio δ th /η Damköhler number Da = l T S L /u δ th , Karlovitz number Ka = (ρ 0 S L δ th /μ 0 ) 0.5 (u /S L ) 1.5 (l T /δ th ) −0.5 and turbulent Reynolds number Re t = ρ 0 u l T /μ 0 for all cases are presented in Table 1 where ρ 0 and μ 0 are the unburned gas density and viscosity, respectively, δ th = (T ad − T 0 )/ max |∇T | L is the thermal flame thickness (with 'L' denoting unstrained laminar flame quantities). Cases A-C are representative of the corrugated flamelets (Ka < 1), thin reaction zones (1 < Ka < 100) and broken reaction zones (Ka > 100) regimes [38] of premixed turbulent combustion respectively. The domain size is 20mm × 10mm × 10mm (8mm × 2mm × 2mm) in cases A and B (case C) and the domain has been discretised by a uniform mesh of dimension 512 × 256 × 256 (1280 × 320 × 320). The mean inlet velocity has been adjusted to match the turbulent flame speed and the temporal evolution of the flame area has been monitored until a quasi-steady state is reached. The statistical stationary state has been achived at a time corresponding to 1.0l T /u , 6.8l T /u 6.7l T /u for cases A-C respectively and this simulation time remains comparable to several previous analyses [47][48][49]. For statistically planar flames, the directions normal to the mean direction of flame propagation (i.e. x 1 -direction) are statistically homogeneous. All the Reynolds/Favre averaging operations have been conducted by ensemble averaging the variables in the directions normal to the mean direction of flame propagation.

Flame-turbulence interaction
The isosurfaces of non-dimensional temperature c T = (T −T 0 )/(T ad −T 0 ) for cases A-C are shown in Fig. 1, which indicates that the flame morphologies in these cases are significantly different from each other and interested readers are referred to Refs. [43,50] for further discussion in this regard. The flame wrinkling can be quantified based on the normalised flame surface area A/A 0 where the flame surface area A is evaluated as: A = V |∇c|dV and A 0 is the initial value of flame surface area based on the one-dimensional steady state laminar flame solution. The quasi-steady state values of A/A 0 are 3.25, 5.0 and 3.25 for cases A, B and C respectively. The inlet turbulence intensity u /S L increases from case A to B, which leads to a greater extent of flame wrinkling in case B than in case A. However, l T /δ th values in cases A and C are considerably different and thus, the extent of flame wrinkling in cases A and C remain comparable despite large differences in u /S L values. It should be noted that Ka has not been modified here in isolation and thus the differences in behaviours of T 1 − T 8 originate not only due to the changes in the variation of Ka but also due to the changes in Da.

Statistical behaviour of turbulent scalar flux
In the case of statistically planar flames, ρu 1 Y α is the only non-zero component of turbulent scalar flux. The variations of ρu 1 Y α × 1/ρ 0 S L |Y αu − Y αb | (where subscripts u and b refer to the values in unburned and burned gases respectively and the multiplier is used for normalisation) withc T are shown in Fig. 2. The signs of ρu 1 Y α for α = H 2 and O 2 are the same but are different to ρu 1 Y α for α = H 2 O. Moreover, it can be seen from Fig. 2 that the qualitative behaviour of ρu 1 Y α for α = H 2 , O 2 and H 2 O does not change in cases A and B but the behaviour in case C is completely different from that in cases A and B. In order to understand this observation the variation of is indicative of a gradient (counter-gradient) transport. Furthermore, for a gradient transport ρu 1 Y α /(−∂Ỹ α /∂x 1 ) provides the density-weighted eddy diffusivity ρD t , whereas an unphysical negative ρD t is obtained for a counter-gradient transport. Figure 3 indicates that a predominantly counter-gradient behaviour has been observed for turbulent scalar fluxes of H 2 , O 2 and H 2 O for cases A and B, whereas a gradient type transport is obtained for these species in case C.
It is important to note that the results in Fig. 3 do not necessarily mean that countergradient (gradient) transport will always be obtained for the flames within the corrugated flamelets and thin reaction zones regimes (broken reaction zones regime). According to Veynante et al. [12], a gradient transport is obtained for N B < 1, whereas a counter-gradient transport is favoured for N B > 1, where N B = τ S L /2αu is the Bray number with α being an efficiency function. The value of α is about 0.5 according to Veynante et al. [12] and thus τ S L /u can be taken to provide a measure of the Bray number N B . The Bray number criterion proposed by Veynante et al. [12] is valid in an order of magnitude sense. The velocity jump due to flame normal acceleration arising from chemical heat release is taken to scale with τ S L , whereas the velocity jump due to turbulence scales with u . A countergradient transport is obtained when τ S L dominates over u , whereas a gradient transport is obtained when u dominates over τ S L . In cases A and B, τ S L /u remains greater than unity (= 8.16 for case A and 1.14 for case B) whereas it is smaller than unity (= 0.41) for case C and accordingly a counter-gradient behaviour has been observed for cases A and B whilst a gradient transport is obtained for case C.

Statistical behaviour of the terms of the turbulent scalar flux transport equation
The variations of T 1 − T 8 withc T in cases A-C are shown in Fig. 4 for α = H 2 , O 2 and H 2 O so that the relative magnitudes of these terms for different species can be compared for flames belonging to the different combustion regimes. The modelling of the terms of the turbulent scalar flux transport equation will be discussed in the subsequent sub-sections. The signs of these terms for H 2 O are different from those of H 2 and O 2 for all cases because H 2 O is a product species, whereas H 2 and O 2 are the reactants. It is worth noting that a positive (negative) contribution of T 1 − T 8 acts to produce a counter-gradient (gradient) transport of turbulent scalar flux ρu 1 Y H 2 O . By contrast, a negative (positive) value of T 1 − T 8 promotes a counter-gradient (gradient) transport of turbulent scalar fluxes of H 2 and O 2 (i.e. ρu 1 Y H 2 and ρu 1 Y O 2 ). The mean and fluctuating pressure gradient terms T 4 and T 5 play leading order roles in cases A and B. For case C, T 4 and T 5 remain comparable to the magnitudes of the contributions of T 1 , T 2 , T 6 and T 8 . The flame normal acceleration sets up a negative mean pressure gradient (i.e. ∂P /∂x 1 < 0 when the mean direction of flame propagation is in the negative x 1 -direction), which gives rise to negative (positive) values of T 4 for reactants (products) because of negative (positive) values of Y α . A negative (positive) fluctuation of Y α for reactants (products) tends to induce a more negative pressure gradient and thus T 5 = −Y α (∂P /∂x 1 ) assumes mostly positive (negative) values for α = H 2 and O 2 (α = H 2 O). The terms due to mean species and velocity gradients T 2 and The molecular diffusion terms T 6 and T 7 play marginal roles in cases A and B but T 6 plays a leading order role in case C, especially for the turbulent scalar flux transport of H 2 . It can be seen from Fig. 4 that the combined action of T 6 and T 7 acts to reduce the extent of counter-gradient transport in cases A and B, whereas it acts to reduce the extent of gradient transport in case C, and this behaviour remains unchanged for all species considered here. The combined contribution of T 6 and T 7 can be split into a molecular diffusion contribution (∝ ∇·(ρD∇(u 1 Y α )) and a dissipation contribution (∝ −2ρD∇Y α · ∇u 1 ). For high values of Re t the molecular dissipation contribution dominates over the diffusion contribution (i.e. |∇·(ρD∇(u 1 Y α )))|<<| − 2ρD∇Y α · ∇u 1 |). For a counter-gradient transport, the quantity 2ρD∇Y α · ∇u 1 assumes positive (negative) values for products (reactants) and thus the contribution of (T 6 + T 7 ) acts to reduce the extent of counter-gradient transport in cases A and B. By contrast, 2ρD∇Y α · ∇u 1 assumes negative (positive) values for products (reactants) for a gradient transport, and thus it acts to reduce the extent of gradient transport in case C.
The reaction rate contribution T 8 assumes positive (negative) values towards the unburned gas side of the flame brush, before becoming negative (positive) on the burned gas side of the flame brush for H 2 and O 2 (H 2 O) for all cases, but for cases A and B this term remains small in comparison to the magnitudes of T 4 and T 5 , whereas its magnitude is comparable to T 1 , T 2 , T 4 , T 5 and T 6 in case C. An increase (decrease) in reaction rate magnitude on the unburned (burned) gas side tends to induce a negative u 1 due to the decay of turbulence under the enhanced viscous action, and this leads to positive (negative) values towards the unburned gas side of the flame brush and negative (positive) values of T 8 on the burned gas side of the flame brush for reactants such as α = H 2 and O 2 (for products such as α = H 2 O) for all cases. However, this action is relatively stronger in case C than in cases A and B because flame-generated turbulence effects are stronger in these cases. The flame-generated turbulence acts to locally enhance the turbulence level in cases A and B, which counters the decay of turbulence with increased chemical activity and thus the relative magnitude of T 8 is smaller than the leading order contributions of T 4 and T 5 in these cases in comparison to case C.
It is worth noting that the observed behaviours from Fig. 4 are consistent with the scaling analyses presented in Chakraborty and Cant [27] and Lai et al. [39], which indicate that T 2 and T 4 are expected to play dominant roles in the turbulent scalar flux transport and T 3 is likely to be of marginal importance for high Reynolds number flames with small values of Da (i.e. Da < 1).
The terms T 2 and T 3 are closed terms in the context of second-moment closure because ρu i u j and ρu j Y α are already modelled in this framework. The terms T 1 , T 4 , T 5 , T 6 , T 7 and T 8 are unclosed and need closures in order to solve Eq. 2. This is often achieved by modelling ρu j u i Y α , (T 4 + T 5 ), (T 6 + T 7 ) and T 8 , which will be discussed next in this paper.

Modelling of turbulent flux of scalar flux ρu
The quantity ρu 1 u 1 Y α is the only non-zero component of ρu j u i Y α for a statistically planar flame. The turbulent flux of ρu i Y α is usually modelled for non-reacting flows in the following manner utilising the gradient hypothesis [34]: The first term on the right-hand side represents the reacting contribution to ρu 1 u 1 c , whereas the combined action of second and third terms represent the effects of turbulence on ρu 1 u 1 c . The last term O(1/Da) originates from the interior of the flame, and this contribution becomes negligible for Da 1. Chakraborty and Cant [24,27] proposed an alternative model of ρu 1 u 1 Y α by incorporating the reacting contribution according to the Bray-Moss-Libby (BML) analysis and the turbulent contribution in Eq. 5i is accounted for by the TDH model. The model by Chakraborty and Cant [24,27] is given by: Equation 5ii will henceforth be referred to as the CC model in this paper. The second term on the right hand side of Eq. 5ii is capable of predicting counter-gradient transport and includes the effects of flame normal acceleration due to chemical heat release (i.e..

−ρu
5ii accounts for the velocity jump across the flame brush with the first term representing the effects of heat release, whereas the second term represents the effects of turbulent velocity fluctuations).
The predictions of the TDH and CC models are compared to ρu 1 u 1 Y α extracted from DNS data in Fig. 5  O, which implies that ρu 1 u 1 Y α exhibits counter-gradient behaviour in these cases. In cases A and B, the CC model has been found to accurately predict the quantitative behaviour of ρu 1 u 1 Y H 2 and ρu 1 u 1 Y O 2 on the unburned gas side of the flame brush, but does not predict the correct trend on the burned gas side of the flame brush and its prediction becomes comparable to the TDH model. This implies that the magnitude of the second term on the right hand side of Eq. 5i-5iii becomes It is worth noting that the TDH model was originally proposed for non-reacting flows using the gradient hypothesis, and thus does not include the effects of flame normal acceleration due to chemical heat release and therefore is not equipped to predict the counter-gradient behaviour. The effects of heat release remain significant even for case C in spite of large values of Ka and thus, it is perhaps not surprising that this model does not perform satisfactorily for all species in all cases considered here. This behaviour has been found to be consistent with previous analyses based on simple chemistry DNS results [24,27,39]. This implies that the choices of species and chemical mechanism do not have significant influence on the agreement of the TDH model with ρu 1 u 1 Y α extracted from DNS data.
The functional form of the second term on the right hand side of Eq. 5ii is derived based on a presumed bi-modal distribution of Y α with peaks at Y αu and Y αb , which is strictly valid for Da > 1 and Ka < 1. Thus, the CC model works relatively satisfactorily for case A (where Da > 1 and Ka < 1) but the probability density function of Y α does not remain bimodal for Ka > 1 and thus its performance worsens progressively with increasing Ka (i.e. from case B to case C). Some recent analyses [51,52] focussed on modelling of the PDF of c where non bi-modal distribution is realised.

Modelling of the pressure gradient terms (T 4 +T 5 )
The variations of the normalised values of (T 4 + T 5 ) withc T in cases A-C are shown in Fig. 6 for α = H 2 , O 2 and H 2 O. It can be seen from Fig. 6 that the net contribution of   [9] and even though this expression is not strictly valid for Da < 1 it can still be used in a qualitative sense. As Y αu > Y αb (Y αu < Y αb ) for reactants (products) and , the quantity Y α assumes negative (positive) values for reactants (products). Chakraborty and Cant [24] proposed Y α = −ρỸ 2 α τ/ρ 0 (Y αu − Y αb ) for both Da > 1 and Da < 1 combustion, and this expression remains valid for all cases considered here for α = H 2 , O 2 and H 2 O. Thus, the negative (positive) values of Y α for the major reactant (product) species lead to negative (positive) values of T 4 = −Y α (∂P /∂x 1 ) (see Fig. 4). A comparison between Figs. 4 and 6 reveals that although the fluctuating pressure gradient term T 5 locally assumes values with a different sign to that of T 4 , the net contribution of (T 4 + T 5 ) follows the sign of the mean pressure gradient term T 4 .
The mean and fluctuating pressure gradient terms are often modelled together because of their similar origin [2]. Several models are available in the existing literature for (T 4 + T 5 ) (see Refs. [23,27,39]). Some of the models for (T 4 + T 5 ), which were originally proposed for non-reacting flows, take the following form [2]: where C 1c , C 2c , C 3c and C 4c are the model parameters. Launder [35] suggested that C 1c = 3.0, C 2c = 0, C 3c = 0 and C 4c = 0.4, and this model will henceforth be referred to as the PL model. Craft [36] adopted a similar model (referred to as the PC model) with C 1c = 3.0, C 2c = 0.5, C 3c = 0 and C 4c = 0. An alternative model (PD model) was suggested by Durbin [37] where C 1c = 2.5, C 2c = 0, C 3c = 0 and C 4c = 0.45. Jones [53] and Bradley et al. [54] modelled (T 4 + T 5 ) in the following manner: where C φ1 = 3.0 andC φ2 = 0.5 are taken for the model by Jones [53] (PJ model) , whereas C φ1 = 3.0 and C φ2 = 0 are considered for the model by Bradley et al. [54] (PB model). Lindstedt and Vaos [31] proposed another alternative model (PLV model) as: where C As = 1/3 and G il is the generalised Langevin coefficient which is a function of Reynolds stress ρu i u j and the mean velocity gradient ∂ũ i /∂x j [31]. It is worth noting that Y α was evaluated using the BML relation: in Refs. [31,54,55] but for the current a-priori analysis Y α is extracted from DNS data. In all cases considered here, Y α can be modelled as: Nishiki et al. [20] also proposed a model (PN model) based on a-priori simple chemistry DNS analysis for flames in the corrugated flamelets regime in the following manner: (9) where C D = 0.8, C E1 = 0.38 and C E2 = 0.66 are the model constants. The first term on the right hand side of the PN model originates from the BML relation for T 4 [9].
The predictions of all the aforementioned models are compared to (T 4 + T 5 ) extracted from DNS data in Fig. 6 for α = H 2 , O 2 and H 2 O. Figure 6 shows that the PL, PC and PD models fail to capture both qualitative and quantitative behaviours of (T 4 + T 5 ) in cases A and B irrespective of the choice of α. By contrast, these models predict the qualitative behaviour of (T 4 + T 5 ) satisfactorily in case C. These models overpredict the magnitude of (T 4 + T 5 ) towards the unburned gas side of the flame brush, whereas quantitative agreement with DNS data remains reasonable towards the burned gas side of the flame brush. The PL, PC and PD models were originally proposed for incompressible non-reacting flows [2,[34][35][36] where the contribution of T 4 = −Y α ∂P /∂x i is identically zero and thus its contribution was ignored. It can be seen from Fig. 4 that the contribution of T 4 remains significant for the cases considered here irrespective of the choice of α. A comparison of cases A-C in Fig. 4 shows that the contribution of T 4 in comparison to T 5 diminishes from case A to case C, and thus the PL, PC and PD models are found to be more successful in capturing the statistical behaviour of (T 4 + T 5 ) in case C. The case C belongs to the broken reaction zones regime and thus it shows some attributes of non-reacting flows, which is also reflected in the more satisfactory performance of the PL, PC and PD models than in cases A and B.
The PJ, PB and PLV models exhibit very similar behaviour. All three models capture the qualitative behaviour of (T 4 +T 5 ) in case A for α=H 2 , O 2 and H 2 O, but these models underpredict (overpredict) the magnitude of (T 4 +T 5 ) towards the unburned (burned) gas side of the flame brush. In case B, these models fail to predict the behaviour of (T 4 +T 5 ) on the unburned gas side of the flame brush, but capture both quantitative and the qualitative behaviours on the burned gas side of the flame brush. The PJ and PB models perform satisfactorily in case C for α=H 2 , O 2 and H 2 O but the PLV model does not capture the qualitative and quantitative behaviours of (T 4 +T 5 ) in this case for all the major species considered here.
The PN model captures the qualitative behaviour of (T 4 +T 5 ) better than the other alternative models in case A but overpredicts the magnitude. In case B, the PN model captures the qualitative behaviour, but underpredicts the magnitude on the unburned gas side and overpredicts the magnitude on the burned gas side of the flame brush. However, in case C, the PN model fails to predict the qualitative behaviour and the magnitude of (T 4 +T 5 ) on the unburned gas side, whereas it accurately captures the behaviour on the burned gas side of the flame brush but overpredicts the magnitude of (T 4 +T 5 ). The PN model was originally proposed for strict flamelet combustion (i.e. Da > 1 and Ka < 1) and these assumptions are rendered invalid for the broken reaction zones regime in case C and thus this model does not perform satisfactorily in this case for all species considered here.
It is worth noting that the PL, PC, PD models, which do not account for the leading order contribution of T 4 = −Y α ∂P /∂x i , are not successful in capturing (T 4 +T 5 ) extracted from DNS data. However, the PJ, PB and PN models, which include T 4 = −Y α ∂P /∂x i are more successful in capturing the behaviour of (T 4 + T 5 ) extracted from DNS data than the PL, PC, PD models. The model parameter and the model expression for the PLV model have been calibrated for the flamelet regime of combustion and thus it is perhaps unsurprising that this model does not perform satisfactorily in case C representing the broken reaction zones regime. The first term on the right hand side of the PN model (9) assumes a bimodal distribution of c, which is not realised in cases B and C and thus the prediction of the PN model exhibits some inaccuracies especially in these cases. It is possible that the methodologies, which parameterise the pdf of c when the presumed bi-modal pdf with impulses at c = 0 and c = 1.0 is not realised, can be more successful in deriving improved models for (T 4 +T 5 ) for non-flamelet regime of combustion.

Modelling of the molecular diffusion terms (T 6 +T 7 )
The variations of normalised (T 6 + T 7 ) withc T in cases A-C are shown in Fig. 7 for α = H 2 , O 2 and H 2 O. It can be seen from a comparison between Figs. 2 and 7 that (T 6 + T 7 ) assumes a sign which is opposite to ρu 1 Y α for all species α in all cases. These terms tend to oppose the dominant behaviour of the turbulent scalar flux ρu 1 Y α . In cases A and B, (T 6 + T 7 ) assumes positive (negative) values for α = H 2 and O 2 (α = H 2 O), whereas ρu 1 Y α exhibits negative (positive) values throughout the flame brush [9,12,24,27,39]. By contrast, in case C, (T 6 + T 7 ) assumes negative (positive) values for α = H 2 and O 2 (α = H 2 O), whereas ρu 1 Y α exhibits a positive (negative) value throughout the flame brush.
According to Bray et al. [9] (T 6 + T 7 ) is modelled as: whereω Y α is the chemical reaction rate of α, and K 1 = 0.85 is the model constant. The model given by Eq. 10 will henceforth be referred to as the DBML model. An alternative model (i.e. DN model) was proposed by Nishiki et al. [20] in the following manner: where C F = 0.4 is a model constant. It is worthwhile to note that both Eqs. 10 and 11 are directly proportional to the mean reaction rateω Y α and thus these models predict non-zero values of (T 6 +T 7 ) within the flame brush. The predictions of the DBML, DN and DC models are compared to (T 6 + T 7 ) extracted from DNS data in Fig. 7. It can be seen from Fig. 7 that both DBML and DN models satisfactorily capture the qualitative behaviours of (T 6 + T 7 ) for α = H 2 , O 2 and H 2 O. However, both DBML and DN models overpredict the magnitudes of (T 6 + T 7 ) in cases A and B but the extent of overprediction is relatively smaller in the case of the DN model.
The DBML analysis was originally proposed for the strict flamelet regime (i.e. Da > 1 and Ka < 1) but even then this model significantly overpredicts (T 6 + T 7 ) in case A where the conditions in terms of Da and Ka for which the model was proposed are satisfied (using K 1 = 0.2 instead of 0.85 provides satisfactory quantitative agreement with DNS data). The same can be said for the DN model where C F = 0.2 instead of 0.4 yields good quantitative agreement with DNS data. Thus, the application of these models for case C (where Da < 1 and > 1) is beyond the scope of their validity. In spite of the above limitation, the DBML model predicts the correct sign of (T 6 + T 7 ) for α = H 2 , O 2 and H 2 O in case C and the quantitative agreement with DNS data also remains reasonable. The choice of characteristic ] in the DBML model may not be appropriate for low-Damköhler number (i.e. Da < 1) combustion in case C and this might be one of the reasons behind the overprediction of the DBML model in this case. The DN model is strictly valid only for counter-gradient transport and it predicts wrong sign for gradient transport [24,27,39]. Thus, the DN model fails to predict the correct qualitative behaviour of (T 6 + T 7 ) for α = H 2 , O 2 and H 2 O in case C where a gradient transport is obtained (see Fig. 4).
The optimum values of K 1 for different cases for α = H 2 , O 2 and H 2 O are listed in Table 3. As the DN model cannot predict the correct sign in the case of gradient transport, the optimum values C F for cases A-C have not been reported in Table 3 Table 3 suggests that it is likely to be dependent on Karlovitz and Lewis numbers. Based on this, the optimum value of K 1 has been parameterised as: According to this parameterisation, K 1 reaches an asymptotic for large values of Ka L . The predictions of the DBML model with K 1 according to Eq. 12 (i.e. DBML-M model) are shown in Fig. 7, which shows that Eq. 12 sigfnificantly improves the model performance and the DBML-M model satisfactorily captures both qualitative and quantitative behaviours of (T 6 + T 7 ). It is worth noting that Eq. 12 provides a possible parameterisation of K 1 and an alternative parameterisation exhibiting similar behaviour is also possible.

Modelling of the reaction rate term T 8
The reaction rate contribution to the turbulent scalar flux transport was modelled by Bray et al. [9] for strict flamelet combustion (i.e. Da > 1 and Ka < 1) in the following manner (i.e. RB model): where C R = 1.0 and φ m = 0.5 are the model constants. The variations of normalised T 8 withc T for cases A-C are shown in Fig. 8 for α = H 2 , O 2 and H 2 O along with the predictions of the RB model. Figure 8 shows that in cases A and B, the term T 8 assumes negative values towards the unburned gas side and positive values on the burned gas side of the flame  . This implies that the correlation between the fluctuations of velocity and reaction rate is fundamentally different in case C than in cases A and B. In the case of counter-gradient transport, a positive fluctuation of velocity induced by an increase in reaction rate magnitude tends to produce negative (positive) values of T 8 for reactant (product) species such as H 2 and O 2 (H 2 O) and this is predominantly responsible for the observed behaviours of T 8 towards the unburned gas side and middle of the flame brush in cases A and B. However, the reaction rate magnitude tends to decrease towards the burned gas side whereas temperature continues to rise which acts to increase the positive fluctuations of velocity due to thermal expansion in the case of predominant counter-gradient transport. This leads to negative (positive) values of T 8 for reactant (product) species such as H 2 and O 2 (H 2 O) and this is predominantly responsible for the observed behaviours of T 8 towards the burned gas side of the flame brush in cases A and B. Due to a predominantly gradient transport in case C, an increase in reaction rate magnitude tends to induce negative fluctuations of velocity, which leads to positive (negative) values of T 8 for reactant (product) species such as H 2 and O 2 (H 2 O) towards the unburned gas side and middle of the flame brush. In case C, the reaction rate magnitude tends to decrease towards the burned gas side whereas the rising temperature acts to decrease velocity due to augmented viscous damping, which is reflected in the negative (positive) values on the burned gas side of the flame brush for α = H 2 and O 2 (α = H 2 O). It can be seen from Fig. 8 that the RB model captures the correct qualitative behaviour of T 8 for α = H 2 , O 2 and H 2 O in cases A and B. Although the quantitative agreement of the RB model with DNS data is not perfect, it is reasonable for these cases and this level of agreement has been found to be consistent with previous findings based on simple chemistry DNS data [24,27,39]. The qualitative and quantitative agreement between the RB model and DNS data is comparatively less satisfactory in case C in comparison to cases A and B. It is worth noting that the RB model was originally proposed for Da > 1 and Ka < 1, and thus the performance of this model progressively worsens with increasing (decreasing) Ka (Da).
It has been found that the performance of the RB model is not significantly dependent on φ m but on C R . The optimum values of C R for different cases for α = H 2 , O 2 and H 2 O are listed in Table 4, which shows that C R decreases from case A to case C and its value is different for different species. Moreover, optimum values of C R for α = H 2 are greater than the corresponding optimum values for α = O 2 and H 2 O. The optimum values of C R for α = O 2 and H 2 O are similar in magnitude. Thus, it can be inferred that the optimum values of C R is dependent on Karlovitz and Lewis numbers, which can be parameterised as: The predictions of the RB model with φ m = 0.5 and C R according to Eq. 14 (i.e. RB-M model) are also shown in Fig. 8, which shows that Eq. 14 significantly improves the model performance and the RB-M model satisfactorily captures both qualitative and quantitative behaviours of T 8 . It is worth noting that Eq. 14 provides a possible parameterisation of C R and an alternative parameterisation exhibiting similar behaviour is also possible.

Final Remarks on the Turbulent Scalar Flux Transport Equation
A summary of the models considered for this analysis is presented in Table 5. Based on the foregoing discussion, the optimal combinations of the closure models for the unclosed terms of the turbulent scalar flux transport equation for different combustion regimes for different species (e.g. Eqs. 5iii, 12 and 14 are dependent on Le α ) are summarised in Table 6 for the convenience of readers and future users of the turbulent scalar flux transport equation. Table 6 provides an idea about the appropriate models for major species in different combustion regimes. It can be appreciated from the information provided in Table 6 that there is a huge scope for improvement in the modelling of the turbulent scalar flux transport in the broken reaction zones regime and this aspect needs further investigation. Further, it becomes clear that for most unclosed terms there is not a single existing model that performs satisfactorily in all regimes of combustion and for all species (see Tables 2-4). The suggested empirical relations accounting for some of these effects for CC, DBML and RB models will require further investigation.

Conclusions
The statistical behaviours of the turbulent scalar flux and the terms of its transport equation for major species have been analysed in the context of RANS using a detailed chemistry DNS database of freely-propagating H 2 −air flames with an equivalence ratio of 0.7 spanning the corrugated flamelets, thin reaction zones and broken reaction zones regimes of combustion. The turbulent scalar flux statistics and its transport have been analysed in detail for the major reactants and products (i.e. H 2 , O 2 and H 2 O). A counter-gradient transport has been observed for the cases considered here representing the corrugated flamelets and thin reaction zones regimes of combustion, whereas a gradient transport is observed for the case representing the broken reaction zones regime. Accordingly, the qualitative behaviours of the terms of the turbulent scalar flux transport equation remain similar for the flames representing the corrugated flamelets and thin reaction zones regimes but these behaviours are different to that observed for the broken reaction zones regime flame considered here. It has been found that the performances of the existing closures for turbulent transport, pressure gradient, molecular diffusion and reaction rate terms show some dependence on the choice of the major species. The models for the unclosed terms of the turbulent scalar flux transport equation which perform satisfactorily in the corrugated flamelets and thin reaction zones regimes of premixed combustion have been identified based on a-priori DNS analysis. However, there is no existing modelling methodology which was originally developed for the broken reaction zones combustion and thus the existing closures for the unclosed terms of the turbulent scalar flux transport equation, which were originally proposed for the strict flamelet combustion, have been found to be mostly inadequate for the broken reaction zones regime. Detailed explanations have been provided for the observed performances of sub-models for the unclosed terms of the turbulent scalar flux transport equation for different regimes of premixed turbulent combustion. The present analysis indicates that there is ample scope for improvement in the modelling of the unclosed terms in the turbulent scalar flux transport equation and this improvement is especially needed for the broken reaction zone regime combustion.