1. Introduction
The unresolved variance of a scalar 
$\phi$ is defined as
\begin{equation} \sigma ^2(\phi;\varDelta)=\widetilde{\phi ^2}-\tilde{\phi}\tilde{\phi}, \end{equation}
where 
$\tilde {\phi }=\overline {\rho \phi }/ \bar {\rho }$ denotes standard Favre filtering using a filter 
$G$ with a characteristic length scale 
$\varDelta >0$. In numerical simulations the scalar variance returns information on the level of resolution of the scalar field by the mesh, and provides a direct measure of the energy contained in the unresolved scalar fluctuations. The mechanisms of production/decay of the variance both for reactive and non-reactive scalars due to turbulent micro-mixing have been at the core of model development since the early days (Moss & Bray Reference Moss and Bray1977; Dopazo Reference Dopazo1979). The companion of the scalar variance is the scalar dissipation rate (Borghi Reference Borghi1988; Bilger, Saetran & Krishnamoorthy Reference Bilger, Saetran and Krishnamoorthy1991; Bray Reference Bray1996; Peters Reference Peters2000; Veynante & Vervisch Reference Veynante and Vervisch2002; Pantano, Sarkar & Williams Reference Pantano, Sarkar and Williams2003) that calibrates the decay rate of the scalar fluctuations much like turbulent dissipation is the companion of the turbulent kinetic energy in the popular two-equation approach to turbulence modelling (Jones & Launder Reference Jones and Launder1972). Therefore, the scalar variance and the scalar dissipation rate are among the fundamental ingredients of turbulent reactive flow modelling providing valuable information that can be used to infer the shape of the subgrid-scale (SGS) probability density function (p.d.f.) of the scalar (progress variable/mixture fraction). The concept of using a presumed p.d.f. for the SGS scalar p.d.f. was then introduced to average the response versus equivalence ratio (or mixture fraction) of reactive–diffusive layers under thermochemical equilibrium (Lockwood & Naguib Reference Lockwood and Naguib1975). This opened up many perspectives for combustion modelling using structural flame approaches in which the conditional response of the flame structure versus a given set of control parameters (mixture fraction, progress of reaction, scalar dissipation rate etc.) is averaged from the presumed p.d.f. Considering the asymptotic bimodal limit for the progress-variable p.d.f., the Bray–Moss–Libby theory was developed in which the scalar variance was key for investigating turbulent flame dynamics (Bray & Moss Reference Bray and Moss1977; Bray, Libby & Moss Reference Bray, Libby and Moss1985). These concepts were later extended to large-eddy simulation (LES) for predicting the distribution of chemical species including finite rate chemistry effects (Jimenez et al. Reference Jimenez, Linan, Rogers and Higuera1997; Pierce & Moin Reference Pierce and Moin2004; Domingo, Vervisch & Veynante Reference Domingo, Vervisch and Veynante2008; Lecocq et al. Reference Lecocq, Richard, Colin and Vervisch2011; Nambully et al. Reference Nambully, Domingo, Moureau and Vervisch2014; Mesquita, Mastorakos & Zedda Reference Mesquita, Mastorakos and Zedda2023). Therefore, accurate estimates of the scalar variance are required at various levels in the simulation of turbulent combustion, but also more generally in reactive flow modelling (Fox Reference Fox2003).
Classic models for the scalar variance in LES can roughly be divided into two distinct categories: transport-equation-based models (Jimenez et al. Reference Jimenez, Ducros, Cuenot and Bedat2001; Pera et al. Reference Pera, Reveillon, Vervisch and Domingo2006; Kolla et al. Reference Kolla, Rogerson, Chakraborty and Swaminathan2009; Keil, Klein & Chakraborty Reference Keil, Klein and Chakraborty2021) and algebraic models (Cook & Riley Reference Cook and Riley1994; Girimaji & Zhou Reference Girimaji and Zhou1995; Pierce & Moin Reference Pierce and Moin1998, Reference Pierce and Moin2004; Balarac, Pitsch & Raman Reference Balarac, Pitsch and Raman2008). The consensus in the literature is that models based on a transport equation explicitly include the effects of reaction and transport. Such models however require a secondary modelling layer in order to close a number of unresolved terms that appear in the transport equation for the scalar variance. These include models for the scalar fluxes, unresolved source terms and a model for the scalar dissipation rate (Bilger et al. Reference Bilger, Saetran and Krishnamoorthy1991; Kolla et al. Reference Kolla, Rogerson, Chakraborty and Swaminathan2009; Keil et al. Reference Keil, Klein and Chakraborty2021). Modelling all of these terms introduces further physical assumptions, additional levels of complexity, but also uncertainty with regards to the exact effects of each submodel's parameters on the evolution of the variance transport equation. In fact, quantifying the effects of numerical dissipation and model uncertainty on the evolution of the transport equation is a daunting task. Algebraic models on the other hand are computationally less expensive and straightforward to implement in LES. Notable examples include the scale-similarity model by Cook and Riley (Cook & Riley Reference Cook and Riley1994) first proposed more than two decades ago, but also gradient-based models (Girimaji & Zhou Reference Girimaji and Zhou1995; Pierce & Moin Reference Pierce and Moin1998). Even though algebraic models do not explicitly account for reaction/transport, such effects may still in principle be modelled implicitly. For instance, the classic dynamic model by Pierce & Moin (Reference Pierce and Moin1998) originally developed for a passive scalar (by employing equilibrium assumptions) may still in principle incorporate the effects of reaction since reaction itself is a source of scalar gradients in the flow. Therefore, it is not surprising that gradient-based models have been used to successfully model the progress-variable variance in LES of turbulent premixed flames with overall good results (Vreman et al. Reference Vreman, Albrecht, van Oijen, de Goey and Bastiaans2008; Vreman, Bastiaans & Geurts Reference Vreman, Bastiaans and Geurts2009a; Vreman et al. Reference Vreman, van Oijen, de Goey and Bastiaans2009b; Mukhopadhyay, van Oijen & de Goey Reference Mukhopadhyay, van Oijen and de Goey2015). Nevertheless, even the development of algebraic models requires the introduction of a number of limiting assumptions/hypotheses that may not always be valid or well justified. Consequently, both of the above modelling approaches do not lead to a general and consistent modelling framework whereby models of arbitrary accuracy can readily be developed.
An alternative modelling framework with the potential to address the above issues is based on the concept of reconstruction/deconvolution. Reconstruction methods aim to recover quantitatively accurate estimates of the unfiltered fields 
$\phi ^*$ from their filtered counterparts 
$\tilde {\phi }$ on the LES mesh. These reconstructions may then be used to estimate unresolved variables (scalar fluxes, Reynolds stresses, etc.) by explicit filtering. Probably the first successful implementation of reconstruction-modelling methods in LES dates to the works of Stolz and Adams who successfully simulated decaying homogeneous turbulence (Stolz & Adams Reference Stolz and Adams1999), boundary-layer flows (Stolz & Adams Reference Stolz and Adams2001) and flows with shocks (Adams & Stolz Reference Adams and Stolz2002). This LES modelling approach called approximate deconvolution modelling (ADM) consisted of a filtering step, followed by reconstruction and regularisation steps. In later works by Mathew et al. (Reference Mathew, Lechner, Foysi, Sesterhenn and Friedrich2003) and Mathew, Foysi & Friedrich (Reference Mathew, Foysi and Friedrich2006), ADM was shown to be equivalent to using a single filtering step using a composite filter. Since then, many non-reacting LES studies have also employed reconstruction-modelling approaches (Schlatter, Stolz & Kleiser Reference Schlatter, Stolz and Kleiser2004; Loginov, Adams & Zheltovodov Reference Loginov, Adams and Zheltovodov2006; Boguslawski et al. Reference Boguslawski, Wawrzak, Paluszewska and Geurts2021). In reacting flows, reconstruction-based models have been successfully applied probably first by Mathew (Reference Mathew2002), and in a number of later studies (Domingo & Vervisch Reference Domingo and Vervisch2015; Wang & Ihme Reference Wang and Ihme2019; Domingo et al. Reference Domingo, Nikolaou, Seltz and Vervisch2020; Datta, Mathew & Hemchandra Reference Datta, Mathew and Hemchandra2022). In the majority of these studies, approximate reconstruction methods were employed while in more recent a priori studies iterative and constrained reconstruction algorithms were developed. Such iterative algorithms were used to successfully model a variety of different unresolved terms in the LES transport equations including the scalar fluxes, and the Reynolds stresses (Domingo & Vervisch Reference Domingo and Vervisch2017; Wang & Ihme Reference Wang and Ihme2017; Nikolaou & Vervisch Reference Nikolaou and Vervisch2018; Nikolaou, Vervisch & Cant Reference Nikolaou, Vervisch and Cant2018; Nikolaou, Minamoto & Vervisch Reference Nikolaou, Minamoto and Vervisch2019). Reconstruction-modelling methods use information from the resolved scales in LES without invoking explicit assumptions on the underlying flow field and reaction mode. Even though in principle reconstruction cannot recover scales below the LES mesh, SGS effects can still be incorporated in a variety of approaches. For instance, Stolz & Adams (Reference Stolz and Adams2001) added an additional relaxation term to the velocity transport equation proportional to the difference between the filtered and filtered-reconstructed velocity components 
$(\bar {u}_i-\overline {u^*_i} )$. Wang & Ihme (Reference Wang and Ihme2019) combined reconstruction with interpolation on a finer grid in an effort to enrich the information provided by the resolved scales on the LES mesh while Gullbrand & Chow (Reference Gullbrand and Chow2003) employed mixed models. In yet another approach, Bull & Jameson (Reference Bull and Jameson2016) assumed SGS dissipation to be sufficiently provided by the numerical scheme while reconstruction modelling was employed only for the unresolved LES scales. This modelling approach was applied to LES of turbulent channel flow, and improved the predictions over LES using the classic Smagorinsky model. In more recent works, data-based methods employing super-resolution reconstruction have also been considered (Fukami, Fukagata & Taira Reference Fukami, Fukagata and Taira2020; Kim et al. Reference Kim, Kim, Won and Lee2021). Therefore, the evidence in the literature seems to suggest that reconstruction modelling is a natural framework for modelling SGS effects as well.
With regards to modelling the scalar variance, a successful modelling framework based on iterative and constrained reconstruction was proposed and tested a priori by Nikolaou & Vervisch (Reference Nikolaou and Vervisch2018). This framework provided quantitatively accurate predictions of the variance at a modest computational cost. Nevertheless, iterative reconstruction can be computationally expensive to apply in realistic-geometry LES. Data-based methods on the other hand (Maulik & San Reference Maulik and San2017; Nikolaou & Vervisch Reference Nikolaou and Vervisch2018) require high-quality and relatively large amounts of data. In addition, their computational efficiency depends strongly on the structure (number of layers, number of neurons) of the network, and there exist in fact bounds above which a neural network will always underperform (more floating-point operations) in comparison to using a simple tabulation approach (Nikolaou, Vervisch & Domingo Reference Nikolaou, Vervisch and Domingo2022). One solution to this problem is to use approximate reconstruction operators derived from truncated Taylor-series expansions of the inverse filter operator. Such operators are computationally efficient, and lead to relatively simple analytic expressions that can be used to derive simple algebraic models with the deconvolution truncation order determining the model accuracy. Such approximate reconstruction methods have been considered before for modelling the scalar variance, for instance, by Balarac et al. (Reference Balarac, Pitsch and Raman2008) and Knudsen, Kim & Pitsch (Reference Knudsen, Kim and Pitsch2010). Another approach that was recently proposed by Nikolaou, Vervisch & Domingo (Reference Nikolaou, Vervisch and Domingo2023) is to use inverse explicit discrete filters derived by solving an optimisation problem. This approach avoids the need to iterate as in most reconstruction algorithms and, as a result, reduces the reconstruction time by orders of magnitude (Nikolaou et al. Reference Nikolaou, Vervisch and Domingo2023). Irrespective however of the reconstruction operator used (approximate, inverse filter, iterative), there does not appear to exist a structured and general reconstruction-based methodology for model development. For instance, it is unclear how the reconstructed fields should be employed, what conditions must be satisfied for the model to be physically realisable and how SGS effects can be incorporated in a general reconstruction-based model. The aim of this work is to address these issues by developing a general modelling framework based on reconstruction.
We begin in § 2 by discussing the most important properties of the scalar variance. In § 3, approximate reconstruction operators based on truncated Taylor-series expansions are derived. These operators are used to derive static models all the while maintaining compatibility with lower reconstruction-order models, but also with classic models in the literature. In § 4 a brief overview of the procedure for developing discrete inverse filters is presented, and in § 5 the classic dynamic models are revisited, and the dynamic procedure is reinterpreted in the context of reconstruction-based modelling. Finally, in § 6 the developed models are tested a priori using direct numerical simulation (DNS) data.
2. Variance properties
It is straightforward to show that the variance has the following properties:
\begin{equation} ({\rm i})\quad \lim_{\varDelta \to 0} \sigma^2(\phi;\varDelta)=0. \end{equation}
For 
$a,b$ constant,
\begin{equation} ({\rm ii}) \quad \sigma^2(a\phi+b;\varDelta)=a^2\sigma^2(\phi;\varDelta). \end{equation}
Provided the filter 
$G \ge 0$ and 
$\phi \in [0,1]$,
\begin{equation} ({\rm iii})\quad 0 \leq \sigma^2(\phi;\varDelta) \leq \tilde{\phi}(1-\tilde{\phi})\leq \frac{1}{4}, \end{equation}which we now prove formally. We begin from the definition of the filtering operation in a single dimension,
\begin{equation} \bar{\phi}(x,t)=\int_{-\infty}^{\infty}G(x-s)\phi(s,t)\,{\rm d}s, \end{equation}
where 
$G(s) \ge 0$, 
$G(-s)=G(s)$ (symmetric) and 
$\int _{-\infty }^{\infty }G(x)\,{{\rm d}\kern0.7pt x}=1$, which ensures that constant scalars are unaffected by the filtering (consistency condition). It is important to note at this point that we assume 
$\varDelta$ to be constant, and that the computational mesh is homogeneous. In addition, a Cartesian coordinate system will be used throughout the analysis that follows. Since three-dimensional (3-D) filters are typically constructed using products of one-dimensional (1-D) filters, i.e. 
$G(\boldsymbol {x})=G(x)G(y)G(z)$, and filtering is conducted using dimensional splitting, it is sufficient to derive the variance bounds in the 1-D case. Favre-filtered variables are defined as usual, 
$\tilde {\phi }(x,t)=\overline {\rho (x,t) \phi (x,t)}/ \bar {\rho }(x,t)$, and the equivalent Favre-filtering operation is
\begin{equation} \tilde{\phi}(x,t)=\int_{-\infty}^{\infty}\tilde{G}(s;x,t)\phi(x-s,t)\,{\rm d}s, \end{equation}
where 
$\tilde {G}(s;x,t)=G(s)\rho (x-s,t) / \bar {\rho }(x,t)$ is the corresponding Favre filter. Clearly, 
$\tilde {G} \ge 0$ for any realisable (positive) value of the density, 
$\int _{-\infty }^{\infty }\tilde {G}(s;x,t)\,{\rm d}s=1$ for all 
$x,t \in \mathbb {R}$, and for constant-density flows, 
$\tilde {G}=G$. The lower bound for the variance can be derived by making use of the Cauchy–Schwarz inequality that bounds from above the product 
$\tilde {\phi }\tilde {\phi }$,
\begin{align} \tilde{\phi}(x,t)\tilde{\phi}(x,t) & =\left( \int_{-\infty}^{\infty}\tilde{G}(s;x,t)\phi (x-s,t)\,{\rm d}s \right) ^2\nonumber\\ &=\left( \int_{-\infty}^{\infty}\sqrt{\tilde{G}(s;x,t)}\sqrt{\tilde{G}(s;x,t)}\phi (x-s,t)\,{\rm d}s \right)^2\nonumber\\ & \leq \underbrace{\int_{-\infty}^{\infty}{\tilde{G}(s;x,t)}\,{\rm d}s}_{1} \ \boldsymbol{\cdot}\ \int_{-\infty}^{\infty}{\tilde{G}(s;x,t)}\phi ^2(x-s,t)\,{\rm d}s\nonumber\\ & =\widetilde{\phi ^2}(x,t). \end{align}
Now suppose 
$\phi (x,t)$ is a bounded scalar specifically 
$\phi \in [0,1]$. Such scalars of interest may, for instance, be the progress variable in premixed flames, the mixture fraction in non-premixed flames, the mass fraction, etc. Since 
$\phi \in [0,1]$, 
$\phi ^2 \le \phi$, hence, 
$\widetilde {\phi ^2} \le \tilde {\phi }$ provided 
$\tilde {G}(s;x,t) \ge 0$. Therefore, we have overall 
$\tilde {\phi }\tilde {\phi }\leq \widetilde {\phi ^2} \le \tilde {\phi }$ and 
$0\leq \sigma ^2=\widetilde {\phi ^2}-\tilde {\phi }\tilde {\phi } \leq \tilde {\phi }(1-\tilde {\phi })\leq \tfrac {1}{4}$.
It is important to note that, for positive filters, if 
$\phi \in [0,1]$ then 
$\tilde {\phi } \in [0,1]$ as well irrespective of the distribution of 
$\phi$, but this is not necessarily the case for non-positive filters. A short proof of this is given in Appendix A, and as we will see later on this has some important modelling consequences. In principle, properties (i)–(iii) must be satisfied by any variance model.
3. Approximate reconstruction using truncated Taylor series
Starting from the 3-D filtering operation that, for convenience, we now write as (using the filter symmetry condition)
\begin{equation} \bar{\phi}(\boldsymbol{x})=\int_{-\infty}^{\infty}G(\boldsymbol{r})\phi(\boldsymbol{x}+\boldsymbol{r})\,{\rm d}\boldsymbol{r} \end{equation}
and expanding 
$\phi (\boldsymbol {x}+\boldsymbol {r})$ in Taylor series around 
$\boldsymbol {r}=\boldsymbol {x}$, it is straightforward to show after some algebraic manipulation that (Nikolaou et al. Reference Nikolaou, Vervisch and Domingo2023)
\begin{equation} \bar{\phi}(\boldsymbol{x})=\phi(\boldsymbol{x})+a_2\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\nabla} \phi+O(\varDelta ^4)=\phi(\boldsymbol{x})+a_2\frac{\partial{}}{\partial{x_i}}\left( \frac{\partial{\phi}}{\partial{x_i}} \right)+O(\varDelta ^4), \end{equation}
where 
$a_2=M_{2}/2$ is half the filter's second moment. For a Gaussian filter, 
$a_2=\varDelta ^2/24$ while for a Helmholtz filter, 
$a_2=\varDelta ^2$; table 1 lists two of the most popular filters (Gaussian/Helmholtz) along with their transfer functions, moments and other important parameters. Taking the Fourier transform (here 
$^{\wedge} $ denotes Fourier transform) of the above expansion we obtain
\begin{equation} \hat{\bar{\phi}}(\boldsymbol{k})=\hat{\phi}(\boldsymbol{k})+a_2(\,j\boldsymbol{k}) (\,j\boldsymbol{k})\hat{\phi}(\boldsymbol{k})+O(k ^4), \end{equation}which we can invert approximately to obtain
\begin{equation} \hat{\phi}(\boldsymbol{k})=\hat{\bar{\phi}}(\boldsymbol{k})\left( 1-a_2 \boldsymbol{k} \boldsymbol{k} +O(k^4) \right)^{-1}=\hat{\bar{\phi}}(\boldsymbol{k})+a_2\hat{\bar{\phi}}(\boldsymbol{k})\boldsymbol{k}\boldsymbol{k}+O(k^4). \end{equation}
The above equation is an approximation of the Fourier transform of the inverse filter – an alternative approach to obtain this for a specific filter is to use the exact transform of the inverse filter used. For a Gaussian filter, for instance, 
$\hat {G}^{-1}=e^{\varDelta ^2 k^2 /24}$ that can be expanded in Taylor series resulting in the same expression as above (Balarac et al. Reference Balarac, Pitsch and Raman2008). Taking an inverse Fourier transform of the above we obtain an approximation 
$\phi ^*$ of the original field 
$\phi$,
\begin{equation} \phi^*(\boldsymbol{x})=\bar{\phi}(\boldsymbol{x})-a_2 \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\nabla} \bar{\phi}(\boldsymbol{x})+O(\varDelta ^4). \end{equation}
For a forward filter having a filter width 
$\varDelta$, the above equation is an approximation of the unfiltered field with a truncation order of 
$O(\varDelta ^4)$. As 
$\varDelta \rightarrow 0$, 
$a_2 \rightarrow 0$ while 
$\bar {\phi }(\boldsymbol {x}) \rightarrow \phi (\boldsymbol {x})$, hence, 
$\phi ^*(\boldsymbol {x}) \rightarrow \phi (\boldsymbol {x})$. The above approximation holds for any filter, and higher-order approximations/reconstructions can be obtained by including higher-order derivatives (Leonard Reference Leonard1997; Carati, Winckelmans & Jeanmart Reference Carati, Winckelmans and Jeanmart2001); however, such expressions become substantially more complex for deriving (simple) algebraic models. The reconstruction operation in (3.5) can also be written in terms of a reconstruction operator 
$R_{\varDelta }=(I-a_2\nabla ^2)$ acting on the filtered scalar 
$\phi$, i.e. 
$R_{\varDelta }[\bar {\phi }]$. This notation will prove useful in the sections that follow.
Table 1. Common explicit positive filters with their transfer functions and even moments 
$M_{2r}$.
3.1. Constant-density flows
The approximate reconstruction operator in (3.5) opens up a number of possibilities for modelling the scalar variance. We begin with constant-density flows where 
$\sigma ^2=\overline {\phi ^2}-\bar {\phi }\bar {\phi }$. One option would be to model the variance using 
$\overline {\phi ^{*2}}-\bar {\phi }\bar {\phi }$ since 
$\overline {\phi ^2}$ is the unresolved part. However, this approach involves two distinct fields (
$\phi, \phi ^*$) and properties (ii)–(iii) of § 2 may not always be satisfied. An alternative option is to use
\begin{equation} \sigma^2=\overline{\phi ^{*2}}-\overline{\phi ^*}\ \overline{\phi ^*}, \end{equation}
where 
$\phi ^*$ is the reconstructed scalar on the LES mesh. Clearly, as 
$\varDelta \rightarrow 0$, the filter must tend to a 
$\delta$ function, and property (i) is satisfied. Property (ii) is always satisfied inherently – this follows from the linearity of the filtering operation and the filter's consistency condition. With regards to property (iii), provided 
$\phi ^* \in [0,1]$ then, as shown in Appendix A, 
$\overline {\phi ^*} \in [0,1]$ as well. As a result, the predictions of the above model will be bounded in a similar manner as the actual variance specifically in the region 
$[0,\overline {\phi ^*}(1-\overline {\phi ^*})]$. For a perfect reconstruction, 
$\phi ^* \rightarrow \phi$, and the model bounds approach the actual variance bounds. For an imperfect reconstruction, the model bounds may not necessarily coincide with the actual variance bounds but guarantee that the theoretical upper bound of 1/4 is never exceeded. This is in contrast with models of the gradient type that are not intrinsically bounded, and whose upper bound depends solely on the magnitude of the scalar gradient.
With the above modelling approach, and an 
$O(\varDelta ^2)$ approximation of 
$\phi$, i.e. 
$\phi ^*=\bar {\phi }+O(\varDelta ^2)$, we obtain the model
\begin{equation} \sigma ^2(\phi)=\overline{\bar{\phi}\bar{\phi}}-\bar{\bar{\phi}} \bar{\bar{\phi}}+O(\varDelta ^2), \end{equation}
which is essentially the scale-similarity model of Bardina, Ferzinger & Reynolds (Reference Bardina, Ferzinger and Reynolds1980) applied to a scalar. The bounds of this model are 
$[0,\bar {\bar {\phi }}(1-\bar {\bar {\phi }})]$. We will refer to this model as the SM2 model for constant-density flows. For an 
$O(\varDelta ^4)$ approximation on the other hand, i.e. 
$\phi ^*=\bar {\phi }-a_2 \nabla ^2 \bar {\phi }$, we obtain
\begin{align} \sigma ^2 &=\overline{\left( \bar{\phi}-a_2\nabla^2\bar{\phi} \right) ^2}-\overline{\left( \bar{\phi}-a_2\nabla^2\bar{\phi} \right) } \boldsymbol{\cdot} \overline{\left( \bar{\phi}-a_2\nabla^2\bar{\phi} \right) }+O(\varDelta ^4)\nonumber\\ & =\underbrace{\overline{\bar{\phi}\bar{\phi}}-\bar{\bar{\phi}}\bar{\bar{\phi}}}_{{\rm SM2}}+2a_2\left( \bar{\bar{\phi}} \nabla ^2 \bar{\bar{\phi}}-\overline{\bar{\phi} \nabla ^2 \bar{\phi}} \right)+O(\varDelta ^4). \end{align}
The above model is composed of the SM2 model term in (3.7) ensuring compatibility, and an additional correction term that arises from the higher-order reconstruction terms retained in the expression for 
$\phi ^*$. We will refer to the above model as the SM4 model. Following this modelling procedure it is clear that one may develop higher-order models by retaining higher-order terms in the Taylor-series expansion all the while ensuring backwards compatibility with lower-order models. In contrast with the SM2 model however, the variance bounds may not be satisfied for the SM4 model since the 
$O(\varDelta ^4)$ reconstruction that involves gradients of 
$\tilde {\phi }$ may cause 
$\phi ^*$ to lie outside its physical limits and, as a result, so will the variance predictions. A remedy to this problem will be presented in the sections that follow. For now, we focus our attention on the SM4 model. If we further expand all the double-filtering operations in (3.8) in Taylor series and truncate these to 
$O(\varDelta ^4)$, we obtain
\begin{align} \sigma ^2&= \overline{\bar{\phi}\bar{\phi}}-\bar{\bar{\phi}}\bar{\bar{\phi}}+2a_2\left( \bar{\bar{\phi}} \nabla ^2 \bar{\bar{\phi}}-\overline{\bar{\phi} \nabla ^2 \bar{\phi}} \right)\nonumber\\ & =\bar{\phi}\bar{\phi}+a_2\nabla ^2 \bar{\phi} \bar{\phi}+O(\varDelta ^4)-\left( \bar{\phi}+a_2\nabla^2\bar{\phi}+O(\varDelta ^4) \right) ^2+ \left( \bar{\phi}\nabla^2 \bar{\phi}-\bar{\phi}\nabla^2 \bar{\phi}+O(\varDelta ^4)\right)\nonumber\\ & =a_2\nabla ^2 \bar{\phi} \bar{\phi}-2a_2\bar{\phi}\nabla ^2 \bar{\phi}+O(\varDelta ^4)\nonumber\\ & =2a_2\frac{\partial{\bar{\phi}}}{\partial{x_i}}\frac{\partial{\bar{\phi}}}{\partial{x_i}}+O(\varDelta^4). \end{align}
The above is none other than the gradient model but having a well-defined model parameter given by 
$2a_2=M_{2}$ that, for a Gaussian filter, is 
$\varDelta ^2/12$ – a similar expression was derived for the Reynolds stresses by Clark, Ferzinger & Reynolds (Reference Clark, Ferzinger and Reynolds1979). Therefore, based on the above procedure, we see that the gradient model is a subset model of the more general scale-similarity model in (3.8) (SM4).
3.2. Variable-density flows
For variable-density flows, the modelling procedure becomes somewhat more complicated but the key idea remains the same: we use approximate reconstructions 
$\rho ^*=R_{\varDelta }[\bar {\rho }]$ and 
$(\rho \phi )^*=R_{\varDelta }[\overline {\rho \phi }]$ in order to obtain an expression for the variance. For variable-density flows, in the same spirit as for constant-density flows, we model the variance using
\begin{equation} \sigma^2=\frac{\overline{\rho^* \phi^* \phi^*}}{{\overline{\rho ^*}}}-\left( \frac{{\overline{\rho^* \phi^*}}}{{\overline{\rho ^*}}} \right)^2=\widetilde{\phi^*\phi^*}-\widetilde{\phi^*}\widetilde{\phi^*}, \end{equation}
where 
$\phi ^*=(\rho \phi )^*/\rho ^*$; table 2 lists all the different operations/notations that will be used throughout the rest of the paper for clarity. The above model satisfies conditions (i) and (ii) always and provided 
$\phi ^* \in [0,1]$ the model bounds are 
$[0,\widetilde {\phi ^*}(1-\widetilde {\phi ^*})]$. Therefore, as 
$\phi ^* \rightarrow \phi$, the bounds approach the actual variance bounds. For an imperfect reconstruction on the other hand, these bounds ensure that the upper bound of 1/4 is never exceeded.
Table 2. Operations involving the normal filter 
$\varDelta$ and the test filter 
$\hat {\varDelta }$.
For an 
$O(\varDelta ^2)$ approximation, 
$\rho ^*=\bar {\rho }+O(\varDelta ^2)$ and 
$(\rho \phi )^*=\overline {\rho \phi }+O(\varDelta ^2)$. Inserting these in (3.10) we obtain the equivalent SM2 model for variable-density flows:
\begin{equation} \sigma^2=\frac{\overline{\bar{\rho}\tilde{\phi}\tilde{\phi}}}{{\bar{\bar{\rho}}}}-\left( \frac{{\overline{(\bar{\rho} \tilde{\phi})}}}{{\bar{\bar{\rho}}}}\right)^2=\left( \tilde{\phi}\tilde{\phi} \right) ^{\breve{ }} - \breve{\tilde{\phi}}\breve{\tilde{\phi}}. \end{equation}
Here 
$\breve {\tilde {\phi }}=\overline {\bar {\rho }\tilde {\phi }}/\bar {\bar {\rho }}$ (notation that is extended in the form 
$({\cdot })^{\breve { }} = \overline {\bar {\rho }\widetilde {{\cdot }}}/\bar {\bar {\rho }}$ in (3.11)). The bounds of the above model are 
$[0,\breve {\tilde {\phi }}(1-\breve {\tilde {\phi }})]$. Note that, for constant-density flows, 
$\breve {\tilde {\phi }}=\bar {\bar {\phi }}$ and (3.11) reduces to (3.7) thus ensuring compatibility with the constant-density case. For an 
$O(\varDelta ^4)$ approximation, the reconstruction operator is given by
\begin{equation} R_{\varDelta}[\bar{u}]=\bar{u}-a_2\nabla ^2 \bar{u}; \end{equation}
therefore, 
$\rho ^*=\bar {\rho }-a_2\nabla ^2 \bar {\rho }+O(\varDelta ^4)$, 
$(\rho \phi )^*=\overline {\rho \phi }-a_2\nabla ^2 \overline {\rho \phi }+O(\varDelta ^4)$, and it is straightforward to show after some algebraic manipulation that
\begin{equation} \left( \frac{\overline{(\rho \phi)^*}}{\overline{\rho^*}}\right)^2=\breve{\tilde{\phi}}\breve{\tilde{\phi}}+2\frac{a_2}{\bar{\bar{\rho}}}\breve{\tilde{\phi}}\breve{\tilde{\phi}}\nabla^2\bar{\bar{\rho}}-2\frac{a_2}{\bar{\bar{\rho}}}\breve{\tilde{\phi}}\nabla^2\overline{\bar{\rho}\tilde{\phi}} +O(\varDelta^4) \end{equation}and
\begin{align} (\rho \phi^2)^* & =\frac{(\rho \phi)^* (\rho \phi)^*}{\rho ^*} \nonumber\\ & =\frac{ \left( \overline{\rho \phi}-a_2\nabla ^2 \overline{\rho \phi}+O(\varDelta ^4) \right) ^2}{\bar{\rho}-a_2\nabla ^2 \bar{\rho}+O(\varDelta ^4)}\nonumber\\ & =\frac{1}{\bar{\rho}}\left( \overline{\rho \phi}{\cdot} \overline{\rho \phi}-2a_2\overline{\rho \phi}\nabla ^2 \overline{\rho \phi}+O(\varDelta ^4) \right) \left(1-\frac{a_2}{\bar{\rho}}\nabla ^2\bar{\rho}+O(\varDelta ^4) \right)^{-1}\nonumber\\ & =\frac{1}{\bar{\rho}}\left( \overline{\rho \phi}{\cdot} \overline{\rho \phi}-2a_2\overline{\rho \phi}\nabla ^2 \overline{\rho \phi}+O(\varDelta ^4) \right) \left(1+\frac{a_2}{\bar{\rho}}\nabla ^2\bar{\rho}+O(\varDelta ^4) \right)\nonumber\\ & =\bar{\rho}\tilde{\phi}\tilde{\phi}+a_2\tilde{\phi}\tilde{\phi} \nabla ^2\bar{\rho}-2a_2\tilde{\phi}\nabla ^2\bar{\rho}\tilde{\phi}+O(\varDelta ^4), \end{align}
which constitutes an 
$O(\varDelta ^4)$ approximation of 
$\rho \phi ^2$. Inserting the above in (3.10) we eventually obtain the equivalent SM4 model for variable-density flows:
\begin{align} \sigma ^2=\underbrace{\left( \tilde{\phi}\tilde{\phi} \right) ^{\breve{ }} - \breve{\tilde{\phi}}\breve{\tilde{\phi}}}_{{\rm SM2}}+\frac{2a_2}{\bar{\bar{\rho}}}\left( \breve{\tilde{\phi}} \nabla ^2 \overline{\bar{\rho}\tilde{\phi}}-\overline{\tilde{\phi}\nabla^2\bar{\rho}\tilde{\phi}} \right) +\frac{a_2}{\bar{\bar{\rho}}}\Big( \overline{\tilde{\phi}\tilde{\phi}\nabla^2\bar{\rho}}+\Big(\tilde{\phi}\tilde{\phi}\Big)^{\breve{}} \nabla^2\bar{\bar{\rho}}-2\breve{\tilde{\phi}}\breve{\tilde{\phi}} \nabla^2\bar{\bar{\rho}}\Big)\!. \end{align}
The above model is composed of the scale-similarity term of the SM2 model in (3.11) augmented by two more correction terms. For constant-density flows, the entire third term in the parentheses disappears, in which case we recover (3.8). If we further expand all filtering operations up to 
$O(\varDelta ^4)$ it is straightforward to show that
\begin{gather} \breve{\tilde{\phi}}=\tilde{\phi}+a_2\nabla^2\tilde{\phi}+\frac{2a_2}{\bar{\rho}}\frac{\partial{\bar{\rho}}}{\partial{x_i}}\frac{\partial{\tilde{\phi}}}{\partial{x_i}}+O(\varDelta^4), \end{gather}
\begin{gather} \breve{\tilde{\phi}}\breve{\tilde{\phi}}=\tilde{\phi}\tilde{\phi}+2a_2\tilde{\phi}\nabla^2\tilde{\phi}+\frac{4a_2}{\bar{\rho}}\tilde{\phi}\frac{\partial{\bar{\rho}}}{\partial{x_i}}\frac{\partial{\tilde{\phi}}}{\partial{x_i}}+O(\varDelta^4) \end{gather}and
\begin{equation} \left( \tilde{\phi} \tilde{\phi} \right)^{\breve{}}=\tilde{\phi}\tilde{\phi}+2a_2\tilde{\phi}\nabla^2\tilde{\phi}+2a_2\frac{\partial{\tilde{\phi}}}{\partial{x_i}}\frac{\partial{\tilde{\phi}}}{\partial{x_i}} +\frac{4a_2}{\bar{\rho}}\tilde{\phi}\frac{\partial{\bar{\rho}}}{\partial{x_i}}\frac{\partial{\tilde{\phi}}}{\partial{x_i}}+O(\varDelta^4). \end{equation}
Since the second and third terms in (3.15) already include a 
$\varDelta ^2$ term in the parameter 
$a_2$ outside of the parentheses, all terms inside the parentheses must be expanded up to 
$O(\varDelta ^2)$ only to keep the overall expansion order to 
$O(\varDelta ^4)$, i.e. 
$\breve {\tilde {\phi }}=\tilde {\phi }$, 
$\bar {\bar {\rho }}=\bar {\rho }$, etc., causing them to vanish. As a result, after some algebraic manipulation the approximation of (3.15) eventually reduces to
\begin{equation} \sigma^2=2a_2\frac{\partial{\tilde{\phi}}}{\partial{x_i}}\frac{\partial{\tilde{\phi}}}{\partial{x_i}}+O(\varDelta ^4), \end{equation}which is the gradient model for variable-density flows. Therefore, the gradient model is a subset model of the more general reconstruction-based SM4 model of (3.15). It is important to note that this derivation of the gradient model is based on a purely mathematical approach using approximate reconstruction operators. No explicit assumptions on the nature/structure of the scalar (passive, reactive) and/or the flow field were made in this process. This implies that the gradient model alone should in principle be sufficient to model both non-reactive and reactive scalars as evidenced by previous LES of turbulent premixed flames (Vreman et al. Reference Vreman, Albrecht, van Oijen, de Goey and Bastiaans2008, Reference Vreman, Bastiaans and Geurts2009a,Reference Vreman, van Oijen, de Goey and Bastiaansb; Mukhopadhyay et al. Reference Mukhopadhyay, van Oijen and de Goey2015) that employed either the static or the dynamic gradient model.
3.3. Post-reconstruction bounding
For a general reconstruction operator 
$R_{\varDelta }$, (3.10) may not necessarily satisfy the variance bounds because 
$\phi ^*=(\rho \phi )^*/ \rho ^*$ may not lie in the range 
$[0,1]$. Possible reasons for this may include 
$\rho ^*$, 
$(\rho \phi )^*$ lying outside their physical minimum/maximum limits because of approximate reconstruction, but also due to numerical errors. A straightforward remedy to this problem is to bound the reconstructed fields. Therefore, (3.10) is still employed, but using the bounded reconstructed fields instead such that 
$\phi ^* \in [0,1]$ by construction. Bounding serves two purposes: (i) that 
$\widetilde {\phi ^*} \in [0,1]$, and (ii) that the bounds for (3.10) are 
$[0,\widetilde {\phi ^*}(1-\widetilde {\phi ^*})]$ always.
One approach to implement bounding is directly in the reconstruction operator 
$R_{\varDelta }$, i.e. use a constrained reconstruction algorithm, and such algorithms were employed in previous a priori (Wang & Ihme Reference Wang and Ihme2017; Nikolaou & Vervisch Reference Nikolaou and Vervisch2018) and a posteriori (Wang & Ihme Reference Wang and Ihme2019) studies. Constrained reconstruction guarantees that 
$\phi ^* \in [0,1]$, but is generally computationally more expensive than unconstrained reconstruction, and may also hinder the convergence of the deconvolution process. A computationally more efficient approach is to use post-reconstruction bounding. Provided the reconstruction operator is sufficiently accurate, post-reconstruction bounding helps to bound the minority of points lying outside their physical boundaries. A general approach for doing so is to use two levels of bounding: a physical one and a numerical one. For a general reconstruction operator 
$R_{\varDelta }$, physical bounding may be applied in a computationally efficient approach using
\begin{gather} \rho^*=\max \left( \min \left( R_{\varDelta}[\bar{\rho}],\rho_h \right)\!,\rho_l \right)\!, \end{gather}
\begin{gather} (\rho\phi)^*=\max \left( \min \left( R_{\varDelta}[\overline{\rho \phi}],(\rho\phi)_h \right)\!,(\rho\phi)_l \right)\!, \end{gather}
where the subscripts 
$l,h$ denote the corresponding low and high values. Alternatively, a smoother function other than maxmin can also be employed. Physical bounding ensures that the reconstructed fields lie within their physical limits. For instance, in the case of premixed flames 
$\rho _l=\rho _p$ the density of the hot products, and 
$\rho _h=\rho _r$ the density of the cold reactants while 
$(\rho \phi )_l=0$ and 
$(\rho \phi )_h=\rho _r$. Following physical bounding, numerical bounding can also be applied directly on 
$\phi ^*$ using
\begin{equation} \phi^*=\max\left(\min \left( (\rho \phi)^*/\rho^*,1 \right)\!,0 \right), \end{equation}
which ensures that any numerical errors in evaluating 
$\phi ^*$ from the physically bounded reconstructions do not cause it to exceed the limits 
$[0,1]$. As we show later on, the above approach substantially improves the reconstructions, hence, the variance predictions. To illustrate the procedure, we consider an 
$O(\varDelta ^4)$ reconstruction operator, 
$R_{\varDelta }=(I-a_2 \nabla ^2)$, in which case we have
\begin{gather} \rho^*=\max \left( \min \left( (I-a_2 \nabla^2)[\bar{\rho}],\rho_h \right)\!,\rho_l \right)\!, \end{gather}
\begin{gather} (\rho\phi)^*=\max \left( \min \left( (I-a_2 \nabla^2)[\overline{\rho \phi}],(\rho\phi)_h \right),(\rho\phi)_l \right)\!. \end{gather}
Calculating 
$\phi ^*$ using (3.22)–(3.24), and inserting all in (3.10) we obtain a scalar variance model bounded exactly in the range 
$[0,\widetilde {\phi ^*}(1-\widetilde {\phi ^*})]$ – we will refer to this particular model as AD4. Note however that the above procedure is general, and can be applied for any general reconstruction operator 
$R_{\varDelta }$ that may well be based on any reconstruction algorithm.
4. Inverse filter reconstruction
Higher-order reconstruction operators 
$R_{\varDelta }$ can be obtained using iterative algorithms such as van Cittert iterations (van Cittert Reference van Cittert1931). A constrained version of this algorithm was introduced in previous work (Nikolaou & Vervisch Reference Nikolaou and Vervisch2018) where it was used to successfully model the progress-variable variance in turbulent premixed flames. In iterative algorithms the iteration count 
$N$ determines the extent of reconstruction with higher 
$N$ corresponding to higher-order reconstruction (larger wavenumbers are recovered). Compared with approximate reconstruction methods such as the one presented in the previous section, iterative methods have superior performance albeit at an increased computational cost that depends on 
$N$, the stencil size of the filter, but also the problem dimension – an expression quantifying this has been developed by Nikolaou et al. (Reference Nikolaou, Vervisch and Domingo2023). An alternative computationally efficient yet high-order reconstruction process recently proposed by Nikolaou et al. (Reference Nikolaou, Vervisch and Domingo2023) is based on the concept of inverse discrete filters. In this approach, the same way a discrete forward explicit filter is used to filter a signal, an inverse discrete explicit filter is applied on the filtered signal in order to recover a reconstruction estimate. This approach was demonstrated to be computationally more efficient (by orders of magnitude) than iterative reconstruction while maintaining the same level of accuracy (Nikolaou et al. Reference Nikolaou, Vervisch and Domingo2023).
The procedure for obtaining an optimised inverse filter begins by first obtaining an optimised forward filter. We begin from the 1-D discrete and explicit filtering operation defined as
\begin{equation} \bar{u} _i = G_d \circledast u_i =\sum_{l=-M}^{M} g_l u _{i+l}, \end{equation}
where 
$M$ is the half-stencil size of the filter and 
$g_l$ are the discrete filter coefficients for each point on the stencil (we will be dealing with symmetric filters). With a mesh spacing 
$h$, 
$x_i=ih$, 
$i \in [0,N_x-1]$, the transfer function 
$\hat {G}_d$ of the above discrete filtering operation is given by
\begin{equation} \hat{G}_d(k_r)=\sum_{l=-M}^{M}g_l{\rm e}^{j k_rh l}=g_0+2\sum_{l=1}^{M}g_l \cos \left( k_rh l \right)\!, \end{equation}
where 
$k_r=2{\rm \pi} r/ (N_xh)$. The coefficients 
$g_l$ in (4.2) completely determine 
$\hat {G}_d$ that ideally should match 
$\hat {G}$ (the actual filter transfer function) for all 
$k_rh \in [0,{\rm \pi} ]$. Standard methods for obtaining the filter coefficients include Newton–Cotes rules, using filters derived from truncated Taylor expansions, etc. As discussed in Nikolaou et al. (Reference Nikolaou, Vervisch and Domingo2023), both of these approaches have some important drawbacks. For instance, filters based on Newton–Cotes rules generally have poorer responses at higher wavenumbers, and are not consistent (
$\hat {G}_d(0) \neq 1$). Filters based on Taylor expansions on the other hand are consistent, but their predictive ability deteriorates above a threshold value of 
$\gamma =\varDelta /h$. In addition, in both of these approaches there is no explicit error control. In order to circumvent such issues, but at the same time to automate the process for developing discrete filters, an optimisation procedure was proposed to obtain the coefficients 
$g_l$. The optimisation problem reads
\begin{equation} \underset{g_l}{\text{arg min}} \left( \int_{kh=0}^{\rm \pi}{\left( \hat{G}_d(kh;g_l)-\hat{G}(kh; \gamma) \right)}^2 d(kh) \right). \end{equation}
To ensure consistency (
$\hat {G}_d(0)=1$), the following constraint is applied:
\begin{equation} g_0+2\sum_{l=1}^{M}g_l=1. \end{equation}
To ensure 
$0 \leq \hat {G}_d \le 1$ (positive transfer function), an additional constraint is introduced,
\begin{equation} \hat{G}({\rm \pi}) \leq g_0+2\sum_{l=1}^{M}g_l\cos(k_rh l) \leq 1 \quad \text{for all } k_rh \in (0,{\rm \pi}]. \end{equation}
The optimisation problem can be solved by providing the actual transfer function 
$\hat {G}$ and half-stencil size 
$M$, but can also be solved adaptively. In the adaptive version, which is the approach used here, 
$\hat {G}$ and the target optimisation error are provided. An error controller then adjusts 
$M$ until the target optimisation error is obtained. Once the optimised 1-D filter coefficients are obtained, for the 3-D case the filtering operation is calculated using dimensional splitting, i.e. the optimised 1-D filter coefficients are used to filter in each coordinate direction.
In a similar fashion, we can also define a discrete reconstruction operation (de-filtering) as
\begin{equation} {R_{\varDelta}[\bar{u}_i]=u^* _i = V_d^N \circledast \bar{u}_i =\sum_{l=-M_{IF}}^{M_{IF}} \beta _l \bar{u} _{i+l}}, \end{equation}
where 
$V_d ^N$ is an approximate discrete inverse filter now acting on the filtered field. The superscript 
$N$ denotes that this filter corresponds to the same operator obtained by applying 
$N$ van Cittert iterations on the filtered field while 
$M_{IF}$ is the stencil size of the inverse filter. It is straightforward to show that the transfer function of 
$V_d^N$ is given by (Nikolaou et al. Reference Nikolaou, Vervisch and Domingo2023)
\begin{equation} \hat{V}_d(k_r)=\sum_{l=-M_{IF}}^{M_{IF}}\beta _l{\rm e}^{j k_rh l}= g_0+2\sum_{l=1}^{M_{IF}}\beta _l \cos \left( k_rh l \right), \end{equation}
and is completely determined by the coefficients 
$\beta _l$. In order to obtain these coefficients we solve the optimisation problem
\begin{equation} \underset{\beta_l}{\text{arg min}} \left( \int_{kh=0}^{\rm \pi}{\left( \hat{V}^N_d(kh;\beta_l)\hat{G}_d(kh;g_l)-\hat{Q}_d^N(kh;g_l) \right)}^2 d(kh) \right)\!, \end{equation}
where 
$\hat {Q}_d^N(kh;g_l)$ is the discrete transfer function of the reconstructed signal obtained using 
$N$ van Cittert iterations with a forward discrete filter 
$G_d$,
\begin{equation} \hat{Q}^N_d(k_r)=1-\left( 1-b\hat{G}_d(k_r) \right)^N \left( 1-\hat{G}_d(k_r) \right)\!. \end{equation}
Note that 
$b$ is a constant typically taken to equal unity. We further impose the constraints,
\begin{gather} \beta_0+2\sum_{l=1}^{M_{IF}}\beta_l=1, \end{gather}
\begin{gather} \beta_0+2\sum_{l=1}^{M_{IF}}\beta_l\cos(k_rh l) < (N+1)\quad \text{for all}, k_rh \in (0,{\rm \pi}]. \end{gather}
The first constraint ensures that 
$\hat {V}_d^{N}(0)=1$ and the second ensures that the upper bound is not exceeded (Nikolaou et al. Reference Nikolaou, Vervisch and Domingo2023). It is important to note that the half-stencil size of the inverse filter 
$M_{IF}$ need not be equal to the half-stencil size of the forward filter 
$M$. In practice, we take the following steps in order to compute optimised inverse filter coefficients 
$\beta _l$.
(i) For a given actual filter transfer function
$\hat {G}$, solve the optimisation problem in (4.3) to obtain an optimised discrete forward filter with coefficients $g_l$ and a transfer function $\hat {G}_d$.(ii) Specify the number of van Cittert iterations
$N$ and calculate $\hat {Q}_d^{N}(k_rh;gl)$.(iii) Specify the optimisation error and solve the optimisation problem in (4.8) to obtain
$\beta _l$, each time adjusting the stencil size $M_{IF}$ until the target optimisation error is obtained.
Reconstruction in the 3-D case is performed by applying the 1-D inverse filter operator in each coordinate direction. The reconstruction operator can be applied on 
$\bar {\rho }$ and 
$\overline {\rho \phi }$ to obtain higher-order reconstructions simply by increasing 
$N$ and solving a new optimisation problem to obtain the corresponding coefficients, and the variance modelled following the procedure outlined in § 3. In other words, the modelling procedure remains the same, all that changes is the reconstruction operator 
$R_{\varDelta }$ for which we now have an automated procedure through which its accuracy can be made arbitrarily high/low.
4.1. Optimised inverse filters
The number of iterations 
$N$ is an input parameter for the above optimisation problem. In practical LES, and using the classic iterative van Cittert algorithm, 3–5 iterations were used in previous works (Stolz & Adams Reference Stolz and Adams1999, Reference Stolz and Adams2001). In fact, 
$N=5$ is sufficient to recover the filter's cutoff wavenumber 
$k_f$ defined as 
$\hat {G}(k_f)=0.5$. For larger 
$N$, wavenumbers that are damped by more than 50 % are recovered, which may not be desirable in practical LES. For instance, the performance of a second-order accurate centred numerical scheme (typical of most LES solvers) drops substantially after about 
$k h>1$ (Lele Reference Lele1992) and these wavenumbers should not be recovered. In practical LES, 
$\varDelta /h$ that determines 
$k_f$ should in principle be chosen such that wavenumbers beyond which the accuracy of the numerical scheme deteriorates are not recovered. In this work, 
$\varDelta /h=4$ is used that results in 
$k_fh \simeq 1$, and which ensures that the filtered fields are well resolved as per the criterion developed in earlier work (Nikolaou & Vervisch Reference Nikolaou and Vervisch2018). As a result, we also use 
$N=5$ to solve the optimisation problem in an effort to mimic what would be done in actual LES. Note that in principle, we can increase 
$N$ and obtain ever improved models; however, in practice, as 
$N$ increases, so does the range of higher wavenumbers being recovered. As a result, an ever-increasing stencil size 
$M_{IF}$ would be required that would make reconstruction more expensive – the exact 
$M_{IF}$ threshold of when this process becomes more expensive than using iterative reconstruction has been derived in our earlier work (Nikolaou et al. Reference Nikolaou, Vervisch and Domingo2023). The choice 
$N=5$ results in an explicit inverse filter for the 
$\gamma =4$ case with 
$M_{IF}=M=5$, and the inverse filter reconstruction being about 5 times faster than iterative reconstruction according to the criterion derived (and validated) by Nikolaou et al. (Reference Nikolaou, Vervisch and Domingo2023).
Tables 3 and 4 list the important parameters for the optimised forward and inverse discrete Gaussian filters obtained by solving the above two optimisation problems for 
$\gamma =4$ and 
$\gamma =8$, respectively, while figure 1 shows the corresponding transfer functions of the actual/discrete forward/inverse filters. Note how the filter cutoff wavenumber increases when applying the inverse filter (almost doubles) that signifies the reconstruction of the damped wavenumbers. Clearly, both at the filtering and reconstruction levels the optimised discrete filters closely match the actual corresponding transfer functions. Application of the inverse filter to 
$\bar {\rho }$, 
$\overline {\rho \phi }$, and then following the procedure outlined in § 3 leads to a new base static model for the scalar variance known as the deconvolution-inverse filter (DEIF) model. For more details on the numerical implementation, the effect of using different filters (Helmholtz, implicit, etc.), different iteration counts 
$N$ and so on, we refer the reader to our earlier work (Nikolaou et al. Reference Nikolaou, Vervisch and Domingo2023).
Table 3. Optimised discrete explicit forward and inverse filter coefficients and properties for 
$\gamma =4$. Note that the filters are symmetric, hence, 
$g_{-l}=g_l$ and 
$\beta _{-l}=\beta _l$.
Table 4. Optimised discrete explicit forward and inverse filter coefficients and properties for 
$\gamma =8$. Note that the filters are symmetric, hence, 
$g_{-l}=g_l$ and 
$\beta _{-l}=\beta _l$.
Figure 1. Actual transfer function of a Gaussian filter 
$\hat {G}$, actual reconstructed transfer function 
$\hat {Q}^N$ obtained using 
$N=5$ van Cittert iterations, discrete transfer function of optimised filter 
$\hat {G}_d$, discrete reconstructed transfer function 
$\hat{Q}_d^N$ obtained using van Cittert iterations and the product between the discrete forward filter transfer function and the optimised inverse discrete filter transfer function 
$\hat {G}_d\hat{V}_d^N$: (a) 
$\gamma =4$ and (b) 
$\gamma =8$.
5. Dynamic models
5.1. The classic dynamic gradient model revisited
The purpose of this section is to discuss the well-known dynamic gradient model. This will serve as a prelude to § 5.2 where the dynamic procedure will be reinterpreted in the context of reconstruction modelling. In turn, this will enable us to derive a structured framework for developing general dynamic models for arbitrary reconstruction operators.
The dynamic procedure was originally introduced by Germano et al. (Reference Germano, Piomelli, Moin and Cabot1991) for modelling the unresolved Reynolds stresses, specifically to correct the predictions of the static Smagorinsky model (Smagorinsky Reference Smagorinsky1963) that was found to be too dissipative in the near-wall region. The procedure was later used to develop a dynamic variance model by Pierce & Moin (Reference Pierce and Moin1998). In this procedure, the base static gradient model is corrected using a dynamic parameter 
$C_d$,
\begin{equation} \sigma ^2 = C_d \varDelta ^2 \frac{\partial{\tilde{\phi}}}{\partial{x_i}}\frac{\partial{\tilde{\phi}}}{\partial{x_i}}, \end{equation}
which must be determined. The original dynamic evaluation of 
$C_d$ essentially involved the following two assumptions.
(i) The dynamic parameter
$C_d$ remains constant between the two filtering levels, at least locally and over the region spanned by the test filter $\hat {\varDelta }$.(ii) The filter itself is idempotent, i.e.
$\hat {\bar {\phi }}=\hat {\phi }$.
Here 
$\hat {\ }$ denotes filtering with a filter having a characteristic length scale 
$\hat {\varDelta }> \varDelta$ (typically 
$\hat {\varDelta }=2\varDelta$). We now derive the original expression for 
$C_d$ formally using the above two assumptions, and next we derive an alternative expression for the model's dynamic parameter irrespective of the filter properties. To proceed, we write (5.1) as
\begin{equation} {\overline{\rho \phi ^2}}-{\bar{\rho} \tilde{\phi} \tilde{\phi}}= C_d \varDelta ^2 \bar{\rho} \frac{\partial{\tilde{\phi}}}{\partial{x_i}}\frac{\partial{\tilde{\phi}}}{\partial{x_i}}, \end{equation}which we test filter to obtain
\begin{equation} \widehat{\overline{\rho \phi ^2}}-\widehat{\bar{\rho} \tilde{\phi} \tilde{\phi}}= C_d \varDelta ^2 { \left( \bar{\rho} \frac{\partial{\tilde{\phi}}}{\partial{x_i}}\frac{\partial{\tilde{\phi}}}{\partial{x_i}} \right)}^{\hat{ }}. \end{equation}
Defining 
$\check {\tilde {\phi }}=\widehat {\bar {\rho }\tilde {\phi }}/ \hat {\bar {\rho }}$ (for filtering operations, refer to table 2), the above can be neatly written as
\begin{equation} \check{\widetilde{\phi ^2}}-\left( \tilde{\phi}\tilde{\phi} \right) ^{\check{}}= C_d \varDelta ^2 { \left( \frac{\partial{\tilde{\phi}}}{\partial{x_i}}\frac{\partial{\tilde{\phi}}}{\partial{x_i}} \right)}^{\check{ }}, \end{equation}
which involves the unresolved quantity 
$\check {\widetilde {\phi ^2}}$. The objective is to obtain a second expression involving 
$\check {\widetilde {\phi ^2}}$ and solve for 
$C_d$. Provided the filter is idempotent we may immediately write an expression for the variance at the test-filter level, i.e.
\begin{equation} \frac{\widehat{\rho \phi ^2}}{\hat{\rho}}-\frac{\widehat{\rho \phi}}{\hat{\rho}}\frac{\widehat{\rho \phi}}{\hat{\rho}}=C_d\hat{\varDelta}^2\frac{\partial{(\widehat{\rho \phi} / \hat{\rho})}}{\partial{x_i}}\frac{\partial{(\widehat{\rho \phi} / \hat{\rho})}}{\partial{x_i}}, \end{equation}which because of the filter idempotency assumption is equivalent to
\begin{equation} \frac{\widehat{\overline{\rho \phi ^2}}}{\hat{\bar{\rho}}}-\frac{\widehat{\overline{\rho \phi}}}{\hat{\bar{\rho}}}\frac{\widehat{\overline{\rho \phi}}}{\hat{\bar{\rho}}}=C_d\hat{\varDelta}^2\frac{\partial{(\widehat{\overline{\rho \phi}} / \hat{\bar{\rho}})}}{\partial{x_i}}\frac{\partial{(\widehat{\overline{\rho \phi}} / \hat{\bar{\rho}})}}{\partial{x_i}} \end{equation}or
\begin{equation} \check{\widetilde{\phi ^2}}-\check{\tilde{\phi}}\check{\tilde{\phi}}=C_d\hat{\varDelta}^2\frac{\partial{\check{\tilde{\phi}}}}{\partial{x_i}}\frac{\partial{\check{\tilde{\phi}}}}{\partial{x_i}}. \end{equation}
It is important to note that, for non-idempotent filters, the equivalent filter incorporating the effects of both the original and test filters should be used in order to write down an expression for the variance at the test-filter level – in the original works by Germano et al. (Reference Germano, Piomelli, Moin and Cabot1991) this was denoted using 
$\hat {\bar {\varDelta }}$, where 
$\bar {\varDelta }$ stood for the original filter (denoted here as 
$\varDelta$).
Equations (5.4) and (5.7) can now be used to solve for 
$C_d$. In practice however, using (5.4) and (5.7) to calculate 
$C_d$ may lead to instabilities (Pierce & Moin Reference Pierce and Moin1998). The alternative approach used by Pierce & Moin (Reference Pierce and Moin1998) is to calculate 
$C_d$ using
\begin{equation} C^M_d=\frac{\langle LM^P \rangle}{\langle M^PM^P \rangle} \end{equation}
instead. In the expression above, 
$L=\widehat {\bar {\rho }\tilde {\phi }\tilde {\phi }}-\hat {\bar {\rho }} \check {\tilde {\phi }}\check {\tilde {\phi }}$ is the Leonard term and 
$M^P={\hat {\varDelta }}^2{\hat {\bar {\rho }}|\boldsymbol {\nabla } \check {\tilde {\phi }}|^2}-\varDelta ^2 \widehat {\bar {\rho }|\boldsymbol {\nabla } \tilde {\phi }|^2}$ is the model term. This approach is analogous to the least-squares minimisation approach proposed by Lilly (Reference Lilly1992) for calculating the dynamic parameter in the Smagorinsky model for the Reynolds stresses (excluding the averaging in the nominator and denominator). The brackets 
$\langle \rangle$ typically denote averaging in statistically homogeneous directions, and serve to further stabilise 
$C_d$ (Pierce & Moin Reference Pierce and Moin1998). If no homogeneous directions exist, an alternative integral-based formulation can be used instead proposed by Ghosal et al. (Reference Ghosal, Lund, Moin and Akselvoll1995).
For non-idempotent filters, we can use the results of § 3 to derive an approximate expression of the unresolved quantity 
$\check {\widetilde {\phi ^2}}$ without invoking the idempotency assumption. A similar procedure was employed by Balarac et al. (Reference Balarac, Pitsch and Raman2008) in order to derive an improved expression for the dynamic parameter of a gradient-based model for the scalar variance in the case of constant-density flows. Extending the procedure to variable-density flows we have, from (3.19),
\begin{equation} \bar{\rho}\widetilde{\phi ^2}=\bar{\rho} \tilde{\phi}\tilde{\phi}+2a_2\bar{\rho} \frac{\partial{\tilde{\phi}}}{\partial{x_i}}\frac{\partial{\tilde{\phi}}}{\partial{x_i}}, \end{equation}which we test filter to obtain
\begin{equation} \widehat{\bar{\rho}\widetilde{\phi ^2}}=\widehat{\bar{\rho} \tilde{\phi}\tilde{\phi}}+2a_2{\left( \bar{\rho}\frac{\partial{\tilde{\phi}}}{\partial{x_i}}\frac{\partial{\tilde{\phi}}}{\partial{x_i}} \right) ^{\hat{}}} \end{equation}or
\begin{equation} \check{\widetilde{\phi ^2}}={ \left( \tilde{\phi}\tilde{\phi} \right) ^{\check{}}}+2a_2{\left( \frac{\partial{\tilde{\phi}}}{\partial{x_i}}\frac{\partial{\tilde{\phi}}}{\partial{x_i}} \right) ^{\check{}}}=\check{\tilde{\phi}}\check{\tilde{\phi}}+2a'_2\frac{\partial{\check{\tilde{\phi}}}}{\partial{x_i}}\frac{\partial{\check{\tilde{\phi}}}}{\partial{x_i}}+2a_2{\left( \frac{\partial{\tilde{\phi}}}{\partial{x_i}}\frac{\partial{\tilde{\phi}}}{\partial{x_i}} \right) ^{\check{}}}, \end{equation}
where the second step follows by expanding 
${ ( \tilde {\phi }\tilde {\phi } ) ^{\check {}}}$ using (3.19) at the test-filter level (
$a_2'$ corresponding to the test-filter second moment using 
$\hat {\varDelta }$ instead of 
$\varDelta$). In order to correct the above base static model by taking into account SGS effects through a dynamic parameter 
$C_d$, the dynamic model should be proportional to both 
$a_2{( \partial {\tilde {\phi }}/\partial {x_i} \partial {\tilde {\phi }}/\partial {x_i} )^{\check {}}}$ and 
$(a'_2\partial {\check {\tilde {\phi }}}/\partial {x_i}\partial {\check {\tilde {\phi }}}/\partial {x_i})$. Since 
$a_2$ is proportional to 
$\varDelta ^2$, and 
$a'_2$ is proportional to 
$\hat {\varDelta }^2$ the dynamic formulation should read
\begin{equation} \check{\widetilde{\phi ^2}}-\check{\tilde{\phi}}\check{\tilde{\phi}}=C_d \hat{\varDelta} ^2 \frac{\partial{\check{\tilde{\phi}}}}{\partial{x_i}}\frac{\partial{\check{\tilde{\phi}}}}{\partial{x_i}}+C_d\varDelta ^2{\left( \frac{\partial{\tilde{\phi}}}{\partial{x_i}}\frac{\partial{\tilde{\phi}}}{\partial{x_i}} \right) ^{\check{}}}. \end{equation}
In contrast to (5.7), (5.12) contains an additional term that is the filtered gradient product. Solving for 
$C_d$ using (5.4) and (5.12) and regularising, we obtain an alternative expression for the dynamic parameter,
\begin{equation} C^B_d=\frac{\langle LM^B \rangle}{\langle M^BM^B \rangle}. \end{equation}
In the above expression, the Leonard term remains the same, but the model term is now 
$M^B={\hat {\varDelta }}^2{\hat {\bar {\rho }}|\boldsymbol {\nabla } \check {\tilde {\phi }}|^2}$. For constant-density flows, (5.13) reduces to that of Balarac et al. (Reference Balarac, Pitsch and Raman2008), which was shown to provide improved predictions for the scalar variance in the case of a non-reactive scalar under forced homogeneous isotropic turbulence. The assumptions involved in deriving (5.8) were shown by Balarac et al. (Reference Balarac, Pitsch and Raman2008) not to be verified, at least when using non-idempotent filters as was done in the original derivation of the method (Pierce & Moin Reference Pierce and Moin1998). This leads to the p.d.f. of 
$C_d$ having substantial probabilities/contributions over a wide range of negative values that of course leads to negative, i.e. unrealisable variance predictions (Balarac et al. Reference Balarac, Pitsch and Raman2008). In contrast, 
$C^B_d$ is always positive by construction. Furthermore, the evaluation of the dynamic parameter using (5.13) was shown to be more robust (Balarac et al. Reference Balarac, Pitsch and Raman2008). For reacting scalars, a thorough a priori testing of both the 
$C^M_d$ and 
$C^B_d$ formulations was conducted by Knudsen et al. (Reference Knudsen, Kim and Pitsch2010) using DNS data of statistically planar turbulent premixed flames. The testing was conducted for a wide range of 
$\varDelta$ values. For small 
$\varDelta$, (5.13) improved the predictions substantially in line with the results of Balarac et al. (Reference Balarac, Pitsch and Raman2008). For larger 
$\varDelta$, it was found that both formulations, but also even the static gradient model, all resulted in similar performances. On this point, it is interesting to note that Mukhopadhyay et al. (Reference Mukhopadhyay, van Oijen and de Goey2015) conducted both a priori and a posteriori studies of turbulent premixed Bunsen flames using both the static and the dynamic gradient model (using 
$C^M_d$), and it was concluded that the dynamic model did not improve the predictions over its static counterpart. Therefore, the overall evidence appears to suggest that the improvements obtained by a dynamic procedure might be limited by the filter size.
5.2. A reconstruction-based interpretation of the dynamic procedure
The role of a dynamic procedure is to correct a base static model, call this 
$B_m$. However, in LES the unresolved variance is unknown, and the only information available to correct 
$B_m$ are the LES variables 
$\bar {\rho }$, 
$\overline {\rho \phi }$, etc. Nevertheless, from the resolved variables one can calculate a corresponding resolved variance 
$\sigma ^2_r$ – an equivalent resolved stress tensor was used for developing a dynamic model for the unresolved Reynolds stresses by Liu, Meneveau & Katz (Reference Liu, Meneveau and Katz1994). If the resolved variance is modelled using a model 
$T_m$, a natural approach is to try to correct the base model of the unresolved variance using 
$\sigma ^2_r$ and 
$T_m$. In fact, the Leonard term in the expression for 
$C^B_d$ in (5.13) (numerator) divided by 
$\hat {\bar {\rho }}$ is the resolved variance 
$\sigma ^2_r$ at the test-filter level,
\begin{equation} \sigma^2_r=\frac{{\widehat{\bar{\rho}\tilde{\phi}\tilde{\phi}}}}{\hat{\bar{\rho}}}-\left( \frac{\widehat{\bar{\rho}\tilde{\phi}} }{ \hat{\bar{\rho}} }\right) ^2=\left( \tilde{\phi}\tilde{\phi} \right)^{\check{}}-\check{\tilde{\phi}}\check{\tilde{\phi}}=\frac{L}{\hat{\bar{\rho}}}. \end{equation}
Also, the remaining denominator term in (5.13) 
$\hat {\varDelta }^2|\boldsymbol {\nabla } \check {\tilde {\phi }}|^2$ is in fact a gradient model 
$T_m$ for the resolved variance. Therefore, the expression for 
$C^B_d$ in (5.13) can also be written as 
$\sigma ^2_r/T_m$, which is the ratio between the resolved and the modelled resolved variance at the test-filter level. This ratio can be interpreted as a measure of the error in modelling the resolved variance 
$\sigma ^2_r$ using 
$T_m$. We may therefore express the scale-similarity assumption as
\begin{equation} \frac{\sigma^2}{B_m} \simeq \frac{\sigma^2_r}{T_m} \simeq C^N_d, \end{equation}
essentially stating that the error in approximating the unresolved variance 
$\sigma ^2$ using a base model 
$B_m$ scales as the error in modelling the resolved variance 
$\sigma ^2_r$ using a base model 
$T_m$. Regularisation (
$\langle \ \rangle$) of the dynamic parameter 
$C^N_d$ then leads to the following expression of the dynamic model:
\begin{equation} \sigma^2=B_m\frac{\langle \sigma^2_R T_m \rangle}{\langle T_mT_m \rangle} =B_mC^N_d. \end{equation}
As shown in § 3, the gradient model is a subset model of more general reconstruction base models 
$B_m$ such as SM4 and AD4. Therefore, if 
$B_m$ is a general reconstruction-based model, a natural approach to obtain a model 
$T_m$ for the resolved variance is to use reconstruction as well but at the test-filter level. Denoting 
$\rho _L=\bar {\rho }$ and 
$\phi _L=\tilde {\phi }$ we can write the resolved variance as
\begin{equation} \sigma^2_r=\frac{{\widehat{\rho_L \phi_L \phi_L }}}{\hat{\rho}_L}-\left( \frac{\widehat{\rho_L\phi_L} }{ \hat{\rho}_L }\right) ^2=\left( \phi_L \phi_L \right)^{\check{}}-\check{\phi}_L \check{\phi}_L. \end{equation}
Clearly, the above equation is analogous to the definition of the unresolved variance, the only difference being the filter size and the definition of the Favre averaging. We are now in a position to write the reconstruction model at the test-filter level for 
$\sigma ^2_r$ as
\begin{equation} T_m=\frac{{\widehat{{\rho^*_L} {\phi^*_L} {\phi^*_L} }}}{\widehat{{\rho^*_L}}}-\left( \frac{\widehat{{\rho^*_L}{\phi^*_L}} }{ \widehat{{\rho^*_L}}} \right) ^2=\left( \phi^*_L \phi^*_L \right)^{\check{}}-\left(\check{\phi^*_L} \right)^2, \end{equation}
where 
$\rho ^*_L$ and 
$\phi ^*_L=(\rho \phi )^*_L/ \rho ^*_L$ are the reconstruction estimates at the test-filter level. Specifically, denoting the reconstruction operator at the test-filter level as 
$R_{\hat {\varDelta }}$, the reconstructions at the test-filter level are given by 
$\rho ^*_L=R_{\hat {\varDelta }}[\hat {\bar {\rho }}]$ and 
$(\rho \phi )^*_L=R_{\hat {\varDelta }}[\widehat {\overline{\rho u}}]$. Provided 
${\phi _L^*} \in [0,1]$ then 
$\check {\phi _L^*} \in [0,1]$ as well, and 
$T_m \in [0, \check {\phi _L^*}(1- \check {\phi _L^*})]$ ensuring that the modelled resolved variance is bounded. Also, since 
$\sigma ^2_r \in [0, \check {\phi }_L(1-\check {\phi }_L)]$, 
$C^N_d \ge 0$ and the variance predictions are always positive. Furthermore, in the limit 
$R_{\varDelta }[\bar {\phi }] \rightarrow \phi$ and 
$R_{\hat {\varDelta }}[\hat {\bar {\phi }}] \rightarrow \bar {\phi }$, 
$B_m \rightarrow \sigma ^2$ while 
$T_m \rightarrow \sigma ^2_r$, therefore, 
$C^N_d\rightarrow 1$, and the dynamic model recovers the actual variance.
The reconstruction operator 
$R_{\hat {\varDelta }}$ acting on the test-filtered variables can be obtained using the same procedure as in § 3 but at the test-filter level. For example, using truncated Taylor-series expansions we have, for a general variable 
$\hat {u}_L=\hat {\bar {u}}$, the reconstruction estimate
\begin{equation} R_{\hat{\varDelta}}[\hat{\bar{u}}]=u^*_L=\hat{u}_L-a_2'\nabla^2 \hat{u}_L+O(\hat{\varDelta}^4)=\hat{\bar{u}}-a_2'\nabla^2 \hat{\bar{u}}+O(\hat{\varDelta}^4), \end{equation}
where 
$a'_2$ is the filter's second moment but using 
$\hat {\varDelta }$ instead of 
$\varDelta$ (see table 1). As an example, for an 
$O(\hat {\varDelta }^2)$ reconstruction at the test-filter level: 
$\rho ^*_L=\hat {\rho }_L=\hat {\bar {\rho }}$, 
$(\rho \phi )^*_L=\widehat {\rho _L \phi _L}=\widehat {\bar {\rho } \tilde {\phi }}$ and 
$\phi ^*_L=\widehat {\bar {\rho }\tilde {\phi }}/\hat {\bar {\rho }}=\check {\tilde {\phi }}$. Inserting these reconstructions in (5.18) we obtain an equivalent SM2 model for the resolved variance,
\begin{equation} T_m=\frac{\widehat{ \hat{\bar{\rho}} \check{\tilde{\phi}}\check{\tilde{\phi}} }}{\hat{\hat{\bar{\rho}}}}-\left(\frac{ \widehat{\hat{\bar{\rho}}\check{\tilde{\phi}}}}{\hat{\bar{\rho}} }\right)^2=\left( \check{\tilde{\phi}}\check{\tilde{\phi}} \right)^{{\unicode{x0311}}{}}-{\unicode{x0311}}{\check{\tilde{\phi}}} {\unicode{x0311}}{\check{\tilde{\phi}}}, \end{equation}and the dynamic parameter for the SM2 model (with regularisation) becomes
\begin{equation} C^N_d=\frac{\left\langle \left( \left( \tilde{\phi}\tilde{\phi} \right)^{\check{}}-\check{\tilde{\phi}}\check{\tilde{\phi}} \right) \left( \left( \check{\tilde{\phi}}\check{\tilde{\phi}} \right)^{{\unicode{x0311}}{}}-{\unicode{x0311}}{\check{\tilde{\phi}}}{\unicode{x0311}}{\check{\tilde{\phi}}} \right) \right\rangle }{\left\langle \left( \left( \check{\tilde{\phi}}\check{\tilde{\phi}} \right)^{{\unicode{x0311}}{}}-{\unicode{x0311}}{\check{\tilde{\phi}}}{\unicode{x0311}}{\check{\tilde{\phi}}} \right) \left( \left( \check{\tilde{\phi}}\check{\tilde{\phi}} \right)^{{\unicode{x0311}}{}}-{\unicode{x0311}}{\check{\tilde{\phi}}}{\unicode{x0311}}{\check{\tilde{\phi}}} \right) \right\rangle}. \end{equation}The overall dynamic SM2 model, call this DSM2-N, is just the base static SM2 model of (3.11) corrected using the above expression for the dynamic parameter,
\begin{equation} \sigma^2=C^N_d \left( \left( \tilde{\phi}\tilde{\phi} \right) ^{\breve{ }} - \breve{\tilde{\phi}}\breve{\tilde{\phi}} \right)\!. \end{equation}
For an 
$O(\hat {\varDelta }^4)$ reconstruction on the other hand: 
$\rho ^*_L=\hat {\rho }_L-a'_2\nabla ^2\hat {\rho }_L+O(\hat {\varDelta }^4)=\,$ 
$\hat {\bar {\rho }}-a'_2\nabla ^2\hat {\bar {\rho }}+O(\hat {\varDelta }^4)$ and 
$(\rho \phi )^*_L=\widehat {({\rho \phi })_L}-a'_2\nabla ^2\widehat {({\rho \phi })_L}+O(\hat {\varDelta }^4)=\widehat {\overline {\rho \phi }}-a'_2\nabla ^2(\widehat {\overline {\rho \phi }}) + O(\hat {\varDelta }^4)$. As before, 
$\phi ^*_L=(\rho \phi )^*_L/\rho ^*_L$. Inserting these reconstructions in (5.18) and following the same procedure as in § 3, we eventually obtain an equivalent SM4 model for the resolved variance,
\begin{align} T_m=\left( \check{\tilde{\phi}}\check{\tilde{\phi}} \right)^{{\unicode{x0311}}{}}-{\unicode{x0311}}{\check{\tilde{\phi}}}{\unicode{x0311}}{\check{\tilde{\phi}}}+\frac{2a'_2}{\hat{\hat{\bar{\rho}}}}\left( {\unicode{x0311}}{\check{\tilde{\phi}}} \nabla ^2 \widehat{\hat{\bar{\rho}}\check{\tilde{\phi}}}-\widehat{\check{\tilde{\phi}}\nabla^2\hat{\bar{\rho}}\check{\tilde{\phi}}} \right) +\frac{a'_2}{\hat{\hat{\bar{\rho}}}}\left( \widehat{\check{\tilde{\phi}}\check{\tilde{\phi}}\nabla^2\hat{\bar{\rho}}}+\left(\check{\tilde{\phi}}\check{\tilde{\phi}}\right)^{{\unicode{x0311}}{}} \nabla^2\hat{\hat{\bar{\rho}}}-2{\unicode{x0311}}{\check{\tilde{\phi}}}{\unicode{x0311}}{\check{\tilde{\phi}}} \nabla^2\hat{\hat{\bar{\rho}}}\right)\!. \end{align}
Expanding the test-filtering operations up to 
$O(\hat {\varDelta }^4)$ in the above SM4 model, as done in § 3, it is straightforward to show after some algebraic manipulation that the above expression reduces to
\begin{equation} T_m=2a_2'\frac{\partial{\check{\tilde{\phi}}}}{\partial{x_i}}\frac{\partial{\check{\tilde{\phi}}}}{\partial{x_i}}, \end{equation}
which is the equivalent gradient model for the resolved variance at the test-filter level. For a Gaussian filter, 
$a_2'=\hat {\varDelta }^2/24$ in which case we recover the gradient model for the resolved variance employed in the classic dynamic models. Therefore, the expression for the model term in (5.13) is a subset of the more general model in (5.23). As a result, if the expression in (5.23) is used to develop a dynamic model, we obtain (with regularisation)
\begin{align} C^N_d=\frac{\left\langle \left( \left( \tilde{\phi}\tilde{\phi} \right)^{\check{}}-\check{\tilde{\phi}}\check{\tilde{\phi}} \right) \left( \left( \check{\tilde{\phi}}\check{\tilde{\phi}} \right)^{{\unicode{x0311}}{}}-{\unicode{x0311}}{\check{\tilde{\phi}}}{\unicode{x0311}}{\check{\tilde{\phi}}}+\dfrac{2a'_2}{\hat{\hat{\bar{\rho}}}}\left( {\unicode{x0311}}{\check{\tilde{\phi}}} \nabla ^2 \widehat{\hat{\bar{\rho}}\check{\tilde{\phi}}}-\widehat{\check{\tilde{\phi}}\nabla^2\hat{\bar{\rho}}\check{\tilde{\phi}}} \right) +\dfrac{a'_2}{\hat{\hat{\bar{\rho}}}}\left( \widehat{\check{\tilde{\phi}}\check{\tilde{\phi}}\nabla^2\hat{\bar{\rho}}}+\left(\check{\tilde{\phi}}\check{\tilde{\phi}}\right)^{{\unicode{x0311}}{}} \nabla^2\hat{\hat{\bar{\rho}}}-2{\unicode{x0311}}{\check{\tilde{\phi}}}{\unicode{x0311}}{\check{\tilde{\phi}}} \nabla^2\hat{\hat{\bar{\rho}}}\right) \right) \right\rangle}{\left\langle \left( \left( \check{\tilde{\phi}}\check{\tilde{\phi}} \right)^{{\unicode{x0311}}{}}-{\unicode{x0311}}{\check{\tilde{\phi}}}{\unicode{x0311}}{\check{\tilde{\phi}}}+\dfrac{2a'_2}{\hat{\hat{\bar{\rho}}}}\left( {\unicode{x0311}}{\check{\tilde{\phi}}} \nabla ^2 \widehat{\hat{\bar{\rho}}\check{\tilde{\phi}}}-\widehat{\check{\tilde{\phi}}\nabla^2\hat{\bar{\rho}}\check{\tilde{\phi}}} \right) +\dfrac{a'_2}{\hat{\hat{\bar{\rho}}}}\left( \widehat{\check{\tilde{\phi}}\check{\tilde{\phi}}\nabla^2\hat{\bar{\rho}}}+\left(\check{\tilde{\phi}}\check{\tilde{\phi}}\right)^{{\unicode{x0311}}{}} \nabla^2\hat{\hat{\bar{\rho}}}-2{\unicode{x0311}}{\check{\tilde{\phi}}}{\unicode{x0311}}{\check{\tilde{\phi}}} \nabla^2\hat{\hat{\bar{\rho}}}\right) \right)^2 \right\rangle}. \end{align}
It is important to note that in contrast to the dynamic parameter for the SM2 model, the above expression is not necessarily always positive as a result of truncating the filtering operations to obtain the expression for 
$T_m$. The dynamic SM4 model DSM4-N then becomes
\begin{align} \sigma ^2 &=C^N_d\left(
\left( \tilde{\phi}\tilde{\phi} \right) ^{\breve{ }} -
\breve{\tilde{\phi}}\breve{\tilde{\phi}}+\frac{2a_2}{\bar{\bar{\rho}}}\left(
\breve{\tilde{\phi}} \nabla ^2
\overline{\bar{\rho}\tilde{\phi}}-\overline{\tilde{\phi}\nabla^2\bar{\rho}\tilde{\phi}}
\right)\right.\notag\\ &\quad + \left. \frac{a_2}{\bar{\bar{\rho}}}\left(
\overline{\tilde{\phi}\tilde{\phi}\nabla^2\bar{\rho}}+\left(\tilde{\phi}\tilde{\phi}\right)^{\breve{}}
\nabla^2\bar{\bar{\rho}}-2\breve{\tilde{\phi}}\breve{\tilde{\phi}}
\nabla^2\bar{\bar{\rho}}\right) \right)\!.
\end{align}
As shown above, the dynamic procedure employed in classic algebraic dynamic models is essentially a two-level approximate reconstruction-based modelling approach. At the normal filter level 
$\varDelta$ reconstruction serves to obtain a base model 
$B_m$ for the unresolved variance 
$\sigma ^2$ while at the test-filter level 
$\hat {\varDelta }$ it serves to obtain a model 
$T_m$ for the resolved variance 
$\sigma ^2_r$. Then, by employing the scale-similarity assumption through (5.16) the ratio of the resolved variance to the modelled resolved variance is used to correct the variance prediction of the base static model.
5.3. Generalized dynamic reconstruction models
In this section we consolidate all the results thus far in an effort to present a general methodology for obtaining a dynamic model for any given general reconstruction operator 
$R_{\varDelta }$ acting on the LES fields, and for any reconstruction operator 
$R_{\hat {\varDelta }}$ acting on the test-filtered LES fields. This also includes bounding of the reconstructed fields. The overall model reads
\begin{equation} \sigma^2=C^N_d \left( \frac{\overline{\rho^* \phi ^* \phi ^*}}{ \overline{\rho ^*}}- \left(\frac{ \overline{\rho^* \phi^*} }{ \overline{\rho^*}} \right)^2 \right)\!. \end{equation}
The bounded reconstructions at 
$\varDelta$ are obtained using
\begin{equation} \left. \begin{gathered} \rho^* =\max(\min(R_{\varDelta}[\bar{\rho}],\rho_h),\rho_l), \\ (\rho \phi)^* =\max(\min(R_{\varDelta}[\overline{\rho\phi}],(\rho\phi)_h),(\rho\phi)_l),\\ \phi^*=\max(\min((\rho\phi)^*/\rho^*,1),0), \end{gathered} \right\} \end{equation}
and at 
$\hat {\varDelta }$ by
\begin{equation} \left. \begin{gathered} \rho^*_L =\max(\min(R_{\hat{\varDelta}}[\hat{\bar{\rho}}],\rho_h),\rho_l), \\ (\rho \phi)^*_L =\max(\min(R_{\hat{\varDelta}}[\widehat{\overline{\rho\phi}}],(\rho \phi)_h),(\rho\phi)_l),\\ \phi^*_L=\max(\min((\rho\phi)^*_L/\rho^*_L,1),0). \end{gathered} \right\} \end{equation}
The dynamic parameter 
$C^N_d$ is then calculated using
\begin{equation} C^N_d=\frac{\left\langle \left( \left( \tilde{\phi}\tilde{\phi} \right)^{\check{}}-\check{\tilde{\phi}}\check{\tilde{\phi}} \right) \left( \left( \phi^*_L \phi^*_L \right)^{\check{}}-\check{\phi^*_L} \check{\phi^*_L} \right) \right\rangle}{\left\langle \left( \left( \phi^*_L \phi^*_L \right)^{\check{}}-\check{\phi^*_L} \check{\phi^*_L} \right) \left( \left( \phi^*_L \phi^*_L \right)^{\check{}}-\check{\phi^*_L} \check{\phi^*_L} \right) \right\rangle}. \end{equation}
Clearly, the above expression is always positive ensuring that we recover realisable predictions of the variance. As we will see in the sections that follow, in the case where 
$R_{\varDelta }=(I-a_2\nabla ^2)$ and 
$R_{\hat {\varDelta }}=(I-a'_2\nabla ^2)$ the above procedure results in a dynamic model far superior than the DSM4-N model developed in the previous section. We will refer to this model as DAD4-N. Furthermore, in the case where the reconstruction operator is given by the optimised inverse discrete filter we will refer to this model as DEIF-N. Note that in this case optimised inverse filters are used both at the normal filter level 
$\varDelta$, and at the test filter 
$\hat {\varDelta }$. At the normal filter level the inverse filter coefficients are listed in table 3 (
$\gamma =4$), and at the test-filter level the inverse filter coefficients are listed in table 4 (
$\gamma =8$).
6. Results
A summary of all the different models tested is given in table 5. These include five static models (SM2, SM4, GR, AD4, DEIF) and their dynamic versions (DSM2-N, DSM4-N, DGR-B, DAD4-N, DEIF-N). Here ‘B’ in the dynamic models indicates that the dynamic parameter is calculated using 
$C^B_d$ from (5.13) (the dynamic gradient model with 
$C^M_d$ is not considered due to its poorer performance in comparison), while 
$N$ indicates using 
$C^N_d$. Note that the 
$C^N_d$ formulation depends on the base static model used and the reconstruction operator used as outlined in § 5.2. The performance of each model is quantified by calculating the mean-squared error between the modelled variance 
$\sigma ^2_m$ and the actual reference variance 
$\sigma ^2$ as obtained on the LES mesh. This is defined using
\begin{equation} e=\frac{1}{N_s}\sum_{\boldsymbol{x},t} {(\sigma^2_m(\boldsymbol{x},t)-\sigma^2(\boldsymbol{x},t))}^2 |_{\tilde{c} \in [\epsilon,1-\epsilon]}, \end{equation}
where 
$N_s=N_xN_yN_zN_t$ is the total sample size and 
$\epsilon =0.05$. The conditioning on 
$\tilde {c} \in [\epsilon,1-\epsilon ]$ ensures that no-reaction regions are not sampled thereby biasing the results. This is because in no-reaction regions the scalar fields are more or less constant and filtering/reconstruction would be nearly perfect. The performance of each model is also evaluated by calculating the conditional variance based on 
$\tilde {c}$ denoted by 
$\sigma ^2|\tilde {c}$. The conditional averages are also time averaged 
$\langle \ \rangle$ over several samples in time to increase the statistical accuracy of the results. In the following section a brief description of the DNS databases used for testing is presented.
Table 5. Summary of the various models tested and the relevant equations used. A Gaussian filter is used for all models for which 
$a_2=\varDelta ^2 /24$. Note that the parameter 
$C^N_d$ is different for each model and calculated using the procedure outlined in § 5.2.
6.1. The DNS databases
The DNS databases correspond to two turbulent premixed flames simulated using detailed chemistry. Database PB corresponds to a statistically planar multi-component-fuel/air flame propagating in an inflow–outflow configuration, and database V97 is a statistically stationary hydrogen/air V-flame anchored on a fixed rod. In both cases, turbulent reactants enter from one end of the computational domain forming a flame that interacts with the incoming turbulence. Figure 2 shows the flow configuration for each flame. Flame PB was simulated using the in-house fully compressible solver SENGA2 (Cant Reference Cant2012) and flame V97 using the TTX solver (Minamoto et al. Reference Minamoto, Fukushima, Tanahashi, Miyauchi and Dunstan2011). The solver SENGA2 employs a tenth-order centred finite-difference scheme for the spatial derivatives and a fourth-order Runge–Kutta scheme for the time stepping, while TTX employs a fourth-order finite-difference scheme and a third-order Runge–Kutta scheme. Boundary conditions in both codes are implemented using the Navier–Stokes characteristics boundary conditions formalism (Thompson Reference Thompson1987; Poinsot & Lele Reference Poinsot and Lele1992; Sutherland & Kennedy Reference Sutherland and Kennedy2003).
Figure 2. Progress-variable isosurface of the leading edge of the flame (
$c=0.1$) for (a) the planar flame and (b) the V-flame. Dimensions are in mm.
Chemistry is modelled using a 49 reaction 15 species chemical mechanism for the planar flame (Nikolaou & Swaminathan Reference Nikolaou and Swaminathan2013), and a 27 reaction 12 species chemical mechanism for the hydrogen/air flame (Gutheil, Balakrishnan & Williams Reference Gutheil, Balakrishnan and Williams1993). Table 6 lists the important parameters for each database: 
$u_{rms}$ is the root-mean-square value of the fluctuating component of the incoming turbulent velocity field, 
$l_{T}$ is the turbulence integral length scale in the reactant side, the turbulence Reynolds number is 
$Re_{{T}} =u_{rms}l_{T}/ \nu _r$, the Damkohler number is 
$Da=(l_{T}/u_{rms})/(\delta / s_{L})$ and the Karlovitz number is 
$Ka=(\delta / \eta _k)^2$, where 
$s_{{L}}$ is the laminar flame speed. The (diffusive) thickness is 
$\delta =\nu _r/s_{{L}}$. The laminar flame thickness is defined as 
$\delta _{{L}}=1/\max (\textrm {d}c/{\textrm {d} x})$, where 
$c$ is the progress variable. The progress variable is defined using species mass fractions, i.e. 
$c=(Y-Y_r)/(Y_p-Y_r)$. For the planar flame, 
$c$ is based on the CO mass fraction, and for the V-flame it is based on the 
$\mathrm {H}_2$ mass fraction. Further details of the simulations, numerical schemes, etc. can be found in Nikolaou & Swaminathan (Reference Nikolaou and Swaminathan2015) and Minamoto et al. (Reference Minamoto, Fukushima, Tanahashi, Miyauchi and Dunstan2011). As noted in table 6, the databases differ substantially in their thermo-chemical properties (laminar flame speed, flame thickness, etc.) and flow configuration (freely propagating and with mean shear). Note that even though the 
$Da$ number for the planar flame is lower, the V-flame experiences smaller-scale wrinkling due to the presence of the rod that ‘re-energizes’ the flow. In the planar flame, turbulence decays with downstream distance 
$x$ and, as a result, the flame experiences larger-scale wrinkling as evidenced in figure 2(a). This has some important modelling consequences since reconstruction for the V-flame must recover a wider range of wavenumbers for an accurate variance estimation. Therefore, for the same reconstruction order, we would expect the variance predictions for the V-flame to be of lower accuracy.
Table 6. Important flame parameters for the planar (PB) and V (V97) flames.
The DNS data are filtered using a Gaussian filter: 
$G(\boldsymbol {x})=( 6 / ( {\rm \pi}{\varDelta } ) ^{2} )^{3/2}$ 
$\exp ( -6\boldsymbol {x} {\cdot } \boldsymbol {x}/ \varDelta ^{2} )$. The analytic form of the filter is approximated in physical space using a 3-D discrete filtering operator that is constructed by the successive application of 1-D discrete filters in each coordinate direction. This implies that the discrete filter operator does not preserve the isotropy of the analytic filter. However, in practice, the filter width is sufficiently larger than the mesh spacing that ensures that any anisotropic effects are negligible. On the boundaries, the discrete filter is applied as usual taking into account the corresponding boundary conditions for each case. In order to simulate an LES, the filtered DNS data are then sampled from the DNS mesh onto a coarser LES mesh. The choice of the LES mesh size 
$h$ is based on the criterion derived by Nikolaou & Vervisch (Reference Nikolaou and Vervisch2018), specifically 
$h= \varDelta /4$ is used. This criterion ensures that the filtered flame thickness is adequately resolved on the LES mesh. Table 7 shows for each filter size the corresponding LES mesh. Note that 
$\varDelta ^+=\varDelta / \delta _L$, and the range of 
$\varDelta ^+$ values considered corresponds to well-resolved LES conditions – all three LES meshes are substantially coarser than the DNS mesh.
Table 7. The DNS and LES meshes for 
$\varDelta /h =4$.
6.2. Static models
In order to test the validity of the modelling framework developed in § 3 we begin by examining the performance of the base (static) models. Table 8 lists the errors for each model, and figure 3 shows the corresponding conditional averages. Instantaneous predictions for each model are shown for the V-flame in figure 4 (similar scatter plots were obtained for case PB). The results in table 8 indicate that the SM4 model substantially improves the predictions over the SM2 model for all three filter widths and for both databases. The SM4 model also performs better than the GR model that, as shown in § 3, is a subset of SM4. This behaviour is also evident in the instantaneous predictions shown in figure 4: the GR model has an inferior correlation with the reference variance, and a similar instantaneous distribution for this model was also obtained by Knudsen et al. (Reference Knudsen, Kim and Pitsch2010). Even though the SM4 model has some negative predictions these are very few, which contributes to the lower overall error for this model. The AD4 model, which is of the same reconstruction order (
$O(\varDelta ^4)$) as the SM4 model, has substantially improved performance over the SM4 model for both databases and all three filter widths. This implies that post-reconstruction bounding has a substantial effect on the model performance; however, part of the improvement also originates from the direct use of the approximately reconstructed fields. This is in contrast to the SM4 model where the approximate reconstructions are used to develop an algbebraic model that entails further truncations of the filtering operations. Furthermore, we observe from figure 3 that the predictions of the AD4 model improve consistently over the entire range of 
$\tilde {c}$ values. Overall, the DEIF model has the best performance as evidenced by the lowest error, and the corresponding substantially improved conditional averages – figure 4(e) shows that this is a direct result of the much improved correlation with the reference variance. Between the two databases, DEIF performs better for the planar flame that exhibits larger-scale wrinkling. This is not surprising since the inverse filter was designed to correspond to 
$N=5$ iterations of the van Cittert algorithm that is sufficient to recover the largest scales (smaller wavenumbers) for this flame. The V-flame on the other hand contains contributions from smaller scales and higher-order reconstruction is required to recover these, which leads to a slight under-prediction of the peak variance. Despite this, the predictions of the DEIF model are still substantially improved in comparison to the rest of the models.
Figure 3. Conditionally averaged variance predictions of the static models against the reference variance at 
$\varDelta ^+=3$ for (a) the planar flame and (b) the V-flame. Note that these are normalised using the maximum of the reference conditional variance. Results are shown for (a) 
$\textrm {PB} - \varDelta ^+=3$, (b) 
$\textrm {V}97 - \varDelta ^+=3$.
Figure 4. Instantaneous variance predictions of the static models against the reference variance for case V97 and 
$\varDelta ^+=3$. The predictions are normalised using the maximum reference variance. Results are shown for (a) SM2, (b) GR, (c) SM4, (d) AD4, (e) DEIF.
Table 8. Mean-squared errors obtained using (6.1) for each of the static models.
Overall, these results combined serve to confirm the validity of the reconstruction-based modelling framework proposed in § 3, and one could employ higher-order reconstructions leading to further improvements. In the following section we examine whether the corresponding dynamic formulations of each of the above static models improve their predictions and to what extent.
6.3. Comparison of $C_d$ formulations
Before quantifying the performance of the dynamic models, it is imperative to examine the distribution of each model's dynamic parameter. The p.d.f.s for each one are shown in figure 5 for case V97 (similar distributions were obtained for the planar flame). In accordance with the results obtained for a passive scalar by Balarac et al. (Reference Balarac, Pitsch and Raman2008) but also for reacting scalars by Knudsen et al. (Reference Knudsen, Kim and Pitsch2010), 
$C^M_d$ can become negative as evidenced by the finite non-zero probability for negative 
$C^M_d$ values in figure 5(a). This leads to negative variance predictions that are not realisable. Perhaps more important however is that for increasing filter widths 
$\varDelta$, the probability of 
$C^M_d$ becoming negative increases whereas the probability of becoming positive decreases. This trend, which is also in accordance with the results obtained by Balarac et al. (Reference Balarac, Pitsch and Raman2008), implies that the predictions for this model deteriorate for increasing 
$\varDelta$. For these reasons, the variance predictions using this dynamic model were found to be inferior to the rest of the models. The formulation in (5.13) (
$C^B_d$) on the other hand is always positive as expected, and more robust as it is less affected by variations in 
$\varDelta$. In the case of the DSM2-N model, the dynamic parameter is always positive by construction that is also reflected in the corresponding p.d.f. The distribution of 
$C^N_d$ for the DSM4-N model on the other hand has some negative contributions owing to truncating the filtering operations to obtain the expression for 
$T_m$ in (5.25). In the case of the DAD4-N model, 
$C^N_d$ is also always positive by construction. It is also important to note that all four 
$C^N_d$ formulations are much less sensitive to variations in filter width in comparison to the formulation for 
$C^M_d$. Also, since the base static models AD4 and DEIF already provide adequately good predictions of the variance, the corresponding p.d.f.s of their dynamic parameters approach a delta function centred around 1 in accordance with the limiting behaviour of the dynamic reconstruction models presented in § 5.
Figure 5. Probability density functions of the different dynamic parameters: (a) 
$C^M_d$, (5.8); (b) 
$C^B_d$, (5.13); (c) 
$C^N_d$ for the SM2 model, (5.21); (d) 
$C^N_d$ for the SM4 model, (5.25); (e) 
$C^N_d$ for the AD4 model, (5.30) and (f) 
$C^N_d$ for the DEIF model, (5.30). The p.d.f.s are for case V97.
6.4. Dynamic models
Table 9 shows the errors for the five dynamic models. The conditional averages are shown in figure 6, and instantaneous predictions are shown in figure 7 for case V97 at 
$\varDelta ^+=3$. It is important to note at this point that all dynamic models employ a test filter 
$\hat {\varDelta }=2\varDelta$, and regularisation of the dynamic parameter is conducted by averaging in the homogeneous directions in accordance with standard practice in the literature. Firstly, by comparing tables 8 and 9 we see that the dynamic procedure improves the predictions of the static models. This results in improved predictions of the peak variance as evidenced from the instantaneous predictions shown in figure 7 despite the generally more distributed predictions observed for the dynamic models. This increased scatter is the result of local variations in the dynamic parameter particularly in regions of intense heat release. In such regions, the standard averaging procedure in homogeneous directions may not be sufficient to smooth out any large variations. This issue has also been identified in past studies on turbulent premixed flames modelling, for instance, when calculating the dynamic Smagorinsky parameter (Klein, Kasten & Chakraborty Reference Klein, Kasten and Chakraborty2015). An alternative averaging procedure that involves conditioning on 
$\tilde {c}$ may help to reduce the scatter as suggested by Klein et al. (Reference Klein, Kasten and Chakraborty2015). Nevertheless, the conditional averages are also improved for all 
$\tilde {c}$ values as we observe from figure 6. The dynamic procedure proposed in § 5 even improves the results of the static SM2 model that is the lowest-order reconstruction model (
$O(\varDelta ^2)$), and which has the largest error of all the static models. In fact, for the largest filter width, the DSM2-N model compares favourably with the dynamic gradient model that is derived from an 
$O(\varDelta ^4)$ reconstruction approximation. These results serve to validate the dynamic modelling procedure that was reinterpreted as a two-level reconstruction-modelling approach in § 5. The DSM4-N model improves the results considerably over the DSM2-N model as a result of the higher-order Taylor expansion for the reconstructed fields. Some negative predictions do occur for reasons discussed in the previous section, however, the overall error is lower than the dynamic gradient model. Furthermore, from the instantaneous predictions in figure 7(d) we see that the overall correlation of the DSM4-N model is improved in comparison to the dynamic gradient model. The DAD4-N model on the other hand has the lowest error, highest correlation and best overall performance of all the algebraic models. The dynamic inverse filter model DEIF-N has the best overall performance and improves the predictions particularly in the leading edge (
$\tilde {c}<0.4$) of the flame surface as evidenced in figure 6(b). It is interesting to note that, for the V-flame, all models show similar responses when comparing their conditional averages, however, they all have substantially different errors and correlations with respect to the reference variance. This serves to illustrate that comparing conditional averages alone is insufficient for model comparison/evaluation.
Figure 6. Conditionally averaged variance predictions of the dynamic models against the reference variance at 
$\varDelta ^+=3$ for (a) the planar flame and (b) the V-flame. Note that these are normalised using the maximum of the reference conditional variance. Results are shown for (a) 
$\textrm {PB} - \varDelta ^+=3$ and (b) 
$\textrm {V}97 - \varDelta ^+=3$.
Figure 7. Instantaneous variance predictions of the dynamic models against the reference variance for case V97 and 
$\varDelta ^+=3$. The predictions are normalised using the maximum reference variance. Results are shown for (a) DSM2-N, (b) DGR-B, (c) SM4, (d) AD4, (e) DEIF.
Table 9. Mean-squared errors obtained using (6.1) for each of the dynamic models.
7. Conclusions
A generalized modelling framework for the scalar variance in LES is presented. This new modelling framework is based on reconstruction of the primary variables with the degree of reconstruction determining the accuracy of the model. It is shown that the classic scale-similarity model of Cook & Riley (Reference Cook and Riley1994) but also the gradient model by Pierce & Moin (Reference Pierce and Moin1998) are subsets of more general models derived using the proposed framework. Another novelty of the proposed framework is that the base static models are physically realisable by construction since the variance predictions are bounded, thus satisfying the relevant properties of the scalar variance.
The methodology for the development of dynamic models stems from a new interpretation of the classic scale-similarity assumption as a two-level reconstruction process whose aim is to correct the prediction of the base static model using information from the resolved scales. This involves a reconstruction step at the normal filter level 
$\varDelta$ that serves to obtain a base static model, and a reconstruction step at the test-filter level 
$\hat {\varDelta }$ that serves to obtain a model for the resolved variance. This naturally leads to an alternative interpretation of the model dynamic parameter 
$C_d$, which is intrinsic to the base model but also to the specific reconstruction operator used.
By employing approximate reconstruction operators obtained using truncated Taylor-series expansions, a range of static and dynamic algebraic models are developed. In addition, a novel and computationally efficient reconstruction operator based on the concept of optimised inverse explicit filters is used to develop a new high-order static and dynamic reconstruction model. All models are eventually tested a priori using high-fidelity DNS data of two turbulent premixed flames. The highest-order reconstruction models show a substantial improvement in comparison to the classic models serving to justify the proposed modelling framework. At the same time, their performance is consistent with variations in filter width despite the large differences in thermo-chemical properties and flow configuration between the two databases.
Even though approximate reconstruction operators were employed, the proposed modelling framework allows for any arbitrary higher-order approximate or iterative reconstruction operator to be used, therefore leading to higher-order variance models. In addition, the proposed dynamic procedure allows SGS effects to be incorporated in a base reconstruction model, thereby extending the capabilities of current reconstruction-based modelling methods. Furthermore, in contrast to classic modelling approaches, the modelling framework presented in this work is general. As such, it can be used to model a wide range of unresolved terms in the governing equations for LES, thereby providing a single, coherent and unified modelling framework for the simulation of turbulent flows.
Acknowledgement
Z.N. acknowledges Dr Y. Minamoto for providing the hydrogen V-flame DNS data.
Funding
Z.N. acknowledges funding received under both the ‘Hydrogen and Energy Decarbonisation’ program by Agence Nationale de la Recherche (ANR) and the REDAFLOW project under the European Union's Horizon 2020 research and innovation programme-Marie Skłodowska-Curie grant agreement no 101019855.
Declaration of interests
The authors report no conflict of interest.
Appendix A
Let 
$\phi (x,t) \in [0,1]$ and 
$\tilde {G}(s;x,t)=G(s)\rho (x-s,t) / \bar {\rho }(x,t) \geq 0$ for all 
$s,x,t \in \mathbb {R}$ with the consistency condition 
$\int _{-\infty }^{\infty }\tilde {G}(s;x,t)\,\textrm {d}s=1$. We have
\begin{equation} \tilde{\phi}(x,t)-1=\int_{-\infty}^{\infty}\tilde{G}(s;x,t)\left( \phi(x-s,t)-1 \right)\,{\rm d}s, \end{equation}
where 
$( \phi (x-s,t)-1 ) \in [-1,0]$ and 
$\tilde {G}(s;x,t)( \phi (x-s,t)-1 ) \le 0$. As a result, 
$\tilde {\phi }(x,t) \leq 1$ for all 
$x,t \in \mathbb {R}$. Also, 
$\tilde {G}(s;x,t) \phi (x-s,t) \ge 0$, hence, 
$\tilde {\phi }(x,t) \ge 0$. Therefore, overall if 
$\phi \in [0,1]$ then 
$\tilde {\phi } \in [0,1]$ as well.
For non-positive filters, 
$\tilde {\phi }$ may not be bounded in 
$[0,1]$ even if 
$\phi \in [0,1]$. For such filters, there always exists by definition at least one region in space and time such that 
$\tilde {G}(s;x,t) < 0$ for all 
$s \in [a,b]$ where we will take 
$b>a$ without loss of generality. In such a case we may define, for all 
$t \in \mathbb {R}$, the function
\begin{equation} \phi(x-s,t)= \begin{cases} v \in [0,1] \ \text{constant}, & \text{for all }(x-s) \in [x-b,x-a] ,\\ 0, & \text{otherwise}. \end{cases} \end{equation}For the above distribution, the filtering integral reduces to
\begin{equation} \tilde{\phi}(x,t)-1=\int_{a}^{b}\tilde{G}(s;x,t)\left(v-1\right)\,{\rm d}s \end{equation}
for all 
$(x-s) \in [x-b,x-a]$ and where 
$\tilde {G}(s;x,t)(v-1) \geq 0$. Therefore, 
$\tilde {\phi }(x,t) \geq 1$ for all 
$(x-s) \in [x-b,x-a]$. Similarly, the lower bound is not satisfied either since for all 
$(x-s) \in [x-b,x-a]$ 
$\tilde {G}(s;x,t)v \leq 0$. In short, we have shown that, for non-positive filters and a bounded scalar 
$\phi \in [0,1]$, there always exists a bounded distribution in 
$[0,1]$ such that 
$\tilde {\phi }(x,t)$ is not bounded in 
$[0,1]$.
