Comparison of Flame Propagation Statistics Based on Direct Numerical Simulation of Simple and Detailed Chemistry. Part 2: Influence of Choice of Reaction Progress Variable

Flame propagation statistics for turbulent, statistically planar premixed flames obtained from 3D Direct Numerical Simulations using both simple and detailed chemistry have been evaluated and compared to each other. To achieve this, a new database has been established encompassing five different conditions on the turbulent combustion regime diagram, using nearly identical numerical methods and the same initial and boundary conditions. The discussion includes interdependencies of displacement speed and its individual components as well as surface density function (i.e., magnitude of the reaction progress variable) with tangential strain rate and curvature. For the analysis of detailed chemistry Direct Numerical Simulation data, three different definitions of reaction progress variable, based on CH4,H2O and O2 mass fractions will be used. While the displacement speed statistics remain qualitatively and to a large extent quantitatively similar for simple chemistry and detailed chemistry, there are pronounced differences for its individual contributions which to a large extent depend on the definition of reaction progress variable as well as on the chosen isosurface level. It is concluded that, while detailed chemistry simulations provide more detailed information about the flame structure, the choice of the reaction progress variable definition and the choice of the resulting isosurface give rise to considerable uncertainty in the interpretation of displacement speed statistics, sometimes even showing opposing trends. Simple chemistry simulations are shown to provide (a) the global flame propagation statistics which are qualitatively similar to the corresponding results from detailed chemistry simulations, (b) remove the uncertainties with respect to the choice of reaction progress variable, and (c) are more straightforward to compare with theoretical analysis or model assumptions that are mostly based on simple chemistry assumptions.


Introduction
The use of fossil resources for the supply of primary energy leads to increasing concerns about emissions and their environmental, and health impact [1]. However, the challenges laid out by the Paris Agreement are daunting [2] and the time scales to phase out the existing infrastructure are substantial [1]. Therefore, combustion remains a topic of considerable importance and Direct Numerical Simulations (DNS) offer valuable insights into the physical and thermochemical phenomena, as reviewed for the early DNS work by Poinsot et al. [3], and more recently by Chen [4]. The ever-growing computing power is an enabler for two trends in recent combustion DNS research: Firstly (not discussed here but underlining the importance of combustion DNS), the availability of high-fidelity datasets has given a significant boost to the development of turbulence and combustion models using machine learning methods [5] and secondly detailed chemical kinetic mod-els for larger and larger hydrocarbon and oxygenated hydrocarbon molecules are being produced [6]. While practical fuels for transportation often contain thousands of distinct chemical compounds [7], there is also an increasing interest in simple fuels such as hydrogen [8,9] or syngas [10,11] that could play an important role during the transition of the fuel landscape. Their main components such as H 2 , CO and CH 4 are conceptually at the base of detailed chemical reaction mechanisms that describe much more complex hydrocarbon combustion [12]. Early work on two-dimensional fully resolved simulations of methane-air flames has been reported for example in [13][14][15]. Peters et al. [13] focused on the three different mechanisms arising from the chemical reaction, flame normal diffusion, and flame curvature, which dictate the statistical behavior of flame displacement speed. A similar analysis was performed by Echekki and Chen [14] to identify the contributions of curvature to the displacement speed by using the governing transport equation for the deficient reactant. Chen and Im [15] analysed the curvature stretch correlation indicating two distinct branches, one for positive and one for negative displacement speeds. The authors also indicated the large importance of curvature stretch which results in a non-linear flame speed correlation with stretch rate. Methane flames in three dimensions have been discussed by Bell et al. [16].
The complexity for parameterising chemical kinetics, even for simple fuels, ranges from a few species and reactions [17] up to hundreds of species and thousands of reactions [18]. Such mechanisms often become unmanageable for combustion DNS in three dimensions, and the increased complexity underlines the importance of accurate thermodynamic parameters for all species together with collision efficiencies of different third bodies involved in pressure-dependent reactions [6]. While detailed chemical mechanisms provide a wealth of information about the flame structure, they also have disadvantages in the context of analysis, and modelling turbulence-chemistry interaction. The choice of reaction progress variable definition together with the value of a suitable isolevel gives rise to uncertainty in the interpretation of results and it becomes more difficult to compare with theoretical predictions or models which are often based on simple chemistry assumptions. Klein et al. [19] discussed the modelling of turbulent scalar fluxes for different species in hydrogen air flames. The statistics of the turbulent scalar flux transport equation [20] and the flame surface density transport equation [21] revealed remarkable differences when performed for different major species. Finally, the effect of choosing different species on the determination of flame area has been analysed in [22].
In the following discussion, calculations using one irreversible reaction will be referred to as simple chemistry (SC) simulations, while calculations using a detailed chemical mechanism with detailed transport will be abbreviated as detailed chemistry (DC) simulations.
It becomes clear from the foregoing discussion that it is not straightforward to identify one, correct physical system representation for combustion DNS and it should be selfevident, that the methods chosen for the scientific analysis should be commensurate with the requirements of a specific scientific purpose, as argued in [23] because highperformance computing is itself responsible for non-negligible carbon emissions [24,25]. Another tradeoff that should be considered is that SC simulations allow for large parametric variations and at least in some respects, it might be possible to extract more physical insight from a larger matrix of simulations compared to a few very detailed DNS runs. The present work focuses on the fluid-dynamic aspects of flame-turbulence interactions and does not deal with emissions or ignition which would, without any doubt, require a more complex treatment of chemical kinetics. Flame propagation statistics in terms of displacement speed and its stretch rate dependence are of fundamental importance in the level-set and Flame Surface Density-based modelling methodologies [26]. The flame stretch, in turn, can be decomposed into tangential strain rate and curvature stretch contributions. The reader is referred to the textbook of Poinsot and Veynante [26] for a general introduction. The definition of displacement speed typically relies on a suitably defined reaction progress variable and the same holds true for mean flame curvature. Moreover, displacement speed statistics are typically taken on a specific isosurface within the flame front, which again requires the definition of a reaction progress variable together with the identification of the isolevel of reaction progress variable. It has been shown in the first part of this work [23] that displacement speed statistics and their interrelation with curvature and tangential strain rate are in very good qualitative and reasonably good quantitative agreement between simple and detailed chemistry DNS. The present work goes into more detail by analysing the individual components of displacement speed as well as the surface density function (i.e., magnitude of the gradient of reaction progress variable) and their interdependencies with curvature and strain. In particular, the influence of the definition of reaction progress variable will be discussed in detail. Different species participating in the chemical reaction are characterised by different individual Lewis numbers and the global Lewis number is known to influence flame propagation statistics. Such analysis has been pioneered by Rutland and Trouvé [27] for generic planar flame simulations showing a pronounced Lewis number dependent correlation between surface curvature and the local flame speed. It was also shown that the local flame structure is altered in curved regions. The effects of Lewis number on strain rate and curvature dependences of displacement speed have been analysed individually by Chakraborty and Cant in references [28] and [29], respectively. In these analyses, once more, statistically planar turbulent premixed flames have been considered and the same holds true for the analysis of the evolution of flame surface density [30] as well as surface density function [31] which are highly relevant quantities for modelling turbulent premixed combustion. While some effects of curvature might be annihilated in a mean sense for statistically planar flames, the mean curvature plays an important role for Lewis number effects on flame speed statistics, as demonstrated by Ozel Erol et al. [32] for spherical flames and by Rasool et al. [33] for Bunsen flames. Because of the Lewis number dependence, it might be expected that displacement speed statistics will be different for reaction progress variable definitions based on different species mass fractions. Hence, the main objectives of the present analysis are:

1.
To identify the influence of the definition of reaction progress variable on the statistics of flame propagation and its correlations with curvature and strain.

2.
To compare results from SC and DC simulations in order to provide information to balance their largely different computational costs as well as their value in terms of degree of detail compared to a larger parametric analysis or larger scale separation.
The same DNS database that was used in reference [23] has been used for this analysis. Although detailed discussions on this database and its numerical implementation were provided in [23], some of that information is repeated here for ensuring the self-contained nature of this paper in Section 2. This will be followed by a detailed discussion of the results and finally, conclusions will be drawn.

Direct Numerical Simulation Database
Two different compressible high order DNS codes SENGA [34] and SENGA2 [35] have been used for the present analysis which are described in detail in part 1 of this paper [23] together with all the necessary information related to this database. Both codes solve the conservation equations of mass, momentum, energy and species (see [34,35]) which can also be found in the textbook of Poinsot and Veynante [26]. For self-consistency, a short description summarizing the main features is provided here in tabular form for the numerical methods and physical models (see Table 1) and a schematic diagram of the computational configuration (see Table 2) is presented in Figure 1a, and the species Lewis numbers are presented in Table 3. Snapshots of the flame contours at the time when statistics are extracted are exemplarily shown for one case only in Figure 1b. For more details and all flame contours the reader is referred to the first part of the paper [23].  [40] are used at all non-periodic boundaries Initial data Homogeneous isotropic turbulent flow field [41] with superimposed laminar flame data Resolution ∆x, y, z = 3 × 10 −5 m which resolves the Kolmogorov scale and as well the thermal flame thickness with 13.7 grid points Duration Simulations were run for one chemical time scale which corresponds to at least 3 (2.14 ) eddy turn over times for cases A, B, C, E (case D)  All flames are initialized with a precomputed laminar flame profile where, in the case of SC, the laminar flame speed and the thermal flame thickness δ th were matched to the DC DNS case by adjusting the viscosity and pre-exponential factor, where δ th is defined as with T being the instantaneous, dimensional temperature and the subscript L refers to laminar flame quantities.
The initial values of normalised root-mean-square turbulent velocity fluctuation u /S L , turbulent length scale to flame thickness ratio l/δ th , Damköhler number Da = l T S L /u δ th and Karlovitz number Ka = (u /S L ) 3/2 (l/δ th ) −1/2 for cases A-E are presented in Table 4 where µ 0 is the unburned gas viscosity. All flames in this analysis belong to the thin reaction zones regime, as defined by Peters [42], and their position on the regime diagram is shown in Figure 1 of [23]. The case names will henceforth be referred to in this paper in such a manner that the second letter (S for SC and D for DC) distinguishes 10 simulations. A reaction progress variable can be defined as where Y is the mass fraction of the chosen species, and subscripts 0 and ∞ indicate the values in the unburned and fully burned gases, respectively. Any major product or reactant with a monotonic behaviour can be chosen for the definition of c and for the purpose of this analysis three major species CH 4 , H 2 O and O 2 have been selected. Based on the reaction progress variable c, an alternative flame thickness similar to the thermal flame thickness can be defined in the following manner: which can vary considerably as shown in Table 5. If not mentioned otherwise, results for each respective species are normalised with the corresponding flame thickness, denoted δ L for the ease of notation. A reaction progress variable can also be defined based on temperature as c T = (T − T 0 )/(T ad − T 0 ). Due to the possibility of superadiabatic temperatures (in particular for non-unity Lewis number flames) the use of c T can potentially be problematic. Nevertheless, in the case of SC DNS c T and c become identical for unity Lewis number, adiabatic, low Mach number conditions, and under these conditions we have δ th = δ L and c T = c. However, for detailed chemistry, this flame thickness can be notably smaller than the thermal flame thickness (see Table 5) and depends on the choice of reaction progress variable. For the present thermo-chemistry, δ L = 0.27 mm and δ th = 0.43 mm are obtained for the stoichiometric methane-air laminar premixed flame. An alternative and frequently used definition of flame thickness is the Zel'dovich flame thickness, given by: where D th,0 is the thermal conductivity in the unburned gas. As explained in detail in [23], under the assumptions detailed in Table 1, it is not possible to simultaneously match δ z and δ th in DC and SC simulations and the thermal flame thickness has been chosen as the reference value in this work. This implies that viscosity is larger in the SC simulations and hence turbulence decays faster and flame wrinkling is slightly reduced for SC when compared to DC at the time when statistics are taken. Following the suggestion of Peters [42] and previous DNS studies on displacement speed statistics [13][14][15]43,44], results in this paper will be presented for a progress variable isosurfaces in the reaction layer, close to the location of maximum reaction rate c ≈ 0.8 for SC simulations. As the use of c T has been abandoned due to the possibility of super-adiabatic temperatures (which results in c T > 1 and leads to negative reaction rates in the Arrhenius expression), there are ambiguities regarding the appropriate c values for alternative definitions of reaction progress variable. Figure 2a shows the mass fractions of selected important species versus axial position from a laminar stoichiometric methane-air flame while Figure 2b shows the corresponding reaction progress variable profiles for the mass fractions of CH 4 , H 2 O, and O 2 . It becomes obvious that c = 0.8 refers to different positions in the flame for different reaction progress variables with the CH 4 (O 2 ) isosurfaces closest to the reactant (product) side. Figure 3 shows the reaction progress variable based on CH 4 , H 2 O, O 2 mass fractions versus reaction progress variable based on temperature c T together with and h 0 f,i being the reaction rate, mixture specific heat at constant pressure in the unburned gas and enthalpy of formation of species i, respectively, and N is the total number of species). It becomes clear from Figure 3a that the value of c = 0.8 corresponds roughly to the location of maximum heat release for reaction progress variable based on H 2 O mass fraction while the corresponding locations for CH 4 (O 2 ) based reaction progress variables are slightly shifted to the left (right) but still close to the location of maximum heat release. For consistency and simplicity, results are always reported for a reaction progress variable value of c = 0.8 for the respective species unless mentioned otherwise. Finally, for completeness, Figure 3b shows the normalised reaction rates for CH 4 , H 2 O, and O 2 versus c T from laminar flame data.
Another important point to note from Figures 2 and 3 is that the mass fraction and reaction progress variable profiles for H 2 O are not strictly monotonic and in fact, they assume a local maximum around c T = 0.7 or c H 2 O ≈ 0.85. Strictly speaking, this invalidates the use of H 2 O for the definition of a reaction progress variable. However, as the effect is small and happens close to the burned gas side, it was decided to keep results for H 2 O for comparison of the definitions of reaction progress variable reported in previous studies [20][21][22].  Premixed combustion in general and displacement speed statistics are known to depend on differential molecular transport effects [26][27][28][29][30][31][32][45][46][47] characterized by the Lewis number Le = α/D with α and D being the molecular thermal diffusion and mass diffusion D of the deficient reactant respectively. In lean mixtures, the Markstein length depends on the Lewis number of the fuel and in rich mixtures on the Lewis number of the oxidizer [45][46][47]. Predictions based on this simple theory are clearly not valid at near stoichiometric conditions and Bechthold and Matalon [48] suggested the following expression: with Φ = φ (Φ = 1/φ) for fuel-rich (fuel-lean) mixtures, with φ and β being the equivalence ratio and Zel'dovich number respectively. For stoichiometric conditions, A Le = 1 and the Lewis number becomes the average of Le F and Le O . Using the values of Table 3 for CH 4 and O 2 gives a Lewis number of Le = 1.04. While expression (1) can be applied to a single fuel it is not applicable to fuel blends and it neglects many species participating in the chemical reaction. Furthermore, while the species Lewis numbers can be taken to be approximately constant the mixture changes from the unburned (i.e., c T ≈ 0) to the burned gas side (i.e., c T ≈ 1) which gives rise to a local effective Lewis number, that might play a role in the context of local flame propagation analysis. While there is still no consensus on an adequate formulation for an effective Lewis number [46], the following local effective Lewis number could be defined following the ideas in [45,46] and extending it to a local analysis including all species excluding nitrogen. The volumetric definition reads: where x i are renormalised mole fractions for all species but excluding N 2 , which does not participate in the reaction. Alternatively, a diffusion-based effective Lewis number can be given by While these definitions are not meant to resolve the ambiguity in defining effective Lewis numbers, they provide a good impression of the changing molecular transport effects within the flame. Figure 4 shows the effective Lewis number according to Equations (6) and (7) as a function of the distance in the direction of flame propagation for the laminar flame data. There are quantitative differences for the definitions given by both expressions, but they qualitatively indicate the same general trend: the Lewis number is larger than unity on the reactant side, drops considerably within the flame and relaxes towards a more positive value on the burned gas side. Unburned gas flame properties are often used in order to correlate flame properties such as turbulent flame speed. In this regard, it is noted that the stoichiometric methane flame has an effective Lewis number slightly larger than unity, while often it is simply referred to as a unity Lewis number flame.

Results
Displacement speed S d is an important quantity for the modelling of turbulent premixed flames using, e.g., level-set or flame surface density [30,42] approaches. The displacement speed is defined as the speed at which the flame surface moves normal to itself with respect to an initially coincident material surface. The displacement speed can be decomposed into a reactive (S r ) normal diffusion (S n ) and tangential diffusion component (S t ) following [13,14].
where ρ, . ω c and D c are the density, chemical reaction rate and diffusivity of reaction progress variable, → N denotes the flame normal vector and κ m the mean flame curvature.
The reaction rates . ω α (where the index α is omitted if not necessarily required) for fuel mass fraction Y α of species α and for reaction progress variable c based on species α (denoted (8)) differ by the factor −1/(Y α,0 − Y α,∞ ). Hence, while the mass fraction related reaction rate will assume positive values for major products and negative values for major reactants, . ω c is always a positive quantity in the context of reaction progress variable.

Mean Variation and Statistics of Displacement Speed and Its Components
Mean variations of normalised displacement speed S d /S L across the flame for all cases are shown in Figure 5 for SC and DC using three different definitions for reaction progress variable based on CH 4 , H 2 O and O 2 mass fractions. The variations between the different cases can be explained by stretch effects [23] and are consistent for all cases shown in Figure 5a  An advantage of the CH 4 mass fraction-based reaction progress variable definition is that the displacement speed can be evaluated throughout the flame without any problems, while the H 2 O mass fraction-based definition becomes singular due to the nonmonotonicity of H 2 O-based reaction progress variable, as discussed earlier. Displacement speed defined by O 2 mass fraction-based reaction progress variable becomes difficult to evaluate towards the burned gas side because |∇c| assumes very small values before c reaches a value of unity. It has to be admitted that there are quantitative differences between displacement speed statistics based on SC and DC simulations, however, they are of the same magnitude as the differences based on DC simulations using different definitions of reaction progress variable.
The behaviour observed from Figure 5 is consistent with the probability density functions (PDF) of normalised displacement speed S d /S L on the c = 0.8 progress-variable isosurfaces, which are shown in Figure 6 for cases A-E and the quantitative evaluation of mean E and standard deviation σ of S d /S L for c = 0.8 are shown in Table 6. Table 6 indicates that the values of σ(S d /S L ) can be split into three groups, cases A, B, cases C, D and case E, and this behaviour is consistent for all cases and in particular when comparing SC and DC.  Similar trends can be observed for the PDF of normalised combined reaction and normal diffusion components of displacement speed (S r + S n )/S L shown for the c = 0.8 isosurface in Figure 7. It is noted that the DC PDFs for reaction progress variable based on H 2 O and O 2 mass fractions tend to be wider than those for SC and DC based on CH 4 . This behaviour will be analysed in more detail by looking next at the individual contributions S r , S n , S t , and the corresponding PDFs for c = 0.8 are shown in Figures 7-9, while the corresponding mean values and standard deviations for S r and S n are shown in Tables 7 and 8. The mean value of S t is close to zero for a statistically planar flame and therefore these numbers are not shown in the form of a table. Figure 8 and Table 7 indicate that the mean values of S n /S L are close to each other for all cases but there are differences based on the different definitions of reaction progress variable. This indicates that the differences in the PDFs of normal diffusion components are driven by the significant differences in the laminar profiles of reaction progress variable across the flame brush when using different species mass fractions for its definition as shown in Figures 2b and 3a. It is also noted that the width of the PDFs of S n /S L tends to be larger for the DC cases compared to the SC simulations. A possible explanation is that ρD c can be taken as a constant for SC simulations while ρD c is temperature-dependent in the DC cases, and hence is likely to introduce additional variations for the DC setup due to the decoupling of reaction progress variable and temperature (see Figure 3). Furthermore, as noted earlier, the SC simulations have higher viscosity and hence slightly smaller turbulence intensity at the time when statistics are taken.
The PDFs of the normalised reaction component of displacement speed S r /S L for the c = 0.8 isosurface are shown in Figure 9 and the first and second moments for the c = 0.8 reaction progress variable value are given in ω c is a unique function of c. In contrast, for the DC case, the reaction rate of an individual species depends not only on temperature but also on the reaction progress of several other species. Furthermore, the standard enthalpy of formation of a pure element is zero which shows that . ω c (and equally . ω α ) and . ω T are not always directly related. Hence, a change of reaction progress variable in the DC case is less directly interlinked with a change in temperature. In all cases (including DC), the turbulent motion of the fluid causes variations of the surface density function |∇c| (which scales with the inverse flame thickness) in such a manner that the variance of S r /S L increases from case A to case E.    Despite the remarkable differences of the PDFs of S r and S n for the different definitions of reaction progress variable, the reaction-diffusion balance in the expression (S r + S n ) remains qualitatively much more similar compared to the individual contributions as discussed before, because a more positive value of S r is balanced by a more negative value of S n and vice versa. Finally, the PDFs of normalised tangential diffusion component of displacement speed S t /S L for the c = 0.8 isosurface are shown in Figure 10. Due to the statistically planar nature of the flames considered here, their mean value is close to zero in all cases and the standard deviation has the same order of magnitude for all definitions of reaction progress variable but tends to increase with increasing turbulence intensity.
The statistics of S d , (S r + S n ), S t remain qualitatively, and to a large extent quantitatively similar for SC and DC simulations and also for different definitions of reaction progress variable. However, there are marked qualitative differences for individual contributions of displacement speed S r and S n . To explain this behaviour, it is important to understand their interrelation with curvature and strain rate as well as the influence of the definition of reaction progress variable, which will be discussed in the next subsection.

Curvature and Tangential Strain Rate Effects on Local Displacement Speed and Its Components
Theoretical studies suggest that flame stretch is the controlling parameter of the flame structure in the limit of weak turbulence and weak flame wrinkling (for an overview see [26,47]). Flame stretch rate is, in turn, a function of flame curvature κ m and tangential strain rate a T = δ ij − N i N j ∂u i /∂x j . The PDFs of a T and κ m have been discussed in part 1 of this paper and do not provide additional insight here. Their joint PDFs are shown in Figure 11 and the correlation coefficients between a T and κ m are reported in Table 9 for c = 0.8.
It can be seen from Figure 11 that there is a negative correlation between a T and κ m in all cases which is the strongest for cases A and B and weakest for cases D and E. The correlation strength is weakly dependent on the definition of reaction progress variable. These results are in very good agreement with an earlier analysis reported in [49]. As a T and κ m are not independent parameters, the dependence of S d and its components with a T can be deduced from the interrelation of S d and its components with κ m and vice versa. Therefore, the following discussion is limited to correlations of displacement speed and its components with mean curvature.
The joint PDFs of S d and (S r + S n ) with κ m are shown in Figures 12 and 13, respectively for the c = 0.8 isosurface and the corresponding correlation coefficients are reported in Table 10 for c = 0.8. It can be seen from Figure 12 that displacement speed is negatively correlated with mean curvature. While the correlation coefficients are close to −1.0 for cases A and B, the correlation becomes non-linear for increasing turbulence intensity, as can be discerned from the considerably smaller values of corr(S d , κ m ) in comparison to −1.0. Table 10 also shows that the strength of linearity or non-linearity depends to some extent on the physio-chemical model and also on the definition of reaction progress variable in the case of DC. Tangential diffusion component of displacement speed S t is deterministically negatively correlated with curvature (see Equation (8)) with corr(S t , κ m ) < −0.99 in all cases, and for all definitions of reaction progress variable, as illustrated in Figure 14. This shows that the non-linear curvature dependence of S d originates from the non-linear curvature dependence of (S r + S n ). For c = 0.8, the correlation strength between (S r + S n ) is weak in all cases. Figure 13 and Table 10 show that there is a weak positive correlation in the case of both SC and DC simulations based on the reaction progress variable defined in terms of H 2 O while the correlation turns out to be mildly negative for DC and reaction progress variable based on CH 4 and O 2 mass fractions. Again, for SC, the behaviour is in good agreement with the DC results and detailed explanations are also provided for the change of correlation strength with the value of reaction progress variable [49], while the present analysis focuses on the differences between SC and DC and also on the differences resulting from different definitions of reaction progress variable.      So far, results have revealed a reasonable qualitative agreement of displacement speed statistics with some quantitative variations when comparing DC with SC, which are of the same order of magnitude, as compared to different definitions of reaction progress variable in the case of DC. As illustrated in Figures 15 and 16, this behaviour changes drastically when looking at the individual contributions of displacement speed S r and S n with mean curvature. For example, there is a weak positive correlation corr(S r , κ m ) in the case of SC and DC based on H 2 O mass fraction, whereas the correlation is clearly negative in the case of DC with reaction progress variable based on CH 4 and O 2 mass fractions for the reaction progress variable isosurface c = 0.8. In other words, there is a qualitative change in corr(S r , κ m ) with the variation of the definition of reaction progress variable, which warrants further explanation. Similar qualitative differences can be seen from the S n − κ m joint PDFs shown in Figure 16. The values of corr(S r , κ m ) and corr(S n , κ m ) are reported in Table 11.   ω c is independent of curvature and thus, according to Equation (8) the correlation between S r and κ m is governed by the correlation between |∇c| and κ m . However, this statement is rendered invalid in the case of DC, where there might be a significant correlation between . ω c and κ m , which also depends on the choice of reaction progress variable. Figures 17 and 18 show the joint PDFs of normalised reaction rate − with mean curvature κ m and the joint PDFs for the surface density function |∇c| with κ m on the c = 0.8 isosurface, while the correlation coefficients are reported in Table 12. For the ease of notation, the dimensionless reaction rate − following. Indeed Figure 17 shows that corr    For DC, Figure 17 clearly indicates a negative correlation between reaction rate and curvature for all definitions of reaction progress variable on the c = 0.8 isosurface. In the case of SC, this would correspond to a scenario with Lewis number larger than unity [31,32] while the opposite would be expected for Le < 1. Indeed, it was discussed in the introduction that a methane flame is characterised by a Lewis number slightly larger than 1.0 according to the parameterisations given by Equations (5)-(7) [45][46][47] which is qualitatively consistent with the trends shown in Figure 17. However, the species CH 4 , H 2 O and O 2 have considerably different Lewis numbers (Le = 0.97, 0.83, 1.11 respectively, according to Table 3) and therefore their behaviours warrant a more detailed discussion.
In Le < 1 SC simulations [31,32], the focusing of reactants takes place at a faster rate than the rate of defocusing of heat at the positively curved zones, which leads to the simultaneous presence of high temperature and reactant concentrations and conversely for negatively curved regions. Hence, high (low) temperature values are associated with positive (negative) curvatures for Le < 1 and the opposite is valid for Le > 1 [27][28][29][30][31][32]. For the SC simulations, the reaction rate depends on temperature according to where B * is the normalised pre-exponential factor, E act is the activation energy, such that a negative (positive) correlation between T and κ m induces a negative (positive) correlation between . ω and κ m and the same holds true for non-dimensional temperature which will be denoted T + in the following discussion and in fact is identical to c T . For a reaction mechanism involving M steps and N species given by: the situation is much more complex and one obtains where the specific forward ( f ) and backward (b) reaction rate coefficients k f /b,β are given by an Arrhenius expression similar to Equation (9). Here, the ν α,β , ν α,β are called stoichiometric coefficients, X α denotes the N different species, M α are the molar masses, c α the molar densities. Equation (11) shows that the reaction rate for species α in the context of DC depends on many more factors compared to the SC 1-step case shown in Equation (9) such that a negative (positive) correlation between T + and κ m does not necessarily induce a correlation of same sign between ω T because the enthalpy of formation of species α can be very small or even zero. Finally, all species share the identical temperature field, such that explanations holding true for SC cannot directly be applied to DC without caution. In fact, the joint PDFs between temperature and curvature for the c = 0.8 isosurface are shown in Figure 19 and the correlation coefficients are reported in Table 13 for the c = 0.8 isosurface. Indeed, by comparing Figures 17 and 19, it becomes obvious that corr(T + , κ m ) > 0 while corr for reaction progress variable based on H 2 O mass fraction, indicating that the value of temperature and rate of formation of species α can be partially decoupled.  Water vapour is characterised by a higher diffusivity and Lewis number Le = 0.83 (see Table 3). For SC the thermo-diffusive mechanisms explained before lead to higher reaction rates in positively curved regions and higher flame wrinkling for Le < 1 compared to Le = 1 flames. As a result of this, positively (negatively) curved regions propagate faster (slower) and finger-like structures start to develop. A transfer of this mechanism to isosurfaces of individual Le < 1 species in DC simulations is obviously not possible because all species transport equations are coupled. By contrast, as discussed next, it appears that the flame wrinkling of the H 2 O mass fraction-based reaction progress variable isosurface is smaller compared to those based on CH 4 and O 2 isosurfaces. This is a result of preferential diffusion and its influence on intermediate reaction steps and can on one hand be attributed to higher diffusivity of H 2 O which smoothens high curvature magnitudes. On the other hand, this effect is further attenuated by higher (reduced) reactivity in the positively (negatively) curved regions. The turbulent flame areas A T normalised by the cross-section of the computational domain A 0 for CH 4 , H 2 O and O 2 are exemplarily given by A T /A 0 = 1.553, 1.487 and 1.557 respectively for case A. The values for temperature isosurfaces (representing the laminar flame temperature for c = 0.8) for the three definitions of reaction progress variable are given by A T /A 0 = 1.541, 1.516, 1.476.
The reduced wrinkling of the H 2 O based reaction progress variable isosurface compared to O 2 based reaction progress variable for case AD can be confirmed from Figure 20a, which shows a simultaneous plot of isosurfaces of reaction progress variable c = 0.8 for H 2 O (red) and O 2 (grey) with the O 2 based isosurface overlapping the H 2 O-based isosurface. Similar observations in terms of flame wrinkling have been reported for detailed chemistry hydrogen-air simulations [22]. Figure 20 also shows isosurfaces of c = 0.8 coloured with non-dimensional temperature for case AD for reaction progress variable based on CH 4 , H 2 O and O 2 mass fractions. The opposite correlation of temperature and curvature (observed in Figure 19) between CH 4 ,O 2 and H 2 O can be clearly seen. While positively curved regions are characterised by lower temperature in the case of CH 4 and O 2 based reaction progress variable isosurfaces, they show higher temperature for H 2 O. In order to understand this apparent contradiction, it is also important to recall that the different isosurfaces are characterised by different positions in the flame (see Figures 2 and 3) and to note that the colour bars have different temperature ranges. As the H 2 O isosurface is less wrinkled than the corresponding temperature isosurface (because of higher mass diffusivity than thermal diffusivity), positively curved regions reach out less into to unburned gas side and hence they are characterised by higher temperatures. Similarly, negatively curved regions reach out less into the burned gas side and are consequently characterised by lower temperatures. The effect can be clearly seen in Figure 21a and is consistent with the positive correlation between T and κ m for the H 2 O-based reaction progress variable isosurface.  In order to finally understand the joint PDFs of S r and κ m , besides reaction rate, it is important to look at the joint PDFs of |∇c| and κ m . An initial observation from Figure 18 is that the joint PDFs for the O 2 -based reaction progress variable behave to some extent similar to the SC case (i.e., they show a small correlation magnitude with a positive and a negative correlating branch for negative and positive curvatures, respectively) while the joint PDFs for reaction progress variable based on CH 4 and H 2 O mass fractions are negatively correlated and the positively correlating branch cannot be observed. It has been explained in part 1 of this paper that (for constant dilatation rate) a T and a n (where a n is the normal strain rate) are negatively correlated. Further, a compressive (negative) normal strain rate causes local thinning which acts to increase the scalar gradient |∇c|. The negative correlation between a T and κ m leads to the principally negative correlation between |∇c| and κ m . However, |∇c| can also assume small values at locations of high negative curvature because of secondary thickening effects induced by high values of dilation rate, which locally overcomes a T to induce extensive normal strain rate (i.e., a N = (∂u i /∂x i − a T ) > 0) [44,50,51] leading to a positively correlating branch for SC and DC based on O 2 -based reaction progress variable definition.
Based on the correlation behaviour of numerator and denominator in the definition of S r and S n , it is possible to infer the correlations of S r and S n with curvature. However, in contrast to low Mach number SC simulations where ρ can be assumed to be constant on a given isosurface, the same assumption does not hold any longer in the case of DC simulations, such that the denominator ρ|∇c| potentially shows a more complex behaviour than |∇c| on its own (it is recalled that ρ deterministically depends on T such that different signs of corr(ρ, κ m ) can be observed for different reaction progress variable definitions, see Figure 19). The Pearson correlation coefficient between two random variables x, y is given by the covariance of both variables divided by the product of their standard deviations: corr(x, y) = cov(x, y)/σ(x)σ(y). For the sign of the correlation, it is sufficient to look at the covariance between the random variables. Following Bohrnstedt and Goldberger [52], the covariance of the product of two random variables x, y with a third random variable z is given by the expression cov(xy, z) = E(x)cov(y, z) + E(y)cov(x, z) + TOM, where E denotes the expected value of a variable and TOM is the third-order moment of the three variables. It can be shown that the third-order moment vanishes under multivariate normality. In the present case, it is not negligible. The above identity shows that In other words, the sign of the correlation S r − κ m depends on the expected values of . ω c and (ρ|∇c|) −1 and the covariances of these variables with κ m . Tables 14 and 15 show the individual contributions of cov(S r , κ m ) and cov(S n , κ m ), respectively, according to Equation (12), exemplarily for reaction progress variable based on CH 4 and H 2 O. These species have been selected because they show similar correlations trends for . ω c − κ m and |∇c| − κ m but corr(S r , κ m ) shows the opposite sign for CH 4 and H 2 O, which requires additional explanation and a similar statement holds true for corr(S n , κ m ). Table 14. Statistical quantities that influence the value of cov(S r , κ m ) (see Equation (12)) for the c = 0.8 isosurfaces for all cases and reaction progress variable based on CH 4  Secondly, the positive correlation T − κ m induces a negative correlation ρ − κ m due to the ideal gas law. This, in turn, reinforces the negative correlation |∇c| − κ m for H 2 O-based reaction progress variable, while it results in a change of sign and weakening of correlation strength for CH 4 -based reaction progress variable because of opposing correlations of ρ and |∇c| with κ m . A similar explanation provides the explanations for the correlation behaviour between S n and κ m (see Table 15).
The foregoing discussion has shown that the statistics of S r , S n and |∇c| with curvature are sensitive to the definition of reaction progress variable and detailed explanations for this behaviour have been provided. Besides that, the statistics are also sensitive to the choice of isosurface level, which is also evident from the correlation coefficients reported in Tables 11-13. Important qualitative differences have also been found for the joint PDFs of normalised reaction rate, surface density function and non-dimensional temperature with mean curvature. These differences have been explained based on preferential diffusion and its influence on intermediate reaction steps which give rise to different degrees of wrinkling of different isosurfaces for different definitions of reaction progress variable.
For further illustration, the joint PDFs of normalised reaction rate, |∇c|, S r , and S n for different definitions of reaction progress variable are exemplarily shown in Figures A1-A4 in Appendix A, but in contrast to earlier figures, the statistics are now taken on the same isosurface corresponding to the O 2 based reaction progress variable value of c = 0.8. It can be seen from Figures A1-A4 in Appendix A that some of the qualitative differences observed earlier diminish if statistics are taken on the same O 2 based reaction progress variable isosurface, and similar qualitative trends have been observed for other reaction progress variable definitions. For example, corr(S r , κ m ) becomes now consistently negative and, in addition, the negative correlation corr(|∇c|, κ m ) vanishes.

Conclusions
Flame propagation statistics from three-dimensional simple chemistry and transport DNS have been compared to results obtained from detailed chemistry and transport simulations using a new database consisting of five simple chemistry (SC) and detailed chemistry and transport (DC) simulations of statistically planar turbulent premixed, stoichiometric methane-air flames for a range of different Damköhler, Karlovitz and turbulent Reynolds number values. Both methodologies use nearly identical numerical methods and the identical initial and boundary conditions. The discussion includes interdependencies of displacement speed and its individual components as well as surface density function and reaction rate with tangential strain rate and curvature using three different reaction progress variable definitions based on CH 4 , H 2 O and O 2 mass fractions. The findings can be summarised as follows:

•
The statistics of S d , (S r + S n ), a T , κ m and their interrelation with each other remain qualitatively and to a large extent quantitatively similar for SC and DC and as well for different definitions of reaction progress variable. • However, there are marked qualitative differences for individual contributions of displacement speed S r and S n and their interrelation with the curvature such that depending on the definition of reaction progress variable different signs of the correlation coefficient have been observed. • Similar qualitative differences have also been found for the joint PDFs of normalised reaction rate . ω + , surface density function |∇c| and non-dimensional temperature T + with mean curvature κ m . • These differences have been explained based on preferential diffusion and its influence on intermediate reaction steps which give rise to different extents of wrinkling of reaction progress variable isosurfaces for its different definitions. This leads to different correlations between temperature and mean curvature, which also induces different dependencies between density and mean curvature.

•
In addition, the different variations of reaction progress variable c across the flame for different definitions of reaction progress variable contribute to differences in |∇c| variation. This may lead to small values of |∇c| towards the burned gas side for some definitions of c, which results in a poorly conditioned definition of displacement speed and its components towards the product side.

•
The reaction rate of some individual species α can be decoupled from temperature or ω T while their interrelation is much more straightforward for a 1-step simplified chemical reaction. While for SC, reaction rate and density can be considered independent of mean curvature for unity Lewis number, adiabatic low Mach number simulations, this simplification does not hold true for detailed chemistry and transport.
Without any doubt, DC simulations provide more detailed information about the flame structure than SC simulations. However, the choice of reaction progress variable together with a suitable isosurface value gives rise to considerable uncertainty in the interpretation of displacement speed statistics, sometimes even showing opposing trends. Single-step chemistry simulations are shown to capture the correct qualitative trends of flame propagation statistics and do not have uncertainties with respect to the choice of reaction progress variable. Therefore, SC simulation results are more straightforward to compare with theoretical analysis or model assumptions, which are usually based on SC assumptions. The present results illustrate the advantages and disadvantages of both SC and DC simulations and aim to help combustion researchers to balance computational cost and associated carbon dioxide emissions with accuracy requirements. Moreover, the advantages of having lots of detailed information in DC needs to be contrasted with the possibility of having a large parametric study or a larger scale separation in SC.
Finally, this work is limited to the flames with an effective Lewis number close to unity. Extending the comparison to different Lewis numbers, together with the identification of suitable definitions for global Lewis number of more complex or dual fuel mixtures is part of the future work.

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