Analysis of the Combined Modelling of Sub-grid Transport and Filtered Flame Propagation for Premixed Turbulent Combustion

Flame surface density (FSD) based reaction rate closure is an important methodology of turbulent premixed flame modelling in the context of Large Eddy Simulations (LES). The transport equation for the Favre-filtered reaction progress variable needs closure of the filtered reaction diffusion imbalance (FRDI) term (i.e. filtered value of combined reaction rate and molecular diffusion rate) and the sub-grid scalar flux (SGSF). A-priori analysis of the FRDI and SGSF terms has in the past revealed advantages and disadvantages of the specific modelling attempts. However, it is important to understand the interaction of the FRDI and SGSF closures for a successful implementation of the FSD based closure. Furthermore, it is not known a-priori if the combination of the best SGSF model with the best FRDI model results in the most suitable overall modelling strategy. In order to address this question, a variety of SGSF models is analysed in this work together with one well established and one recent FRDI closure based on a-priori analysis. It is found that the success of the combined FRDI and SGSF closures depends on subtle details like the co-variances of the FRDI and SGSF terms. It is demonstrated that the gradient hypothesis model is not very successful in representing the SGSF term. However the gradient hypothesis provides satisfactory performance in combination of a recently proposed FRDI closure, whereas unsatisfactory results are obtained when used in combination with another existing closure, which was shown to predict the FRDI term satisfactorily in several previous analyses.


Introduction
The complexity of the system of partial differential equations describing turbulent reactive flows can be simplified assuming single-step chemistry and a unity Lewis number (=thermal diffusivity/mass diffusivity). The mass fractions of the reactive species and the non-dimensional temperature in a premixed flame can be utilised to define a reaction progress variable c assuming the values c= 0 on the reactant side and c = 1 in the fully burned products [1]. The Large Eddy Simulation (LES) approach is considered a promising technique in order to improve the prediction quality of computational fluid dynamic simulations. The LES filtering operation of a quantity Q with a Gaussian filter kernel G (r) is given as: In Eq. 1, is the filter width andQ = ρQ/ρ refers to the well-known Favre filtering operation. In this framework, the transport equation for the filtered reaction progress variable becomes: ∂ (ρc) /∂t+∇· (ρũc) = −∇ · (ρuc −ρũc) + ∇ · (ρD∇c) +ω c . ( Here ρ, u, D andω c denote density, velocity, progress variable diffusivity and reaction rate respectively. The numerical solution requires closure for both terms on the right hand side of Eq. 2, i.e. for the sub-grid scalar flux (SGSF) denoted as: T SGSF i := ρu i c − ρũ ic (where u i is the i th component of velocity) as well as the filtered reaction diffusion imbalance (FRDI) term T F RDI = ∇ · (ρD∇c) +ω c . In the context of LES of premixed flames, turbulent scalar flux modelling has been investigated by a number of authors in the past [2][3][4][5][6][7][8]. A variety of concepts exist for turbulent premixed combustion modelling (e.g. Gequation [9], artificially thickened flame concept [10], Flame Surface Density (FSD) closure [2,[11][12][13]). In this work we focus on the FSD concept. In this context T F RDI is expressed as: where S d is the displacement speed of a given c isosurface and =|∇c|/ |∇c| is the wrinkling factor. The generalised FSD is defined as gen = |∇c| and the surface weighted filtering operation (·) S is given by (Q) S = Q |∇c|/|∇c| (see [2]). Interested readers are referred to [11,12] for examples of LES implementations of FSD based reaction rate closures. A-priori analysis of both FRDI and SGSF terms individually was reported by the authors in the past [13][14][15]. Recent a-posteriori analyses [16,17] indicated complex interaction of the closures of T F RDI and T SGSF i , which is the main focus of this work. The rest of the paper is organised in the following manner. The details related to the DNS database and its filtering will be explained in the next section. Following this, the closures of T F RDI and T SGSF i used in this work will be briefly summarized. This will be followed by the analysis of the models on term by term basis in order to understand how they interact. Finally conclusions will be drawn.

DNS Database
A DNS database of turbulent premixed flames under decaying turbulence with single step Arrhenius type irreversible chemistry has been considered for the current analysis, consisting of five flames with global Lewis number Le = 0.34 (case A), 0.6 (case B), 0.8 (case C), 1.0 (case D) and 1.2 (case E). The initial values of the normalised turbulent rootmean-square (rms) velocity fluctuation u /S L = 7.5, integral length scale to thermal flame thickness ratio l/δ th =2.45, Damköhler number Da = lS L /δ th u =0.33, and Karlovitz number Ka = u /S L 3/2 (l/δ th ) −1/2 = 13.0 are taken to be the same for all cases considered here where S L is the unstrained laminar burning velocity, δ th = (T ad −T 0 )/ max |∇T | L is the thermal flame thickness with T ad and T 0 being the adiabatic flame temperature and the reactant temperature respectively. Note that the subscript 'L' refers to the unstrained laminar flame quantities. The heat release parameter τ = (T ad − T 0 ) /T 0 and the Zel'dovich number β = T ac (T ad − T 0 ) /T 2 ad are taken to be 4.5 and 6.0 respectively where T ac is the activation temperature. Standard values of Prandtl number (P r = 0.7) and ratio of specific heats (γ g = 1.4) have been used.
The simulation domain is taken to be a cube of 24.1δ th × 24.1δ th × 24.1δ th which is discretised using a uniform Cartesian grid of 230×230×230 points ensuring 10 grid points are kept within δ th . Spatial derivatives for all internal grid points are evaluated using 10th order central differences but the order of discretization gradually drops to an one-sided 2 nd order scheme at the non-periodic boundaries. Time integration is carried out using an explicit 3rd order low storage Runge-Kutta scheme. The boundary conditions in the mean flame propagation direction (aligned with negative x 1 -direction) are taken to be partially non-reflecting, whereas boundaries in transverse directions are taken to be periodic. The turbulent velocity fluctuations are initialised using a homogeneous isotropic incompressible velocity field. The reacting flow field is initialised by a steady planar unstrained premixed laminar flame solution. In all cases flame-turbulence interaction takes place under decaying turbulence and all non-dimensional numbers mentioned earlier have to be understood as initial values here and in the remainder of the text. The simulations have been carried out for one chemical time scale (i.e. t chem = δ th /S L ), which corresponds to 3.34 initial eddy turn over times (i.e. 3.34t e = 3.34l/u ). By this stage, u /S L decayed by 50 %, whereas l/δ th increased by a factor of 1.7 in the unburned gas ahead of the flame. The DNS database considered here was used previously in several studies and interested readers are referred to references [18][19][20][21][22][23][24] for more details.
The values of u /S L , l/δ th and Re t used for the current analysis remain comparable to several previous analyses [2,10,[25][26][27] which were used for a-priori DNS modelling. Furthermore, FSD and SDR models proposed based on a-priori DNS analyses using this database [28][29][30] have been found to be in good agreement with a-posteriori assessments based on actual LES simulations [16,[31][32][33]. The aforementioned facts provide the confidence in the findings of the present analysis.
For the purpose of the a-priori analysis carried out in this work, the DNS data has been explicitly filtered using the Gaussian filter kernel given by Eq. 1. Results will be presented from ≈ 0.4 δ th where the flame is almost resolved, up to ≈ 2.8 δ th where the flame becomes fully unresolved and is comparable to the integral length scale l. The result of the filtering operation is a dataset with the same dimensions as the original DNS database i.e. 230×230×230 points. For the purpose of a-priori analysis, one has to decide if gradients of a filtered variable should be evaluated numerically based on the DNS grid (filter) size DNS or based on the size of the convolution filter conv , which corresponds in our case to conv = n DNS with n = 4, 8,12,16,20,24,28. Liu et al. [34] noted that for a meaningful a-priori test model expressions should be based solely on variables sampled on the coarse grid. The same logic is followed in this work and all model expressions are evaluated on the coarse grid. Consequently, one obtains truncation errors due to finite difference formulas in addition to the modelling errors in this approach. Nevertheless, it is worth mentioning that sampling the filtered data on the fine grid is useful for examining the potential of a model formulation. The difference between both approaches will be further illustrated in Section 6.

Closures for Sub-grid Scalar Flux (SGSF)
The most conventional model for (ρuc −ρũc) is the gradient hypothesis model given by: where μ t is the eddy viscosity and Sc t is the turbulent Schmidt number, which are taken to be C s = 0.18, Sc t = 1.0. Equation 4 will henceforth be referred to as the GH model. Richard et al. [6] proposed a model (RF) for T SGSF i in the following manner (C L = 0.12): The first termρC L u ∂c/∂x i on the right hand side of Eq. 5 is responsible for gradient transport (GT), whereas the second term −ρ 0 S L M i (c −c) accounts for counter-gradient transport (CGT). Variants of this model have been discussed in Chakraborty and Klein [14,15] and only one representative model will be considered here for the sake of brevity. Nevertheless it is worth mentioning that Lecocq et al. [8] suggested that replacingc on the right hand side of Eq. 3 byc conveniently includes the CGT modelling along with the FRDI term. Accordingly, a model for ∂ T SGSF i /∂x i can be formulated and will henceforth be denoted as T I M because it implicitly accounts for CGT if |∇c| in the LES closure of the unclosed term in the Favre-filtered reaction progress variable equation (refer to Eq. 3) is replaced by |∇c|.
It is worth noting that this model can only be formulated in terms of the divergence of SGSF subject to the assumption that (ρS d ) S ≈ ρ 0 S L . For the purpose of this work, the wrinkling factor is modelled according to Eq. 10, i.e. in agreement with one of the FSD models introduced below.
Clark et al. [35] proposed a model in the context of momentum transport closure for incompressible flows which takes the following form for the SGSF closure in this particular context: Equation 7 will henceforth be referred to as the CM (Clarks model) model in this paper. The CM model allows for both GT and CGT. This can be explained based on the following scaling analysis, which shows that Eq. 7 is indeed related to Eq. 5. Under strict flamelet assumption one obtains [36]:ũ i ≈c (u i ) P + (1−c) (u i ) R , which upon using Eq. 7 results in T CM i ∼ρ 2 (u i ) P −(u i ) R |∇c| 2 subject to the scaling |∇c| 2 ∼c (1−c) / 2 . This shows that the CM model scales as [36]. Hence the CGM model inherently accounts for the possibility of CGT (GT) when (u i ) P >(u i ) R (u i ) P <(u i ) R and ∂c/∂x i > 0. It was demonstrated in [14,15] that the CM model performs very well in a-priori analysis despite the fact that it was originally proposed for momentum transport in the context of incompressible flows.

Closures for Filtered Reaction Diffusion Imbalance (FRDI) Term
The analysis of Ma et al. [16] demonstrates that the Fureby model [11] seems to be a promising representative of an algebraic FSD model with the following modification (the unity is not present in the original model formulation) where it has been again assumed that (ρS d ) S ≈ ρ 0 S L and the modelled wrinkling factor [11]. The fractal dimension D f , evaluated using an empirical parameterisation, and the efficiency function are given by [11]: where δ L = D/S L is the Zel'dovich flame thickness. It is worth noting that the original model formulation by Fureby [11] reads: has neglected the resolved contribution of wrinkling and so, in the presence of no turbulence (or under fully resolved condition) the wrinkling factor becomes zero. This is the reason for suggesting the modified version given by Eq. 8 where the resolved component ρ 0 S L |∇c| is explicitly accounted for. The underestimation of the original model formulation has also been demonstrated in Ref. [13]. The quantity T FU will henceforth refer to Eq. 8 in this paper. An alternative FRDI model particularly developed to represent high pressure flames has recently been proposed by Keppeler et al. [37] and will be considered here, where O = 2.2 and I = max (δ L Ka −1/2 , 2δ L ) denote outer and inner cutoff scale respectively: Here Ka , the subgrid Karlovitz number, and the fractal dimension D f are given by [37]: The function F (c) can be approximated as 4 . For further details the reader is referred to the original papers [11,37], where these models have been proposed. A comparison between Eqs. 8 and 10 reveals that T K1 can also be considered as an algebraic FSD model. Noting that the term 4.5c(1 −c)F (c) −1 assumes a value of the order of unity [37], and the term k can be understood as the wrinkling factor. In addition to the differences between F and k , and the expressions for the fractal dimension D f , the resolved FSD |∇c| in Eq. 8 is replaced by |∇c| in Eq. 10. A comparison of Eq. 10 with Eq. 3 indicates that it is useful to analyse the following variants of Eq. 10: In other words, the models T K1, * , T K2, * (see Eqs. 12ii and 12iii) correspond to T K1 , T K2 but the term 4.5c (1−c) F (c) −1 is omitted. In this paper, the term T F U, * represents the Fureby model where |∇c| is replaced with |∇c| (Eq. 12iv) and T F U, * * equals T F U multiplied with 4.5c (1−c) F (c) −1 (Eq. 12v). These variants will be only considered once during the correlation analysis in order to separate the effects of the various terms from each other. Finally, it is important to note that the approximation (ρS d ) S ≈ρ 0 S L is not valid even for planar flames with Lewis numbers considerably smaller than unity. Hence, Chakraborty and Cant [22] suggested the following scaling in the context of Reynolds Averaged Navier-Stokes (RANS) simulations: It will be explicitly mentioned in the text and table captions if the terms T F U , T K1 , T K2 are used in conjunction with a modelling assumption other than (ρS d ) S ≈ρ 0 S L .

Postprocessing Methodology
LES models are supposed to accurately capture the local behavior of the corresponding term. Hence the Pearson correlation coefficient [38] is in this work considered an important indicator for the quality of a LES closure model and it should be as close to unity as possible. However, it is only a measure of linear dependence between two quantities and hence invariant under multiplication of the model with a constant. The second step in the analysis is therefore to compare the magnitude of the modelled expressions with the magnitude of the corresponding DNS counterpart. The approach adopted in this work is to first prepare conditional plots of each quantity. However, instead of showing all line plots for all cases, all filter widths and all models, a normalized mean deviation between closure and exact value of an unclosed term T is computed in the following manner: The divergence of SGSF (ρuc−ρũc) conditionally averaged onc for the unity Lewis number flame D is considered here to illustrate the definition given by Eq. 14 . The SGSF obtained from DNS data (T DNS ) for a filter width of /δ th = 0.8 is shown in Fig. 1 along with its absolute value which is used for the normalization shown in Eq. 14. The GH model (see Eq. 4) prediction T GH and the absolute value of the modelling error |T DNS − T GH | are shown as well. The normalized deviation is then given by the integral value of |T DNS − T GH | divided by the integral value of T DNS , which in this particular example is given by: ε = 0.0283/0.0222 = 1.275. This implies that the error is larger than 100 % which is due to the wrong sign of the GH model prediction.
The average error over all Lewis numbers and all filter widths Le, , and the average error over all filter widths for one specific Lewis number case for a given model T Model are defined in the following manner to condense the data further, Similar abbreviations with obvious meaning will be used for the correlation coefficients. As an example c Le, (T Model ) denotes the Pearson correlation coefficient between T Model and T DNS averaged over all Lewis numbers and filter widths.

Analysis of Sub-grid Scalar Flux Models
Before analysing the SGSF models (given by Eqs. 4-7) it is essential to provide some physical insight into the database used in this work. Veynante et al. [39] have pointed out that the competition between scalar transport due to turbulent velocity fluctuations and flame Fig. 1 Magnitude of the SGSF model expression ∇ · (ρuc−ρũc) × δ th /ρ 0 s L for the GH model T GH compared with its DNS counterpart T DNS for flame D and a filter width of /δ th = 0.8. The relative error is given by ε = 1 0 |T DNS − T GH |dc / 1 0 |T DNS |dc (see Eq. 14). Both integrands are shown as well in the figure normal acceleration determines the statistical behaviour of turbulent scalar flux: when the effects of flame normal acceleration are stronger than the effects of turbulent velocity fluctuation a counter gradient transport (CGT) is observed and vice versa. For our database this means that a strong CGT is observed when the Lewis number is considerably smaller than unity (Le = 0.34, 0.6, cases A, B). This is due to the effect that all Le < 1 flames are thermo-diffusively unstable because reactants diffuse into the reaction zone at a faster rate than the rate at which heat is conducted out, which leads to simultaneous occurrence of high temperature and reactants concentration. For flames with Le ≈ 1 turbulence effects become relatively strong and the amount of CGT is reduced for cases C, D, E with Le = 0.8, 1.0, 1.2. Details are presented in [14,15] and are not repeated here. In contrast to [14,15], the analysis of the SGSF models is done here in terms of the divergence of the scalar flux T SGSF i rather than the individual vector components. This is due to the fact that only the scalar expression ∂ T SGSF i /∂x i can be combined with T F RDI . Consequently, the evaluation of ∂ T SGSF i /∂x i includes truncation errors due to finite difference formulas in addition to modelling errors. In particular, the magnitude of gradients is under-predicted if the numerical differentiation is done on the LES grid size. This effect is illustrated in Fig. 2 where the divergence of the SGSF term extracted from explicitly filtered DNS data is numerically evaluated based on both the DNS grid (DN S) and the LES grid (LES) for filter widths /δ th = 0.4 and /δ th = 2.8 in cases A and D. It can be seen from Fig. 2 that evaluating the finite differences on the coarse grid results in a considerable underprediction of the magnitude of the gradients and hence of the divergence of SGSF.
As a result, it is difficult to differentiate between the (theoretical) performance of a model T SGSF i and the performance of the full expression ∂ T SGSF i /∂x i contained in the transport equation. It is also important to note that the expression T SGSF i may also depend on the gradients of flow variables. In that case the numerical errors will deteriorate the model performance further. In order to discuss these effects further, Fig. 3 shows the correlation coefficients for the different model expressions as well as the relative error as defined by Eqs. 14 and 15 for different Lewis number cases considered here. The results in Fig. 3a and b rely on the numerical differentiation performed on the DNS grid for discretising the divergence operator, whereas Fig. 3c and d are obtained by using the finite difference formulas based on the LES grid size.
As described in much more detail in [14,15] it can be seen from Fig. 3a and c that the gradient hypothesis model (GH) is in all cases negatively correlated with SGSF. Due to the increasing amount of CGT the correlation becomes more negative with decreasing Le. Richard's model (RF, Eq. 5) and the implicit model IM (Eq. 6) show a relatively high positive correlation for flames with small Le but the correlation strength decreases when turbulence effects become more dominant for high Le cases. Essentially it can be seen that the expression ρ 0 S L (|∇c| − |∇c|) changes the strong negative correlation of gradient type models (see Eq. 4 and first part of Eqs. 5 and 6) into a positive correlation, which shows that this term is of CGT type. Note, that assuming M i ≈ const the second term in Eq.5 can be written as ρ 0 S L (|∇c| − |∇c|) when the divergence of the SGSF term is calculated, as pointed out by Ma et al. [16]. Only Clarks model (Eq.7) shows a reasonably high positive correlation for all cases considered here. A comparison between Fig. 3a and 3c reveals that the correlation strength decreases when finite difference expressions are calculated on the coarse LES grid. The effect is most pronounced for the CM model. Considering next the deviation in Fig. 3b it is clear that the CM model has the smallest overall error Le, . In contrast to Fig. 3d the error Le, is of comparable magnitude for the models RF, CM and IM when the differentiation is done on the coarse grid. In particular, it can be seen that the IM model provides a slightly smaller error than the CM model. This behaviour can be explained in the following manner. Let us consider two models for approximating ρu i c − ρũ ic where model 1 properly predicts SGSF but model 2 over-predicts SGSF. Taking the divergence on a coarse grid gives rise to an under-prediction of ∇ i (ρu i c−ρũ ic ) magnitude as discussed before, and this can potentially result in a situation where model 2 performs better than model 1 on a particular grid, because the overprediction is compensated by numerical differentiation on a coarse grid.

Analysis of Generalised FSD Models
In this section, we focus on the modelling of the generalized FSD. In order to avoid introducing additional nomenclature, the models T F U , T K1 , T K2 have to be understood in this section as the models for gen only (rather than (ρS d ) S gen ). Table 1 shows the averaged correlation coefficients c Le, for T F U , T K1 , T K2 together with the variants T K1, * , T K2, * , T F U, * , T F U, * * . Furthermore, Table 1 shows the averaged individual case specific correlation coefficients c for the FRDI models T F U , T K1 , T K2 considered in detail in this work. It is clear from Table 1 that the correlation strength decreases consistently with decreasing Le. The individual case-specific correlation coefficients c for the model variants T K1, * , T K2, * , T F U, * , T F U, * * are not shown in Table 1 because this data does not provide any additional information in comparison to the averaged correlation coefficients c Le, , and also because these model variants (i.e. T K1, * , T K2, * , T F U, * , T F U, * * ) are not discussed further in this paper. It can be seen that T K1 , T K2 have slightly higher positive correlation coefficients than that of T F U . This effect is entirely due to the inclusion of the additional factor 4.5c (1−c) F (c) −1 as can be seen from the correlation coefficients of the model variants T K1, * , T K2, * (where 4.5c (1−c) F (c) −1 has been removed) and T F U, * * (where the term 4.5c (1−c) F (c) −1 has been included). Furthermore, Table 1 indicates that the models where gen ∝ |∇c| (e.g. T K1 , T K1, * , T F U, * ) show comparable correlation coefficients to the models where gen ∝ |∇c| (e.g. T K2 , T K2, * , T F U ).
A detailed investigation reveals that the amplitude spectrum ofc andc differs mainly at high frequencies. To illustrate this effect we consider a planar laminar back to back flame with Le = 1.0, because this configuration allows to perform a Fourier transformation in a natural manner due to periodicity. Figure 4a show the corresponding profiles of reaction progress variable c and the normalized density ρ. In Fig. 4b the distribution ofc andc  Fig. 4 Sketch of a planar laminar back to back flame with Le = 1.0: Profiles of a normalised density ρ/ρ 0 and reaction progress variable c. b filteredc and Favre filteredc reaction progress variable, assuming a filter size of /δ th = 2.8; c magnitudes of ∂c/∂x and ∂c/∂x normalized with δ th . The numerical differentiation is done with respect to either the DNS grid size or the LES filter width, as indicated in the subfigure; The x-coordinate is normalised with the thermal flame thickness δ th in subfigures a-c; d Semi-logarithmic plot of the one-sided amplitude spectra ofc andc as functions of the non-dimensionalised wave number are shown for /δ th = 2.8. The magnitude of the gradients of these two quantities (see Fig. 4c) is numerically evaluated using both DNS grid size as well as the LES filter width /δ th = 2.8. Finally Fig. 4d displays the one sided amplitude spectra ofc andc. Two effects can be seen from Fig. 4c and d, which are discussed in the following manner. The spectral content ofc andc differs mainly at high wavenumbers, i.e. the part which is in the sub-filter range. Ifc andc are sampled and differentiated on an increasingly coarser grid the Nyquist theorem shows that only the lowest frequencies can play a role and this is exactly the frequency range where both quantities are similar in terms of their frequency content. Secondly, it can be seen (refer to Fig. 4c) that the difference between |∇c| and |∇c| evaluated numerically on the LES grid is very small in comparison to the exact value of |∇c| obtained using the DNS grid. This explains why the correlation coefficients of the models proportional to |∇c| and the models proportional to |∇c| approach each other for large filter sizes, because they have to be compared to exact quantities evaluated on the DNS grid.
In order to assess the performances of the three different models further, it is important to understand the variation of the volume averaged value of the modelled FSD with filter size. Ideally this quantity should not vary with the filter size as the total flame surface area remains unchanged. Figure 5 shows the relative error of total model and real flame surface area as a function of filter width . It can be seen from Fig. 5 that Fureby's model slightly under-predicts the flame surface area but it has the correct behaviour in the limit of small filter size. In contrast, Keppeler's model tends to over-predict the flame surface area for large filter width and in addition it under-predicts when the flame is well resolved. The amount of under-prediction is even larger than the flame wrinkling observed for this filter width, which can be explained by the fact that the modelled wrinkling factor k assumes values smaller than unity for <δ L . This could be easily fixed by taking the maximum of unity and modelled k . The same behaviour can also occur for Fureby's model when the unity in Eq. 8 is omitted [16]. Despite the fact that T K1 and T K2 tend to overpredict the flame wrinkling in comparison to this DNS data for large filter widths, it is worth noting that this model was calibrated and shown to work satisfactorily in a-posteriori analysis for a variety of flames.

Analysis of Filtered Reaction Diffusion Imbalance Models
The next step in the analysis is to consider the modelling of (ρS d ) S because this quantity needs to be modelled in addition to gen in order to model the FRDI term (ρS d ) S gen . Table 2 shows the correlation coefficients using the standard assumption (ρS d ) S ≈ ρ 0 S L in comparison with the exact expression obtained from DNS as well as combined with the different models for gen (i.e. T F U , T K1 , T K2 ). The results for the approximation given by Eq. 13, i.e. the Lewis number correction, are not shown in Table 2 because the correlation coefficient is invariant under this multiplication. A comparison between Tables 2 and 1 reveals two effects. Firstly, neglecting the variation of (ρS d ) S by assuming (ρS d ) S = ρ 0 S L results in a considerable reduction of the correlation coefficient. Secondly, including the exact expression for (ρS d ) S results in an increased correlation between     Figure 6 shows the effects of (ρS d ) S estimation on the error incurred due to FSD based closure. As pointed out by Chakraborty and Cant [22] the approximation (ρS d ) S ≈ ρ 0 S L does not hold true for non-unity Lewis number flames and consequently there is an increasing trend of when Le deviates from unity. This can especially be seen for cases A and B. It is remarkable that using the correction given by Eq. 13 brings the errors down nearly to the level one gets using the exact numbers from DNS for (ρS d ) S .

Analysis of Combined Modelling of FRDI and SGSF
In order to explain the effect of combined SGSF and FRDI modelling on the evaluation of −∇ · (ρuc −ρũc) + ∇ · (ρD∇c) +ω c , it is useful to note that the correlation coef- The FRDI term and its variance are larger in magnitude than the SGSF term. If the random variables representing expressions for SGSF and FRDI are considered to be  (18) It turns out that the sum of the three covariances in Eq. 18 is on average of the order of 5 % of Cov T F RDI DNS , T F RDI LES and consequently one can expect that the correlation coefficients of the combined model expressions (SGSF+FRDI) are within a few percentages of those for FRDI listed in Table 2 where the small differences observed are to a large extent due to the interaction of the 3 co-variances shown in Eq. 18. In fact, Table 3 shows that this indeed is the case. It is interesting to observe that the Fureby model performs slightly better when combined with the RF or the CM model which can be expected from a-priori analysis of sub-grid scalar flux models. The K2 model is rather insensitive to combination with a flux model. However, the Keppeler model in its original formulation (K1) performs best when combined with the worst performing scalar flux model i.e. the GH model. This can be understood from the fact that using |∇c| instead of |∇c| acts as an implicit CGT model which in combination with the GH model indeed behaves similar to the RF model as can be seen from Fig. 3. Indeed the terms ρ 0 S L (|∇c| − |∇c|) and ρ 0 S L (|∇c| − |∇c|) seem to be particularly good in representing CGT. Therefore adding another flux model which is good in representing CGT ultimately results in a worse overall behaviour because CGT effects are taken twice into account. The highest correlation coefficients are achieved for the Keppeler  Table 4 which shows the errors based on the conditional quantities. It is worth noting again, that the Lewis number correction (Eq. 13) considerably improves the model behavior in terms of the magnitude of the model expression but it does not affect the correlation coefficients.

Conclusions
The performance of several models for subgrid-scale scalar flux and filtered reaction diffusion imbalance (FDRI) in the context of turbulent premixed flames has been assessed on an individual basis and in terms of their combined contribution. The analysis is based on apriori investigations of DNS data for a range of different Lewis numbers. The main focus of this analysis is to assess a new FRDI closure and in particular the effect of combining SGSF and FRDI models. The main findings of the paper can be summarized as follows: • The gradient hypothesis model (GH) is negatively correlated with the DNS SGSF term in all cases. The correlation becomes increasingly negative with decreasing Lewis number. By contrast, Richard's flux model (RF) shows a clear positive correlation in particular when CGT is dominant, i.e. for small Le cases. Clark's model (CM) shows overall the best performance in representing SGSF. However, it is important to note that the performance of the models is affected by numerical errors because the transport equation contains the divergence of the SGSF term, rather than its model. • Keppeler's model KP1 exhibits a slightly higher positive correlation than the Fureby's model. It has been shown that this is due to the appearance of the additional term 4.5c(1 −c)F (c) −1 in this model expression. However, Fureby's model, with the modification used in this work, is more successful in predicting the flame surface area for the cases considered in this work, in particular in response to the variation of filter size. • The magnitude of the FRDI term dominates over the SGSF term. However, an ideal combination of FRDI and SGSF models does not necessarily consist of the best individual models. However, the resulting closure depends on subtle details, such as the co-variances of the FRDI and SGSF terms. • The approximation (ρS d ) S ≈ ρ 0 S L is in the context of non-unity Lewis number flames not at all applicable even for statistically planar flames. The scaling (ρS d ) S ≈ ρ 0 S L /Le improves the magnitude of the surface-weighted filtered values of density weighted displacement speed considerably but the correlation remains low. The results show clearly that a more accurate modelling of (ρS d ) S is necessary in the context of LES.
The a-priori analysis performed in this work helped to explain some of the observations from a-posteriori analysis reported by Allaudin et al. [17]. For example, it is possible to explain that results obtained with the Keppeler model (KP1) will not improve when a more advanced SGSF closure than GH is used. Regarding the expectation that a FRDI model proportional to |∇c| should perform worse than a model proportional to |∇c|, it has been demonstrated that both versions yield similar results when the gradients are evaluated on a coarse LES grid. However, a recent a-posteriori analysis [16] showed that using |∇c| can result in undesirable flame thickening and that using |∇c| is sometimes required to counter this effect. To explain this phenomenon it is certainly important to include numerical discretization errors in the analysis, in particular because transport ofc is often achieved with the help of flux-limiter schemes that are prone to numerical diffusion. Understanding and explaining these effects in detail is beyond the scope of this analysis and the state of the art of a-priori analysis.