Heng Xiao
email: hengxiao@vt.edu

Paola Cinnella
email: paola.cinnella@ensam.eu

Kevin T Crofton

Quantification of model uncertainty in RANS simulations: A review

Keywords:

In computational fluid dynamics simulations of industrial flows, models based on the Reynolds-averaged Navier-Stokes (RANS) equations are expected to play an important role in decades to come. However, model uncertainties are still a major obstacle for the predictive capability of RANS simulations. This review examines both the parametric and structural uncertainties in turbulence models. We review recent literature on data-free (uncertainty propagation) and data-driven (statistical inference) approaches for quantifying and reducing model uncertainties in RANS simulations. Moreover, the fundamentals of uncertainty propagation and Bayesian inference are introduced in the context of RANS model uncertainty quantification. Finally, the literature on uncertainties in scale-resolving simulations is briefly reviewed with particular emphasis on large eddy simulations.

Introduction

Turbulence affects natural and engineered systems from sub-meter to planetary scales yet it is among the last unsolved problems in classical physics. Accurate predictions of turbulent flows are of vital importance for the design, analysis, and operation of many critical systems in aerospace engineering such as aircraft, spacecraft, and gas turbine engines. The dynamics of fluid flows is described by the where u i , p, x i and t are, respectively, the flow velocity, pressure, and spatial and temporal coordinates. Although simpler in form than the partial differential equations governing the above-mentioned problems, incompressible NS equations cover a very wide variety of flow configurations and bear the key difficulty that leads to the turbulence modeling dilemma, i.e., the nonlinear convective term in Equation (1b). Equation ( 1 )/ i j j , the NS equations admit chaotic solutions when the Reynolds number is beyond some flow-dependent critical value. As the Reynolds number increases, eventually the flow reaches a state of motion characterized by strong three-dimensional and unsteady chaotic fluctuations of the velocity and pressure fields, which is referred to as the turbulent regime.

Landscape of turbulence modeling

Turbulent flows are characterized by a wide range of spatial and temporal scales. Consequently, performing direct numerical simulations (DNS) by solving the NS equations and resolving all the turbulence scales are prohibitively expensive, particularly for high Reynolds number flows. Practically used turbulence modeling strategies range from DNS with the highest fidelity, where all physics of spatial and temporal scales are resolved and no modeling is involved, to Reynolds averaged Navier-Stokes (RANS) simulations with the lowest fidelity, where the entire range of turbulent flow scales is modeled. This model hierarchy is illustrated in Fig. 1, with the top represented by the most physics-resolving and computationally expensive approach (DNS) and the bottom by the most empirical and computationally affordable approach (RANS). Lower fidelity models toward the bottom of the hierarchy involve more flow-dependent, uncertain closures than the higherfidelity, scale-resolving approaches towards the top of the hierarchy. On the other hand, high-fidelity, scale-resolving models are more susceptible to influences from numerical uncertainties as well as initial and boundary conditions.

A compromise between DNS and RANS simulations at two ends of the spectrum is large eddy simulation (LES), in which only the larger, more energetic scales are resolved, while scales below a cutoff threshold are filtered out. The filtered Navier-Stokes equations contain a subgrid-scale (SGS) stress that is unclosed and needs to be modeled. The SGS stress term represents the interactions between the filtered and resolved scales, which result from the nonlinear convection term [START_REF] Sagaut | Large Eddy Simulations for Incompressible Flows: An Introduction[END_REF]. Large eddy simulations have significantly reduced computational costs compared to DNS for shear flows far removed from wall boundaries. Unfortunately, they remain prohibitively expensive for wall bounded flows at high Reynolds number due to the small yet energetic scales dominating the dynamics in the near-wall regions [START_REF] Spalart | Detached-eddy simulation[END_REF]. This challenge has led to the development of methods combining LES in free shear regions with RANS models or other simplified models (e.g., boundary layer equation or law of the wall) in the under-resolved near-wall regions. Such approaches include hybrid RANS/LES models [START_REF] Fröhlich | Hybrid LES/RANS methods for the simulation of turbulent flows[END_REF][START_REF] Chaouat | The state of the art of hybrid rans/les modeling for the simulation of turbulent flows[END_REF] and wall-modeled LES [START_REF] Cabot | Approximate wall boundary conditions in the large-eddy simulation of high Reynolds number flow[END_REF][START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF][START_REF] Kawai | Wall-modeling in large eddy simulation: length scales, grid resolution, and accuracy[END_REF][START_REF] Yang | Integral wall model for large eddy simulations of wall-bounded turbulent flows[END_REF], among others.

While scale-resolving simulations such as DNS, LES, and hybrid RANS/LES provide more insights of fluid flow physics, in many simulations of engineering turbulent flows such as those for aerodynamic design and optimization, the quantities of interest depend on the mean flow only, and the instantaneous flow fields are not of concern. In these cases it is desirable to solve for the mean flow more efficiently. For that purpose, the instantaneous velocity u i and pressure p are decomposed into the sum of the mean2 components U i and P and the fluctuations u i and p , respectively. Substituting the decomposition into the Navier-- Stokes equations and taking the ensemble-average leads to the RANS equations:

U x 0 i i = (2a) U t U U x P x Re U x x u u x ( ) 1 . i i j j i i j j i j j 2 + = + (2b)
The RANS equations are similar in form to the Navier-Stokes equations except for the term involving the tensor u u i j . As with the SGS stress term in the filtered NS equations for LES, this term stems from the nonlinear convection term in the NS equation and represents the crosscomponent covariance among the velocity fluctuations. It is often referred to as Reynolds stress due to its formal similarity to the viscous stresses and is denoted as u u .

ij i j = ( 3
)
Since the velocity fluctuations are not available in RANS simulations, one must resort to closure models to supply Reynolds stresses, which lies at the root of most efforts of turbulence modeling. The choice of the appropriate modeling level remains a matter of expert judgment. In particular, it inevitably involves a compromise between computational cost and predictive accuracy. Even after a given fidelity level is selected (e.g., RANS or LES), several possible closure models may be designed for relating the unclosed terms to the resolved variables. These closure models differ both by their mathematical structure and by the associated model parameters. The common practice in turbulence modeling is to leave the choice of a specific closure model to user judgment and to treat model parameters as adjustable coefficients that are generally calibrated to reproduce simple, canonical flows. Both of the preceding aspects, however, represent sources of uncertainty in the prediction of new flows. Recent development of turbulence modeling in RANS, LES, and hybrid approaches has been reviewed by Durbin [START_REF] Durbin | Some recent developments in turbulence closure modeling[END_REF]. Despite considerable progress recently made in LES and hybrid RANS/LES models (e.g., [START_REF] Sagaut | Large Eddy Simulations for Incompressible Flows: An Introduction[END_REF][START_REF] Spalart | Detached-eddy simulation[END_REF][START_REF] Fröhlich | Hybrid LES/RANS methods for the simulation of turbulent flows[END_REF][START_REF] Girimaji | Partially-averaged Navier-Stokes model for turbulence: a Reynoldsaveraged Navier-Stokes to direct numerical simulation bridging method[END_REF][START_REF] Xiao | A consistent dual-mesh framework for hybrid LES/RANS modeling[END_REF]), RANS models are expected to remain the workhorse in engineering practice for decades to come, due to their much lower computational costs and superior robustness. For this reason, our review mainly focuses on the quantification and reduction of uncertainties in RANS models.

The landscape of RANS-based turbulence modeling has not changed for decades. The stagnation is evident from two observations as illustrated in Fig. 2. First, the number of wind tunnel tests performed in a typical design cycle of a commercial airplane was reduced from 75 in the 1970s to 10 in the 1990s, but this number has been stagnant since then, with turbulence models being the major bottleneck in predictive accuracy [START_REF] Johnson | Thirty years of development and application of CFD at boeing commercial airplanes, seattle[END_REF]. Second, most of the currently used turbulence models were developed decades ago and provide unsatisfactory performance for many flows. Generations of researchers have labored for many decades on dozens of turbulence models, yet none of them achieved predictive generality. Flow-specific tuning and fudge functions are still an indispensable part of RANS simulations [START_REF] Spalart | Philosophies and fallacies in turbulence modeling[END_REF]. Current development of improved turbulence models faces the dilemma of conserving the low computational costs and high robustness of RANS approaches while incorporating as much physics as possible.

Origin of uncertainties in RANS models

A recent review on data-driven turbulence modeling strategies [START_REF] Duraisamy | Turbulence modeling in the age of data[END_REF] classified the model uncertainties in RANS simulations into four levels, including uncertainties due to information loss in the Reynolds-averaging process, uncertainties in representing the Reynolds stress as a functional form of the mean field, uncertainties in the choice of the specific function, and uncertainties in the parameters of a given model. In this review, we will focus on the uncertainties due to the choice of functional forms and parameters in the turbulence models. Fig. 3 shows a graphical representation of different sources of model uncertainties in typical RANS models.

The following observations about the Reynolds stress tensor have profound implications for turbulence modeling and RANS model uncertainty quantification. First, it is a covariance tensor of velocity fluctuations as pointed out above, and mathematically any covariance tensor must be symmetric positive semi-definite. This is referred to as realizability requirement. Second, it appears in the RANS momentum equation through its divergence . While the Reynolds stress as a symmetric rank-two tensor has six independent components, the divergence as a forcing term only has three components. The majority of existing turbulence models use the Reynolds stress as the target of modeling (Fig. 3). The rationale behind this choice is that the divergence form makes it easier to ensure conservation of momentum. That is, in this form the momentum is introduced into the system by the modeled Reynolds stress only through the boundaries and not within the volume. In contrast, directly constructing such a conservative forcing term is not straightforward [START_REF] Perot | A New Approach to Turbulence Modeling[END_REF]. Therefore, in the remainder of this paper we discuss only turbulence models based on the Reynolds stress .

Reynolds stress based turbulence models require prescribing a constitutive relation for as a function of the mean flow fields. The most widely used class of models, generally known as linear eddy viscosity models, relies on the Boussinesq analogy (see, e.g., [START_REF] Wilcox | Turbulence Modeling for CFD[END_REF]). This assumption states that the anisotropic part of behaves similarly to the viscous stress tensor of a Newtonian fluid, i.e. it is a linear function of the local mean flow rate-of-strain S ij :

k S 2 3 2 ij ij t ij + = (4a) S U x U x with 1 2 , ij i j j i = + (4b)
where ij k ij

2 3 + is the Reynolds stress anisotropy, k u u i i ii 1 2 1 2 = =
is the turbulent kinetic energy with a summation over index i implied, ij is the Kronecker delta (or the second order identity tensor in its vector form I), and the eddy viscosity t is the proportionality scalar.

The limitations of the Boussinesq assumption have been widely recognized in the literature, particularly for flows with separation, streamline curvature, or strong pressure gradients (see, e.g., [START_REF] Wilcox | Turbulence Modeling for CFD[END_REF] for a review). Since it is often not possible to know beforehand if one or more of such flow features will be present in a new flow configuration, predictions based on the RANS equations are flawed by a structural (i.e. model-form) uncertainty [START_REF] Draper | Assessment and propagation of model uncertainty[END_REF][START_REF] Kennedy | Bayesian calibration of computer models[END_REF]. Several attempts have been made to overcome the weaknesses of linear eddy-viscosity models, e.g., by developing nonlinear eddy viscosity models [START_REF] Speziale | On nonlinear k-l and k-ε models of turbulence[END_REF], explicit algebraic Reynolds stress models (EARSM) [START_REF] Gatski | On explicit algebraic stress models for complex turbulent flows[END_REF], and Reynolds stress transport models (RSTM) [START_REF] Wilcox | Turbulence Modeling for CFD[END_REF][START_REF] Launder | Progress in the development of a Reynoldsstress turbulence closure[END_REF]. All such models rely on more sophisticated constitutive relations than Equation (4). Nevertheless, such sophisticated models lack the robustness of the simple linear eddy viscosity models. For example, cubic eddy viscosity models involve many more parameters, which are difficult to calibrate with available data [START_REF] Ray | Learning an eddy viscosity model using shrinkage and Bayesian calibration: a jet-in-crossflow case study[END_REF]. As another example, the Reynolds stress transport equations have a pressure-strain-rate that needs to be modeled, and the predictive performance of RSTM are highly sensitive to its modeling. Consequently, the lack of robustness restricts these advanced models to a small fraction of practical turbulent flows despite their theoretical superiority [START_REF] Spalart | Philosophies and fallacies in turbulence modeling[END_REF], and no turbulence models are able to accurately predict the flow physics in all circumstances. The importance of model uncertainty is clearly illustrated in Fig. 4a, which shows the predicted pressure distribution on the wing section of a Common Research Model (CRM) predicted by a number of turbulence models. A large scattering of the predictions is observed, particularly downstream of the shock wave generated at the upper wing surface.

In addition to the structural uncertainties, parametric uncertainties arising from the coefficients closure models also have to be accounted for. Such coefficients are usually calibrated against experimental data for a set of simple flows (e.g., the decay of homogeneous and isotropic turbulence, flat plate boundary layers, and simple shear flows), which are generally far from practical applications. Moreover, the calibration data suffer from measurement errors, which inevitably impair the credibility of the calibrated parameters. Finally, many of the nominal coefficients found in the RANS modeling literature may not correspond to best-fit of calibration dataset, but were chosen based on numerical considerations. In practice, the closure coefficients are often empirically re-tuned by using heuristic and trial-and-error approaches in order to fit available data for a target class of flows. Fig. 4b illustrates the effect of varying only one of the seven parameters in the algebraic model of Baldwin and Lomax [START_REF] Baldwin | Thin layer approximation and algebraic model for separated turbulent flows[END_REF]. In particular, the location of the shock wave at the airfoil upper surface and the post-shock pressure are very sensitive to the varied coefficient [START_REF] Cinnella | Review of uncertainty quantification in turbulence modelling to date[END_REF].

Both the parametric and the structural uncertainties mentioned above are of epistemic nature, i.e. theoretically they could be reduced when better knowledge of turbulent flow physics and/or more abundant or more accurate data become available. This is in contrast to aleatory uncertainties, which arise from intrinsic variability of a process, e.g., uncertainties in manufactured geometries [START_REF] Parussini | Fictitious domain with least-squares spectral element method to explore geometric uncertainties by non-intrusive polynomial chaos method[END_REF][START_REF] Liu | Quantification of airfoil geometryinduced aerodynamic uncertainties-comparison of approaches[END_REF], operation conditions of turbines or aircraft [START_REF] Avdonin | Quantification of the impact of uncertainties in operating conditions on the flame transfer function with non-intrusive polynomial chaos expansion[END_REF] or inflow conditions [START_REF] Gorlé | Quantifying inflow and RANS turbulence model form uncertainties for wind engineering flows[END_REF][START_REF] Mariotti | Stochastic analysis of the impact of freestream conditions on the aerodynamics of a rectangular 5:1 cylinder[END_REF]. In practice, reducing epistemic uncertainties by leveraging additional knowledge (e.g., by developing more advanced models to incorporate such knowledge) is far from straightforward. Additionally, sophisticated models may lack numerical robustness or incur excessive computational costs. Except for a few canonical examples, it is challenging, if not impossible, to identify the dominant source of uncertainty with definitive evidence, even for a given flow and a specific turbulence model. For instance, in many cases it is possible to improve the results of a model flawed by structural inadequacy by over-tuning its closure parameters. However, such over-tuning typically leads to poor predictions when applying the model to different flows from the calibration flows. Such a phenomenon is referred to as over-fitting in statistics and machine learning [START_REF] James | An Introduction to Statistical Learning[END_REF].

Approaches for quantifying uncertainties in turbulence models

Empirical assessment of uncertainties in turbulence models dates back to the early days of turbulence modeling, but rigorous treatments of such uncertainties in a statistical framework is only a recent development. While it is a consensus that aleatory uncertainties are best represented in a probabilistic framework, different approaches have been pursued for epistemic uncertainties. Because epistemic uncertainties come from lack of knowledge, it is a philosophical question whether to treat such uncertainties in probabilistic framework. In the Bayesian framework, all sources of uncertainty are represented as subjective beliefs and assigned a measure of probability. This review primarily focuses on Bayesian approaches. However, many other non-Bayesian or non-probabilistic approaches for treating epistemic uncertainties exist. Examples include imprecise probability theory [START_REF] Klir | Uncertainty and Information: Foundations of Generalized Information Theory[END_REF], probability bounds analysis [START_REF] Ferson | What Monte Carlo methods cannot do[END_REF][START_REF] Ferson | Different methods are needed to propagate ignorance and variability[END_REF][START_REF] Liu | Arithmetic with uncertain numbers: rigorous and (often) best possible answers[END_REF], Dempster-Shafer evidence theory [START_REF] Shafer | A Mathematical Theory of Evidence[END_REF], fuzzy sets [START_REF] Zadeh | Fuzzy sets, Fuzzy Sets, Fuzzy Logic, and Fuzzy Systems: Selected Papers by Lotfi a Zadeh[END_REF], and credal sets [START_REF] Bernardini | Bounding Uncertainty in Civil Engineering[END_REF]. For an overview and applications of some of these approaches, see Refs. [START_REF] Oberkampf | Verification and Validation in Scientific Computing[END_REF][START_REF] Roy | A comprehensive framework for verification, validation, and uncertainty quantification in scientific computing[END_REF].

Current approaches for quantifying the model-form uncertainties associated with RANS simulations can be classified into parametric and non-parametric approaches 3 depending on where the uncertainties are introduced. In parametric approaches, uncertainties are introduced to the closure coefficients of chosen turbulence models, based on which the overall prediction uncertainties are assessed. Although neglecting uncertainties in the model forms and constrained by the baseline models, the parametric approach has the advantage of being non-intrusive and thus readily available to CFD practitioners. On the other hand, non-parametric approaches directly investigate the uncertainties on modeled terms (fields in RANS solvers), e.g., the eddy viscosity [START_REF] Dow | Quantification of structural uncertainties in the k-ω turbulence model[END_REF], source terms in the turbulent transport equations [START_REF] Singh | Using field inversion to quantify functional errors in turbulence closures[END_REF], or the Reynolds stress itself [START_REF] Xiao | Quantifying and reducing model-form uncertainties in Reynolds-averaged Navier-Stokes simulations: a data-driven, physics-informed Bayesian approach[END_REF][START_REF] Ling | Reynolds averaged turbulence modelling using deep neural networks with embedded invariance[END_REF]. An advantage of these approaches is that the uncertainties of modeled terms reveal more physical insights than the uncertainties of the model coefficients, e.g., allowing the flow regions more prone to model inaccuracies to be identified. However, nonparametric approaches also introduce new challenges, since the uncertainties are now quantified for spatial fields, which theoretically have infinite degrees of freedom. The dimensionality (and thus the cost of the uncertainty quantification) increases with the size of mesh used to discretize the RANS equations. Additionally, such methods are intrusive by nature and thus are less friendly to industrial practitioners who are limited to black-box CFD solvers.

In addition to the parametric/non-parametric classification, it is possible to distinguish forward and backward methods, also referred to as data-free and data-driven approaches as illustrated in Fig. 6. Forward (data-free) methods consist in propagating some pre-specified probability distributions on the closure coefficients (or on the modeled terms) through the RANS equations and investigating the uncertainty distribution of the solution (Fig. 6a). On the other hand, backward (data-driven) methods consist in assimilating available data to infer the coefficient distributions or model errors (Fig. 6b). Such inferred distributions then become available for propagation through the RANS equations in a subsequent prediction step as in the forward analysis. The applicability of the calibrated RANS models to new flows remains as a main concern for both parametric and non-parametric approaches. Table 1 shows a classification of the literature based on their parametric/non-parametric and forward/backward characteristics. Note that the classification omitted data-driven methods that primarily focused on developing turbulence models [e.g., [START_REF] Weatheritt | A novel evolutionary algorithm applied to algebraic modifications of the RANS stress-strain relationship[END_REF][START_REF] Weatheritt | The development of algebraic stress models using a novel evolutionary algorithm[END_REF] rather than quantifying their uncertainties. A roadmap is provided in Fig. 5 to help the reader navigate through this review.

The rest of the paper is organized as follows. A brief review of available techniques for uncertainty propagation, data assimilation and statistical inference is presented in Section 2. In Section 3 we review parametric and multi-model approaches, the latter of which partly accounts for model-form uncertainties. Section 4 is dedicated to nonparametric approaches, which target model-form uncertainties. For completeness, an overview of uncertainties in scale-resolving approaches, and more specifically LES, are briefly reviewed in Sections 5. Finally, conclusions, future research, and perspectives are presented in Section 6.

Fundamentals of probability and statistics for uncertainty quantification

Probability and statistics lie at the core of most of the work reviewed in this work. Therefore, we provide a brief overview of the relevant methods in this section in the context of quantifying and reducing RANS model uncertainties. Based on these foundations, we briefly introduce the algorithms used for uncertainty propagation (forward analysis) and Bayesian inference (backward analysis). In particular, we discuss some commonly used methods for exact and approximate Bayesian inferences.

Representation, sampling, and propagation of model uncertainties

In the probabilistic approach, the uncertain quantities of concern in the RANS model, such as the model coefficients, can be represented as random variables. A random variable Z is a scalar function that may take a range of possible values z, referred to as realizations. A vector of random variables

Z Z Z [ , , ] n 1 =
, indexed by integers, is a random vector. An example is the combination of coefficients in a RANS model. A random field Z x ( ) is a field of random variables indexed by the spatial coordinate x. It is also referred to as stochastic process when the index is time coordinate t. Random field is a generalization of random vectors to the continuous limit. A continuous random variable can be characterized by its probabilistic distributions such as cumulative distribution function or probability density function p z ( ). Common quantities of interest in uncertainty quantification are statistical moments of the random variables such as expectation Z [ ] and variance z ar[ ], which can be ob- tained via integration over all possible outcomes of Z, e.g.,

Z z p z dz [ ] ( ) , = (5a) Z z Z p z dz ar[ ] ( [ ]) ( ) . 2 = (5b)
The expectations and variances of random vectors and random fields can be obtained by applying Equation ( 5) to each component thereof, recalling that random vectors and random fields are collections of random variables indexed by integers and real numbers, respectively. Moreover, a random vector is further characterized by its covariance matrix K Z Z ov( , )

ij i j =
, which represents the correlation among the components of Z . A generalization of the covariance matrix of random vectors to random fields leads to covariance kernel x x K ( , ), which indicates the pair-wise covariance between the random variables x Z ( ) and x Z ( ) corresponding to locations x and x . The most commonly used covariance kernel for the random fields representing model discrepancies is the squared exponential kernel:

x x x x K l ( , ) exp | | 2 , 2 2 2 = ( 6
)
with σ and l indicating variance and length scale, respectively. Such a kernel implies that the correlation between two random variables depends on their corresponding indexing locations. The farther apart the two locations x and x are, the smaller the correlation between x Z ( )

and x Z ( ) is.
In this work, we consider a RANS-based CFD model M y : ( ; ) , which is parameterized by and maps the latent variables (e.g., geometry, boundary conditions) to an observable output y. The multidimensional random variable can be a vector of model coefficients in parametric approaches or a spatial field in non-parametric approaches, e.g., Reynolds stress field x ( ) or eddy viscosity field x ( )

t
. Two types of analyses can be performed:

• Uncertainty propagation (forward analysis): When the probability distribution p ( ) of the model parameters is known, the prob- ability distribution p y ( ) of the output can be obtained by (i) sam- pling the specified distribution p ( ), e.g., by using a Monte Carlo method, (ii) evaluating the model M, and (iii) aggregating the propagated samples. A typical algorithm for plain Monte Carlo sampling is presented in Appendix A.1. The probability distribution p ( ) that is known on the parameters is referred to as the prior distribution.

• Bayesian inference (backward analysis): When data D is available on the output y, which may be noisy, biased, or incomplete, the Fig. 6. Illustration of uncertainty propagation (forward analysis) and statistical inference (backward analysis) in the context of RANS simulations. Uncertainty propagation (forward analysis) involves propagating specified prior distributions on the input (e.g., angle of attack/AoA, Reynolds number, model coefficients, or modeled terms such as Reynolds stresses) through a RANS simulation code and investigate the uncertainties in the solutions (quantities of interests/QoIs, e.g., lift and drag coefficients). Statistical inference (backward analysis) involves assimilating available measurement data to reduce uncertainties in the aforementioned input (e.g., AoA or Reynolds number). The inferred input distributions can be subsequently propagated to make predictions on the QoIs.

Table 1

Classification of literature of RANS model uncertainty quantification based on parametric/non-parametric approaches and data-free (forward)/data-driven (backward) approaches. Works in multi-model approaches are listed along with parametric approaches. input probability distribution of can be inferred. The result is the posterior distribution p ( | ) D of given data D , representing an updated distribution from the prior distribution p ( ) after observing the data.

Uncertainty propagation (forward analysis)

Techniques to propagate uncertainties can be classified into two categories [see ref. [START_REF] Le Maître | Spectral Methods for Uncertainty Quantification: With Applications to Computational Fluid Dynamics[END_REF], Chapter 1.4]: spectral methods [START_REF] Ghanem | Stochastic Finite Elements: a Spectral Approach, revised edition[END_REF] and Monte Carlo (MC) methods [START_REF] Glasserman | Monte Carlo Methods in Financial Engineering[END_REF]. Spectral methods discretize the uncertainty space of the random variables by using orthogonal basis functions. This is done in a similar way in which orthogonal basis functions (e.g., Fourier functions or orthogonal polynomials) are used for the spatial discretization of deterministic PDEs. In uncertainty quantification, spectral methods have faster statistical convergence but they depend on the smoothness of the prior and the function that maps the inputs to outputs. Another barrier for spectral methods is the "curse of dimensionality": the number of function evaluations needed to accurately describe the statistics increases exponentially with the cardinality of the parameter space. Monte Carlo methods, on the other hand, approximate the solution by using random samples from the input uncertainty space and are not adversely affected by its dimensionality. However, the convergence rate is uniformly slow at a rate of N ( )

1/2

O

, where N is the number of samples [START_REF] Glasserman | Monte Carlo Methods in Financial Engineering[END_REF].

While the Monte Carlo based uncertainty quantification seems straightforward, the slow convergence rate poses a major challenge in applications where the computational cost of propagating each sample is high, as is the case for CFD simulations. Accelerating the statistical convergence of Monte Carlo methods has been a topic of intensive research, and numerous techniques for variance reduction have been proposed. Examples include stratified sampling, Latin hypercube sampling [START_REF] Helton | Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems[END_REF], importance sampling, and control variate [START_REF] Glasserman | Monte Carlo Methods in Financial Engineering[END_REF]. A recent development is multilevel Monte Carlo (MLMC) methods [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF][START_REF] Barth | Multi-level Monte Carlo finite element method for elliptic PDEs with stochastic coefficients[END_REF][START_REF] Mishra | Multi-level Monte Carlo finite volume methods for nonlinear systems of conservation laws in multi-dimensions[END_REF], where simulations on coarser meshes are used as control variate of those on fine meshes to reduce the variances. A generalization of MLMC has led to multi-fidelity Monte Carlos methods [START_REF] Müller | Multilevel Monte Carlo for two phase flow and Buckley-Leverett transport in random heterogeneous porous media[END_REF][START_REF] Müller | Solver-based vs. grid-based multilevel Monte Carlo for two phase flow and transport in random heterogeneous porous media[END_REF][START_REF] Peherstorfer | Survey of multifidelity methods in uncertainty propagation, inference, and optimization[END_REF], where a sequence of models with ascending fidelities (e.g., empirical formulas, panel methods, RANS, LES) are combined for input uncertainty propagation, with lower-fidelity models used as control variate of higher fidelity models as in the MLMC methods. However, so far these methods have been primarily used for propagating input uncertainties and not model uncertainties. One difficulty associated with multi-level and multi-fidelity methods is the possible non-trivial interactions between model uncertainties and numerical discretization uncertainties.

Another approach for overcoming the difficulty of expensive model simulations are surrogate models or response surface methods. In these methods, a surrogate of the original model, e.g., in the form of splines, polynomial chaos, or neural networks, are first constructed based on data obtained by evaluating the original model M at a number of design points. The surrogate models provide an approximate functional mapping M y :

that replaces the true mapping M for use in the subsequent sample propagation. Once constructed, the surrogate models can be evaluated at negligible computational costs. However, as with spectral methods, a main difficulty for the surrogate model approach is the curse of dimensionality, which makes it impractical for high dimensional input space.

Statistical inference (backward analysis)

Most of the works on inference of model uncertainties (referred to as backward analysis above) are based on Bayes' theorem: D is the total probability of observing the data, which normalizes the posterior probability.

p p p p ( | ) ( | ) ( ) ( ) , = D D D ( 7

Bayesian inference based on Markov chain Monte Carlo sampling

Theoretically, evaluating the posterior can be straightforward using the following procedure similar to the plain Monte Carlo sampling: (i) draw samples from the prior, (ii) evaluate the likelihood for each sample, and (iii) aggregate the samples to estimate the posterior. However, this is much more challenging than in the forward analysis above. In the forward analysis the probability distribution is known, and thus one can draw more samples from the high probability regions, e.g., by using stratified sampling [START_REF] Glasserman | Monte Carlo Methods in Financial Engineering[END_REF]. In contrast, Bayesian inference involves sampling from the posterior, the high probability regions of which is not known a priori. For example, samples drawn from regions with high prior probability may turn out to have very small likelihood after an expensive model evaluation, which may lead to very small posterior probability (see Equation ( 7)). Therefore, plain Monte Carlo methods are rarely used due to its difficulty in efficiently targeting the high posterior regions. Instead, Markov chain Monte Carlo (MCMC) methods are commonly used, which are a class of sequential sampling strategies in which the next sampled state only depends on the current state. Such a strategy allows the sampling to focus on high probability regions with occasional excursion to low probability regions (tails). Given a target distribution, the MCMC algorithm samples from that distribution by constructing a Markov chain whose stationary distribution coincides with the target distribution. A typical MCMC algorithm with Metropolis-Hastings sampling is detailed in Appendix A.2 and illustrated graphically in Fig. 7.

While the MCMC is the golden standard of Bayesian inference and posterior sampling, a major challenge of its application is that it requires a large number of samples to achieve statistical convergence. Typically the required number of samples range from [START_REF] Wilcox | Turbulence Modeling for CFD[END_REF] 5

O to (10 ) 6 O
, with the specific number depending on the shape of the posterior distribution and the effectiveness of the sampling. In CFD applications, each evaluation involves a simulation that takes hours or even weeks to run depending on the complexity of the flow configuration. For example, RANS simulations of a jet in crossflow, which is a geometrically simple yet industrially relevant case, needed [START_REF] Wilcox | Turbulence Modeling for CFD[END_REF] 7

O
grid points and

O CPU hours to run on a high performance computing cluster [START_REF] Ray | Learning an eddy viscosity model using shrinkage and Bayesian calibration: a jet-in-crossflow case study[END_REF][START_REF] Ray | Bayesian parameter estimation of a k-ε model for accurate jet-in-crossflow simulations[END_REF]. Clearly, it is impractical to perform a full RANS simulation for each evaluation of likelihood in the MCMC sampling. This is not only due to the large number of required samples but also because of the sequential nature of traditional MCMC algorithms -the next proposed sample depends on the evaluated posterior at the current state.

As in the uncertainty propagation discussed above, surrogate models are commonly used for likelihood evaluation in MCMC-based model uncertainty quantification to alleviate the high computational cost of RANS simulations [START_REF] Ray | Learning an eddy viscosity model using shrinkage and Bayesian calibration: a jet-in-crossflow case study[END_REF][START_REF] Ray | Bayesian parameter estimation of a k-ε model for accurate jet-in-crossflow simulations[END_REF][START_REF] Edeling | Bayesian estimates of parameter variability in the k-ε turbulence model[END_REF]. Efficient sampling of high dimensional spaces with MCMC is a topic of active research, with many methods proposed in the past few years, e.g., by adaptively constructing local approximations during the sampling and by using the likelihood to inform the sampling [see, e.g., [START_REF] Moselhy | Bayesian inference with optimal maps[END_REF][START_REF] Parno | Transport map accelerated Markov chain Monte Carlo[END_REF].

Another difficulty arises from the physical constraints among the state variables (e.g., parameters in closure models or Reynolds stresses at different spatial locations), which is particularly relevant for RANS model uncertainty quantification. For example, in the parametric approach such constraints on the parameters dictates that points in some regions in the state space may yield nonphysical solutions or fail to converge at all. Consequently, such regions should be excluded when using MCMC to sample the posterior. Again, this can be done by building surrogate models from simulation data [START_REF] Ray | Learning an eddy viscosity model using shrinkage and Bayesian calibration: a jet-in-crossflow case study[END_REF][START_REF] Ray | Bayesian parameter estimation of a k-ε model for accurate jet-in-crossflow simulations[END_REF][START_REF] Guillas | Bayesian calibration of the constants of the k-ε turbulence model for a CFD model of street canyon flow[END_REF]. The fraction of excluded regions increases exponentially with the dimension of the sample space. Finally, it is noted that such a surrogate approach is also restricted to state spaces with low dimensions.

Approximate Bayesian inference based on MAP estimation

The MCMC method provides the most accurate sampling of the posterior but requires a large number of samples. When the exact probability is not critical and only the low order moments such as the mean and the variance are important, various approximate Bayesian inference methods can be used [e.g., [START_REF] Mons | Reconstruction of unsteady viscous flows using data assimilation schemes[END_REF][START_REF] Iglesias | Ensemble Kalman methods for inverse problems[END_REF]. These methods use maximum a posteriori (MAP) probability estimate to obtain the mode (peak) of the posterior and not the full posterior distribution.

The MAP estimate can be computed in several ways, among which the most commonly used are variational methods and ensemble methods. Both methods are used in data assimilation with a wide range of applications ranging from numerical weather forecasting to subsurface flow characterization. Both variational methods and ensemble methods have been adopted for parameter inferences. To this end, the system state is first augmented to include both the observable, physical state t ( ) i (e.g., velocities, pressure, and/or turbulent kinetic energy) and parameters (e.g., model coefficients or Reynolds stress discrepancies, which are not observable and need to be inferred). Specifically, z is written as a vector formed by stacking the unknown parameters and the physical states i :

z [ , , ; ] , n 1 = ( 8
)
where indicates vector transpose, and [ ,,, ]

r 1 2
= is a vector of r parameters. When computing the MAP estimate, the following objective function is to be minimized:

z z z J y [ ] [ ] = + H P 2 R 2 1 1 (9)
where P and R are the covariance matrices of the state z and the ob- servation errors, respectively, with A A A similarly defined; H is the observation matrix, which maps the state space to the observation space, typically reducing the dimension dramatically. Its interpretation will be further detailed in the context of the ensemble Kalman filtering algorithm (see Appendix A.3).

Obtaining the MAP estimate is equivalent to minimizing the cost function J in Equation ( 9) under the constraint imposed by the models describing the physical system (i.e., RANS equations in case of turbulent flows), during which the set of parameters minimizing the discrepancies between the prediction and the observation data is sought.

In variational methods the minimization problem is often solved by using gradient descent methods, with the gradient obtained with adjoint methods. In contrast, ensemble methods use samples to estimate the covariance of the state vector, which is further used to solve the optimization problem. Variational methods have been the standard in data assimilation and still dominate the field, while ensemble methods such as ensemble Kalman filtering have matured in the past decades and are making their way to operational weather forecasting. Hybrid approaches combining both approaches are an area of intense research and have been explored in CFD applications [START_REF] Mons | Reconstruction of unsteady viscous flows using data assimilation schemes[END_REF].

Recently, ensemble Kalman filtering (EnKF) [START_REF] Evensen | The ensemble Kalman filter: theoretical formulation and practical implementation[END_REF][START_REF] Evensen | Data Assimilation: The Ensemble Kalman Filter[END_REF] has been widely used in inverse modeling to estimate model uncertainties [START_REF] Xiao | Quantifying and reducing model-form uncertainties in Reynolds-averaged Navier-Stokes simulations: a data-driven, physics-informed Bayesian approach[END_REF][START_REF] Iglesias | Ensemble Kalman methods for inverse problems[END_REF]. In EnKF-based inverse modeling, one starts with an ensemble of model parameter values drawn from their prior distribution. The filtering algorithm uses a Bayesian approach to assimilate observation data (e.g., data from experiments and high-fidelity simulations) and produces a new ensemble that represents the posterior distribution. In parametric or field inference of concern here, the EnKF method is used in an iterative manner to find the states that optimally fits the model and data with uncertainties of both accounted for, which is essentially a derivative-free optimization. As such, it is referred to as the iterative ensemble Kalman method. This is in contrast to the EnKF-based data assimilation as used in numerical weather forecasting, where the observation data arrive sequentially. The algorithm for the iterative ensemble Kalman method is presented in Appendix A.3.

EnKF has some well known limitations due to its assumptions of linear models and Gaussian distributions, and theoretically they would perform poorly for non-Gaussian priors and highly nonlinear forward models. However, despite the above-mentioned limitations, EnKF methods have been successfully used in a wide range of applications. Mathematicians have performed analyses to shed light on why they have worked well in view of their theoretical limitations [START_REF] Ernst | Analysis of the ensemble and polynomial chaos Kalman filters in Bayesian inverse problems[END_REF][START_REF] Schillings | Analysis of the ensemble Kalman filter for inverse problems[END_REF].

Parametric and multi-model approaches

In this review we use "parametric approaches" to refer to methods that quantify the uncertainty associated with RANS simulations by investigating primarily the sensitivity of the results to the closure coefficients. As mentioned in Section 1, we will use "forward approaches" to refer the methods that consist of perturbing the closure coefficient according to some probability distribution function and quantifying the output uncertainty on the computed solution. This is in contrast to "backward approaches", in which observations are used to infer the model coefficients. In both the forward and backward approaches, the model structure is fixed and only the uncertainty on the coefficients is quantified. This can nevertheless be used to learn about structural inadequacy of the model, as will be shown later. In some cases, one of the outcomes of the inference process is an estimate of the plausibility of a given model based on the available observations, i.e. the inference may also provide some guidelines for model selection. Finally, we will discuss multi-model approaches in which the uncertainty on the model choice is tackled by considering a set of alternative model structures.

Uncertainties in RANS model parameters

All RANS models have some closure coefficients. A typical example is provided by the well-known k-ε model, initially proposed by Jones and Launder [START_REF] Jones | The prediction of laminarization with a two-equation model of turbulence[END_REF]. In this model, the Reynolds stress tensor is modeled by using the Boussinesq approximation in Equation ( 4), and the turbulent viscosity t is computed by solving additional transport equations for the turbulent kinetic energy k and the turbulent dissipation ε:

C k t µ 2 = (10a) k t U k x x k x i i k i t k i + = + + P (10b) t U x C k C k x x , i i k i t i 2 1 2 + = + + P (10c)
where k P is the production of turbulent kinetic energy through energy extraction from the mean flow gradient:

S S : , k ij ij = P ( 11
)
and the double dot ":" indicates tensor contraction.

The k-ε model above involves coefficients C µ , C 1 , C 2 , k , and . The nature of these coefficients leads to ambiguity regarding their values, and a set of flow-independent optimal values are unlikely to exist [START_REF] Edeling | Bayesian estimates of parameter variability in the k-ε turbulence model[END_REF]. The above-mentioned coefficients are traditionally calibrated to reproduce results of a few canonical flows. One of such canonical flows is the decaying homogeneous isotropic turbulence. In this flow the k and ε equations (10b) and (10c) simplify to

dk dt , = (12)
d dt C k and . 2 2 = ( 13
)
These equations can be solved analytically to give

k t k t t ( ) , n 0 0 = ( 14
)
with reference time t nk

/ 0 0 0 = and the exponent n C 1/( 1) 2 =
, the latter of which leads to:

C n n 1 . 2 = + ( 15
)
The standard value for n is such that C 1.92 2 =

. However, this is by no means a hard requirement and other models do use different values for The coefficient C µ is calibrated by considering the approximate balance between production and dissipation which occurs in free shear flows or in the inertial part of turbulent boundary layers. This balance can be expressed as

U x C k U x . k t µ 1 2 2 2 1 2 2 = = = P ( 16
)
Here, indices 1 and 2 indicate coordinates in streamwise and wallnormal directions, respectively. Equation ( 16) can be manipulated together with the turbulent-viscosity hypothesis

U x / t 12 1 2 = to yield U x ( / ) 12 1 2 1 = , which in turn yields C k . µ 12 2 = ( 17
)
DNS data [START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF] were used to show that k 0.30 [START_REF] Girimaji | Partially-averaged Navier-Stokes model for turbulence: a Reynoldsaveraged Navier-Stokes to direct numerical simulation bridging method[END_REF] (except close to the wall), and thus C 0.09 µ = is the recommended value [START_REF] Pope | Turbulent Flows[END_REF]. Again, however, different models use different values for C µ . For example, it was found that C 0.085 µ in the case of the RNG k-ε model. Another fundamental flow to be considered is the fully developed plane channel flow, which implies that Dk Dt D Dt / / 0 = = . The re- sulting simplified governing equations leads to the following constraint among several parameters [START_REF] Pope | Turbulent Flows[END_REF]:

C C C ( ), µ 2 1/2 2 1 = ( 18
)
where κ is the von Karman constant. It should be noted that the nominal coefficients in the k-ε model satisfy this constraint only approximately, leading to 0.43, instead of the standard value of 0.41 = . However, even the "standard values" has been questioned recently, with κ determined to fall in the range [0.33, 0.45] based on experimental data in the literature [START_REF] Zanoun | Evaluating the law of the wall in two-dimensional fully developed turbulent channel flows[END_REF].

The following constraint between C 1 and C 2 can be found by manipulating the governing equations of uniform (i.e., U x / constant 1 2 = ) shear flows [START_REF] Pope | Turbulent Flows[END_REF]:

C C 1 1 . k 2 1 = P ( 19
)
Tavoularis and Karnik [START_REF] Tavoularis | Further experiments on the evolution of turbulent stresses and scales in uniformly sheared turbulence[END_REF]

<

, would lead to contracting jet, which is nonphysical [START_REF] Durbin | Personal Communication[END_REF].

The parameter k can be considered a turbulent Prandtl number and represents the ratio of the momentum eddy diffusivity and the TKE diffusivity. These quantities are usually close to unity, which is why the standard value for k is assumed to be 1.0. However, no experimental data can be found to justify this assumption [START_REF] Platteeuw | Uncertainty quantification applied to the k-ε model of turbulence using the probabilistic collocation method[END_REF], leading to a range of recommended values among the different variations of the k-ε model. For instance, the RNG k-ε model uses 0.72 k = [START_REF] Wilcox | Turbulence Modeling for CFD[END_REF]. The parameter controls the diffusion rate of ε, and its value can be determined by using the constraint imposed by Equation ( 18), i.e. Similar uncertainties affect the coefficients of other turbulence models. Margheri et al. [START_REF] Margheri | Epistemic uncertainties in RANS model free coefficients[END_REF] discuss in further detail the uncertainties in the coefficients of the k-ε model and Menter's SST k-ω model [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF] and characterized their probability distributions by using generalized polynomial chaos approximations of extensive literature databases. Recently Shaefer et al. [START_REF] Schaefer | Uncertainty Quantification and Sensitivity Analysis of SA Turbulence Model Coefficients in Two and Three Dimensions[END_REF] also investigated the uncertainties in the coefficients of the SA model [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF], Wilcox' k-ω model, and the SST k-ω model, pointing out the large epistemic intervals on their values.

C C C ( ) .

Parametric uncertainty in RANS models: forward approaches

In light of the scattering in closure coefficients of RANS models as reviewed above, several uncertainty quantification (UQ) analyses have focused on quantifying the effect of such uncertainties on the output quantities of interest (QoI). Among the pionneering efforts on forward sensitivity analysis of RANS models is the work of Turgeon et al. [START_REF] Turgeon | Application of a sensitivity equation method to the k-ε model of turbulence[END_REF]. They investigated the effect of uncertainty in the C µ , C 1 , C 2 , k and of the standard k-ε turbulence model (combined with wall functions) on the solution output. The uncertainty analysis was based on a generalized sensitivity equation method [START_REF] Turgeon | A general continuous sensitivity equation formulation for the k-ε model of turbulence[END_REF], i.e. using sensitivity derivatives to propagate uncertainties in the turbulence model coefficients to the solution. In these papers, the uncertainty intervals of the turbulence coefficients are taken arbitrarily, since finding information about the range of uncertainty in the coefficients is not straightforward. The results presented for the flow past a flat plate and the flow over a backward-facing step, which is a severe configuration for RANS models, show that the uncertainty in the model coefficients is not sufficient to account for the observed discrepancies between the predictions and the measurements. An interesting by-product was the identification of the most influential parameters based on the scaled sensitivities. For the flow over a backward-facing step, parameters C

1 and C 2 are found to exert the strongest influence on the wall friction coefficient C f and thus on the reattachment point location. Fig. 8 shows the nominal prediction and the uncertainty range for the distribution of C f downstream of the step (panel a) and of its scaled sensitivities (panel b), defined as

C C , f f j j nom,

=

where j is the j-th model parameter and j nom, is its nominal value. The method was finally applied to an airfoil flow, showing the increasing sensitivity of the solution to the RANS coefficient for larger angles of attack.

Sensitivity-based analyses provide only an uncertainty band around the nominal solution. To obtain more information about the uncertainty of the solution, and specifically its full probability distribution given some input joint probability of the model parameters, UQ techniques (e.g., the MC method presented in Appendix A.1) can be used to propagate an assigned joint distribution on the closure coefficients across the model. For instance, Platteeuw et al. [START_REF] Platteeuw | Uncertainty quantification applied to the k-ε model of turbulence using the probabilistic collocation method[END_REF] used experimental databases and DNS results, along with physical constraints on the coefficients to construct realistic a priori approximations of the input distributions for the different coefficients of the standard k-ε with wall functions [START_REF] Launder | Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc[END_REF]. Their final set of uncertain coefficients includes the model parameters C µ , C 2 , k , the wall function parameters κ (i.e. the von Karman constant), and the log-law constant, as well as the turbulence intensity imposed at the free-stream boundary. A probabilistic collocation method was used to efficiently propagate the input joint distribution through a zero-gradient flat plate flow configuration. Mean flow variations as a consequence of the turbulence model uncertainty were found to be large enough (at least compared to numerical errors) to encompass the experimental data available for the friction coefficient distribution along a flat plate. They also carried out a sensitivity analysis of the output QoI, showing that the solution was more sensitive to the wall function parameters than to other model parameters. Fig. 9 shows the uncertainty range obtained by assigning a normal probability density to the von Karman constant, (0.417, 0.0127) N

, while keeping other parameters constant. The most probable solution is in slightly better agreement with the experimental data. On the other hand, the predicted uncertainty interval encompasses the data.

Forward UQ for the k-ε turbulence model with wall functions was also carried out by using the Latin hypercube sampling method [START_REF] Dunn | Uncertainty quantification of turbulence model coefficients via Latin hypercube sampling method[END_REF]. This was used to propagate distributions of the input coefficients estimated from the data from Pope [START_REF] Pope | Turbulent Flows[END_REF] for the flow past a backwardfacing step, and the mean values were reported for the flow output parameters of interest along with their associated uncertainties. The results showed that model coefficient variability had significant effects on the streamwise velocity component in the recirculation region near the reattachment point and turbulence intensity along the free shear layer. The reattachment point location, pressure, and wall shear were also significantly affected.

In the above-mentioned works, the uncertainty distributions of the input parameters were all obtained in a largely subjective manner. The specification of such prior distribution has an impact on the output probability distributions. To reduce such uncertainties it is possible to use analytical relationships allowing to express the closure coefficients in terms of basic properties of canonical flows (e.g., the power-law exponent of the free decay of turbulent kinetic energy in isotropic turbulence). Following this idea, Margheri et al. [START_REF] Margheri | Epistemic uncertainties in RANS model free coefficients[END_REF] carried out an extensive literature survey and collected a large amount of experimental and numerical data characterizing the input coefficient distributions for the Launder-Sharma low-Reynolds number k-ε and Wilcox k-ω models. The collected data exhibited a significant scattering, which confirmed the hypothesis that the uncertainties in the measured or computed basic flow properties leads to uncertainties in the RANS model coefficients. location of the turbulent kinetic energy. The observed inaccuracies were ascribed to structural uncertainties of turbulence models, which are not accounted for by the parametric data-free approaches.

Parametric uncertainty in RANS models: backward approaches

Statistical inference of model parameters

Forward parametric approaches strongly rely on the availability of reliable data for constructing the coefficient probability intervals or joint distributions. Unfortunately this information is inevitably incomplete and subject to errors. Additionally, it remains restricted to rather simple flow configurations, and it is difficult to extend such data for robust predictions of different flows. Finally, data are only available for observable quantities (e.g., pressures and velocities) and not for the closure coefficients themselves. However, an inverse statistical problem can be solved to infer the input coefficients and possibly their uncertainties. Once obtained, this information can be propagated back through the model to estimate uncertainty intervals on the output QoIs.

The inverse statistical problem can be solved by using a deterministic or probabilistic approach. In the deterministic approach, a set of optimal closure coefficients is obtained by minimizing the model error with respect to some reference data. For instance, Margheri et al. [START_REF] Margheri | Epistemic uncertainties in RANS model free coefficients[END_REF] utilized the gPC response surfaces generated for their forward UQ analyses to find optimal combinations of model coefficients that lead to minimum global error on the mean and friction velocities with respect to DNS data for the turbulent channel flow case. Their findings suggest that the values of the model coefficients recommended in literature, which are generally set as default in commercial and open-source CFD codes, do not fall within the best-fit range. Note however that such deterministic estimates do not provide information on the variability of the optimal coefficients or their validity for a different flow case.

In order to quantify and reduce the uncertainties on model coefficients while simultaneously providing an estimate of model-form uncertainties, it is possible to use Bayesian inference techniques as in Section 2.3. In such an approach, a priori knowledge or assumptions about the coefficients is updated by using available data. When data are highly uncertain or sparse, the updated information will exhibit little difference from the prior distribution. As more data arrive, it is possible to further update the model, thus refining the initial estimate. In the Bayesian calibration process, a key ingredient is the likelihood function in Equation ( 7), which may carry information about observational noise on the data and model-form uncertainty. The latter being the gap between the average model predictions and the "truth", as will be discussed later in Section 3.3.2.

Cheung et al. [START_REF] Cheung | Bayesian uncertainty analysis with applications to turbulence modeling[END_REF] performed the first application of Bayesian uncertainty quantification techniques for calibrating turbulence models and making probabilistic predictions for new flows. They used MCMC sampling to carry out Bayesian calibration of the Spalart-Allmaras model from velocity and skin friction data for three boundary layers with zero, adverse, and favorable pressure gradients. This effort enabled the estimation of the whole posterior joint probability distribution of the coefficients (instead of deterministic values) as well as a comparison of competing models for the likelihood function (noted M 1 , M 2 , and M 3 ) relating the observed data to the model output. As an ex- ample, Fig. 11 shows the marginal posterior distributions obtained for the von Karman constant κ and the coefficient c ,1 , along with their joint scatter plot when using the stochastic model M 3 . Bayesian calibration is able to discover a posterior correlation between these two parameters, showing the importance of calibrating all parameters simultaneously. The MCMC-based calibration process involved a large number of boundary layer calculations (32,768 samples), each based on a full Navier-Stokes incompressible flow solver. Ray and co-workers [START_REF] Ray | Learning an eddy viscosity model using shrinkage and Bayesian calibration: a jet-in-crossflow case study[END_REF][START_REF] Ray | Bayesian parameter estimation of a k-ε model for accurate jet-in-crossflow simulations[END_REF][START_REF] Lefantzi | Estimation of K-ε Parameters Using Surrogate Models and Jet-in-crossflow Data[END_REF][START_REF] Ray | Robust Bayesian calibration of k-ε model for compressible jet-in-crossflow simulations[END_REF] used a similar approach to infer the model coefficients for a more complex configuration, namely, a jet-in-cross-flow. For example, experimental data were used to calibrate the parameters in a nonlinear eddy viscosity model [START_REF] Ray | Learning an eddy viscosity model using shrinkage and Bayesian calibration: a jet-in-crossflow case study[END_REF], where surrogate models were used to reduce the computational burden of the MCMC sampling.

Kato and Obayashi [START_REF] Kato | Approach for uncertainty of turbulence modeling based on data assimilation technique[END_REF] used ensemble Kalman filtering [START_REF] Evensen | The ensemble Kalman filter: theoretical formulation and practical implementation[END_REF][START_REF] Evensen | Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics[END_REF] to determine the values of the parameters of the Spalart-Allmaras turbulence model for the zero-pressure gradient flat plate boundary layer at M 0.2 = and Re 5 10 6 = × . The data were velocity profiles and wall pressures generated by the same model using a known set of coefficients (equal to the nominal ones). An advantage of using synthetic data is to remove structural uncertainty, since the trained model is the same used to generate the data. The results show the ability of the EnKF method to identify the correct model parameters for a relatively low computational cost (ensembles of 100 function evaluations, i.e. CFD calculations). The approach has been extended to more complex flows around airfoils [START_REF] Kato | A data assimilation methodology for reconstructing turbulent flows around aircraft[END_REF], establishing a general framework for combining experimental fluid dynamics and CFD for predictions.

An even more efficient way of finding the optimal coefficients is to maximize the likelihood function by using gradient-based methods. This corresponds to finding the set of closure coefficients corresponding to the maximum probability of observing the data. The main drawback of this approach is that only deterministic sets of coefficients are obtained as an outcome of the calibration. Papadimitriou and Papadimitriou [START_REF] Papadimitriou | Bayesian uncertainty quantification of turbulence models based on high-order adjoint[END_REF] obtained variance estimates of the optimal coefficients by using the Hessian of the likelihood function with respect to the parameters θ. They found that the posterior variance due to the overall observational uncertainty (e.g. to the discrepancy between the model output and the data) plays a dominant role. This indicates that coefficient calibration alone is not sufficient to match the data, and that the bias introduced by the model structure is mostly responsible for the discrepancy. Unfortunately, Hessian calculations require computing the second sensitivity derivatives of the model with respect to the parameters, which is a highly intrusive and delicate task and is not compatible with black-box Navier-Stokes solvers.

Bayesian strategies similar to that of Cheung et al. [START_REF] Cheung | Bayesian uncertainty analysis with applications to turbulence modeling[END_REF] can also provide estimates of the uncertainty associated with the model form, grounded in uncertainties in the space of model closure coefficients. This can be achieved by calibrating the model separately against several sets of data. The spread in the posterior estimates of closure coefficients across calibration scenarios provides a measure of the need for readjusting the model coefficients to compensate for the inadequacy. An example of such a sensitivity study is given by Edeling et al. [START_REF] Edeling | Bayesian estimates of parameter variability in the k-ε turbulence model[END_REF], where the Launder-Sharma model was calibrated separately against 13 sets of flat-plate boundary layer profiles from Kline et al. [START_REF] Kline | AFOSR-IFP-Stanford Conference: Proceedings Held at[END_REF]. The results showed a significant variation in the most-likely closure-coefficients values for the different pressure gradients, despite the relatively restricted class of flows (flat plate boundary layers) considered for the calibrations.

The main lessons learned from the preceding exercise are: (i) there are no universal values for the closure parameters of the turbulence models; (ii) the parameters need to adjust continuously when changing the dataset to compensate the intrinsic inadequacy (simplifying modeling assumptions) of the chosen model (see, e.g., the variation of the marginal posterior pdf for κ, reported in Fig. 12a); and (iii) as a result, closure coefficients obtained by calibrating the model against a given boundary layer are generally not valid for the prediction of a different one.

The variability of closure coefficients for the observed flow can however be used as a measure of model inaccuracy when predicting a new flow. In Edeling et al. [START_REF] Edeling | Bayesian estimates of parameter variability in the k-ε turbulence model[END_REF], this is done by summarizing the posterior variability of the parameters within and in between calibration datasets (called hereafter scenarios) by means of probability boxes (pboxes), commonly used in Bayesian statistics to summarize the joint effect of parametric and epistemic model-form uncertainties [START_REF] Oberkampf | Verification and Validation in Scientific Computing[END_REF]. Pboxes are constructed as the envelope of the empirical cumulative distribution functions of the output predicted using different posteriors of the parameters. An example of p-box for the nondimensional velocity u + predicted at a non-dimensional wall distance y 46.2 = + for a pipe flow boundary layer is given in Fig. 12b. Analogous results are obtained at various locations across the boundary layer, thus leading to an estimate of the uncertainty bounds on the predicted velocity profile (Fig. 12c). The p-box prediction is found to encompass the experimental uncertainty intervals, leading however to an overly conservative estimate of the uncertainty bounds.

Accounting for structural uncertainties in RANS models

A delicate step in Bayesian calibration is the construction of a statistical model relating the true (unseen) process to the data via the model, which is directly related to the definition of the likelihood function. This should consider at least the fact that the observed quantities differ from the true ones by the experimental (observational) noise, which may be expressed through the relation:

z , = + (21)
with ζ the true value for z and ϵ a random vector representative of the experimental noise. The experimental data noise ϵ is often assumed to be independently distributed without spatial correlation, and it is modeled as a Gaussian process with diagonal covariance matrix, i.e., (0, ) N [START_REF] Cheung | Bayesian uncertainty analysis with applications to turbulence modeling[END_REF][START_REF] Edeling | Bayesian estimates of parameter variability in the k-ε turbulence model[END_REF]. Theoretically, the true value for ζ could be obtained as an output of the model y, once a suitable set of parameters θ has been identified, i.e.

y ( ) =

. In practice, however, no model is perfect. Even if there is no parameter uncertainty, so that we know the true values of all the inputs required to make a particular prediction of the process being modeled, the predicted value will not equal the true value of the process [START_REF] Brynjarsdóttir | Learning about physical parameters: the importance of model discrepancy[END_REF]. The discrepancy is due to model inadequacy. It is even possible that the physically true value of a calibration parameter gives a worse fit and less accurate future prediction than other values, simply because of the simplifying assumptions upon which the model has been built. Conversely, it is dangerous to interpret calibration results as estimates of the true physical values of those parameters.

A general framework to include the model inadequacy term in the stochastic model was first proposed in Ref. [START_REF] Kennedy | Bayesian calibration of computer models[END_REF]. Model discrepancy can be taken into account by introducing an additional error term to the statistical model as in Equation ( 21), which could be of additive nature, i.e., z y ( )

= + = + + (22)
or of multiplicative nature:

z y ( ) = + = + ( 23
)
where the symbol denotes the Hadamard (element-wise) multiplication. Note that all quantities above, z, y, η, and ϵ, are spatial fields and should be written as x z ( ), x y ( ), x ( ), and x ( ), respectively. The spatial dependence is omitted for brevity. The choice of model-inadequacy formulation largely depends on the nature and prior knowledge about the observed quantity z. In Equations ( 22) and ( 23), η is a random field representative of the model inadequacy, i.e., of the fact that true value is not equal to the code output but with some systematic deviations. For instance, Cheung et al. [START_REF] Cheung | Bayesian uncertainty analysis with applications to turbulence modeling[END_REF] chose a multiplicative error model for relating the measured and computed velocity profiles, so that the no-slip boundary condition at the solid wall is satisfied by any realization of the stochastic model.

When an additive model inadequacy term is used, it becomes difficult to separate its effect from that of the observational error. As a consequence, both terms are often merged together. In all cases, the random variable η may involve additional parameters proper to the statistical model introduced for describing the error behavior, referred to as hyperparameters. Sometimes these are known beforehand or are estimated independently based on likelihood maximization criteria [START_REF] Nitschke | Model-form and predictive uncertainty quantification in linear aeroelasticity[END_REF], but most often they need to be calibrated from the data along with the physical model parameters . Another important point is that η is expected to correlate modeling errors for a QoI evaluated at different locations in the flow field or even for various QoIs and different datasets. For instance, Cheung et al. [START_REF] Cheung | Bayesian uncertainty analysis with applications to turbulence modeling[END_REF] introduced a multiplicative term to calibrate the Spalart-Allmaras model from velocity profiles and skin friction distributions for three boundary layer data sets. In their work, all the competing multiplicative statistical models describe the inadequacy term as Gaussian process, i.e., (1, ) . Consequently, the observations can also be modeled as a Gaussian process, and thus the likelihood function can be written as follows: 1

N z = = = + ( 24
)
where N is the dimension of the data vector d, and the covariance matrix of the true process is

y y diag( ) diag( ) = based on the defi- nition y ( ) =
in Equation ( 23). It can be seen that the covariance matrix of the modeled multiplicative term has impact upon the true process ζ and thus influences the formulation of the likelihood function as shown in Equation [START_REF] Tinoco | Summary data from the sixth AIAA CFD drag prediction workshop: CRM cases[END_REF]. Different statistical models can be obtained for different choices of the covariance matrix, and more complex choices for the inadequacy term have also been investigated [START_REF] Oliver | Bayesian uncertainty quantification applied to RANS turbulence models[END_REF]. Their results showed the choice of spatial correlation structure for the modeling inadequacy played an important role in the Bayesian model selection.

Although the use of model inadequacy terms such as those of Equations ( 22) and ( 23) is helpful in alleviating parameter overfitting problems and in estimating how well the calibrated model is able to fit the data, the approach suffers from several limitations:

• The correction terms are specific to the observed QoI and cannot be re-used for the prediction of a different (unobserved) QoI.

• The correction often depends on the spatial distribution of observed data for the calibration scenario, and can be hardly applied to, e.g., a different geometry.

• Even if the same QoI and geometrical configuration are considered, the inadequacy terms calibrated for a given dataset may not be valid for the prediction of a different scenario (e.g., operating condition).

The non-universality of the inadequacy term is well illustrated by the results of Edeling et al. [START_REF] Edeling | Bayesian estimates of parameter variability in the k-ε turbulence model[END_REF], who used a statistical model involving a multiplicative model-inadequacy term similar to the correlated model M 3 of Cheung et al. [START_REF] Cheung | Bayesian uncertainty analysis with applications to turbulence modeling[END_REF]. It was observed that the expected value of the model inadequacy term, as calibrated from data, varies significantly from case to case. Additionally, for some of the calibration datasets the posterior values taken by the correction term are much higher than for the other cases, indicating that parameter adjustment was not sufficient for the model output to capture the data.

Accounting for multiple models: Bayesian model selection and averaging

An interesting outcome of Bayesian calibration is the possibility of deriving statistical criteria for model selection, i.e., for choosing the best model in some statistical sense among a class of competing models. This consists in providing estimates of the posterior probability of each model in the considered set of models M M M { , , , }

I 1 2 = M
given the observed data. The "model" here should be interpreted in a broader sense, including not only physical models (e.g., k-ε, k-ω, and Reynolds stress models) with associated coefficients but also statistical models, e.g., the covariance kernel used to construct likelihood functions [START_REF] Cheung | Bayesian uncertainty analysis with applications to turbulence modeling[END_REF][START_REF] Oliver | Bayesian uncertainty quantification applied to RANS turbulence models[END_REF] as in Equations [START_REF] Launder | Progress in the development of a Reynoldsstress turbulence closure[END_REF][START_REF] Ray | Learning an eddy viscosity model using shrinkage and Bayesian calibration: a jet-in-crossflow case study[END_REF][START_REF] Tinoco | Summary data from the sixth AIAA CFD drag prediction workshop: CRM cases[END_REF]. Model probabilities are obtained as an outcome of parameter calibration as introduced above. First, each model in the set M is assigned a probability M ( ) j , j I 1, , = , based on prior knowledge (e.g., from expert elicitation) or the lack thereof, in which case a noninformative, uniform distribution is chosen. Additionally, the prior distributions for the closure coefficients θ or statistical hyperparameters associated with each model are also specified. If data D are available, the prior probability mass function (pmf) can be updated according to Bayes' theorem, leading to the posterior pmf of model M j :

M p M M p M M j I ( | ) ( | ) ( ) ( | ) ( ) 1, , j j j i I i i 1 = = = D D D ( 25
)
where p M ( | ) j

D

is the evidence for model M j that normalizes the posterior pdf of the model parameters θ, as in Equation [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF]. The evidence can be computed at the end of the calibration by numerically integrating the numerator of Equation ( 7), using the posterior samples of θ. This can be a challenging process requiring special techniques [e.g., 111]). The estimated pmf of the models can subsequently be used for predictions by choosing the model with maximum posterior probability in the case of model selection, or alternatively by weighting the various posterior predictive distributions for the QoI with the posterior pmf in the case of model averaging.

The approach above has been used for Bayesian model selection and calibration in RANS simulations. It was found difficult to identify a single best model for a range of flows. Consequently, predicting new (unobserved) flow scenarios based on a single closure model calibrated on a limited dataset may lead to biased results, and thus Bayesian model selection is insufficient. Oliver and Moser [START_REF] Oliver | Bayesian uncertainty quantification applied to RANS turbulence models[END_REF] calibrated the combination of four eddy viscosity models and three statistical models by using DNS data of plane channel flows and compared the posterior probabilities and predictive capabilities. The results showed that the considered data slightly favored Chien's k-ε model [START_REF] Chien | Predictions of channel and boundary-layer flows with a low-Reynoldsnumber turbulence model[END_REF] with an inhomogeneous stochastic model for the inadequacy, but no clear winner emerged with a dominantly high posterior probability. Edeling et al. [START_REF] Edeling | Predictive RANS simulations via Bayesian model-scenario averaging[END_REF] systematically demonstrated the difficulty of identifying a single best model without ambiguity. They used Bayesian inference to compute the posterior probabilities of five turbulence models ranging from simple algebraic eddy viscosity models to sophisticated Reynolds stress models by using DNS data of 13 boundary layer flows of various configurations. The posterior pmf for each dataset are presented in Fig. 13, which suggests that none of the models has a consistently higher probability than other models for all datasets, and the probabilities of all models are highly flow-dependent. As a consequence, it was not possible to select a single best model valid for all flow configurations. Moreover, somewhat surprisingly, the Reynolds stress model was not the most plausible one for all flows despite its theoretical superiority; on the other hand, after calibration the algebraic model performed rather well over a wide range of flow configurations.

The difficulty of making predictions with a single calibrated model clearly calls for a framework based on multi-model ensembles. Multimodel approaches have been used in aerodynamics [START_REF] Poroseva | Improving the predictive capability of turbulence models using evidence theory[END_REF] and many other applications [START_REF] Diomede | Discharge prediction based on multi-model precipitation forecasts[END_REF][START_REF] Duan | Multi-model ensemble hydrologic prediction using Bayesian model averaging[END_REF][START_REF] Tebaldi | The use of the multi-model ensemble in probabilistic climate projections[END_REF]. Bayesian modeling averaging is among the most widely used multi-model approaches, where the posterior of the predicted quantity is [START_REF] Draper | Assessment and propagation of model uncertainty[END_REF][START_REF] Hoeting | Bayesian model averaging: a tutorial[END_REF]:

p I p M M ( | , ) ( | ) ( | ), i i i 1 = = D M D ( 26
)
given calibration data D and a set of models M . In this framework the posterior of is an average of I posterior predictive distributions corresponding to I competing models weighted by their respective posterior model probabilities as computed from Equation [START_REF] Baldwin | Thin layer approximation and algebraic model for separated turbulent flows[END_REF]. A significant recent development is the Bayesian model-scenario averaging (BMSA), which is an extension of the classical Bayesian model averaging as shown in Equations ( 25) and ( 26) above. BMSA accounts for uncertainties on the choice of the calibration flow configuration (referred to as scenario). It predicts the QoI for a new scenario S ˜(not used for model calibration) as a weighted average of the predictions provided by a set of models M { } k . It is important to stress here that, for nonlinear systems, averaging the posterior predictive distributions of the QoI obtained by propagating the posterior pdf of the parameters for various scenarios through each model, as in Equation ( 27), is radically different from creating a mixture of the K pdfs for the closure coefficients and propagating it through the model. Specifically, Ray et al. [START_REF] Ray | Robust Bayesian calibration of k-ε model for compressible jet-in-crossflow simulations[END_REF] showed that latter provided unsatisfactory predictions, albeit being less expensive computationally.

In the BMSA prediction, the posterior probability of model M i is the outcome of the multiple calibration process after application of Equation [START_REF] Baldwin | Thin layer approximation and algebraic model for separated turbulent flows[END_REF]. On the other hand, the scenario probability S ( ) k needs to specified a priori and represents the user's belief about the similarity between calibration scenario S k to prediction scenario S ˜as far as the prediction of is concerned. When a physically justified prior is not available, a non-informative, uniform pmf can be used, implying equal probabilities for all scenarios. However, this may overestimate the posterior variance for , which leads to an overly pessimistic estimate of the prediction uncertainty [START_REF] Edeling | Predictive RANS simulations via Bayesian model-scenario averaging[END_REF]. To address this issue, Edeling et al. [START_REF] Edeling | Bayesian estimates of parameter variability in the k-ε turbulence model[END_REF] proposed an empirical scheme for choosing the scenario prior, with S ( ) k being inversely proportional to the scattering of all models trained on scenario S k when predicting the QoI for S ˜. The rationale is that if a calibration scenario S k is similar to the prediction scenario S ˜, the models would give similar predictions of the QoI.

A major drawback of BMSA is its high computational cost, since it requires I K

• stochastic calculations, each requiring forward propaga- tion of a posterior parameter pdf (corresponding to a model/scenario combination) through the CFD model. The computational cost can be drastically reduced to I K • deterministic CFD simulations by propagating though S ˜only the set of parameters with maximum posterior probability for each model and calibration scenario [START_REF] Edeling | Bayesian predictions of Reynolds-averaged Navier-Stokes uncertainties using maximum a posteriori estimates[END_REF], instead of propagating the full pdf. With this simplification, the BMSA approach was applied to complex flow configurations such as the transonic threedimensional flow around the ONERA M6 wing.

As noted by Draper [START_REF] Draper | Assessment and propagation of model uncertainty[END_REF], multi-model approaches still introduce biases in the prediction because of the subjective selection of a finite set of models. However, they play a useful role in reducing the bias compared to predictions based on a single model. An averaged model is a way of obtaining a conservative prediction for an unseen configuration. Indeed, the result will not be as good as the a priori unknown best model but will not be as bad as the worst one. Additionally, BMSA provides an estimate of the solution variance based on the solution variability among the competing models.

Non-parametric approaches

Motivation and overview

The parametric and multi-model approaches introduced in Section 3 explore the uncertainties in the model coefficients and in the model choices. However, it is possible that the true solution lies outside the region in the solution space reachable by the parametric approaches. For example, it is well-known that linear eddy viscosity models are intrinsically not capable of predicting the secondary flows in a square duct. Such a feature is driven by the anisotropy of the Reynolds stresses, but the Boussinesq assumption that is inherent to linear eddy viscosity models excludes this part of the solution space. This intrinsic deficiency cannot be remedied by the calibration of coefficients. An ensemble or averaging of linear eddy viscosity models would not be able to predict such a feature either, because all the models would strongly agree on the wrong solution. A larger portion of the solution space could be covered by introducing a wider variety of models (namely, non-Boussinesq) in the multi-model ensemble. However, the choice of the set of models remains subjective and the selection of a finite set of models prevents the approach from exploring the entire solution space, limiting it to only the portion spanned by the chosen model ensemble. In order to go beyond these limitations, an intriguing possibility is to introduce uncertainties directly into the turbulent transport equations or the modeled terms such as the Reynolds stress or eddy viscosity. Such nonparametric approaches allow for more general estimates of the model inadequacy than the parametric approaches. As illustrated conceptually in Fig. 14, the solution space explored by parametric approaches is a subspace of that explored by nonparametric approaches, and the true solution may lie outside the former space.

We use the wing-body juncture flow as an example to motivate the use of nonparametric approaches in exploring solution spaces for RANS model uncertainty quantification. This configuration consists of an airfoil attached to a flat plate, which is representative of the flows at the wing-fuselage connection of fixed-wing aircraft and blade-hub assembly in turbomachinery. This flow features an abrupt stagnation of the mean flow at the leading edge and a horseshoe vortex around the juncture of the wing and the body as shown in Fig. 15. Due to the high non-equilibrium turbulence, the Reynolds stress and strain rate S at the leading edge region are not aligned with each other, and thus the Boussinesq assumption fails. [START_REF] Duraisamy | Turbulence modeling in the age of data[END_REF] shows clearly the misalignment between orientations of RANS-modeled Reynolds stress (with SST k-ω model [START_REF] Wu | Quantifying model form uncertainty in RANS simulation of wing-body junction flow[END_REF]) and the experimentally measured Reynolds stress [START_REF] Devenport | Time-dependent and time-averaged turbulence structure near the nose of a wing-body junction[END_REF] at two locations, particularly at the near-wall point P 2 (see Fig. 15). Consequently, when exploring uncertainties in the RANS simulations for this flow, the velocity samples obtained with parametric approaches (gray lines in Fig. 17a) based on the Boussinesq assumption, are not able to encompass the truth (× symbols). This is because such a parametric approach is not able to account for the different eigen-directions Fig. 14. A conceptual illustration of the merit of the non-parametric approach RANS model uncertainty quantification, i.e., the ability to explore the solution space more thoroughly. Figure inspired by Soize [START_REF] Soize | A comprehensive overview of a non-parametric probabilistic approach of model uncertainties for predictive models in structural dynamics[END_REF].

15.

Configuration of the wing-body junction flow, illustrating the points where orientations of the Reynolds stress tensors (Fig. 16) are presented and three lines where the velocity profiles (Fig. 17 of the RANS-modeled and the true Reynolds stresses. In contrast, a nonparametric approach that perturbs the RANS-modeled Reynolds stresses, including their eigen-directions (see gray arrows in Fig. 16), can effectively span a range covering the true solution [START_REF] Wu | Quantifying model form uncertainty in RANS simulation of wing-body junction flow[END_REF], which is illustrated in Fig. 17b.

A number of nonparametric approaches have been proposed to quantify model uncertainties in RANS simulations, which can be broadly classified into two categories:

(1) those introducing uncertainties into the model forms, e.g., turbulent transport equations (for fields k, ω, or ) [START_REF] Singh | Using field inversion to quantify functional errors in turbulence closures[END_REF][START_REF] Parish | A paradigm for data-driven predictive modeling using field inversion and machine learning[END_REF], and (2) those introducing uncertainties into the model outputs, e.g., the turbulent viscosity field [START_REF] Dow | Quantification of structural uncertainties in the k-ω turbulence model[END_REF] or the Reynolds stress field [START_REF] Xiao | Quantifying and reducing model-form uncertainties in Reynolds-averaged Navier-Stokes simulations: a data-driven, physics-informed Bayesian approach[END_REF][START_REF] Ling | Reynolds averaged turbulence modelling using deep neural networks with embedded invariance[END_REF][START_REF] Emory | Modeling structural uncertainties in Reynoldsaveraged computations of shock/boundary layer interactions[END_REF][START_REF] Emory | Modeling of structural uncertainties in Reynolds-averaged Navier-Stokes closures[END_REF].

At the algorithmic level, the different parametric and nonparametric approaches outlined above target different stages of the algorithms in turbulence modeling, i.e., the parametric level, the PDE level, and the intermediate field level. This is illustrated in Fig. 3 previously by using linear eddy viscosity models and Reynolds stress transport models as examples. The intermediate fields and PDEs where uncertainties are introduced are highlighted in shaded (orange) boxes in Fig. 3. At a fundamental level, however, they differ from each other in their respective assumptions on where the RANS model uncertainties originate from: the coefficients, the model form of the transport equations, the eddy viscosity field, or the Reynolds stress itself. As reviewed above, even for a specific flow it is difficult to identify the exact source of the model uncertainty (see Section 1.2) due to the coupling among various levels of uncertainties. As such, any such statements on the relative importance of different sources of uncertainties are likely to be not only flow-specific but also weak and inconclusive [START_REF] Ray | Robust Bayesian calibration of k-ε model for compressible jet-in-crossflow simulations[END_REF]. Consequently, the relative advantages of various approaches are far from clear as of now.

Concerning the comparison between parametric and nonparametric approaches, the parametric approaches allow straightforward extrapolation of the calibrated coefficients to additional flow configurations that are not in the calibration dataset. However, naive extrapolation may lead to an over-fitted model with reduced predictive capability, particularly when the generalization of the coefficients cannot be justified. On the other hand, extrapolating a calibrated field from nonparametric approaches (generally dependent on space and time coordinates) is a much more delicate task. As to the comparison between model-form-based and model-output-based UQ approaches, research so far suggests that model-form based approaches are more robust as they involve only mild perturbations of equations in the original models [START_REF] Singh | Using field inversion to quantify functional errors in turbulence closures[END_REF]. On the other hand, model-output perturbation approaches make it easier to utilize benchmark (DNS, LES, or experimental) data for the Reynolds stress or turbulent viscosity, because the quantities being perturbed or inferred have better physical anchoring. Both categories of approaches will be reviewed and compared below.

Introducing uncertainties in turbulent transport equations

The parametric and multi-model approaches are restricted to the chosen baseline models. An immediate extension of these approaches is to perturb the model forms in a non-parametric way, i.e., by modifying the source terms in the turbulent transport equations (e.g., for k, ω, and

). This choice is based on the assumption that errors in the turbulent transport equations rather than the structural uncertainties, e.g., those associated with the Boussinesq assumption, are the dominant source of the prediction errors in RANS simulations. The uncertainties introduced in this approach depend on the specific form of the baseline turbulence model. Taking the k-ω equation for example, a multiplicative discrepancy field x

) is introduced to the source terms of the ω transport equation by Singh and Duraisamy [START_REF] Singh | Using field inversion to quantify functional errors in turbulence closures[END_REF]:

x D Dt k U D k U T k U ( ) ( , , ) ( , , ) ( , , ) i i = + P ( 28
)
where ω is the turbulent frequency; P , D , T indicate production, dissipation, and transport, respectively, of ω. This formulation is equivalent to introducing an additive discrepancy x ( ( ) 1) P but has better conditioning than the latter [START_REF] Singh | Using field inversion to quantify functional errors in turbulence closures[END_REF]. The discrepancy field

x

) can be inferred by using DNS or experimental data of velocities or other quantities of interest, e.g., drag, lift, pressure coefficient, and surface friction. Assuming the velocity is the data to be used, the inference can be cast as the following optimization problem:

J U U argmin , with ( ) L opt dns 2 2 = = ( 29
)
where L 2 indicates L 2 norm. In cases where other derived quantities g (e.g., drag and lift) are used in the optimization, an observation operator H is needed to map the solution to these quantities, i.e., g

U [ ] = H
, and the cost function would be g g

J ( )

L dns 2

=

. The inferred discrepancy x ( ) opt is a correction that allows the baseline k-ω model to agree with the data. The discrepancy field β resides in a space of very high dimensions with a dimension equal to the number of cells in the CFD mesh, and thus the optimal solution is not unique. In the terminology of inverse modeling, this problem is ill-posed and needs to be regularized. The deviation of β from 1 is used as a penalty to regularize the problem, which leads to the following cost function [START_REF] Singh | Using field inversion to quantify functional errors in turbulence closures[END_REF]:

x J U U with ( ) ( ) 1 L L dns 2 2 2 2 = + ( 30
)
where γ is a regularization parameter. The second term, x ( ) 1, prevents the corrected model from deviating too much from the baseline model. With such a regularization, the corrected model is constrained to explore only the vicinity of the baseline solution, which greatly reduces the dimension of the search in the high-dimensional space of possible discrepancy fields β. The inferred discrepancy field can be subsequently used to guide the improvement of the baseline model and to develop data-driven correction schemes. Singh and Duraisamy [START_REF] Singh | Using field inversion to quantify functional errors in turbulence closures[END_REF] used velocities from DNS databases to infer the discrepancy field in plane channel flows at frictional Reynolds numbers ranging from Re 395 = to 4200. The results are shown in Fig. 18. It can be seen that the profiles of discrepancies β for different Reynolds numbers are qualitatively similar. This a priori study suggests that the knowledge gained in one flow can be extended to other flows of similar configurations where data are not available. The end product is a datadriven correction function q ( ) for the baseline model obtained by posing the discrepancy term β as a function of non-dimensionalized mean flow variables (e.g., S and , both properly normalized with local quantities [START_REF] Ling | Evaluation of machine learning algorithms for prediction of regions of high Reynolds averaged Navier Stokes uncertainty[END_REF], as well as the ratio D / P between production and dissipation [START_REF] Singh | Machine-learning-augmented predictive modeling of turbulent separated flows over airfoils[END_REF]). Choosing flow variables q rather than spatial co- ordinates x as the input of the regression enables generalization of the learned function in different flows, possibly at different spatial scales. Singh et al. [START_REF] Singh | Machine-learning-augmented predictive modeling of turbulent separated flows over airfoils[END_REF] showed predicted pressure of the S809 airfoil at Re 2 10 6 = × by using the SA model augmented by the correction function, which was trained with the inferred discrepancy field by using the from the S814 airfoil at Reynolds numbers Re 1 10 6

= ×

and Re 2 10 6 = × . Although the correction scheme is applied on a few specific models (k-ω or SA model), generalization to additional models (e.g., k-ε model or Reynolds stress transport model) is straightforward. On the other hand, since the corrected model is obtained by perturbing the transport equations in the baseline model, it is still constrained by the limitation of the latter. For example, if a linear eddy viscosity model is chosen as baseline, the corrected model would still be limited by the Boussinesq assumption.

Introducing uncertainties in turbulent viscosity

Most of the widely used turbulence models (e.g., k-ε, k-ω, and SA models [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF]) are linear eddy viscosity models, which model the Reynolds stress in the form S 2

k t 2 3 =
as in Equation (4a), where t is the turbulent eddy viscosity. It can thus be assumed that the model uncertainty in RANS simulations can be attributed to the discrepancies in the predicted eddy viscosity field and subsequently introduce uncertainties thereon. For flows with ground truth of mean velocities (e.g., U dns from DNS data), one can define an optimal eddy viscosity field opt that minimizes the discrepancy between the computed velocity u and the ground truth velocity. Finding the optimal viscosity opt amounts to solving the following optimization problem [START_REF] Dow | Quantification of structural uncertainties in the k-ω turbulence model[END_REF]:

J U argmin , with ( ) t L opt dns 2 t 2 = = (31)
where ) t indicates the dependence of the velocity field on the eddy viscosity field through the RANS equations. The optimization is further constrained by the positivity and smoothness of , which can be built into the cost function or enforced in the optimization procedure. For example, Dow and Wang [START_REF] Dow | Quantification of structural uncertainties in the k-ω turbulence model[END_REF] used the following cost function:

J U ( ) t L t L dns 2 2 2 2 = + ( 32
)
where a regularization term with is incorporated to promote smoothness of the viscosity field with parameter γ controlling the desired smoothness. The optimization problem can be solved with gradient descent methods, where the gradient J/ t of the cost function with respect to the control variable x ) t can be obtained very efficiently by using adjoint methods. Alternatively, the optimization problem in Equation ( 31) can also be solved by the iterative ensemble Kalman method [START_REF] Iglesias | Ensemble Kalman methods for inverse problems[END_REF], which can be considered a derivative-free optimization that uses the state covariance, estimated from Monte Carlo samples, instead of the Jacobian. The iterative Ensemble Kalman method has been used to infer the Reynolds stresses discrepancies by using sparse observation data of velocities [START_REF] Xiao | Quantifying and reducing model-form uncertainties in Reynolds-averaged Navier-Stokes simulations: a data-driven, physics-informed Bayesian approach[END_REF] (see Section 4.4).

The viscosity obtained by using optimization methods can be potentially used in two ways for flows whose configurations are similar to that from which data is available:

(1) to improve predictions and reduce uncertainties, or (2) to quantify uncertainties by building statistical models for the discrepancies in the RANS-modeled eddy viscosity.

The first approach would involve building a functional mapping from the mean flow field to the eddy viscosity or its discrepancies. However, as of the writing of this review, the authors are not aware of any published research pursuing this approach. A machine-learning based approach to predict discrepancies of RANS-modeled Reynolds stresses has been investigated [START_REF] Wang | Physics-informed machine learning approach for reconstructing Reynolds stress modeling discrepancies based on DNS data[END_REF][START_REF] Wu | Physics-informed machine learning approach for augmenting turbulence models: a comprehensive framework[END_REF](see Section 4.4), and one can envision a similar approach to be used on the eddy viscosity. On the other hand, the second approach has been pursued by Ref. [START_REF] Dow | Quantification of structural uncertainties in the k-ω turbulence model[END_REF], which is detailed below. Dow and Wang [START_REF] Dow | Quantification of structural uncertainties in the k-ω turbulence model[END_REF] first used DNS data from plane channel flows to infer an optimal eddy viscosity field t opt . They further constructed a zero-mean Gaussian process for the logarithmic discrepancy log( / )

t t log rans = .
Equivalently, the field of true eddy viscosity x ( ) t is modeled as random field as follows:

x x K log log with (0, ( , )) t t rans log log = + GP ( 33
)
where the covariance kernel K was chosen as a squared exponential function with its hyperparameters including variance σ and length scale l (see Equation ( 6)). These hyperparameters were determined by using maximum likelihood estimation with the optimal eddy viscosity field used as data. After the hyperparameters were determined, they sampled the Gaussian processes to obtain realizations of possible eddy viscosity fields (Fig. 19a) in similar yet slightly different geometries, e.g., plane channel wavy walls. This slight extrapolation is based on the assumption that the eddy viscosity discrepancies x ( ) log in a class of similar flows conform to the same statistical distribution. Such realizations of the eddy viscosity obtained from the Gaussian process were used to solve the RANS equations to obtain an ensemble of velocity predictions as shown in Fig. 19b. The obtained ensemble represents the uncertainties in RANS-predicted velocities, which can be further processed to obtain uncertainties for other quantities of interests. This methodology has recently been extended to more complex flows in a Ubend channel [START_REF] Hayek | Adjoint-based optimization of RANS eddy viscosity model for U-bend channel flow[END_REF]. As with the UQ approach based on transport equations [START_REF] Singh | Using field inversion to quantify functional errors in turbulence closures[END_REF], all the predictions in the ensemble are still constrained by the Boussinesq assumption originating from the baseline model.

Introducing uncertainties in Reynolds stresses

Reynolds stress plays a unique and particularly important role in RANS modeling -it is the term through which most turbulence models enter the RANS momentum equations, as can be seen from Fig. 3. In the derivation of the RANS equations there is a closure problem. That is, the Reynolds stress term in the obtained averaged-equations needs to be modeled. While a Reynolds stress transport equation (RSTE) can be derived from the NS equations, the RSTE itself contains even more unclosed terms. On the other hand, if the true Reynolds stress field is supplied to the RANS equation, theoretically the true velocity and all other quantities can also be obtained, provided that the numerical uncertainties are negligible and that the RANS equations are wellconditioned. As the Reynolds stress is the only modeled term in the RANS equations, inaccuracy in its modeling is the source of model-form uncertainty for RANS simulations, at least for single phase, fully turbulent flows without transition [START_REF] Pope | Turbulent Flows[END_REF]. Transition modeling is an important topic [START_REF] Zhang | An efficient Bayesian uncertainty quantification approach with application to k-ω-γ transition modeling[END_REF] but it is beyond the scope of this review. The insight on the importance of Reynolds stress in turbulence modeling was abstracted as composite model theory [START_REF] Oliver | Validating predictions of unobserved quantities[END_REF], which is detailed in Appendix B.

Based on the observations above, it is natural to introduce uncertainties to the Reynolds stresses. So far, two distinct approaches have been proposed to characterize the uncertainties in the Reynolds stresses:

• formulating a stochastic differential equation (SDE) for the Reynolds stress discrepancy tensor driven by a Wiener process (random walk forcing model) [START_REF] Oliver | Representing Uncertainty Due to Inaccurate Models[END_REF], and • using realizability constraints to guide the perturbations of single- point Reynolds stresses.

In both approaches the Reynolds stress discrepancy is considered a random tensor field characterized by physical constraints (e.g., conservation laws or realizability). Both approaches are introduced below.

Stochastic differential equation of Reynolds stress discrepancy

In the first approach, several forms of SDEs were explored for the Reynolds discrepancy field in a plane channel flow, a typical one of which reads as follows [START_REF] Oliver | Representing Uncertainty Due to Inaccurate Models[END_REF]:

C dU dx d dx dx dx dW dx ( ) ( ) pr 2 2 2 5/4 2 7/4 2 + = + (34)
with the three terms indicating production, diffusion, and residual, respectively; x 2 is the wall-normal coordinate and U is the horizontal velocity; W indicates a Wiener process; C pr , C , and C are coefficients to be calibrated from data. The SDE has a form that is similar to, but simpler than, the Reynolds stress transport equations. Specifically, the SDE shares the same convection-diffusion-production form 4 as the RSTE, but the SDE has a stochastic residual term on the right hand side of the SDE in place of the unclosed terms (e.g., triple velocity correlation and pressure-rate-of-strain) in the exact RSTE. The solution to the SDE provides an indication of the uncertainties in the Reynolds stresses, which can be propagated to the velocities and other quantities of interests. The SDE-based approach yields uncertainties for the entire field x ( ), which is in contrast to the single-point realizability constraints examined in Section 4.4.2. If one considers the discrepancy x ( a tensorial random field, the cross-component and spatial correlations are both accounted for through the SDE. Unfortunately, the construction of the SDE heavily relies on physical insights and modeling heuristics. Consequently, it is not straightforward to extend the formulation above to more complex flows beyond plane channel flows.

Estimating uncertainty bounds guided by realizability maps

In the second approach, perturbations are introduced directly to the modeled Reynolds stresses, based on which uncertainty propagation and statistical inferences are performed. A common starting point of these methods is the following decomposition of the Reynolds stress tensor:

k k I a I V V 2 1 3 2 1 3 = + = + ( 35
)
where k is the turbulent kinetic energy, which indicates the magnitude of ; I is the second-order identity tensor; a is the anisotropy tensor;

V v v v [ , , ] 1 2 3 = and diag[ , , ] 1 2 3
= where

1 2 3
+ + = are the orthonormal eigenvectors and eigenvalues of a, respectively, indicating the shape (aspect ratio) and orientation of , if the latter is visualized as an ellipsoid [START_REF] Simonsen | Turbulent stress invariant analysis: clarification of existing terminology[END_REF].

Transformation of the eigenvalues leads to invariants that can be mapped to the well-known Lumley triangle [START_REF] Lumley | Computational modeling of turbulent flows[END_REF] or the recently proposed barycentric triangle [START_REF] Banerjee | Presentation of anisotropy properties of turbulence, invariants versus eigenvalue approaches[END_REF], both of which provide a map for all realizable states of turbulence. Any realizable turbulence state can be mapped to a point within or on the edge of the triangles after the respective transformations. In the case of the barycentric map, the 4 The convection term disappears in mean equations of the plane channel flow.

following linear transformation from eigenvalues ( , , )

1 2
3 of the anisotropy a to the barycentric coordinates c c c ( , , )

1 2
3 is adopted:

c 1 1 2 = (36a) c 2( ) 2 2 3 = (36b) c 3 1 . 3 3 = + (36c)
The barycentric triangle and the mapping above are similar to the Lumley triangle but overcomes several shortcomings of the latter, including (i) the tendency to cluster towards the isotropic state and (ii) the nonlinearity in the mapping from the eigenvalues to tensor invariants. Like the Lumley triangle, the barycentric triangle has clear physical interpretation in that it indicates the componentality of the turbulence [START_REF] Emory | Componentality-based wall-blocking for RANS models[END_REF][START_REF] Emory | Estimating Model-form Uncertainty in Reynolds-averaged Navier-Stokes Closures[END_REF]. For example, the upper corner (c 1 3 = ) corre- sponds to three-component isotropic turbulence while the lower left corner (c 1 2 = ) corresponds to two-component axisymmetric turbu- lence, which occurs in flows close to a solid wall (e.g., point P 2 in Fig. 15).

The realizability requirements on Reynolds stresses have been studies extensively in the early years of turbulence model development. Efforts from Schumann [START_REF] Lumley | Computational modeling of turbulent flows[END_REF], Lumley [START_REF] Schumann | Realizability of Reynolds-stress turbulence models[END_REF], and Pope [START_REF] Pope | PDF methods for turbulent reactive flows[END_REF], among others, have led to a class of realizable Reynolds stress models [START_REF] Speziale | On the realizability of Reynolds stress turbulence closures[END_REF]. However, in the context of quantifying model-form uncertainties in RANS simulations, Emory et al. [START_REF] Emory | Modeling structural uncertainties in Reynoldsaveraged computations of shock/boundary layer interactions[END_REF][START_REF] Emory | Modeling of structural uncertainties in Reynolds-averaged Navier-Stokes closures[END_REF] pioneered the use of realizability maps to guide the exploration of Reynolds stress uncertainties. They proposed introducing separate perturbations to k, , and V re- sulting from the decomposition above to obtain a few representative limiting states:

k I V V 2 1 3 , = + (37)
where indicates perturbed states from the RANS-predicted baseline, e.g., rans = + . The initial focus was placed on the eigenvalues, as the realizability map provides a straightforward and rigorous bound on how they can be perturbed. One possibility of perturbing the anisotropy is to perturb it towards one-component (1C), two-component (2C), and three-component (3C) limiting states of realizable turbulence, represented by the three corresponding vertices of the barycentric triangle (see Fig. 20a).

Nevertheless, the realizability map does not provide a direct bound on the magnitude k and the eigenvectors V. In order to utilize the realizability map to bound k and V, it is important to recognize that k, , and V are not independent but intimately coupled. They are dif- ferent characteristics of the same Reynolds stress tensor, which is governed by a coupled Reynolds stress transport equation (RSTE). In fact, with some algebra the RSTE can be transformed to three individual transport equations for the turbulent kinetic energy(TKE) k, eigenvalues , and eigenvectors V as well as their discrepancies [START_REF] Pope | Turbulent Flows[END_REF][START_REF] Thompson | A strategy for the eigenvector perturbations of the Reynolds stress tensor in the context of uncertainty quantification[END_REF], although only the TKE transport equation (10b) is commonly used in turbulence modeling. The coupling among the three variables can be utilized in many ways. For example, the anisotropy bounds obtained from the realizability map [START_REF] Emory | Modeling of structural uncertainties in Reynolds-averaged Navier-Stokes closures[END_REF] can be used to estimate the bounds on the TKE production S : k = P , which is further substituted into transport equation (10b) to obtain the TKE corresponding to the limiting states [START_REF] Gorlé | Epistemic uncertainty quantification for RANS modeling of the flow over a wavy wall[END_REF]. The obtained TKE fields can be used to estimate their uncertainties. Similarly, Thompson et al. [136] exploited the coupling to estimate the uncertainties in the eigenvectors by using the realizability bounds of the eigenvalues. Unfortunately, a large number of unclosed terms in the Reynolds stress transport equation makes it much more difficult than estimating the uncertainties in the TKE. Recently, the eigenvectors perturbation has also been investigated [START_REF] Iaccarino | Eigenspace perturbations for uncertainty estimation of single-point turbulence closures[END_REF][START_REF] Mishra | Uncertainty estimation for Reynolds-Averaged Navier-Stokes predictions of high-speed aircraft nozzle jets[END_REF]. Two extreme bounding cases of the perturbation for Reynolds stress eigenvectors are considered. In one case the semi-major axis of the Reynolds stress ellipsoid is aligned with the stretching eigen-direction of the mean rate of strain tensor; in another case the semi-major axis of the Reynolds stress ellipsoid is aligned with the compressive eigen-direction of the mean rate-of-strain tensor. The two limiting states are chosen to explore the possible extreme scenarios of turbulent production.

In summary, the above-mentioned studies by Iaccarino and coworkers [START_REF] Emory | Modeling structural uncertainties in Reynoldsaveraged computations of shock/boundary layer interactions[END_REF][START_REF] Emory | Modeling of structural uncertainties in Reynolds-averaged Navier-Stokes closures[END_REF][START_REF] Iaccarino | Eigenspace perturbations for uncertainty estimation of single-point turbulence closures[END_REF][START_REF] Thompson | A strategy for the eigenvector perturbations of the Reynolds stress tensor in the context of uncertainty quantification[END_REF][START_REF] Gorlé | Epistemic uncertainty quantification for RANS modeling of the flow over a wavy wall[END_REF] used barycentric triangle as guide to comprehensively explore the limiting states of Reynolds stresses. They form an efficient, physics-based scheme to estimate RANS model uncertainty by using only five simulations. Moreover, the parameterization scheme of Reynolds stress perturbations becomes the foundation of more sophisticated methods that use statistical inference and machine learning to quantify and reduce the RANS model uncertainties [START_REF] Xiao | Quantifying and reducing model-form uncertainties in Reynolds-averaged Navier-Stokes simulations: a data-driven, physics-informed Bayesian approach[END_REF][START_REF] Wang | Physics-informed machine learning approach for reconstructing Reynolds stress modeling discrepancies based on DNS data[END_REF][START_REF] Wu | Physics-informed machine learning approach for augmenting turbulence models: a comprehensive framework[END_REF].

Systematic sampling of Reynolds stress uncertainty within realizability constraints

In statistical inference for quantifying and reducing model-form uncertainties, it is insufficient to merely perturb Reynolds stresses towards limiting states [START_REF] Xiao | Quantifying and reducing model-form uncertainties in Reynolds-averaged Navier-Stokes simulations: a data-driven, physics-informed Bayesian approach[END_REF]. Rather, statistical sampling and inference requires a systematic scheme for parameterizing the perturbations to the TKE, the eigenvalues, and particularly the eigenvectors. Perturbations on k and can be represented as random fields, albeit with realizability constraints. To this end, the uncertainties associated with RANS modeled TKE k rans can be represented in terms of its logarithmic discrepancy, in a similar way to that for the eddy viscosity [START_REF] Wang | Quantification of uncertainties in turbulence modeling: a comparison of physics-based and random matrix theoretic approaches[END_REF] in Equation [START_REF] Klir | Uncertainty and Information: Foundations of Generalized Information Theory[END_REF]. However, parameterizing the perturbations on the eigenvectors is more challenging due to the need to maintain their orthonormal property, which is necessary to ensure that the perturbed Reynolds stresses remain symmetric positive semidefinite tensors. To this end, it is most convenient to represent the perturbation from V to V as a rigid-body rotation, i.e., Q V V = Q being an ortho- normal rotation matrix representing the perturbation. In fact, a rotation can be represented more compactly by using a set of Euler angles ( , 1

3 ). That is, any rigid-body rotation in a three-dimensional space (with a few rare exceptions) can be achieved by the following three consecutive intrinsic rotations about the axes of the local coordinate system (x-y-z) of the rigid body [START_REF] Goldstein | The Euler angles[END_REF]: (i) a rotation about the z axis by angle 1 , (ii) a rotation about the x axis by 2 , followed by (iii) another rotation about its z axis by 3 . The Euler-angle based representation has been used for quantifying RANS model-form uncertainties [START_REF] Wang | Quantification of uncertainties in turbulence modeling: a comparison of physics-based and random matrix theoretic approaches[END_REF]. Alternatively, the same transformation can be represented as a unit quaternion. Euler's rotation theorem states there exists a unique axis of unit vector n n n n [ , , ] 1 2

3 passing through the origin and an angle ϑ such that V can be obtained via rotating V by ϑ about an axis n, and thus the rigid-body rotation can be represented by a unit quaternion [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF]:

n n n h cos 2 , sin 2 , sin 2 , sin 2
1 2 3 = ( 38
)
where h = 1. In uncertainty quantification and machine learning for RANS modeling, the two representations of Reynolds stress perturbation based on Euler angle and unit quaternion have been compared, and the latter was found to be superior [START_REF] Wu | Representation of Reynolds stress perturbations with application in machine-learning-assisted turbulence modeling[END_REF].

Random matrix approach for quantifying Reynolds stress uncertainty

The realizability constraint of Reynolds stresses plays a critical role in all the RANS model-form uncertainty quantification methods outlined above. However, physics-based decomposition as in Equation ( 35) is only one of the possible ways to guarantee realizability. Xiao et al. [START_REF] Xiao | A random matrix approach for quantifying model-form uncertainties in turbulence modeling[END_REF] proposed an alternative approach where the Reynolds stress tensor is modeled as a 3 3 × random matrix that conforms to a maximum entropy distribution defined on the set of positive semi-definite matrices. Reynolds stress uncertainty can thus be estimated by directly sampling from the defined distribution, with the realizability of all samples guaranteed without using the realizability maps. The validity of the random matrix approach can be clearly seen from the equivalence among the following three interpretations of the Reynolds stress realizability. That is, a Reynolds stress tensor is physically realizable if and only if it satisfies one of the following conditions:

(1) it is the covariance matrix of a real-valued vector (i.e., the velocity), (2) it resides within or on the edge of the barycentric triangle (or Lumley triangle) after transformations (e.g., Equation [START_REF] Liu | Arithmetic with uncertain numbers: rigorous and (often) best possible answers[END_REF] for the former), or (3) it is symmetric and positive semi-definite.

The three conditions are, respectively, (i) the origin of the realizability constraint, (ii) the foundation for the physics-based approach, and (iii) the basis for the random matrix approach. All three conditions above are equivalent [START_REF] Xiao | A random matrix approach for quantifying model-form uncertainties in turbulence modeling[END_REF].

The random approach and the physics-based approach are compared in detail in Wang et al. [START_REF] Wang | Quantification of uncertainties in turbulence modeling: a comparison of physics-based and random matrix theoretic approaches[END_REF]. It was concluded that both approaches yield qualitatively similar results, particularly when the perturbations are small and far away from the limiting states (edges of the barycentric triangle). When the perturbations are large, physicsbased perturbations may result in samples falling outside the barycentric triangle, which must be capped and result in a slightly distorted distribution. In contrast, the random matrix approach does not have this issue. Another important difference is that the physics-based approach perturbs the three components (magnitude, shape, and orientation) separately, while the random matrix approach perturbs all three components simultaneously, with k, , and V implicitly con- strained by the maximum entropy principle.

Quantifying and reducing Reynolds stress uncertainties with data

The works reviewed above all involved forward analysis, i.e., propagation of uncertainties introduced in the Reynolds stresses to velocities or derived quantities of interest. As with the parametric approaches and other non-parametric approaches introducing uncertainties in viscosity (Section 4.3) and turbulent transport equations (Section 4.2), backward analysis (statistical inference) can also be performed on Reynolds stresses to quantify and reduce uncertainties in RANS model predictions. The objective is to find a Reynolds stress field that yields the best agreement with the data (e.g., sparse observations of velocities) accounting for the state covariance and the error covariance. Both variational methods and ensemble methods introduced in Section 2.3 can be used. Xiao et al. [START_REF] Xiao | Quantifying and reducing model-form uncertainties in Reynolds-averaged Navier-Stokes simulations: a data-driven, physics-informed Bayesian approach[END_REF] used the ensemble Kalman method to infer the Reynolds stress and full-field velocities from sparse velocity data. Constraints and empirical prior knowledge about the Reynolds stress field and its discrepancies are built into the inference, specifically including:

(1) realizability at any point, (2) smoothness of the Reynolds stress field and its discrepancy for incompressible flows, and (3) empirical knowledge on the regions where Reynolds stress discrepancies are large. Utilizing these constraints and prior knowledge greatly reduces the dimension of the inverse problem, which has the same effects as the regularization terms in the full-field inversion and optimization problems as in Equations ( 29) and [START_REF] James | An Introduction to Statistical Learning[END_REF]. The realizability is ensured by parameterizing the Reynolds stress in terms of the physics-based decomposition as in Equation [START_REF] Ferson | Different methods are needed to propagate ignorance and variability[END_REF]. The smoothness is ensured by representing the random fields in terms of its Karhunen-Loeve expansion, truncated to the first n terms:

x x w ( ) ( ) n 1 = = ( 39
)
where w are random variables to be inferred, and x { ( )} n 1 = are a set of orthogonal basis functions corresponding to the covariance kernel of the random field . The basis functions can be computed from the Fredholm integral equation by solving an eigenvalue problem for the kernel and embody the empirical knowledge on Reynolds stress discrepancy as mentioned above. Fig. 21 shows representative results from a fully developed square duct flow, presenting the in-plane velocities before and after incorporating the data, i.e., prior and posterior of velocity distributions. It clearly shows that the velocity predictions at all cross-sections are markedly improved, even in locations where velocity observations are not available. The calibrated discrepancy can even be used to correct square duct flows at a higher Reynolds number or flows in different geometries such as a rectangular duct [START_REF] Wu | A Bayesian calibration-prediction method for reducing model-form uncertainties with application in RANS simulations[END_REF]. Furthermore, it was demonstrated that incorporating empirical knowledge is clearly valuable and has similar effects as increasing the amount of observation data [START_REF] Wang | Incorporating prior knowledge for quantifying and reducing model-form uncertainty in RANS simulations[END_REF].

Among the prior knowledge used for the statistical inference, the regions where Reynolds stress discrepancies are large have primarily depended on input from users based on their empirical knowledge. However, the combination of physical and modeling insights with modern data science has opened new opportunities. Gorlé et al. [START_REF] Gorlé | The deviation from parallel shear flow as an indicator of linear eddy-viscosity model inaccuracy[END_REF] proposed an analytical marker function based on the deviation from parallel shear flow and used it to predict discrepancies in RANS-modeled Reynolds stress. Their ideas are based on the insightful observation that commonly used eddy viscosity models were developed and tuned for parallel shear flows (boundary layers). A departure from such flows typically leads to violations of assumptions in these models. Moreover, emerging machine learning techniques have made it possible to provide more accurate maps of where large discrepancies exist. Ling and Templeton [START_REF] Ling | Evaluation of machine learning algorithms for prediction of regions of high Reynolds averaged Navier Stokes uncertainty[END_REF] developed a machine learning method to evaluate potential inadequacy of RANS models by using DNS databases. This approach has been recently applied to more complex flows (e.g., jet-incrossflow [START_REF] Ling | Uncertainty analysis and data-driven model advances for a jet-in-crossflow[END_REF]). The results include several fields of binary labels (whether the specified model assumption is violated), which could be further processed to obtain a variance of Reynolds stress discrepancy that can be incorporated into the covariance kernel field.

The same decomposition scheme in Equation ( 35), which has been used for Reynolds-stress-based uncertainty quantification and statistical inferences, can be used as a parameterization scheme for correcting RANS-predicted Reynolds stresses by using machine learning and training data. Tracey et al. [START_REF] Tracey | Application of supervised learning to quantify uncertainties in turbulence and combustion modeling[END_REF] represented discrepancies in barycentric coordinates as a function of local mean-flow variables and leveraged machine learning to train the function. Wang et al. [START_REF] Wang | Physics-informed machine learning approach for reconstructing Reynolds stress modeling discrepancies based on DNS data[END_REF] and Wu et al. [START_REF] Wu | Physics-informed machine learning approach for augmenting turbulence models: a comprehensive framework[END_REF] developed a more systematic strategy to predict discrepancies in the magnitude, anisotropy, and orientation of the Reynolds stress tensor in terms of an invariant feature set for a set of tensor variables of the mean flow (e.g., S, , p, k), referred to as the in- tegrity basis [START_REF] Ling | Machine learning strategies for systems with invariance properties[END_REF]. They showed improved results in the prediction of Reynolds stresses and mean velocities [START_REF] Wu | Physics-informed machine learning approach for augmenting turbulence models: a comprehensive framework[END_REF] for two canonical flows, i.e., separated flows over periodic hills and secondary flows in a square duct. As an alternative, Ling et al. [START_REF] Ling | Reynolds averaged turbulence modelling using deep neural networks with embedded invariance[END_REF] proposed a neural network architecture with embedded invariance properties to learn and predict the coefficients of an objective formulation for the tensorial function f S ( , ) = . These works illustrated how physical constraints can be embedded in machine learning. Nevertheless, a unique challenge for directly correcting or predicting the Reynolds stress tensors with datadriven models is the possible ill-conditioning of the RANS equations. For example, small errors in the machine-learning-predicted Reynolds stresses can lead to large errors in the propagated velocities [START_REF] Wu | RANS equations with data-driven Reynolds stress closure can be ill-conditioned[END_REF]. In order to overcome this difficulty, Wu et al. [START_REF] Wu | Physics-informed machine learning approach for augmenting turbulence models: a comprehensive framework[END_REF] proposed a scheme of learning the linear and nonlinear parts of the Reynolds stress separately, with the linear part treated implicitly to improve model conditioning. Numerous other approaches have been proposed for augmenting and improving turbulence models based on machine learning [START_REF] Duraisamy | Turbulence modeling in the age of data[END_REF]. Finally, data-driven, machine-learning based methods have also been used in improving CFD models of thermal fluids flow with focus on boiling flows in nuclear reactor thermo hydraulics [e.g., [START_REF] Liu | Data-driven modeling for boiling heat transfer: using deep neural networks and high-fidelity simulation results[END_REF][START_REF] Chang | Reynolds-averaged turbulence modeling using type I and type II machine learning frameworks with deep learning[END_REF][START_REF] Chang | Classification of machine learning frameworks for datadriven thermal fluid models[END_REF][START_REF] Hanna | Coarse-grid computational fluid dynamic (CG-CFD) error prediction using machine learning[END_REF] and in high-Mach number flows [START_REF] Wang | Predicting Reynolds stresses in high-Machnumber hurbulent boundary layers with physics-informed machine learning[END_REF].

Spatial correlations in Reynolds stress discrepancy

Most of the approaches for RANS uncertainty quantification reviewed above have focused on the uncertainty bounds of the Reynolds stress at a single point. This is partly because the mathematical rigor of such a bound can only be maintained by the realizability constraint. However, an equally important source of uncertainty comes from the spatial variation of the Reynolds stress discrepancy. After all, it is the divergence of the Reynolds stress field that appears in the RANS momentum equations. When quantifying uncertainties in RANS-based predictions, Emory et al. [START_REF] Emory | Modeling of structural uncertainties in Reynolds-averaged Navier-Stokes closures[END_REF] specified a spatial field for the eigenvalue perturbations based on the empirical understanding on the performance of the RANS model for the particular problem. Xiao et al. [START_REF] Xiao | Quantifying and reducing model-form uncertainties in Reynolds-averaged Navier-Stokes simulations: a data-driven, physics-informed Bayesian approach[END_REF][START_REF] Xiao | A random matrix approach for quantifying model-form uncertainties in turbulence modeling[END_REF] used the same argument to define a non-stationary Gaussian process model for the Reynolds stress discrepancies. They used Karhunen-Loeve (KL) expansions [START_REF] Le Maître | Spectral Methods for Uncertainty Quantification: With Applications to Computational Fluid Dynamics[END_REF] to approximately represent the perturbation field with leading modes in KL expansion (see Equation ( 39)). Specification of such spatial distribution is probably the weakest link in the entire process of Reynolds stress-based uncertainty quantification.

Since the true Reynolds stress and the RANS modeled counterpart are described by their respective transport equations, the model discrepancies should conform to a transport equation of the same structure. Building upon such insight, Edeling et al. [START_REF] Edeling | Data-free and data-driven RANS predictions with quantified uncertainty[END_REF] proposed a "return-to-eddy-viscosity" model, which is a transport equation with a source term describing the deviation of the turbulence state from equilibrium state assumed by the linear eddy viscosity models:

Dc Dt a k c c x c x ( ) c i t c i 1 1 1 rans 1 1 1 = + + (40)
where c 1 rans is the barycentric coordinates corresponding to the baseline RANS modeled Reynolds stress tensor; a c 1 and c 1 are model coefficients to be calibrated. A similar PDE is formulated for c 2 . These heuristically justified, physics-inspired PDEs provide a bound for the Reynolds stress field. Moreover, the model coefficients in the PDEs above can be calibrated by using data and Bayesian inference, and the calibrated equations are further used for predictions [START_REF] Edeling | Data-free and data-driven RANS predictions with quantified uncertainty[END_REF].

More recently, Wu et al. [START_REF] Wu | Constructing PDE-Informed covariance kernel for quantification of model uncertainties[END_REF] utilized the fundamental connection between PDEs and covariance to provide a physically anchored covariance structure, which has a clear advantages over purely statistical covariance structures previously used for model discrepancies [START_REF] Xiao | Quantifying and reducing model-form uncertainties in Reynolds-averaged Navier-Stokes simulations: a data-driven, physics-informed Bayesian approach[END_REF][START_REF] Xiao | A random matrix approach for quantifying model-form uncertainties in turbulence modeling[END_REF]. Specifically, they constructed an approximate, linearized PDE for the model discrepancy:

D Dt x x or more compactly ( ) i T i + = = S L S ( 41
)
where δ denotes the field of model discrepancy such as the discrepancy in the RANS-modeled Reynolds stresses, and S indicates the unclosed source terms. Equation ( 41) can be generalized as ( ) = L S , where L corresponds to the linearized differential operator on the left-hand side. In previous works [START_REF] Xiao | Quantifying and reducing model-form uncertainties in Reynolds-averaged Navier-Stokes simulations: a data-driven, physics-informed Bayesian approach[END_REF][START_REF] Xiao | A random matrix approach for quantifying model-form uncertainties in turbulence modeling[END_REF], purely statistical covariance structures such as the squared exponential kernel in Gaussian processes were specified for the model discrepancy δ. The physics-inspired transport equation [START_REF] Roy | A comprehensive framework for verification, validation, and uncertainty quantification in scientific computing[END_REF] requires the specification of a kernel of the source term S and provides a physical covariance structure of the error term δ by transforming the covariance with the differential operator as follows:

( ) , S 1 1 = L L ( 42
)
where 1 L is the inverse operator of the linearized PDE [START_REF] Roy | A comprehensive framework for verification, validation, and uncertainty quantification in scientific computing[END_REF], and denotes covariance. Wu et al. [START_REF] Wu | Constructing PDE-Informed covariance kernel for quantification of model uncertainties[END_REF] showed that such a physics-inspired covariance structure better accounts for the spatial correlation of the discrepancy term δ than the squared exponential covariance kernel. Sample results for flow over periodic hills are presented in Fig. 22, which shows the first three modes (i.e., x { ( )} 1 3 = as in Equation ( 39)) obtained by using a squared exponential kernel (Fig. 22a) and a physics-informed kernel (Fig. 22b), e.g., from Equation [START_REF] Dow | Quantification of structural uncertainties in the k-ω turbulence model[END_REF]. In this geometry, the general flow direction is from left to right. The streamlinealigned covariance structure endowed by the convection is evident, while the modes obtained from the squared exponential kernel exhibits nonphysical, spatially isotropic structures.

Uncertainties in large eddy simulations

Large eddy simulation (LES) is a turbulence simulation method that resolves larger scale turbulence and models sub-grid scales [START_REF] Sagaut | Large Eddy Simulations for Incompressible Flows: An Introduction[END_REF]. For LES performed on an adequate mesh, most of the important turbulence scales are resolved except in near-wall regions of wall-bounded flows. Consequently, uncertainties associated with the subgrid scale (SGS) model no longer dominate. Instead, LES are influenced by uncertainties of a number of sources that are of comparable order of magnitude, including:

(1) uncertainties due to SGS models, including their parameters, (2) uncertainties associated with initial and boundary conditions, (3) uncertainties in the numerical discretization (mesh and numerical scheme).

Note that items 2-3 are not model uncertainties but input uncertainties and numerical uncertainties, respectively. This is in stark contrast to RANS simulations, where the model uncertainty clearly dominates other uncertainties. The literature on prediction accuracies of LES is vast and is mostly from deterministic perspectives. In such frameworks, the problem of concern should be more precisely referred to as errors and not uncertainties. Nevertheless, some studies have tackled the problem from an uncertainty quantification point of view, which are shortly reviewed here. Due to the much higher computational cost of LES as compared to RANS simulations, most studies are limited to uncertainty propagation and sensitivity analysis, i.e., propagation of assumed probability distributions through an LES solver to investigate the sensitivity of the output quantities with respect to the input. In order to reduce the number of samples and overall computational costs for uncertainty propagation, many of the studies reviewed below built surrogate models by using different methods, e.g., polynomial chaos expansion [START_REF] Le Maître | Spectral Methods for Uncertainty Quantification: With Applications to Computational Fluid Dynamics[END_REF][START_REF] Najm | Uncertainty quantification and polynomial chaos techniques in computational fluid dynamics[END_REF], probabilistic collocation method [START_REF] Tatang | Direct Incorporation of Uncertainty in Chemical and Environmental Engineering Systems[END_REF], or sparse grid method [START_REF] Bungartz | Sparse grids[END_REF].

A review of recent work about quantification and reduction of uncertainties arising in LES is presented in the following of this section. We point out here that high-fidelity simulations such as LES and DNS generally have smaller uncertainties than RANS simulations. However, even DNS have their own uncertainties, e.g., due to statistical averaging and numerical methods, which must be considered in many situations when using DNS data for RANS model development and calibration. While an in-depth discussion of DNS uncertainties is beyond the scope of the current review, in Appendix C we survey a few aspects that are most relevant for RANS modeling.

Uncertainties in SGS models

Traditional LES computations rely on the explicit introduction of a closure model for the subgrid-scale terms arising from the filtering of the Navier-Stokes equations. A large number of SGS models have been developed over years, almost all of which require specification of model constants, although some (e.g. dynamic Smagorinsky model) allow for a dynamic computation of the parameters from a test filter [START_REF] Germano | A dynamic subgridscale eddy viscosity model[END_REF]. The most widely used SGS model is the algebraic Smagorinsky model, which models the SGS viscosity as:

C S S S S ( ) | ¯| with | ¯| 2 ¯¯, s g ij ij sgs 2 = ( 43
)
where S ¯is the rate-of-strain based on the filtered velocity field, g is the grid size, and C s is a coefficient that needs to be specified (referred to as Smagorinsky constant, usually chosen to be in the range from 0.1 to 0.2) and has the effect of determining the strength of SGS dissipation. Meyers and Sagaut [START_REF] Meyers | On the model coefficients for the standard and the variational multi-scale Smagorinsky model[END_REF] derived the exact expression of the Smagorinsky constant C s by using Pope's formulation for the turbulent kinetic energy spectrum. The derivation demonstrate that C s depends both on the specific flow and on the filter, indicating that it should be treated as an uncertain quantity. The dependence of C s on the filter size deserves special attention for LES with implicit filtering, where the filter size is not explicitly specified but determined by the local grid size (see further discussions on numerical uncertainties in Section 5.3). Lucor et al. [START_REF] Lucor | Sensitivity analysis of large-eddy simulations to subgrid-scale-model parametric uncertainty using polynomial chaos[END_REF] performed LES for decaying homogeneous isotropic turbulence and propagated the uncertainties associated with the Smagorinsky constant. Specifically, they reconstructed accurately the solution statistics with a typical number of 22 samples. They carried out uncertainty propagation corresponding to different grid resolutions and found that an optimal value of the constant can be found for each level of grid refinement [START_REF] Meyers | On the model coefficients for the standard and the variational multi-scale Smagorinsky model[END_REF]. This finding confirmed the close interactions between the SGS model and the numerical discretization.

Meldi et al. [START_REF] Meldi | Is the Smagorinsky coefficient sensitive to uncertainty in the form of the energy spectrum?[END_REF] investigated the sensitivity of the C s constant to the algebraic function and its parameters used to describe the initial energy spectrum. Khalil et al. [START_REF] Khalil | Uncertainty quantification in LES of a turbulent bluff-body stabilized flame[END_REF] performed LES for turbulent bluffbody stabilized flame and studied the uncertainties associated with Smagorinsky constant, Prandtl number, and Schmidt number. Safta et al. [START_REF] Safta | Uncertainty quantification in LES of channel flow[END_REF] investigated LES of channel flow and studied uncertainties associated with model coefficients C µ and C in the k sgs model [START_REF] Yoshizawa | Bridging between eddy-viscosity-type and second-order turbulence models through a two-scale turbulence theory[END_REF]. Unlike previous uncertainty propagation studies for LES, Templeton et al. [START_REF] Templeton | Calibration and Forward Uncertainty Propagation for Large-eddy Simulations of Engineering Flows[END_REF] first used Bayesian inference to calibrate the model coefficients C µ and C based on a DNS database of forced isotropic turbu- lence in a periodic box [START_REF] Li | A public turbulence database cluster and applications to study Lagrangian evolution of velocity increments in turbulence[END_REF]. The quantified uncertainties in the calibrated model coefficients were then propagated to predictions in LES of turbulent channel flows. Tran et al. [START_REF] Tran | A sparse grid method for Bayesian uncertainty quantification with application to large eddy simulation turbulence models[END_REF] also used Bayesian inference to quantify the uncertainties associated with the Smagorinsky constant, the filter length, and the exponent in van Driest damping function.

While the above-mentioned studies used parametric approaches to address uncertainties associated with model coefficients, non-parametric UQ approaches for LES recently started drawing attention. Jofre et al. [START_REF] Jofre | A framework for characterizing structural uncertainty in large-eddy simulation closures[END_REF] estimated the structural uncertainties in the SGS stress model in LES of the canonical plane channel flows. They perturbed the SGS stresses obtained from baseline model in a similar way as in the RANS simulations [START_REF] Emory | Modeling structural uncertainties in Reynoldsaveraged computations of shock/boundary layer interactions[END_REF][START_REF] Emory | Modeling of structural uncertainties in Reynolds-averaged Navier-Stokes closures[END_REF][START_REF] Gorlé | A framework for epistemic uncertainty quantification of turbulent scalar flux models for Reynolds-averaged Navier-Stokes simulations[END_REF]. By directly introducing perturbation into the SGS stresses, the explored uncertainty space is no longer constrained by the baseline SGS model.

Uncertainties in the boundary conditions for LES

Boundary conditions are a crucial ingredient of the overall model in LES, as they may influence the development of shear and boundary layers and transition to turbulence. Congedo et al. [START_REF] Congedo | Numerical prediction of turbulent flows using Reynolds-averaged Navier-Stokes and large-eddy simulation with uncertain inflow conditions[END_REF] investigated the sensitivity of LES to uncertainties in the numerical inlet conditions by studying the turbulent flow in a pipe with an axisymmetric expansion. The study focused on the effect of the inlet bulk velocity, swirl ratio, and turbulent intensity on the resulting mean flow fields. The results were compared with experimental data, for which an estimate of the observational uncertainty was available. On the other hand, Carnevale et al. [START_REF] Carnevale | Uncertainty quantification: a stochastic method for heat transfer prediction using LES[END_REF] studied the heat transfer in a channel with pins and used uncertainty propagation to investigate the flow sensitivity to the Reynolds number, which is representative of the uncertainties associated with inlet mean velocity, fluid density, or geometrical variations. They compared the results to those of RANS simulations, showing that the epistemic uncertainty due to the modeling, i.e., RANS model versus LES model, dominates the aleatoric uncertainties such as the solution sensitivity to the Reynolds number. However, compared to RANS simulations, the LES results are more sensitive to the inlet Reynolds number.

Uncertainties due to the numerical discretization

As has been pointed out above, LES suffer from strong interactions between modeling and numerical errors. This is particularly true for LES with implicit filtering, which is dominant in practical LES, where the mesh is part of the model in that it provides the local filtering bandwidth, as is evident from Equation [START_REF] Singh | Using field inversion to quantify functional errors in turbulence closures[END_REF]. In theory, only the dynamics of the large scales is computed and the smaller scales are modeled. In practice, scale separation in LES is difficult to establish, since the low-pass filtering arises from a complex combination of implicit filtering by the grid and the discretization schemes. Even when explicit filters are applied, the approximations introduced by the discretization methods modify the actual shape of the filter function. The intricate interactions between SGS modeling errors and numerical errors (and the ill-defined filter resulted therefrom) have attracted attention (e.g., [START_REF]An analysis of numerical errors in large-eddy simulations of turbulence[END_REF][START_REF] Vreman | Comparision of numerical schemes in large-eddy simulation of the temporal mixing layer[END_REF][START_REF] Kravchenko | On the effect of numerical errors in large eddy simulations of turbulent flows[END_REF]]). Here we mention a few studies that analyzed the numerical parameters from a probabilistic perspective. For example, Meldi et al. [START_REF] Meldi | Quantification of errors in large-eddy simulations of a spatially evolving mixing layer using polynomial chaos[END_REF] performed LES for a spatially evolving mixing layer and studied the uncertainty propagation for grid stretching ratio in the turbulent and transitional regions. Mariotti et al. [START_REF] Mariotti | Stochastic sensitivity analysis of large-eddy simulation predictions of the flow around a 5:1 rectangular cylinder[END_REF] studied the flow around a 5:1 rectangular cylinder and propagated the uncertainties associated with grid resolution in the spanwise direction and the weight of the explicit low-pass filter.

With the increasing availability of computational resources and the increasing use of LES in industrial simulations, uncertainty quantification in LES is expected follow a similar development path as for RANS but with equal emphasis on all the above-mentioned sources of uncertainties. It will evolve from the current data-free, parametric approaches to more sophisticated, data-driven, non-parametric approaches, and from the current proof-of-concept studies to gradual deployment in industrial simulations.

Conclusions and future research

This review summarized techniques for quantifying uncertainties associated with turbulence models in computational fluid dynamics simulations. We focused on uncertainty quantification in RANS models, because they are expected to remain the workhorse tool for industrial CFD simulations in decades to come, thanks to their lower computational costs and better robustness than scale-resolving methods. Quantifying uncertainties in RANS predictions are of strategic importance towards the goal of certified numerical simulations of fluid flows.

The literature survey shows that RANS uncertainty quantification has been a rapidly evolving field in the past decade. Most of the recent research focused on statistical approaches to estimate prediction uncertainties due to turbulence models and on data-driven methods to reduce such uncertainties. Development of such statistically rigorous techniques for quantifying and reducing RANS model uncertainties has been fostered by:

(1) the considerable increase of computer resources, (2) the ever-increasing mass of high fidelity experimental and numerical data, and (3) the development of statistical sampling and inference methods guided by physical constraints and prior knowledge in turbulence modeling.

This review classifies existing literature of model uncertainty quantification into parametric and non-parametric approaches, which are discussed separately. In the parametric approaches, uncertainties are introduced into the coefficients in RANS closure models. That is, the coefficients are modeled as random variables, whose prior distributions are then propagated to the predictions through RANS simulations or updated by incorporating observation data within the Bayesian inference framework. Extensions of the parametric approaches are multimodel approaches such as Bayesian model averaging and Bayesian model-scenario averaging methods. In these methods, predictions of new flow configurations (scenarios) are formulated as an average of predictions from an ensemble of competing models, weighted by their respective posterior probabilities and the similarity of respective calibration scenarios to the prediction scenario. Parametric and multimodel ensemble methods are robust, non-intrusive, and relatively mature. When combined with surrogate models that replace RANS models to allow for efficient sampling, they can be used in uncertainty quantification involving complex, three-dimensional engineering flows.

A drawback of parametric approaches is that any calibration and inference of the parameters are inevitably based on, and will influence, the entire flow field. However, a turbulent flow may simultaneously contain regions ranging from equilibrium regions that are well predicted by simple models to highly non-equilibrium regions (e.g., separation, shock waves, streamline curvature, rotation) where even advanced models may fail. Non-parametric approaches provide an attractive alternative to tackle turbulence modeling uncertainties while accounting for locality of turbulent flows. These approaches rely on random fields to represent the RANS model discrepancy, which are estimated from physical bounds and further propagated to predictions or from observation data. However, the unique challenge here is that the uncertainty propagation and statistical inference involve random fields of much higher dimensions. Sampling and inference in such a high-dimensional space remain an active field of research with many open challenges and opportunities.

Another thrilling subject for future research is the application of non-parametric approaches for predictions. Data assimilation and machine learning algorithms have been recently applied to extrapolating estimated discrepancy fields to configurations that are relatively close to the ones contained in the training flows [START_REF] Wu | Physics-informed machine learning approach for augmenting turbulence models: a comprehensive framework[END_REF]. However, using the estimated uncertainties to drastically different configurations remains a delicate and possibly dangerous task. Introducing sound physical constraints in the representation of the discrepancy and using physicsbased transport equations for describing its spatial correlation structure seems to be a promising approach [START_REF] Wu | Constructing PDE-Informed covariance kernel for quantification of model uncertainties[END_REF]. Bridging uncertainty quantification and data-driven modeling, such a physics-informed approach has the potential of yielding RANS models that can predict turbulent flows with quantified uncertainties, paving the way toward certified CFD simulations [START_REF] Duraisamy | Turbulence modeling in the age of data[END_REF].

  Navier-Stokes (NS) equations. While many applications in aerospace engineering involve compressible flows, reacting flows, or two-phase flows, for illustration purposes we restrict our attention to the NS equations for incompressible flows of constant-property, Newtonian fluids, shown below:

  ) is normalized with respect to a reference length L ref , a reference velocity U ref , and the density ρ and viscosity μ of the fluid. The parameter Re U L µ / ref ref = is the Reynolds number, a measure of the relative importance of inertia to viscous forces. Because of the nonlinearity of the convection terms u u x (

Fig. 1 .

 1 Fig. 1. A schematic representation of the hierarchy of turbulence modeling approaches based on computational costs and the amounts of resolved versus modeled physics. Figure inspired by Sagaut et al. [1]. Abbreviations: DNS, direct numerical simulations; LES, large eddy simulations; RANS, Reynoldsaveraged Navier-Stokes.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Stagnation of turbulence modeling in the past few decades (shaded regions), showing (a) the number of wind tunnel tests required in the design cycle of commercial aircraft in the past five decades [14] and (b) the time at which commonly used models were developed.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Examples of uncertainties in RANS predictions of pressure coefficient C p distribution on wings and airfoils due to (a) model form and (b) model coefficients. Panel (a) shows the C p profile on a CRM wing-body configuration at 4. 0 angle of attack. Results are from the 6 th AIAA CFD drag prediction workshop based on different RANS models, including k-ε model, k-ω model, SA model, SA with quadratic constitutive relation (QCR), and EARSM. The location of the presented pressure distribution is indicated by the red/solid line on the wing (see inset; showing the port half of the fuselage and the wing only). Figure reprinted with permission from Tinoco et al. [24]. Panel (b) shows the C p profile on a NACA0012 airfoil in a transonic flow with freestream Mach number 0.8 and Reynolds number 9 10 6 × , obtained from RANS simulations with the algebraic model of Baldwin and Lomax [25]. The figure shows the effect of varying C wk , one of the seven model parameters, from 0.25 to 1, adopted from an unpublished report of the second author [26]. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

  The true Reynolds stress field x ( ) and the discrepancies in the RANS-modeled Reynolds stress x ( ) rans are examples of random fields in RANS model uncertainty quantification.

  ) which states that the posterior probability p ( | ) D is proportional to the prior p ( ) and the likelihood p ( | ) D . The prior p ( ) summarizes all available knowledge about θ before observing the data D . The like- lihood function p ( | ) D represents the probability of observing the data from a process described by the model M ( ) parameterized by θ. In the context of RANS uncertainty quantification, evaluating p ( | ) D for a given realization of the model parameters θ involves running the CFD code and is thus a costly operation. Finally, p ( )

Fig. 7 .

 7 Fig. 7. Illustration of Markov chain Monte Carlo sampling of a banana-shaped posterior (shaded contour) in a two-dimensional state space. The sampled distribution is illustrated with the trace of past samples (dots) and the marginal distributions (histograms plotted on the horizontal and vertical axes). Image obtained by using the MCMC demonstration code (https://chi-feng.github.io/ mcmc-demo/) by Chi Feng of MIT.

C 2 .

 2 For instance, the RNG k-ε model uses a modified value C ˜1.68 2 = , and the k-τ model (essentially a k-ε model rewritten in terms of k/ = ) uses C 1.83 2 = [91]. Nevertheless, experimental results suggest that most data agrees with n 1.3 = , which corresponds to C 1.77 2 = [92].

Fig. 8 .

 8 Fig. 8. Flow past a backward facing step at Re 50000 h = . Sensitivity of the k-ε model to the closure coefficients. Plots of the skin friction and its sensitivities versus the longitudinal position behind the step. Figures reproduced with permission from Turgeon et al. [48].

Fig. 9 .

 9 Fig. 9. Distribution of C f for the flow along a semi-infinite flat plate with zero pressure gradient and 99% uncertainty interval. Sensitivity of the k-ε model to the von Karman constant κ. Figure reproduced with permission from Platteeuw et al. [50].

Fig. 10 .

 10 Fig. 10. Normalized probability density function (pdf/max(pdf)) of the Launder-Sharma k-ε model coefficients recovered through gPC. Figures reproduced with permission from Turgeon et al. [51].

Fig. 11 .

 11 Fig. 11. Calibration of the Spalart-Allmaras model from the flat plate flow data, showing (a) the posterior distributions and (b) scatter plots of the inferred parameters κ and C 1 by using different statistical models for the inadequacy term. Figures reproduced with permission from Cheung et al. [62].

Fig. 12 .

 12 Fig. 12. Sample posterior distributions and p-box predictions of a new flow based on 13 separate calibrations of the k-ε model. Figures reproduced with permission from Cheung et al. [66].

D

  BMSA yields the posterior distribution of as follows: , each corresponding to the propagation of the parameter posterior obtained by calibration of model I against scenario K through the new prediction scenario S ˜. The average is weighted by the corresponding posterior model probability M

Fig. 13 .

 13 Fig. 13. Posterior probabilities M ( | ) i k D of five turbulence models for 13 ca- libration datasets (boundary layers of various external pressure gradients). The set of models includes a simple algebraic model (Baldwin and Lomax [25]), one-equation and two-equation eddy viscosity models (SA model [100], k-ε model, and k-ω model [10]), and a Reynolds stress model (stress-ω model [10]). Numbers on the horizontal axis denote identification codes for datasets (flow configurations). Figure reproduced with permission from Edeling et al. [60].

  ) are presented.

Fig. 16 .

 16 Fig. 16. Comparison of orientations (as indicated by two eigenvectors v 1 and v 2 ) of Reynolds stresses tensor from RANS modeling (with SST k-ω model) and experimental measurement of [119] at freestream locations P 1 (panels a and b) and near-wall location P 2 (panels c and d). Grey arrows indicate the perturbations on the eigenvectors for exploring uncertainties in RANS-predicted Reynolds stresses, which is a non-parametric approach. The third eigenvector v 3 of and S can be uniquely determined from v v v 3 1 2 = × and are thus omitted, where × indicates cross product between two vectors. Figures reproduced from Wu et al. [118] (unpublished manuscript).

Fig. 17 .

 17 Fig. 17. Comparison of parametric and nonparametric approaches for model-form uncertainty by using RANS simulations (with SST k-ω model) of a wing-body junction flow as an example. This figure compares the mean velocities at three locations (shown in Fig. 15) in front of the leading edge of a wing-body juncture obtained by (a) perturbing turbulent kinetic energy only and (b) perturbing the full Reynolds stress, corresponding to parametric and nonparametric approaches, respectively. Figures reproduced from Wu et al. [118] (unpublished manuscript).

Fig. 18 .Fig. 19 .

 1819 Fig. 18. Profiles of the inferred correction function β in plane channel flows at Reynolds numbers Re ranging from 180 to 4200. Figure reproduced with permission from Singh and Duraisamy [43].

Fig. 20 .

 20 Fig. 20. Model-form uncertainty cation through perturbation of Reynolds stresses within the physically realizable limit enclosed by the Barycentric triangle. (a) Comparison between the perturbation schemes of Emory et al. [54] and Xiao et al. [44]. (b) Perturbation scheme based on random matrix theory [58], which is compared with the physics-based perturbations Xiao et al. [44] shown in (a). Legend in Panel (b): baseline RANS prediction ; DNS data ; sample mean ; perturbed states . Figure reproduced with permission from Xiao et al. [58].

Fig. 21 .

 21 Fig. 21. Inference full-field in-plane mean velocity of the fully developed turbulent flow in a square duct, showing the lower left quadrant. (a) Prior velocity ensemble and (b) velocity ensemble at four spanwise locations with comparison to baseline and benchmark results. The velocity profiles in the prior ensemble are scaled by a factor of 0.3 for clarity. The upper half of the domain is omitted due to diagonal symmetry. Figure reprinted with permission from Xiao et al. [44].

Fig. 22 .

 22 Fig. 22. Comparison of modes (eigen-functions) obtained from a purely statistical kernel (squared exponential function) and a PDE-informed kernel as in Equation (42). Panels (a) and (b) reproduced from Xiao et al. [44] and Wu et al. [152], respectively.

  shear flows and reported values between 1.33 and 1.75, with a mean around 1.47. Note, however, that Equation (19) becomes 2.09 with the standard values for C 1 and C 2 , which is significantly different from the mentioned experimental values. Note that, regardless of the uncertainties, the coefficients have to satisfy the constraint C C

				2	>	1	as has
	been shown through numerical experiments by Ray et al. [23]. The
	physical reason behind this delineation is that the ratio C C / 2	1 corre-
	sponds to the spreading rate of a free jet. A ratio of C C / 2	1 1 < , or
	equivalently C	2	C	1

measured / k P for several uniform

Note that several definitions exist for the mean or average quantities [see, e.g., 10]. The most general one is the statistical ensemble average, which however is rarely used in current practice due to the large number of independent flow realizations required for convergence. For statistically steady flow, time average is used instead based on an ergodicity hypothesis. The same is also used for unsteady flows, although its validity is still controversial.

Here we have used the terminology ("parametric" and "non-parametric") rather liberally, which is closely related to, but not strictly consistent with, the standard terminology in the statistics literature. In statistics, parametric models refer to those parameterized by a finite set of parameters, while non-parametric models refer to those with infinite degrees of freedom (e.g., spatial random fields).

Acknowledgment

PC acknowledges support from French Agence Nationale de la Recherche (ANR). HX acknowledges support and mentoring from the Department of Aerospace and Ocean Engineering at Virginia Tech and particularly Prof. C.J. Roy and Prof. E.G. Paterson. The authors would like to thank Dr. J.-L. Wu, Dr. J.-X. Wang, Dr. R. P. Dwight, Dr. W. N. Edeling, and M. Schmeltzer for their collaborations and Mr. C. Michelén-Ströfer for useful suggestions on the manuscript. The authors gratefully acknowledge Dr. J.-L. Wu for providing figures, materials, and suggestions for this review. Finally, the authors thank the three anonymous reviewers for their constructive and valuable comments, which helped improving the quality and clarity of this review article.

Appendix A. Algorithms in uncertainty quantification

Appendix A. [START_REF] Sagaut | Multiscale and Multiresolution Approaches in Turbulence[END_REF]

. Plain Monte Carlo sampling

The algorithms for plain Monte Carlo sampling is rather straightforward. Given the probability distribution p ( )of the model parameters, Monte Carlo simulations can be used to obtain the distributions of the output. Specifically, the procedure of uncertainty propagation based on plain Monte Carlo simulation is as follows:

(1) Sampling. Draw a number of samples { , , , }

from the specified prior probability distribution p ( ).

(2) Propagation. For each of the sample, the model is evaluated to obtain the outputs y y y { , , , } This procedure is illustrated pictorially in Fig. 6a earlier.

Appendix A.2. Exact Bayesian inference with Markov chain Monte Carlo sampling

Much like the ergodicity assumption for the ensemble averaging to obtain the RANS equations, the MCMC sampling requires the ergodicity assumption. That is, any set within the state space can be reached from any other set with nonzero probability within finite steps. The MCMC procedure with Metropolis-Hastings sampling algorithm is as follows:

(1) Initialize the state (0) .

(2) Based on the current state z i ( ) , make a proposal of next state (e.g., a random walk), i.e., sample z z z q ( | ) This procedure is illustrated in Fig. 7. Intuitively, the sampler always accepts to go to a more likely state, which increases samples in high posterior probability regions. On the other hand, it also allows for the possibility of going to less likely states which allows for exploring the tails (rare events regions in the state space) and increase mixing (traveling back and forth in different regions).

Appendix A.3. Approximate Bayesian inference with iterative Ensemble Kalman method

In the example below, we assume the velocity at some locations is the observed physical state for notation simplicity. The augmented system state z x ( ) is written as a vector formed by stacking the unknown parameters and the physical states x ( ):

in which indicates vector transpose, and [ , , , ]

= is a vector of r parameters. Given the prior distributions for parameters ( ) to be inferred and the covariance matrix R of the observations y obs , the EnKF based inversion algorithm proceeds as follows:

(1) Sampling of prior distribution. From the prior distributions of the parameters, M samples are drawn. Each sample consists of a combination of values for . (2) Propagation. The output y ˆi are computed by using the updated parameters from the previous analysis step (or from the initial sampling if this is the first propagation step). The propagation is performed until next converged results are obtained. The ˆindicates predicted quantities that will be corrected in the analysis step below. The propagation is performed for each sample in the ensemble, leading to the propagated ensemble

Each sample z ˆj is a vector containing a realization of the velocity field and the parameters (see Equation ( 8)). The mean z ¯and covariance P of the propagated ensemble are estimated from the samples.

(3) Analysis/Correction. The computed physical fields (velocities) ˆi in the whole field are compared and sampled to compare with observations i obs . The ensemble covariance P and the error covariance Rare used to compute the Kalman gain matrix K as follows:

Each sample is then corrected as follows by using the Kalman gain matrix:

where superscript z j is the corrected system state; [ ,, ]

are the velocity, the part of the system state vector that can be observed; His the observation matrix. After the correction, the analyzed state contains updated velocities and parameters. × + , which maps a vector in the m dimensional state space to a vector in the n dimensional observation space. While point measurements of velocities are used as observations, other derived quantities such as lift, drag, pressure coefficients, surface coefficients, or velocities along a line of sight can be also used by choosing appropriate observation operators. For all forms of experimental data, the observation matrix Hin the filtering techniques relates the simulated system states to the observed quantities, i.e., y = z H . It is a mapping from system state space to the observation space. Example of observation operator is shown here. Consider the simple system shown in Fig. 23 to illustrate the principle. The simulation domain is discretized with 6 cells and the quantity of concern is the horizontal velocity only. Hence, the state vector has a dimension of 6 by 1. Three quantities are observed, a volumetric measurement of the velocity at cell 1, a velocity measurement at point B (which is the average of cells 2, 3, 5, 6), and an integrated measurement of the velocity along the line C, with weight factors of 1/2, 1/3, 1/6 for cells 4, 5, 6, respectively. The mapping y = z H can be written as:

Appendix B. Composite model theory and openbox treatment of model inadequacy

Introducing uncertainties into Reynolds stresses, in both parametric and non-parametric approaches, is motivated by a key consensus in the turbulence modeling community: Reynolds stress is the source of uncertainty in the RANS equations. This consensus is formulated officially as "composite model theory" by Oliver et al. [START_REF] Oliver | Validating predictions of unobserved quantities[END_REF]. That is, RANS simulations are based on reliable theories describing the conservation laws of mass, momentum, and energy, but contain approximate embedded models to account for the unresolved or unknown physics, i.e., the Reynolds stresses, leading to a composite model. This theory clearly separates the numerical model (simulator) into two components, i.e., rigorous equations and approximate closure models, and states that uncertainties should be introduced where they originate physically. This insight resulted in open-box approaches for uncertainty quantification and statistical inference. This is a major advance in model uncertainty quantification in RANS simulations compared to the earlier framework of Kennedy and O'Hagan [START_REF] Kennedy | Bayesian calibration of computer models[END_REF], where model inadequacy are introduced directly to the quantities of interest or the observed quantities and the numerical model (simulator) is treated as a blackbox. The open-box and blackbox approaches are compared schematically in Fig. 24. Composite models are ubiquitous in various disciplines of science and engineering. For example, in multiphase flow simulations, models are used to describe interphase mass and momentum exchanges in averaged equations [START_REF] Ma | Using statistical learning to close two-fluid multiphase flow equations for a simple bubbly system[END_REF][START_REF] Ma | Using statistical learning to close two-fluid multiphase flow equations for bubbly flows in vertical channels[END_REF]; in climate and weather modeling, parameterization are used to account for unresolved or unknown physics including radiation, cloud, and boundary layer processes [177][START_REF] Tiedtke | Representation of clouds in large-scale models[END_REF][START_REF] Chevallier | A neural network approach for a fast and accurate computation of a longwave radiative budget[END_REF]. In all these examples, the conservation laws are all expressed in well-grounded PDEs, albeit containing unclosed terms.

Appendix C. Uncertainties in DNS and their impact on RANS modeling

DNS data has long been considered the golden standard for evaluating the merits of turbulence models [START_REF] Mansour | Reynolds-stress and dissipation-rate budgets in a turbulent channel flow[END_REF][START_REF] Jeyapaul | Assessment of higher-order RANS closures in a decelerated planar wall-bounded turbulent flow[END_REF]. Evaluation of turbulence models can be done either a priori by comparing RANS-predicted Reynolds stresses with DNS data or a posteriori by comparing the fields solved by using the model of concern with mean fields from DNS. However, DNS are, like experimental observations, affected by more or less large uncertainties that may affect the comparisons: these mainly consist of sampling errors and discretization errors.

Although it is the instantaneous quantities such as velocities and pressure that are solved for in DNS, their statistical moments (e.g., means and covariance) are usually the quantities of interest, obtained by averaging a sufficiently large number of temporally uncorrelated samples of instantaneous fields. Sampling errors are caused by the fact that the samples may be correlated, and that the number of samples may not be sufficient to achieve statistical convergence, see Hoyas and Jiménez [START_REF] Hoyas | Reynolds number effects on the Reynolds-stress budgets in turbulent channels[END_REF] for a discussion.

For a properly performed DNS, the mesh must resolve all relevant flow scales and the sampling error is generally considered dominant. However, the mesh usually has to be chosen based on empirical judgment and sampling and discretization errors may be coupled. A Bayesian approach to account for sampling errors when estimating discretization errors is proposed in Oliver et al. [START_REF] Oliver | Estimating uncertainties in statistics computed from direct numerical simulation[END_REF].

Recently, evaluations of emerging data-driven turbulence models motivated an exercise that involves solving the RANS equations with specified Reynolds stresses, which is referred to as propagation. It has long been assumed that propagating accurate Reynolds stresses would lead to accurate velocities. However, numerous authors [START_REF] Wu | Physics-informed machine learning approach for augmenting turbulence models: a comprehensive framework[END_REF][START_REF] Thompson | A strategy for the eigenvector perturbations of the Reynolds stress tensor in the context of uncertainty quantification[END_REF][START_REF] Thompson | A methodology to evaluate statistical errors in DNS data of plane channel flows[END_REF][START_REF] Poroseva | On the accuracy of RANS simulations with DNS data[END_REF][START_REF] Wang | A comprehensive physics-informed machine learning framework for predictive turbulence modeling[END_REF] found significant discrepancies between the propagated velocities and the DNS velocities. On the other hand, Wang et al. [START_REF] Wang | A comprehensive physics-informed machine learning framework for predictive turbulence modeling[END_REF] performed the same propagation for fully developed turbulent flows in square ducts at various Reynolds numbers and found that the propagated velocities agree with DNS data satisfactorily. Such apparently conflicting findings were explained by different model conditioning in various flows, i.e., different sensitivity levels of the mean velocities to Reynolds stresses [START_REF] Wu | RANS equations with data-driven Reynolds stress closure can be ill-conditioned[END_REF].

In addition to the preceding sources of uncertainty, DNS also suffers from uncertainties associated with the specification of the boundary conditions in a similar way as LES does. Among studies based on probabilistic approaches we mention Ko et al. [START_REF] Ko | Sensitivity of two-dimensional spatially developing mixing layers with respect to uncertain inflow conditions[END_REF], who examined the sensitivity in the DNS of two-dimensional plane mixing layers to uncertainties in the inflow boundary conditions.