## 2 Numerical Methodology and DNS Database
An increase in CFL number by a factor 2.0 did not significantly affect the values of \(S_{L}\) and \(\delta_{th}\) (i.e., less than 0.1% difference). Time integration is performed using an explicit third-order low-storage Runge-Kutta scheme, and spatial derivatives for all internal grid points are evaluated using a 10 th order central difference scheme, but the order of accuracy gradually drops to a one-sided second-order scheme at the non-periodic boundaries. The SENGA code is well established in the scientific community, and its implementation has been verified several times in the past. Exemplarily, it is mentioned that for a Taylor-Green vortex, the maximum deviation in enstrophy with respect to reference data [26; 27] is 2.5%, while kinetic energy can be considered to be identical with those reference solutions. The boundary conditions in the mean flame propagation direction are taken to be partially nonreflecting, whereas \begin{table}
\begin{tabular}{c c c c c c} \hline \hline
**Case** & \(Re_{t}\) & \(u^{\prime}/S_{L}\) & \(l\delta_{th}\) & \(Da\) & \(Ka\) \\ \hline A & 11.67 & 1.0 & 4.58 & 4.58 & 0.47 \\ B & 87.5 & 7.5 & 4.58 & 0.61 & 9.6 \\ C & 175.0 & 15.0 & 4.58 & 0.31 & 27.16 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Initial parameters for the three turbulent premixed flames considered in this analysis. boundaries in transverse directions are taken to be periodic. The computational cost is of the order of \(10^{5}\) CPU hours (on Intel Xeon E5) for each of the cases considered here. Figures 1-3 show the instantaneous distribution of reaction progress variable \(c\) with superimposed isocontours corresponding to \(c=0.1,\ 0.5,\ 0.9\) for cases A-C in two different sections cutting the computational domain in \(x\)-\(y\) and \(y\)-\(z\) direction. In all figures, the flame propagates from right to left along the negative \(x\)-direction. In this context, the reaction progress variable \(c\) is defined based on the reactant mass fraction \(Y_{R}\) as follows: \(c=\left(Y_{R0}-Y_{R}\right)/\left(Y_{R0}-Y_{R\infty}\right)\) where the subscripts \(0\) and \(\infty\) refer to the values in the unburned reactants and fully burned products, respectively. It can be observed from Figure 1 that \(c\)-isosurfaces are mostly parallel to each other in case A, whereas the \(c=0.1\) isosurface is more distorted than the \(c=0.9\) in cases B and C. The Karlovitz number increases from case A to case C (see Table 1), which results in a larger length scale separation between \(\delta_{th}\) and \(\eta\). As a consequence, energetic turbulent eddies are more likely to perturb the preheat zone for high values of \(Ka\) and local flame thickening can be observed in case C. Furthermore, it can be also observed from Figures 1-3 that flame wrinkling increases from case A to case C. While turbulent structures tend to be relatively isotropic for cases B and C, case A clearly shows larger structures for the velocity fluctuations which are aligned with the mean flame propagation direction (i.e. \(u^{\prime}/S_{L}\)) compared to those which are normal to the \(x\)-direction (e.g., \(v^{\prime}/S_{L}\)). [PERSON] et al. [28] derived the following expression of \(\langle\rho u_{i}^{{}^{\prime\prime}}u_{j}^{{}^{\prime\prime}}\rangle\) based on a presumed bi-modal probability density function (PDF) of reaction progress variable \(c\). \[ \begin{array}{ll}\langle\rho u_{i}^{{}^{\prime\prime}}u_{j}^{{}^{\prime\prime }}\rangle=\langle\rho\rangle\{c\}&(1-\{c\})[\langle u_{i}\rangle_{P}-\langle u _{i}\rangle_{R}]\big{[}\langle u_{j}\rangle_{P}-\langle u_{j}\rangle_{R}]\\ &+\langle\rho\rangle(1-\{c\})\langle u_{i}^{\prime}u_{j}\rangle_{R}+\langle \rho\rangle\{c\}\langle u_{i}^{\prime}u_{j}\rangle_{P}+O(\gamma_{c})\end{array}  \tag{3}\]  
Here, \(\{Q\}\) is the Favre-averaged/filtered value of a quantity \(Q\), defined as \(\{Q\}=\langle\rho Q\rangle/\langle\rho\rangle\) and the corresponding Favre fluctuation is given by \(Q^{\prime\prime}=Q-\{Q\}\).