# Efficient prediction of turbulent flow quantities using a Bayesian
hierarchical multifidelity model

S. Rezaeiravesh<EMAIL_ADDRESS>T. Mukha<EMAIL_ADDRESS>P. Schlatter<EMAIL_ADDRESS>Department of Mechanical, Aerospace and
Civil Engineering, The University of Manchester, UK  SimEx/FLOW, Engineering
Mechanics, KTH Royal Institute of Technology, SE-100 44 Stockholm, Sweden
Swedish e-Science Research Centre (SeRC), Stockholm, Sweden

###### Abstract

High-fidelity scale-resolving simulations of turbulent flows quickly become
prohibitively expensive, especially at high Reynolds numbers. As a remedy, we
may use multifidelity models (MFM) to construct predictive models for flow
quantities of interest (QoIs), with the purpose of uncertainty quantification,
data fusion and optimization. For numerical simulation of turbulence, there is
a hierarchy of methodologies ranked by accuracy and cost, which include
several numerical/modeling parameters that control the predictive accuracy and
robustness of the resulting outputs. Compatible with these specifications, the
present hierarchical MFM strategy allows for simultaneous calibration of the
fidelity-specific parameters in a Bayesian framework as developed by Goh et
al. [11]. The purpose of the multifidelity model is to provide an improved
prediction by combining lower and higher fidelity data in an optimal way for
any number of fidelity levels; even providing confidence intervals for the
resulting QoI. The capabilities of our multifidelity model are first
demonstrated on an illustrative toy problem, and it is then applied to three
realistic cases relevant to engineering turbulent flows. The latter include
the prediction of friction at different Reynolds numbers in turbulent channel
flow, the prediction of aerodynamic coefficients for a range of angles of
attack of a standard airfoil, and the uncertainty propagation and sensitivity
analysis of the separation bubble in the turbulent flow over periodic hills
subject to the geometrical uncertainties. In all cases, based on only a few
high-fidelity data samples (typically direct numerical simulations, DNS), the
multifidelity model leads to accurate predictions of the QoIs accompanied with
an estimate of confidence. The result of the UQ and sensitivity analyses are
also found to be accurate compared to the ground truth in each case. Supported
by the present results, the hierarchical multifidelity model can be proposed
as a viable solution to the outer-loop problems in a range of computational
fluid dynamics (CFD) involving turbulent flows, especially in situations when
there is no direct correlation between the data of low- and high-fidelities;
occurring for instance in wall modeling for turbulent flows. In addition, the
methodology is applicable in other problems in computational physics, as long
as there can be a hierarchy of simulation fidelities defined, including
climate modeling and biomedical applications.

###### keywords:

Multifidelity modeling, Uncertainty quantification, Hierarchical models,
Bayesian inference, Turbulent flows.

##  1 Introduction

In science and engineering, different computational models can be derived to
make realizations of the quantities of interest (QoIs) of a process or an
event happening in reality. High-fidelity (HF) models can result in highly
accurate and robust realizations, but running them is typically
computationally expensive. In contrast, different low-fidelity (LF) models
with lower computational cost can be developed for the same process which,
however, potentially lead to lower accuracy QoIs due to partial or completely-
missing physics captured by the model. On the other hand, in different
applications arising in uncertainty quantification (UQ), data fusion, and
optimization numerous realizations of the QoIs are required associated to the
samples taken from the space of the inputs/parameters in order to make
reliable estimations (these non-intrusive problems are referred to as the
outer-loop problems). In this regard, multifidelity models (MFM) can be
constructed by combining realizations of the HF and LF models such that a
balance between the overall computational cost and predictive accuracy is
achieved. The goal is to provide, by combining HF and LF models, an estimate
of the QoI that is better than any of the models alone.

In the recent years, different types of MFMs have been applied to a wide range
of problems, see e.g. the recent review by Peherstorfer et al. [28]. The use
of the MFMs in studies of turbulent flows can be greatly advantageous,
considering the wide range of engineering applications relying on these flows
and also the high cost generally involved in the HF computations (such as
scale-resolving simulations) and experiments of the turbulent flows. There is
a distinguishable hierarchy in the fidelity of the computational models
utilized for simulation of turbulence, see e.g. Ref. [35]. Let us consider the
wall-bounded turbulent flows where a turbulent boundary layer forms at the
wall boundaries. Direct numerical simulation (DNS) can provide the highest-
fidelity results for a given turbulent flow, however it can become
prohibitively expensive at high Reynolds numbers which are relevant to
practical applications. The computational cost can be reduced by employing
large eddy simulation (LES) which aims at directly resolving the scales larger
than a defined size and modeling the unresolved effects. At the lowest cost
and fidelity level, Reynolds-averaged Navier-Stokes (RANS) simulations can be
performed which avoid directly resolving any flow fluctuations by resorting to
a statistical description of turbulence. Between RANS and wall-resolving LES,
other approaches such as hybrid RANS-LES and wall-modeled LES can be
considered, see [35, 21]. Although this clear hierarchy is extremely
beneficial when constructing MFMs as will be thoroughly discussed and
demonstrated in the present paper, there is a challenge to be dealt with: the
realizations of different turbulence simulation approaches are, in general,
sensitive to various modeling and numerical parameters as well as inputs. At
lower fidelities like RANS, modeling effects are dominant while as moving
towards LES and DNS, numerical factors become more relevant, including grid
resolution and discretization properties. Hereafter, these fidelity-specific
controlling parameters are referred to as tuning or calibration parameters.

Combining training data from different turbulence-simulation approaches, MFMs
are constructed over the space of design and uncertain parameters/inputs. An
appropriate approach to construct MFMs for turbulent flow problems should
systematically allow for simultaneous calibration of the tuning parameters. An
appropriate methodology which is employed in the present study is the
hierarchical multifidelity predictive model developed in [15, 11] in which the
calibration parameters of the involved fidelities are estimated using the data
of the higher fidelity models. This MFM which is hereafter referred to as HC-
MFM, can also incorporate the observational uncertainties. The HC-MFM can be
seen as an extension of the model by Higdon et al. [15] which was employed to
combine experimental (field) and simulation data. A fundamental component of
this class of MFMs is the Bayesian calibration of the computer models as
described in the landmark paper by Kennedy and O’Hagan [20]. At each level of
the MFM, the Gaussian process regression (GPR) [30] is employed to construct
surrogates for the simulators.

The application of the HC-MFM in the field of computational fluid dynamics
(CFD) and turbulent flows is novel, and we make specific adaptations suitable
for turbulence simulations. In this regard, the present paper aims at
assessing the useful potential of the HC-MFM by applying it to three examples
relevant to engineering wall-bounded turbulent flows. To highlight the
contributions of the present work, the existing studies in the literature
devoted to the development and application of the multifidelity models to CFD
and turbulent flows are briefly reviewed here by classifying them according to
their underlying MFM strategy. 1\. A model was originally introduced by
Kennedy and O’Hagan [19] where a QoI at each fidelity is expressed as a first-
order autoregressive model of the same QoI at the immediately lower fidelity.
Co-Kriging using Gaussian process regression to construct surrogates are
classified in this category, see e.g. Ref. [4, 39] for applications in
turbulence simulations. To enhance the computational efficiency of the co-
Kriging for several fidelity levels, recursive algorithms have been proposed
and applied to CFD problems, see [13, 29]. 2\. A class of MFMs has been
developed based on non-intrusive polynomial chaos expansion (PCE) and
stochastic collocation methods, see [25, 26]. Here, an additive or a
multiplicative term is considered to correct the LF model’s predictions
against the HF model. 3\. The multi-level multifidelity Monte Carlo (MLMF-MC)
models [3, 8] are appropriate for the UQ forward problems. These models are
developed by combining multilevel [9] and control-variate [27] MC methods to
improve the rate of convergence of the stochastic moments of the QoIs compared
to the the standard MC method. Jofre et al. [18] applied MLMF-MC models to an
irradiated particle-laden turbulent flow. The HF model was considered to be
DNS and the two LF models were based on a surrogate particle approach and
lower resolutions for flow and particles. 4\. Other models including the
hierarchical Kriging model based on GPR where the predictions of a LF model
are taken as the trend in the HF Kriging, see Han and Görtz [14].

Recently, Voet et al. [39] compared inverse weighted distance-, PCE- and co-
Kriging-based MFMs using the data of RANS and DNS for the turbulent flow over
a periodic hill, and concluded that the co-Kriging model outperforms the the
others in terms of accuracy. This is the first (and to our knowledge only)
study where multifidelity models have been applied to engineering-relevant
RANS and DNS data for the purpose of uncertainty propagation. Voet et al. [39]
also found that the performance of the co-Kriging can deteriorate when there
is no significant correlation between the RANS and DNS data and at the same
time there is a significant deviation between them. Motivated by this
deficiency, we adapt and use the HC-MFM where the discrepancy between the data
(and not their correlation) over the space of design parameters is learned
using independent Gaussian processes. The model is absolutely generative and
can be extended to an arbitrary number of fidelity levels. Besides the
systematic calibration of the fidelity-specific parameters during its training
stage, the HC-MFM is also capable of handling uncertain data, as for instance
happens when QoIs are turbulence statistics computed over a finite time-
averaging interval (recently, a framework is proposed to combine these
observational uncertainties with parametric ones, see Ref. [34]). Relying on
these characteristics, the HC-MFM is suitable to be applied to the data of
various turbulence simulation methodologies to address different types of the
outer-loop problems. In contrast to all the previous studies (at least in
CFD), we adopt a Bayesian inference to construct the HC-MFM, a feature which
results in more accurate models as well as estimating confidence intervals for
the predictions.

The rest of the paper is organized as follows. In Section 2, various elements
of the HC-MFM approach are introduced and explained. Section 3 is devoted to
the application of the HC-MFM to an illustrative example, turbulent channel
flow, polars for an airfoil, and analysis of the geometrical uncertainties in
the turbulent flow over a periodic hill. The summary of the paper along with
the conclusions is presented in Section 4.

##  2 Method

In this section, the hierarchical MFM with calibration (HC-MFM) developed by
Goh et al. [11], which forms the basis for the present study is reviewed. We
will proceed by sequentially going through the aspects of Gaussian process
regression (GPR), model calibration, and eventually the HC-MFM formulation.

###  2.1 Gaussian Process Regression

In general, the Gaussian processes (GPs) provide a way to systematically build
a representation of the QoI as a function of the various inputs to the model.
Eventually, regression can be performed by evaluating the Gaussian process at
new inputs not seen by the model before. Let
$\mathbf{x}\in\mathbb{X}\subset\mathbb{R}^{p_{\mathbf{x}}}$ represent the
controllable inputs and parameters, adopting the notation from Ref. [20]. As a
convention, all boldface letters are hereafter considered to be a vector or a
matrix. The design and uncertain parameters appearing in optimization and UQ
analyses, respectively, can also be classified as $\mathbf{x}$. A Gaussian
process $\hat{f}(\cdot)$, see e.g. Ref. [30], can be employed to map the
inputs $\mathbf{x}$ to the QoIs or outputs $y\subset\mathbb{R}$ of the
computer codes (simulators) or field data. For a finite set of training
samples $\\{\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{n}\\}$ with
corresponding observations $\\{y_{1},y_{2},\ldots,y_{n}\\}$, the collection of
$\\{\hat{f}(\mathbf{x}_{1}),\hat{f}(\mathbf{x}_{2}),\ldots,\hat{f}(\mathbf{x}_{n})\\}$
will have a joint Gaussian (multivariate normal) distribution, [30]. The GP
$\hat{f}(\mathbf{x})$ is written as,

$\hat{f}(\mathbf{x})\sim\mathcal{GP}\left(m(\mathbf{x}),k(\mathbf{x},\mathbf{x}^{\prime})\right)\,,$
(1)

which is fully described by its mean $m(\mathbf{x})$ and covariance function
$k(\mathbf{x},\mathbf{x}^{\prime})$ defined as,

$\displaystyle m(\mathbf{x})$ $\displaystyle=$
$\displaystyle\mathbb{E}[\hat{f}(\mathbf{x})]\,,$ (2) $\displaystyle
k(\mathbf{x},\mathbf{x}^{\prime})$ $\displaystyle=$
$\displaystyle\mathbb{E}[(\hat{f}(\mathbf{x})-m(\mathbf{x}))(\hat{f}(\mathbf{x}^{\prime})-m(\mathbf{x}^{\prime}))]\,.$
(3)

In general, the GPs can be used in the case of having observation noise
$\bm{\varepsilon}$ in the $y$ data. Using an additive error model, we have,

$y(\mathbf{x})=\hat{f}(\mathbf{x})+\bm{\varepsilon}\,,$ (4)

where the noises are assumed to be independent and have Gaussian distribution
$\bm{\varepsilon}\sim\mathcal{N}(0,\sigma^{2})$.

In the Gaussian process regression (GPR), given a set of training data
$\mathcal{D}=\\{\mathbf{x}_{i},y_{i}\\}_{i=1}^{n}$ the posterior and posterior
predictive distributions of $\hat{f}(\cdot)$ and $y$, respectively, at test
inputs $\mathbf{x}^{*}\in\mathbb{X}$ can be inferred in a Bayesian framework,
see e.g. Ref. [30]. To this end, first a prior distribution for
$\hat{f}(\mathbf{x})$, see Eq. (1), is assumed through specifying functions
for the mean and covariance in Eqs. (2) and (3), where there are unknown
hyperparameters $\bm{\beta}$ in the functions. Using the training data, the
posterior distribution of $\bm{\beta}$ is learned. As a main advantage of the
GPR, the predictions at test samples will be accompanied by an estimate of
uncertainty.

###  2.2 Model Calibration

As pointed out in Section 1, the outputs of computational models (simulators)
at a given $\mathbf{x}$ may depend on different tuning or calibration
parameters, $\mathbf{t}\in\mathbb{T}\subset\mathbb{R}^{p_{\mathbf{t}}}$. Given
a set of observations, these parameters can be calibrated through conducting a
UQ inverse problem which can be expressed in a Bayesian framework, see Kennedy
and O’Hagan [20]. The calibrated model can then not only be employed for
prediction, but also for fusion of the field and simulation data, see Higdon
et al. [15]. Consider $n_{1}$ data samples
$\\{(\mathbf{x}_{i},y_{i})\\}_{i=1}^{n_{1}}$ are observed for a physical
process $\zeta(\mathbf{x})$. To statistically model the observations, a
simulator $\hat{f}(\mathbf{x},\bm{\theta})$ can be employed in which the
$\bm{\theta}$ are the true or optimal values of $\mathbf{t}$ and are to be
estimated from the training data. However, in general, it is possible that
even the calibrated simulator $\hat{f}(\mathbf{x},\bm{\theta})$ produces
observations which systematically deviate from reality. To remove such a bias,
a model-discrepancy term $\hat{\delta}(\mathbf{x})$ can be added to the
simulator, [20, 15]. In many applications, particularly in CFD and turbulent
flow simulations, the flow solver can be run only a limited number of times.
In any realization, the adopted values for the tuning parameters $\mathbf{t}$
are not necessarily optimal and hence potentially lead to the outputs which
are systematically different from the QoIs in reality. For the described
calibration problem, the Kennedy and O’Hagan model [20] reads as,

$\begin{cases}y_{i}=\hat{f}(\mathbf{x}_{i},\bm{\theta})+\hat{\delta}(\mathbf{x}_{i})+\varepsilon_{i}\,&,\quad
i=1,2,\cdots,n_{1}\\\ y_{i}=\hat{f}(\mathbf{x}_{i},\mathbf{t}_{i})\,&,\quad
i=1+n_{1},2+n_{1},\cdots,n_{2}+n_{1}\end{cases}\,,$ (5)

where, $\hat{\cdot}$ specifies a GP and $n_{2}$ is the number of simulated
data. Note that the samples $\\{\mathbf{x}_{i}\\}_{i=1+n_{1}}^{n_{2}+n_{1}}$
are not necessarily the same as $\\{\mathbf{x}_{i}\\}_{i=1}^{n_{1}}$ at which
the observations are made. The index $i$ should be seen as a global index
which implies that a different model is used for each of the two subranges of
$i$. Given the $n_{1}+n_{2}$ data, the posterior distribution for the
calibration parameters $\bm{\theta}$ along with that of the hyperparameters in
the GPs, $\bm{\beta}$ is estimated. Further details are provided in the
section below.

###  2.3 The Hierarchical Multifidelity Model with Automatic Calibration (HC-
MFM)

Goh et al. [11] extended the model (5) to an arbitrary number of fidelity
levels which together form a modeling hierarchy for a physical process. As a
main feature of the resulting MFM, each fidelity can, in general, have its own
calibration parameters and also share some calibration parameters with other
fidelities. The basics of the MFM comprising three fidelity levels are
explained below, noting that adapting the formulation to any number of
fidelities with different combinations of parameters is straightforward. We
assume that the fidelity of the models decreases from $M_{1}$ to $M_{3}$, and
in practice due to the budget limitations, the number of training data
decreases with increasing the model fidelity. The HC-MFM for three fidelities
reads as [11],

$\displaystyle\begin{cases}y_{M_{1}}(\mathbf{x}_{i})\,={\hat{f}(\mathbf{x}_{i},\bm{\theta}_{3},\bm{\theta}_{s})}+{\hat{g}(\mathbf{x}_{i},\bm{\theta}_{2},\bm{\theta}_{s})}+{\hat{\delta}(\mathbf{x}_{i})}+\varepsilon_{1_{i}}\,&,i=1,2,\cdots,n_{1}\\\
y_{M_{2}}(\mathbf{x}_{i})\,={\hat{f}(\mathbf{x}_{i},\bm{\theta}_{3},\mathbf{t}_{s_{i}})}+{\hat{g}(x_{i},\mathbf{t}_{2_{i}},\mathbf{t}_{s_{i}})}+\varepsilon_{2_{i}}\,&,{i=1+n_{1},2+n_{1},\cdots,n_{2}+n_{1}}\\\
y_{M_{3}}(\mathbf{x}_{i})={\hat{f}(\mathbf{x}_{i},\mathbf{t}_{3_{i}},\mathbf{t}_{s_{i}})}+\varepsilon_{3_{i}}\,&,i=1+n_{1}+n_{2},2+n_{1}+n_{2},\cdots,n_{3}+n_{1}+n_{2}\end{cases}\,,$
(6)

where subscript $s$ denotes the parameters which are shared between the
models, whereas $\mathbf{t}_{2}$ and $\mathbf{t}_{3}$ are the calibration
parameters specific to fidelities $M_{2}$ and $M_{3}$, respectively. The
noises are assumed to have Gaussian distributions with zero mean. At each
fidelity level, the associated simulator is created by adding a model
discrepancy term to the simulator describing the immediately lower fidelity.
Concatenating all training data, an augmented vector $\mathbf{Y}$ of size
$n_{1}+n_{2}+n_{3}$ is obtained, for which the covariance matrix can be
written in terms of the covariances of $\hat{f}(\cdot)$, $\hat{g}(\cdot)$,
$\hat{\delta}(\cdot)$ and the observational noise,

$\displaystyle\bm{\Sigma}={\bm{\Sigma}_{f}}+\begin{bmatrix}{\bm{\Sigma}_{g}}&\mathbf{0}_{(n_{1}+n_{2})\times
n_{3}}\\\ \mathbf{0}_{n_{3}\times(n_{1}+n_{2})}&\mathbf{0}_{n_{3}\times
n_{3}}\\\
\end{bmatrix}+\begin{bmatrix}{\bm{\Sigma}_{\delta}}&\mathbf{0}_{n_{1}\times(n_{2}+n_{3})}\\\
\mathbf{0}_{(n2+n3)\times
n_{1}}&\mathbf{0}_{(n_{2}+n_{3})\times(n_{2}+n_{3})}\\\
\end{bmatrix}+\begin{bmatrix}\bm{\Sigma}_{\varepsilon_{1}}&\mathbf{0}_{n_{1}\times
n_{2}}&\mathbf{0}_{n_{1}\times n_{3}}\\\ \mathbf{0}_{n_{2}\times
n_{1}}&\bm{\Sigma}_{\varepsilon_{2}}&\mathbf{0}_{n_{2}\times n_{3}}\\\
\mathbf{0}_{n_{3}\times n_{1}}&\mathbf{0}_{n_{3}\times
n_{2}}&\bm{\Sigma}_{\varepsilon_{3}}\\\ \end{bmatrix}\,.$ (7)

Appropriate kernel functions should be chosen to express the structure of the
covariances. Using samples $i$ and $j$ of the inputs and parameters, the
associated element in the covariance matrix $\bm{\Sigma}_{f}$ will be obtained
from,

$\bm{\Sigma}_{f_{ij}}=\text{cov}(\mathbf{x}_{i},\mathbf{t}_{3_{i}},\mathbf{t}_{s_{i}},\mathbf{x}_{j},\mathbf{t}_{3_{j}},\mathbf{t}_{s_{j}})=\lambda_{f}^{2}\,k(\bar{d}_{f_{ij}})\,,$
(8)

where $\lambda_{f}$ is a hyperparameter and $\bar{d}_{f_{ij}}$ is the scaled
Euclidean distance between two samples $i$ and $j$ over the space of
$(\mathbf{x},\mathbf{t}_{3},\mathbf{t}_{s})$:

$\displaystyle\bar{d}^{2}_{f_{ij}}$ $\displaystyle=$
$\displaystyle\bar{d}^{2}(\mathbf{x}_{i},\mathbf{x}_{j})+\bar{d}^{2}(\mathbf{t}_{3_{i}},\mathbf{t}_{3_{j}})+\bar{d}^{2}(\mathbf{t}_{s_{i}},\mathbf{t}_{s_{j}})$
(9) $\displaystyle=$
$\displaystyle\sum_{l=1}^{p_{\mathbf{x}}}\frac{(x_{l_{i}}-x_{l_{j}})^{2}}{\ell^{2}_{f_{x_{l}}}}+\sum_{l=1}^{p_{\mathbf{t}_{3}}}\frac{(t_{3_{l_{i}}}-t_{3_{l_{j}}})^{2}}{\ell^{2}_{f_{{t_{3}}_{l}}}}+\sum_{l=1}^{p_{\mathbf{t}_{s}}}\frac{(t_{s_{l_{i}}}-t_{s_{l_{j}}})^{2}}{\ell^{2}_{f_{{t_{s}}_{l}}}}\,.$
(10)

Here, $p_{\mathbf{x}}$, $p_{\mathbf{t}_{3}}$, and $p_{\mathbf{t}_{s}}$ specify
the dimension of $\mathbf{x}$, $\mathbf{t}_{3}$, and $\mathbf{t}_{s}$,
respectively. Correspondingly, the length-scale over the $l$-th dimension of
each of these spaces is represented by $\ell_{f_{x_{l}}}$,
$\ell_{f_{{t_{3}}_{l}}}$, and $\ell_{f_{{t_{s}}_{l}}}$, respectively. These
length-scales are among the hyperparameters $\bm{\beta}$ to be learned when
constructing the HC-MFM. There are various options for modeling the covariance
kernel function $k(\cdot)$, see e.g. Refs. [30, 12], among which the
exponentiated quadratic and Matern-5/2 [23] functions are used in the examples
in Section 3. These two functions respectively read as,

$k(\bar{d}_{f_{ij}})=\exp(-0.5\,\bar{d}^{2}_{f_{ij}})\,,$ (11)

and,

$k(\bar{d}_{f_{ij}})=\left[1+\sqrt{5}\,\bar{d}_{f_{ij}}+\frac{5}{3}\bar{d}^{2}_{f_{ij}}\right]\exp(-\sqrt{5}\,\bar{d}_{f_{ij}})\,.$
(12)

Similar expressions can be derived for
$\bm{\Sigma}_{g_{ij}}=k_{g}(\mathbf{x}_{i},\mathbf{t}_{2_{i}},\mathbf{t}_{s_{i}},\mathbf{x}_{j},\mathbf{t}_{2_{j}},\mathbf{t}_{s_{j}})$
and ${\bm{\Sigma}_{\delta_{ij}}=k_{\delta}(\mathbf{x}_{i},\mathbf{x}_{j})}$
appearing in Eq. (7). This leads to introducing new hyperparameters associated
to the GPs. Note that, given how the training vector $\mathbf{Y}$ and
associated inputs are assembled, correct combinations of training data for the
inputs and parameters will be used in the kernels. The unknown parameters to
estimate include calibration parameters in different models, $\bm{\theta}$,
and hyperparameters $\bm{\beta}$ appearing in the GPs. Following the Bayes
rule, the posterior distribution of these parameters given the training data
$\mathbf{Y}$ can be inferred from [11, 20, 15],

$\pi(\bm{\theta},\bm{\beta}|\mathbf{Y})\propto\pi(\mathbf{Y}|\bm{\theta},\bm{\beta})\pi_{0}(\bm{\theta})\pi_{0}(\bm{\beta})\,,$
(13)

where $\pi(\mathbf{Y}|\bm{\theta},\bm{\beta})$ specifies the likelihood
function and $\pi_{0}(\cdot)$ represents a prior distribution. Note that all
priors are assumed to be independent. For all GPs in Section 3, the prior
distribution for $\lambda$ appearing in the covariance matrices such as Eq.
(8) is taken to be half-Cauchy whereas the length-scales $\ell$ in Eqs. (11)
and (12) are assumed to have Gamma distributions. The exact definition of the
priors will be provided later for each case in Section 3, and Table 2
summarizes the formulation of the standard distributions used as priors. The
standard-deviation of the noises are assumed to be the same for which a half-
Cauchy prior is adopted. For the calibration parameters $\bm{\theta}$,
Gaussian or uniform priors are considered. In some cases, we may consider a
constant mean function for the GPs, where a Gaussian distribution is used as
the prior. Due to this, the predictions of a trained MFM when it is used to
extrapolate in $\mathbf{x}$ (outside of the range of training samples) should
be used with caution. To avoid potential inaccuracies, in general, more
elaborate mean functions can be used when constructing the HC-MFMs (this,
however, is not the subject of the present study).

Given the training data $\mathbf{Y}$, a Markov Chain Monte Carlo (MCMC)
technique can be used to draw samples from the posterior distributions of
$\bm{\theta}$ and $\bm{\beta}$, and hence construct a HC-MFM. In the present
study, the described HC-MFM (6) has been implemented in Python using the PyMC3
[36] package with the NUTS MCMC sampling approach [16]. As it will be shown in
Section 3.6, the MCMC sampling method may lead to more accurate results
compared to the point estimators.

After being constructed, an HC-MFM can be used for predicting the QoI $y$ for
any new sample taken from the space of inputs $\mathbf{x}$. The accuracy of
the predicted QoIs will be assessed by measuring their deviation from the
validation data of the highest fidelity ${M_{1}}$. As detailed in Ref. [11],
the joint distribution of the training $\mathbf{Y}$ and new $y^{*}$
(associated to a test sample $\mathbf{x}^{*}$) conditioned on
$\bm{\theta},\bm{\beta}$ will have a multivariate normal distribution with a
covariance matrix of the same structure as $\bm{\Sigma}$ in Eq. (7). For any
joint sample drawn from the posterior distribution of
$\pi(\bm{\theta},\bm{\beta}|\mathbf{Y})$, a sample prediction for $y^{*}$ is
made. Repeating this procedure for a large number of times, valid estimations
for the posterior of the predictions $y^{*}$ can be achieved. Therefore,
estimating the confidence in the predictions is straightforward. Note that at
this stage, various UQ analyses can be performed using the HC-MFM as a
surrogate of the physical process over $\mathbf{x}$.

##  3 Results and Discussion

Four examples are considered to which the HC-MFM described in the previous
section is applied. The first example in Section 3.1 is used to validate the
implementation of the MFM, and the next three examples are relevant to
fundamental and engineering analysis of wall-bounded turbulent flows.

###  3.1 An Illustrative Example

Consider the following analytical model taken from Forrester et al. [5] to
generate high- and low-fidelity samples of the QoI $y$ for input $x\in[0,1]$,

$\begin{cases}y_{H}(x)=(\theta x-2)^{2}\sin(2\theta x-4)\\\
y_{L}(x)=y_{H}(x)+B(x-0.5)-C\\\ \end{cases}\,.$ (14)

In Ref. [5], $\theta$ is taken to be fixed and equal to $6$, but here it is
treated as an uncertain calibration parameter that is to be estimated during
the construction of the MFM. Note that the notations of the general model (5)
can be adopted for Eq. (14). For simulator $\hat{f}(x,t)$ and model
discrepancy $\hat{\delta}(x)$ the covariance matrix in Eq. (8) is used with
the exponentiated quadratic kernel (11). The following prior distributions are
considered: $\lambda_{f},\lambda_{\delta}\sim\mathcal{HC}(\alpha=5)$,
$\ell_{f_{x}},\ell_{f_{t}},\ell_{\delta_{x}}\sim\Gamma(\alpha=1,\beta=5)$,
$\varepsilon\sim\mathcal{N}(0,\sigma)$ with
$\sigma~{}\sim\mathcal{HC}(\alpha=5)$, and $\theta\sim\mathcal{U}[5.8,6.2]$.
Here, $\mathcal{HC}$, $\Gamma$, $\mathcal{N}$ and $\mathcal{U}$ denote the
half-Cauchy, Gamma, Gaussian, and uniform distributions, respectively, see
Table 2. The HF training samples are taken at $x=\\{0,0.4,0.6,1\\}$,
therefore, $n_{H}=4$ is fixed. To investigate the effect of $n_{L}$, three
sets of LF samples of size $10$, $15$, and $20$ are considered which are
generated by Latin hypercube sampling from the admissible space
$[0,1]\times[5.8,6.2]$ corresponding to $x$ and $t$ (uncalibrated instance of
$\theta$), respectively.

Using the data, the HC-MFM (6) for problem (14) is constructed. The first row
in Figure 1 shows the predicted $y$ with the associated $95\%$ confidence
interval (CI) along with the training data and reference true data generated
with $\theta=6$. For all $n_{L}$, the predicted $y$ is closer to HF data than
the LF data, however, for $n_{L}=15$ and $20$, the agreement between the mean
of the predicted $y$ and the true data is significantly improved. A better
validation can be made via the plots in the second row of Figure 1, where the
predicted $y$ and true values of $y_{H}$ at $50$ uniformly-spaced test samples
for $x\in[0,1]$ are plotted. Clearly, increasing the number of the LF samples
while keeping $n_{H}=4$ fixed, improves the predictions and reduces the
uncertainty. In the third row of Figure 1, the posterior densities of $\theta$
are presented. In all cases, a uniform (non-informative) prior distribution
over $[5.8,6.2]$ was considered for $\theta$. Only for $n_{L}=20$, the
resulting posterior density of $\theta$ is high near the true value $6$.
Therefore, it is confirmed that, as explained by Goh et al. [11] the main
capability of the HC-MFM (6) is in making accurate predictions for $y$ and
only if a sufficient number of training data is available, accurate
distributions for the calibration parameters are also obtained. This is shown
here by fixing $n_{H}$ and increasing $n_{L}$, which is favorable in practice.
It is also noteworthy that if $\theta$ was known and hence treated as a fixed
parameter, then even with $n_{L}=10$ very accurate predictions for $y$ could
be already achieved (not shown here).

It also should be noted that the mean of the posterior distribution of
$\sigma$, the noise standard deviation, is found to be negligible, as
expected. This is in fact the case for all other examples in this section.

|  |
---|---|---
|  |
|  |
Figure 1: (Top) Predicted QoI $y$ by HC-MFM (6) along with the training and
true data, (Middle) predicted $y$ versus true observations at $50$ test
samples of $x\in[0,1]$ with error bars representing $95\%$ CI, (Bottom)
posterior probability density function (PDF) of $\theta$ based on $10^{4}$
MCMC samples. The $y_{H}$ and $y_{L}$ training data are generated from Eq.
(14) using $B=C=10$. The training data includes $4$ HF samples combined with
(left column) $10$, (middle column) $15$, and (right column) $20$ LF samples.
The true data is generated by Eq. (14) using $\theta=6$.

###  3.2 Turbulent Channel Flow

Turbulent channel flow is one of the most canonical wall-bounded turbulent
flows. The flow develops between two parallel flat walls which are apart by
the distance $2\delta$, and the flow is periodic in the streamwise and
spanwise directions. Channel flow is fully defined by a Reynolds number, for
instance the bulk Reynolds number $\mbox{Re}_{b}=U_{b}\delta/\nu$, where
$U_{b}$ and $\nu$ specify the streamwise bulk velocity and kinematic
viscosity, respectively. Among different possible QoIs, here we only focus on
the friction velocity $\langle u_{\tau}\rangle$, as a function of Reynolds
number. This quantity is defined as $\sqrt{\langle\tau_{w}\rangle/\rho}$,
where $\tau_{w}$ and $\rho$ are the magnitude of the wall-shear stress and
fluid density, respectively, and $\langle\cdot\rangle$ represents averaging
over time and the periodic directions. Three fidelity levels are considered:
DNS ($M_{1}$), WRLES ($M_{2}$), and a reduced-order algebraic model ($M_{3}$),
where the fidelity reduces from the former to the latter. We use the DNS data
of Refs. [17, 22, 44]. The WRLES of channel flow have been performed at
different Reynolds numbers without any explicit subgrid-scale model using
OpenFOAM [40] which is an open-source finite-volume flow solver. For the
details of simulations see Ref. [32], where it was shown that for a fixed
resolution in the wall-normal direction, variation of the grid resolutions in
the wall-parallel directions could significantly impact the accuracy of the
flow QoIs. Therefore, in the context of the HC-MFM, the calibration parameters
for WRLES are taken to be ${\Delta x^{+}}$ and ${\Delta z^{+}}$, which are the
cell dimensions in the streamwise and spanwise directions, respectively,
expressed in wall-units (${\Delta x^{+}}=\Delta x\,u^{\circ}_{\tau}/\nu$ where
$u^{\circ}_{\tau}$ is the reference $u_{\tau}$ from DNS). At the lowest
fidelity, the following reduced-order algebraic model is considered which is
derived by averaging the streamwise momentum equation for the channel flow in
the periodic directions and time,

$\langle
u_{\tau}\rangle^{2}/U_{b}^{2}=\frac{1}{\mbox{Re}_{b}}\frac{\mbox{d}}{\mbox{d}\eta}\left((1+\zeta(\eta))\frac{\mbox{d}\langle
u\rangle/U_{b}}{\mbox{d}\eta}\right)\,,$ (15)

where $\eta$ is the distance from the wall normalized by the channel half-
height $\delta$, and $\zeta(\eta)$ is the normalized eddy viscosity $\nu_{t}$.
Reynolds and Tiederman [31] proposed the following closed form for
$\zeta(\eta)$,

$\displaystyle\zeta(\eta)=\frac{\nu_{t}(\eta)}{\nu}=\frac{1}{2}\left[1+\frac{{\kappa}^{2}\mbox{Re}^{2}_{\tau}}{9}\left(1-(\eta-1)^{2}\right)^{2}\left(1+2(\eta-1)^{2}\right)^{2}\left(1-\exp(\frac{-\eta\mbox{Re}_{\tau}}{A^{+}})\right)^{2}\right]^{1/2}-\frac{1}{2}\,,$
(16)

where $\mbox{Re}_{\tau}=\langle u_{\tau}\rangle\delta/\nu$ is the friction-
based Reynolds number, and $\kappa$ and $A^{+}$ are two modeling parameters.
At any $\mbox{Re}_{b}$ (and given value of $\kappa$ and $A^{+}$), Eq. (15) is
integrated over $\eta\in[0,1]$ and is iteratively solved using Eq. (16) to
estimate $\langle u_{\tau}\rangle$. Expressing the channel flow example in the
terminology of MFM (6), $\langle u_{\tau}\rangle/U_{b}$ is the QoI $y$,
$x=\mbox{Re}_{b}$, $\mathbf{t}_{3}=(\kappa,A^{+})$, and
$\mathbf{t}_{2}=({\Delta x^{+}},{\Delta z^{+}})$. The training data set
consists of the following databases. For DNS, $\langle u_{\tau}\rangle$ is
taken from Refs. [17, 22, 44] at $\mbox{Re}_{b}=5020$, $6962$, $10000$,
$20000$, $125000$ and $200400$. In total, $16$ WRLES $\langle u_{\tau}\rangle$
samples are obtained from a design of experiment based on the prior
distributions ${\Delta x^{+}}\sim\mathcal{U}[15,85]$ and ${\Delta
z^{+}}\sim\mathcal{U}[9.5,22]$ at $\mbox{Re}_{b}=5020,6962,10000$, and
$20000$. Here, we do not consider the observational uncertainty in $\langle
u_{\tau}\rangle$ which could, for instance, be due to finite time-averaging in
DNS and WRLES, but in general the HC-MFM could take such information into
account. The reduced-order model (15) which is computationally cheap is run
for $10$ values of $\mbox{Re}_{b}$ in range $[2000,200200]$. For each
$\mbox{Re}_{b}$, $9$ joint samples of $(\kappa,A^{+})$ are generated assuming
$\kappa\sim\mathcal{U}[0.36,0.43]$ and $A^{+}\sim\mathcal{U}[26.5,29]$ (note
that $\kappa$ is the von Kármán coefficient). For all the Gaussian processes
in the MFM (6), the exponentiated quadratic covariance function (11) is used.
The prior distribution of the hyperparameters are set as the following:
$\lambda_{f},\lambda_{g}\sim\mathcal{HC}(\alpha=5)$,
$\lambda_{\delta}\sim\mathcal{HC}(\alpha=3)$,
$\ell_{f_{x}},\ell_{f_{t_{3}}},\ell_{g_{x}},\ell_{g_{t_{2}}},\ell_{\delta_{x}}\sim\Gamma(\alpha=1,\beta=5)$,
and the noise standard deviation $\sigma\sim\mathcal{HC}(\alpha=5)$ (assumed
to be the same for all fidelities).

(a) |  | (b) |
---|---|---|---
Figure 2: (a) Mean prediction of $\langle u_{\tau}\rangle/U_{b}$ and
associated $95\%$ CI along with the training data and validation data from DNS
of Refs. [17, 22, 44], (b) diagonal: posterior density of parameters $\kappa$,
$A^{+}$, ${\Delta x^{+}}$, and ${\Delta z^{+}}$; off-diagonal: contour lines
of the joint posterior densities of these parameters. The value of the contour
lines increases from the lightest to darkest color.

Using the described training data in the HC-MFM (6) and running the MCMC chain
for $7000$ samples, after an initial $2000$ samples discarded due to burn-in,
the model is constructed. According to Figure 2(a), the predicted mean of
$\langle u_{\tau}\rangle/U_{b}$ follows the trend of the DNS data. This
approximately holds even at high Reynolds numbers, where there is a large
systematic error in the algebraic model and no WRLES data is available. As
expected, in this range due to scarcity of the DNS data, the uncertainty in
the predictions is high. The plot in Figure 2(b) shows the joint posterior
distributions of the calibration parameters $\kappa$, $A^{+}$, ${\Delta
x^{+}}$, and ${\Delta z^{+}}$ along with the histogram of each parameter. As
mentioned above, the prior distributions of all of these parameters were
assumed to be uniform and mutually independent. However, the resulting
posterior densities for $\kappa$ and $A^{+}$ are not uniform and the samples
of these parameters are correlated. More interestingly, the peak of the
posterior density of $\kappa$ is close to the value of $0.4$ that is assumed
to be universal across various flows and Reynolds numbers within the
turbulence community. In contrast, the posteriors of ${\Delta x^{+}}$ and
${\Delta z^{+}}$ are found to be still close to the prior uniform
distributions and no correlation between their samples is observed. This may
be at least partially be due to the fact that the number of the WRLES data is
limited as they are obtained only at 4 Reynolds numbers. Nevertheless, over
this range of the Reynolds number the posterior prediction of the QoI $\langle
u_{\tau}\rangle$ is very accurate and has the lowest uncertainty, which seems
to indicate that the grid spacing indeed does not lead to a bias for the shear
stress computation. This may in fact be an important aspect for building
future wall models.

###  3.3 Polars for the NACA0015 Airfoil

In this section, the HC-MFM (6) is applied to a set of data for the lift and
drag coefficients, $C_{L}$ and $C_{D}$, respectively, of a wing with a
NACA0015 airfoil profile at Reynolds number $1.6\times 10^{6}$. The angle of
attack (AoA) between the wing and the ambient flow is taken to be the design
parameter $x$. The data comprises of the following sources with respective
fidelities in brackets: wind-tunnel experiments by Bertagnolio [2] ($M_{1}$),
detached-eddy simulations (DES) ($M2$) and two-dimensional RANS ($M_{3}$) both
by Gilling et al. [10]. In their numerical study, Gilling et al. [10]
investigated the sensitivity of the DES results to the resolved turbulence
intensity (TI) of the fluctuations imposed at the inlet boundary. The
sensitivity was found to be particularly significant near the stall angle.
Therefore, when constructing an MFM, the calibration parameter $t_{2}$ in
fidelity $M_{2}$ is taken to be the TI.

For the covariance of the Gaussian processes $\hat{f}(x)$ and
$\hat{g}(x,t_{2})$ in HC-MFM (6), the exponentiated quadratic and Matern-5/2
kernel functions (11) and (12) are used, respectively. The following prior
distributions are assumed for the hyperparameters:
$\lambda_{f},\lambda_{g}\sim\mathcal{HC}(\alpha=5)$,
$\ell_{f_{x}}\sim\Gamma(\alpha=1,\beta=5)$,
$\ell_{g_{x}},\ell_{g_{t_{2}}}\sim\Gamma(\alpha=1,\beta=3)$, and the noise
standard deviation $\sigma\sim\mathcal{HC}(\alpha=5)$ (assumed to be the same
for all fidelities). To make the model capable of capturing large-scale
separation, the stall AoA, $x_{\rm stall}$, is included as a new calibration
parameter in the MFM. Our suggestion is to consider a piecewise kernel
function for the covariance of $\hat{\delta}(x)$ where $x_{\rm stall}$ is the
merging point. If the kernel functions for the AoAs smaller and larger than
$x_{\rm stall}$ are denoted by $k_{\delta_{1}}(\cdot)$ and
$k_{\delta_{2}}(\cdot)$, respectively, then the covariance function for
$\hat{\delta}(x)$ may be constructed as,

$\Sigma_{\delta_{ij}}=\lambda_{\delta_{1}}^{2}\varphi(x_{i})k_{\delta_{1}}(\bar{d}_{\delta_{ij}})\varphi(x_{j})+\lambda_{\delta_{2}}^{2}\varphi(x_{i})k_{\delta_{2}}(\bar{d}_{\delta_{ij}})\varphi(x_{j})\,,$
(17)

where, $\bar{d}_{\delta_{ij}}$ is defined similar to those in Eqs. (9) and
(10) but only in $x$, and $\varphi(x)$ is a function to smoothly merge the two
covariance functions. In particular, we use the logistic function,

$\varphi(x)=\left[1+\exp(-\alpha_{\rm stall}(x-x_{\rm stall}))\right]^{-1}\,,$
(18)

where $\alpha_{\rm stall}$ is a new hyperparameter. The kernel functions
$k_{\delta_{1}}(\cdot)$ and $k_{\delta_{2}}(\cdot)$ are both modeled by the
Matern-5/2 function (12). As the prior distributions, we assume
$\lambda_{\delta_{1}}\sim\mathcal{HC}(\alpha=3)$,
$\lambda_{\delta_{2}}\sim\mathcal{HC}(\alpha=5)$,
$\ell_{\delta_{1}}\sim\Gamma(\alpha=1,\beta=5)$,
$\ell_{\delta_{2}}\sim\Gamma(\alpha=1,\beta=0.5)$, $\alpha_{{\rm
stall}}\sim\mathcal{HC}(\alpha=2)$, and $x_{\rm
stall}\sim\mathcal{N}(14,0.2)$.

The admissible range of $x=\text{AoA}$ is assumed to be
$[0^{\circ},20^{\circ}]$ over which the experimental and RANS data are
available, see Refs. [2, 10]. The training HF data are taken to be a subset of
size $7$ of the experimental data of Ref. [2]. The rest of the experimental
data are used to validate the predictions of the MFM. For the purpose of
examining the capability of the MFM in a more challenging situation, the
training HF samples are explicitly selected to exclude the range of AoAs where
the stall happens. The DES data of Gilling et al. [10] are available at $7$
$\text{AoA}\in[8^{\circ},19^{\circ}]$ and $5$ different values of TI.
Employing these training data in the HC-MFM (6) and drawing $10^{4}$ MCMC
samples after excluding an extra $5000$ initial samples for burn-in, the
predictions for $C_{L}$ and $C_{D}$ shown in Figure 3(a,c) are obtained. The
expected value of the predictions has a trend similar to that of the
experimental validation data of Ref. [2] and is not diverging towards either
the physically-invalid RANS data or scattered DES data at $\text{AoA}\gtrsim
10^{\circ}$. A more elaborate comparison is made through scatter plot of the
MFM predictions against the validation data in Figure 3(b,d). For both $C_{L}$
and $C_{D}$, the agreement between the predicted mean values with the
validation data at lower AoAs (before stall) is excellent and for most of the
higher AoAs, even near and after the stall, is very good. Due to the scarcity
of the HF training data and also systematic error in the RANS and DES data,
the error bars at the predicted values can be relatively large, as more
evident in the case of $C_{D}$ in Figure 3(d).

(a) |  | (b) |
---|---|---|---
(c) |  | (d) |
Figure 3: (a) Lift coefficient $C_{L}$ and (c) drag coefficient $C_{D}$
plotted against the angle of attack (AoA): the HC-MFM (6) is trained by the
experimental data of Ref. [2] (yellow circles), as well as the DES (squares)
and RANS (crosses) data by Gilling et al. [10]. The DES are performed in Ref.
[10] with the resolved turbulence intensities $\text{TI}=0\%,0.1\%,0.5\%,1\%$,
and $2\%$ at the inlet. The validation data (red triangles) are also taken
from the experiments of Ref. [2]. The mean prediction by the HC-MFM (6) is
represented by the solid line along with associated $95\%$ confidence
interval. (b) $C_{L}$, (d) $C_{D}$ predictions by HC-MFM plotted against the
validation data. The error bars represent the $95\%$ CI.

Figure 4 shows the posterior densities of different parameters appearing in
the MFM constructed for $C_{L}$ and $C_{D}$. As expected, the distribution of
the kernels’ hyperparameters varies between the two QoIs. But more
importantly, the posterior distributions of $x_{\rm stall}$ and calibrated TI
are also dependent on the QoI. This clearly shows the suitability of the
present class of MFMs in which calibration of the parameters of different
fidelities is performed as a part of constructing the MFM for a given QoI. The
alternative strategy, which is common in practice (e.g. for co-Kriging models
without calibration), is to calibrate the LF models against the HF data of one
of the QoIs and then run the calibrated LF model to make realizations of all
QoIs. However, given the present results, this strategy seems clearly less
efficient and leads to inferior quality of prediction.

|  |  |  |
---|---|---|---|---
(a) |  |  |  |
|  |  |  |
(b) |  |  |  |
Figure 4: Posterior PDFs of the calibration parameters and hyperparameters of
the GPs appearing in the HC-MFM (6) for (a) $C_{L}$ and (b) $C_{D}$.
Associated training data and predictions are shown in Figure 3. In the plots
of $\ell_{g}$, the blue and red histograms are corresponding to the AoA and
TI, respectively.

As a general goal, an MFM constructs a surrogate for the QoIs in the space of
the design/controlled parameters aiming for the surrogate outputs to be as
close as possible to the HF data. In this regard, the MFMs facilitate applying
different types of sample-based UQ techniques and optimization, see e.g. [37].
In connection with the present example, consider a UQ forward problem to
estimate the stochastic moments of $C_{L}$ and $C_{D}$ due to the variation of
the AoA. For instance, assume $\text{AoA}\sim\mathcal{N}(15,0.1)$ degrees.
This results in the following estimations for the expectation and variance of
$C_{L}$ and $C_{D}$ with the associated $95\%$ confidence intervals:
$\mathbb{E}_{x}[C_{L}]=1.1775\pm 0.1037$, $\mathbb{V}_{x}[C_{L}]=2.6675\times
10^{-5}\pm 4.5096\times 10^{-5}$, $\mathbb{E}_{x}[C_{D}]=5.6891\times
10^{-2}\pm 2.5722\times 10^{-2}$, and $\mathbb{V}_{x}[C_{D}]=1.0618\times
10^{-5}\pm 9.5994\times 10^{-6}$. Note that without the HC-MFM, and only based
on the data of RANS or/and DES, such estimations would be at best inaccurate,
but in general impossible to make.

###  3.4 Effect of geometrical uncertainties in the periodic hill flow

In this last flow case, we consider the turbulent flow over periodic hills
with geometrical uncertainties, see the sketch in Figure 5. The outline in
blue corresponds to the configuration studied in several prior works
(hereafter, baseline geometry), for example by Fröhlich et al. [6] and more
recently by Gao et al. [7]. The latter reference can be consulted for a good
overview of other previous efforts. The shape of the hill is defined by six
segments of third-order polynomials, see e.g. [42] for the exact definition.
The flow Reynolds number is $5600$.

In Ref. [42], a parameterization of the geometry was introduced by scaling the
length of the hill. Particularly, the authors performed a series of DNS for
$L_{x}/h=3.858\alpha+5.142$, where $L_{x}$ is the length of the geometry, $h$
is the height of the hill, and $\alpha$ is a parameter. The value $\alpha=1$
corresponds to the baseline geometry. The corresponding DNS data set for
several values of $\alpha$ has been made publicly on Github, which was
extended in 2021 with additional data introducing a new parameter $\gamma$:

$L_{x}/h=3.858\alpha+5.142\gamma\ .$ (19)

The effect of $\alpha$ and $\gamma$ on the geometry is illustrated in Figure 5
with red and black curves, respectively. The purpose of the present example is
to demonstrate how the HC-MFM can be used to economically assess the effect of
uncertain parameters $\alpha$ and $\gamma$ on the flow QoIs. To that end we
combine the DNS data discussed above from Ref. [42] with data from RANS
simulations performed by us using ANSYS Fluent v19.5 [1]. Particularly, for
the DNS we use the data for $\alpha\in\\{0.5,1.0,1.5\\}$ and
$\gamma\in\\{0.4166,1.0,1.5834\\}$. The same values are also used for the
RANS, complemented by two additional samples for both $\alpha$ and $\gamma$
for which DNS results are not available. These are selected based on the
Gauss-Legendre quadrature rule and are equal to $\alpha=\\{0.702,1.297\\}$ and
$\gamma=\\{0.653,1.347\\}$. Therefore, there is a total of $5\times 5$ samples
over the space of $\alpha-\gamma$ corresponding to which RANS simulations are
performed. For the uncertainty propagation and sensitivity analysis, see
below, we assume $\alpha$ and $\gamma$ to be independent, and
$\alpha\sim\mathcal{U}[0.448,1.552]$ and $\gamma\sim\mathcal{U}[0.356,1.644]$.

Figure 5: The geometry of the periodic hill simulations, illustrating the
effect of parameters $\alpha$ and $\gamma$ using three sets of values for
them.

The standard $k$ \- $\omega$ turbulence model is used in the RANS simulations,
as defined in [1] based on the work of Wilcox [41]. The available low-
Reynolds-number correction to the model was not used. The model depends on a
number of parameters which are listed in [1, p. 61]. It can be shown [41, p.
134] that the parameters $\alpha_{\infty}$, $\beta^{*}_{\infty}$, $\sigma$ and
$\beta_{i}$ are coupled to the von Karmán coefficient $\kappa$ through the
following equation:

$\alpha_{\infty}=\beta_{i}/\beta_{\infty}^{*}-\sigma\kappa^{2}/\sqrt{\beta_{\infty}^{*}}.$
(20)

To illustrate the automatic calibration capability of the HC-MFM, we assume
$\kappa$ to be uncertain and perform simulations for 5 sample values of
$\kappa\in\\{0.348,0.367,0.4,0.433,0.452\\}$, which follow the Gauss-Legendre
quadrature rule. To prescribe the desired value of $\kappa$, we set
$\alpha_{\infty}=0.52$, $\sigma=0.5$, $\beta_{i}=0.0708$ (as suggested by
Wilcox [41, p. 135]), and manipulate $\beta_{\infty}^{*}$ according to (20).
Note that the considered training samples include the standard choice
$\beta_{\infty}^{*}=0.09$, corresponding to $\kappa=0.4$. Using $5$ samples
for each of $\alpha$, $\gamma$ and $\kappa$ and using a tensor-product rule, a
total of 125 RANS simulations were performed for this study.

For RANS simulations, quadrilateral cells were used to discretize the
computational domains, with the total grid size ranging from $\approx 150\cdot
10^{3}$ to $\approx 500\cdot 10^{3}$ cells, depending on the domain size as
defined by $\alpha$ and $\gamma$. Since the selected turbulence model requires
accurate resolution of the boundary layer, the selection of the mesh size was
guided by the discretization in the corresponding DNS case [42]. Specifically,
we adopted the same number of cells in the streamwise and wall-normal
directions as in the DNS, and applied size grading in the wall-normal
direction to ensure low values of $y^{+}$. This means that in the streamwise
direction the mesh may be unnecessarily fine for RANS, but since these
simulations are two-dimensional and steady-state, they are still negligibly
cheap compared to the DNS.

Second-order numerical schemes were used to discretize the RANS equations in
ANSYS Fluent [1]. Specifically, for the convective fluxes a second-order
upwind scheme was used. The coupled solver was used for pressure-velocity
coupling, which did an excellent job at converging the simulations.

####  3.4.1 Creating ground truth and verifying the model

Using all the $9$ DNS data sets available from Ref. [42], the response surface
of the QoIs at the space of $\alpha-\gamma$, and associated PDF of the QoIs
can be estimated. These will be used as the ground truth or reference to
evaluate the performance of the HC-MFM in the following analyses. The
interpolation from the DNS samples to an arbitrary mesh covering the whole
admissible range of $\alpha$ and $\gamma$ can be done using polynomial-based
methods such as PCE (used here) or Lagrange interpolation, as well as GP
regression. Based on the data available from Ref. [42] and the performed RANS
simulations, different QoIs can be considered. Hereafter, to demonstrate the
power of the method, we take the normalized height of the separation bubble,
$H_{\rm bubble}/h$ at the streamwise location $x/h=2.5$ as the QoI.
Alternatively other locations $x/h$ as well as different flow quantities could
be considered. The response surface of the QoI and associated PDF are
illustrated in Figure 6. Based on the pattern of the isolines, we can observe
that the parameter $\alpha$ exhibits a stronger influence on the QoI than
$\gamma$. This can be quantitatively confirmed via the values of the total
Sobol indices [38] as reported in Table 1. The resulting PDF is bimodal with
two peaks around the most probable observed values of $H_{\rm bubble}/h$ at
$x/h=2.5$.

|
---|---
Figure 6: (Left) Isolines of the response surface and (Right) PDF of $H_{\rm
bubble}/h$ at $x/h=2.5$ due to the variation of $\alpha$ and $\gamma$ using
all the $9$ DNS data of Ref. [42] (represented by the symbols in the left
plot). These plots are considered as ground truth or reference for evaluate
the performance of the multifidelity model.

The reference posterior distribution of $\kappa$ as the RANS calibration
parameter can now be inferred. To this end, the HC-MFM described in Section
3.4.2 is constructed using $9$ DNS data sets of Ref. [42] and $125$ RANS
simulations. The prior distribution of $\kappa$ is taken to be uniform over
the range of $[0.3,0.5]$. This non-informative prior distribution removes any
bias towards any particular value in the distribution of $\kappa$. Through a
Bayesian inference via an MCMC method, the sample posterior distribution of
$\kappa$ shown in Figure 7 (left) is obtained. Note that this calibration is
in fact a pure UQ inverse problem, see e.g. [37], where all the RANS data are
utilized to construct a surrogate for $\kappa$, and the DNS data are used as
training data to infer the distribution of $\kappa$. The estimated mean and
standard deviation of the posterior distribution of $\kappa$ are $0.47087$ and
$0.05387$, respectively. As compared to the standard value $0.4$ being used in
the literature, the estimated mean is somewhat larger. However, from a
physical point of view, one would not expect an accurate value of $\kappa$ for
this type of flow due to the separated nature and the relatively low Reynolds
number.

Another advantage of using all the available RANS and DNS data in the HC-MFM
is that the implementation (algorithm and coding) of the HC-MFM can be
verified. As shown in Figure 7 (Right), the prediction of the HC-MFM
constructed by combining all DNS and RANS data sets completely agree with the
predictions of the single-fidelity model based on only the DNS data. Note that
the predicted marginal PDF of the QoI in the HC-MFM is the same as the
reference PDF in Figure 6.

|
---|---
Figure 7: (Left) The sample posterior PDF of $\kappa$ and (right) sample joint
PDF of $H_{\rm bubble}/h$ at $x/h=2.5$ obtained from the HC-MFM using all $9$
DNS data sets of Ref. [42] along with the $125$ RANS simulations performed in
the present study. The RANS simulations are performed using $5$ samples of
$\kappa$ equal to $0.348$, $0.367$, $0.4$, $0.433$, and $0.452$. The marginal
PDFs on the top and right axes are found using the kde (kernel density
estimation) method.

####  3.4.2 Application of the HC-MFM

Adopting the general notation of Section 2, the HC-MFM for the present example
can be written as,

$\displaystyle\begin{cases}y_{M_{1}}(\mathbf{x}_{i})\quad\,\,={\hat{f}(\mathbf{x}_{i},\theta_{2_{i}})}+{\hat{\delta}(\mathbf{x}_{i})}+\varepsilon_{1_{i}}\,&,\quad
i=1,2,\cdots,n_{1}\\\
y_{M_{2}}(\mathbf{x}_{i})={\hat{f}(\mathbf{x}_{i},t_{2_{i}})}+\varepsilon_{2_{i}}\,&,\quad
i=1+n_{1},2+n_{1},\cdots,n_{2}+n_{1}\end{cases}\,,$ (21)

where $M_{1}$ and $M_{2}$ denote DNS and RANS, respectively, the design
parameters are $\mathbf{x}_{i}=(\alpha_{i},\gamma_{i})$, and, $t_{2}$ and
$\theta_{2}$ refer to the simulated and calibrated instances of $\kappa$,
respectively. The kernel of $\hat{f}(\mathbf{x},t_{2})$ is taken to be the
exponentiated quadratic function (11), while for $\hat{\delta}(\mathbf{x})$,
the Matern-5/2 function (12) is employed. For the hyperparameters, the
following prior distributions are considered:
$\kappa\sim\mathcal{U}[0.3,0.5]$, $\lambda_{f}\sim\mathcal{HC}(\alpha=5)$,
$\ell_{f_{\mathbf{x}}},\ell_{f_{t_{2}}}\sim\Gamma(\alpha=1,\beta=5)$,
$\lambda_{\delta}\sim\mathcal{HC}(\alpha=1)$, and
$\ell_{\delta_{\mathbf{x}}}\sim\Gamma(\alpha=1,\beta=1)$. In this example, the
uncertainty in the DNS and RANS data is neglected. Despite this, when
implementing the model in PyMC3 [36], the prior of the Gaussian noise standard
deviation is set to be $\sigma\sim\mathcal{HC}(\alpha=5)$ (same for both
fidelities). But as expected, the mean and standard deviation of the posterior
distribution of $\sigma$ are obtained to be approximately zero.

The hyperparameters will be inferred from the combined set of $n_{1}$ DNS and
$n_{2}$ RANS data. In the analyses to follow, we use all the RANS simulations
as low-fidelity (LF) data, therefore $n_{2}=125$. Two subsets of the DNS data
of Ref. [42] with size $n_{1}=4$ and $5$ are taken to be the high-fidelity
(HF) data. Combining these two HF data sets with the LF data, Case-A and
Case-B data sets are obtained for multifidelity modeling. Figure 8 represents
the samples of these two cases in the space of $\alpha-\gamma$ parameters.

|
---|---
Figure 8: Schematic representation of the samples from $\alpha$ and $\gamma$
corresponding to the low-fidelity (LF), $\times$, high-fidelity (HF), $\Box$,
and all available DNS data from Ref. [42], $\circ$. In the text, the left and
right plots are referred to as Case-A and Case-B, respectively. Note that for
both cases, there are 5 samples for $\kappa$ associated to each of the LF
samples represented here.

Before constructing the HC-MFM, it is important to look at the LF data. In
Figure 9, the PDF of the QoI due to the variation of $\alpha$ and $\gamma$ for
different training samples of $\kappa$ is represented. The two expected yet
important observations are that the PDF of the QoI is significantly influenced
by the value of $\kappa$ used in the RANS simulations, and the fact that the
PDFs are much different from the ground truth PDF shown in Figure 6. The
larger influence of $\kappa$ compared to $\alpha$ and $\gamma$ is also
reflected in the associated Sobol indices [38], as reported in Table 1. Note
that the PDF of the QoI considering the simultaneous variation of $\alpha$,
$\gamma$, and $\kappa$ is shown in Figure 9(f).

|  |
---|---|---
(a) | (b) | (c)
|  |
(d) | (e) | (f)
Figure 9: (a-e) PDF of $H_{\rm bubble}/h$ at $x/h=2.5$ due to the variation of
$\alpha$ and $\gamma$ using the RANS data simulated with $\kappa$ equal to
$0.348$, $0.367$, $0.400$, $0.433$, $0.452$, respectively. Note that $5\times
5$ samples are taken from the $\alpha-\gamma$ space at each of these
constant-$\kappa$ simulations. The PDF in (f) is obtained using all the
$5\times 5\times 5$ samples from $\alpha$, $\gamma$ and $\kappa$.

The response surface and PDF of the QoI obtained from only the HF data of
Case-A and Case-B are illustrated in Figure 10. For case-A with $n_{1}=4$ HF
data, the PDF of the QoI has a plateau which makes it clearly different from
the reference PDF in Figure 6. By adding only one more DNS data point and
obtaining Case-B, the response surface becomes more similar to the reference,
however, the associated PDF is still single mode. The improved predictions
through the application of the HC-MFM are shown in Figure 11. Compared to the
HF-data in Figure 10, the PDFs of the QoI clearly exhibit a second peak for
the values between $0.5$ and $0.6$. This peak has been introduced by adding
the LF data, see Figure 9(f), and this is, in fact, the task for the
multifidelity model to adjust the involved hyperparameters such that the
fusion of the data at the two fidelities leads to a PDF similar to the
reference. Clearly, for Case-B with only 5 DNS data samples included, the PDF
and response surface of the QoI are very close to the ground truth in Figure 6
(9 DNS). This can also be confirmed by plotting the associated HC-MFM
predictions against the reference data at all test points in the
$\alpha-\gamma$ plane, see Figure 12. The joint PDF of these two sets of data
for both considered cases is narrow, specifically for Case-B, and hence
implies a low point-to-point deviation of the predictions from the reference
values. It is also interesting to look at the posterior distribution of the
RANS calibration parameter $\kappa$. It is not surprising that the resulting
PDFs from the multifidelity data sets Case-A and Case-B are different in spite
of having the same uniform prior distribution. Two important observations here
are the following: First, in contrast to the previous examples in the present
study, even with a small number of HF-data, i.e. Case-A, a significantly
informative PDF for the LF calibration parameter is obtained. Second, for
Case-B, the posterior distribution of $\kappa$ is very similar to the
reference case where $9$ DNS data sets are used, see Figure 7. Thus, depending
on the case, the HC-MFM is capable of calibrating the model parameters in a
fairly accurate way at the same time of constructing an accurate predictive
model. This somewhat challenges the conclusion which could be drawn from the
previous examples in the present study and also in Goh et al. [11], where the
priority of the HC-MFM is found to make accurate predictions for the QoI
rather than providing accurate calibration of the fidelity-specific
parameters.

|
---|---
|
Figure 10: (Top) Isolines of the response surface and (bottom) PDF of $H_{\rm bubble}/h$ at $x/h=2.5$ due to the variation of $\alpha$ and $\gamma$ using the HF data of (left) Case-A and (right) Case-B. The data are taken from the DNS of Ref. [42] and are specified by dots in the top plots.  |
---|---
|
Figure 11: (Top) Isolines of the response surface and (bottom) PDF of $H_{\rm bubble}/h$ at $x/h=2.5$ due to the variation of $\alpha$ and $\gamma$ obtained from the HC-MFM with the data of (left) Case-A and (right) Case-B. In each of the plots in the bottom row, the PDF resulting from the HC-MFM is compared to the PDFs of the ground truth (Ref., Figure 6), low-fidelity data (LFM, Figure 9), and high-fidelity data (HFM, Figure 10). |
---|---
|
Figure 12: (Top) Joint and marginal PDFs, and (bottom) sample posterior
distribution of $H_{\rm bubble}/h$ at $x/h=2.5$ for (left) Case-A and (right)
Case-B data sets. In the top plot, the contours belong to the joint PDF with
associated values specified in the color bar.

To conclude the periodic hill example, we can quantify stochastic moments of
the QoI as well as the Sobol sensitivity indices [38] due to the uncertainty
in $\alpha$ and $\gamma$. These UQ measures are integral quantities over the
admissible range of the parameters and can be computed using the reference
data (all available DNS), LF, HF, and MF data sets. Note that for the LF data,
the uncertainty in $\kappa$ is also taken into account. Noting the parameters
$\alpha$, $\gamma$, and $\kappa$ are uniformly distributed and independent
from each other, the results summarized in Table 1 are obtained using the
generalized polynomial chaos expansion [43] for the reference and LF data
sets, and the Monte Carlo method for the multifidelity cases. All the UQ
analyses have been performed using UQit [33]. In general, for all cases but
the pure LF data, the prediction of the mean and standard deviation of the QoI
are close to the reference. For the total Sobol indices, the closest estimates
to the reference values is obtained from the HC-MFM applied to Case-B, and on
the second rank, Case-A. Noting the improvement of the Sobol indices accuracy
in each of the multifidelity cases compared to the estimates from the
associated HF data, the effectiveness of the HC-MFM is once again confirmed.
This is an important outcome considering the forward UQ problems and global
sensitivity analyses are of most relevance in CFD applications.

Table 1: Estimated mean, standard deviation, and total Sobol indices of the QoI $R=H_{\rm bubble}/h$ at $x/h=2.5$ due to the uncertainty in $\alpha$ and $\gamma$. For the LF (RANS) data the uncertainty due and sensitivity with respect to parameter $\kappa$ is also included. For the Case-A and Case-B data sets used for multifidelity modeling, see Figure 8.  |  | Moments due to $\alpha$, $\gamma$ | Total Sobol Indices of $R$ w.r.t.
---|---|---|---
| Data set | $\mathbb{E}[R]$ | $\mathbb{S}[R]$ | $\alpha$ | $\gamma$ | $\kappa$
| Reference | 0.48104 | 0.05027 | 0.94551 | 0.08094 | –
| Low-Fidelity (LF) | 0.51999 | 0.01666 | 0.44275 | 0.02689 | 0.60623
Case-A | High-Fidelity (HF) | 0.50047 | 0.05045 | 0.81448 | 0.19144 | –
| Multifidelity (MF) | 0.49209 | 0.05238 | 0.88133 | 0.12680 | –
Case-B | High-Fidelity (HF) | 0.46981 | 0.05198 | 0.82327 | 0.18237 | –
| Multifidelity (MF) | 0.48437 | 0.05215 | 0.91355 | 0.10097 | –

###  3.5 Keeping the RANS Parameter Fixed

In all the examples in the present study, the fidelity-specific calibration
parameters are involved in the multifidelity modeling and posterior
distributions for them are learned during the construction of the MFM. But,
the methodology behind the HC-MFM described in Section 2 is general and
flexible to be directly applicable to the cases where the fidelity-specific
parameters are kept fixed. To demonstrate this, let us apply the HC-MFM to the
example of the periodic hill and use the RANS data at constant values of
$\kappa$, the RANS modeling parameter in Eq. (20). We use the Case-B data sets
shown in Figure 8 which means having $5$ and $25$ samples for the DNS and RANS
simulations, respectively, in the $\alpha-\gamma$ space. The validation of the
PDF of the QoI, $H_{\rm bubble}/h$ at $x/h=2.5$, of the HC-MFM for two values
of $\kappa$ is represented in Figure 13. Adopting $\kappa=0.433$
($\beta_{\infty}^{*}=0.084$) shows a clear improvement in the predictions
compared to the standard value $\kappa=0.4$ ($\beta_{\infty}^{*}=0.09$). This
observation is consistent with the posterior PDF of $\kappa$ in Figure 12,
where the mode of the distribution is higher than $0.4$. From this test, not
only the validity of the HC-MFM for fixed values of the fidelity-specific
parameters is confirmed, but also the fact that such accuracy-controlling
parameters should be actively part of the data generation and hence the
construction of HC-MFM is emphasized.

Another important point is that the good predictive accuracy in Figure 13 is
obtained despite the poor correlation between the DNS and RANS data. This is
because of accurate construction of the model discrepancy term in Eq. (21) and
also accurate estimation of various hyperparameters. It is also noteworthy
that the present example is comparable to what is performed by Voet et al.
[39] using $7$ DNS and $30$ RANS data sets using different multifidelity
modeling approaches while fixing the value of the coefficients in the RANS
closure model. However, a direct comparison between the two studies is not
possible since in Ref. [39], the plots of the MFM predictions versus reference
values of the QoI as in Figure 13 are not provided.

|
---|---
Figure 13: Joint and marginal PDFs of $H_{\rm bubble}/h$ at $x/h=2.5$ for
Case-B data sets using fixed values of $\kappa$ equal to (left) $0.4$ and
(right) $0.433$. Note that the PDF of the QoI due to the variation of $\alpha$
and $\beta$ corresponding to theses $\kappa$ values is plotted in Figure 9 (c)
and (d), respectively.

###  3.6 Impact of Replacing the MCMC by a Point-Estimator

In all the examples presented in this study, the HC-MFM is constructed using
the MCMC method to draw samples from the posterior distribution of the
hyperparameters and calibration parameters, i.e. $\bm{\beta}$ and
$\bm{\theta}$, respectively, in the Bayes formula (13). Similarly, to predict
the sample distribution of the QoI, direct samples from these posterior
distributions are used in the HC-MFM. As an alternative to these sample-based
methods within the Bayesian framework, point estimators such as maximum
a-posteriori probability (MAP) and maximum likelihood estimators (MLE) can be
adopted. The point-estimated values are considered to be the representatives
of the corresponding distribution. Note that the use of the uniform priors in
Eq. (13) makes the MAP estimations identical to MLE’s. Our investigations
showed that using point estimators instead of the MCMC method, could
deteriorate the accuracy of the HC-MFM predictions, independent from how the
LF and HF data are combined.

For instance, according to Figure 14, the PDF of the QoIs predicted by the HC-
MFM using the MAP estimator is significantly worse than what is given by the
MCMC method as shown in Figure 12. This is an important message of the present
study noting that all previous multifidelity studies in the literature
relevant to the fluid flows have been based on using the point estimators, see
e.g. Voet et al. [39] where the MLE is adopted.

|
---|---
Figure 14: (Top) Joint and marginal PDFs, and (bottom) sample posterior
distribution of $H_{\rm bubble}/h$ at $x/h=2.5$ for (left) Case-A and (right)
Case-B data sets. Here a MAP estimator is used to construct the HC-MFM, in
contrast to Figure 12 and the rest of the examples in the present study which
are obtained using an MCMC method.

##  4 Summary and Conclusions

The Bayesian hierarchical multifidelity model with automatic calibration (HC-
MFM) developed by Goh et al. [11] is adapted to several examples relevant to
wall-bounded turbulent flows. The HC-MFM is general, accurate, applicable to
an arbitrary number of fidelity levels and well-suited to the simulations of
turbulent flows since as a part of the MFM construction the fidelity-specific
parameters can be automatically calibrated using the training data of higher
fidelities. This is an important feature noting that in all approaches for
simulation of turbulence, different numerical and modeling uncertain
parameters can influence the accuracy of the QoIs. Because of using the
Gaussian processes, the predictions made by the HC-MFM are accompanied with
confidence intervals. Moreover, it is possible to incorporate the
observational uncertainties in the data at all fidelity levels, and hence
perform various UQ analyses for combinations of different types of
uncertainties, see Ref. [34].

Based on the examples, the following main conclusions can be made. 1\. For a
fixed number of HF training data, the HC-MFM prioritizes the prediction of
QoIs so that they become as close as possible to the HF validation data, while
the posterior distributions of the calibration parameters are found to be
accurate only if sufficiently many LF training data are provided. A similar
conclusion was drawn by Goh et al. [11] by systematically increasing the
amount of both HF and LF data. For the periodic hill subject to geometrical
uncertainties, Section 3.4.2, fixing the number of the LF data and considering
two sets of HF data, the posterior distribution of the RANS (LF) parameter was
found to be close to what would be obtained by using all the available HF
data. This, again, confirms the importance of providing sufficient LF training
data through well exploring the space of the design and calibration
parameters. 2\. When there are more than one QoI, the posterior distribution
of the calibration parameters may depend on the QoI, see Section 3.3.
Therefore, the calibration parameters are more numerical than physical and
hence, predictions by the HC-MFM for a QoI can be more accurate than the case
of a-priori calibrating the low-fidelity models against high-fidelity data of
another QoI (an example is calibrating a RANS closure model by the HF data of
the lift coefficient, and then using the calibrated model in a simulation
aiming for the drag coefficient with optimal accuracy). 3\. As show in Section
3.6, the method for estimating the hyperparameters and parameters in the HC-
MFM can significantly affect the resulting predictive accuracy. In fact, the
MCMC sampling method is shown to result in more accurate predictions compared
to a point estimator like MAP. This important point is usually overlooked in
most of, if not all, the previous studies regarding the multifidelity modeling
in CFD. 4\. If the fidelity-specific calibration parameters are kept fixed,
the HC-MFM is still applicable without any need to modifying its general
formulation. Obviously, the predictive accuracy of the model will depend on
the validity of the value chosen for such parameters when generating the
training data. The success of the HC-MFM relies on the accurate modeling of
the discrepancy terms between different fidelities, and also the use of the
MCMC methods to find optimal values for underlying hyperparameters through
exploration of the parameter space.

The present study may be extended in several directions. For instance, in
addition to the scalar QoIs, spatio–temporal fields can be considered in the
HC-MFM, making it possible to predict full flow fields by combining different
fidelities. Such an approach may be particularly interesting as an alternative
to the more black-box machine-learning tools when it comes to super-resolution
and related methods. Another potential extension is in the combination of the
HC-MFM with a Bayesian optimization for CFD applications and turbulent flow
problems [24]. In this case, the surrogate for the optimizer is based on the
MFM, and is thus potentially cheaper to evaluate during the optimization
process. In particular, applications in flow control for turbulence, where the
main computational time lies in the evaluation of the objective function (i.e.
the CFD solver), may greatly benefit from a well-calibrated multifidelity
model.

## Acknowledgments

This work has been supported by the EXCELLERAT project, which has received
funding from the European Union’s Horizon 2020 research and innovation
programme under grant agreement No 823691. Additional funding was provided by
the Knut and Alice Wallenberg Foundation (KAW). The computations were enabled
by resources provided by the Swedish National Infrastructure for Computing
(SNIC), partially funded by the Swedish Research Council through grant
agreement no. 2018-05973.

##  Appendix A PDF of A Set of Standard Distributions

For the prior distribution of the parameters and hyperparameters of the
multifidelity models in Section 3, a set of standard distributions was used.
Table 2 summarizes the PDF and associated support of such distributions.

Table 2: PDF and associated support of the standard distributions used in Section 3. Distribution | PDF | Support in $x$
---|---|---
Uniform ($\mathcal{U}$) | $\rho(x;\alpha,\beta)=(\beta-\alpha)^{-1}$ | $[\alpha,\beta]$
Gaussian ($\mathcal{N}$) | $\rho(x;\alpha,\beta)=(2\pi\beta^{2})^{-1/2}\exp\left(-\frac{1}{2}(\frac{x-\alpha}{\beta})^{2}\right)$ | $(-\infty,\infty)$
Half-Cauchy ($\mathcal{HC}$) | $\rho(x;\alpha)=2\left[\pi\alpha\left(1+\left(\frac{x}{\alpha}\right)^{2}\right)\right]^{-1}$ | $[0,\infty)$
Gamma ($\Gamma$) | $\rho(x;\alpha,\beta)=\beta^{\alpha}x^{\alpha-1}\exp(-\beta x)/\Gamma(\alpha)$ | $(0,\infty)$

## References

  * [1] ANSYS Inc.  ANSYS Fluent Theory Guide. Release 2019R3, 2019.
  * [2] F. Bertagnolio.  NACA0015 Measurements in LM Wind Tunnel and Turbulence Generated Noise.  Number 1657(EN) in Denmark. Forskningscenter Risoe. Danmarks Tekniske Universitet, 2008.
  * [3] H. R. Fairbanks, A. Doostan, C. Ketelsen, and G. Iaccarino.  A low-rank control variate for multilevel Monte Carlo simulation of high-dimensional uncertain systems.  J. Comput. Phys, 341:121–139, 2017.
  * [4] J. Fatou Gomez.  Multi-fidelity Co-Kriging Optimization using Hybrid Injected RANS and LES.  PhD thesis, Delft University of Technology, 2018.
  * [5] A. I. Forrester, A. Sóbester, and A. J. Keane.  Multi-fidelity optimization via surrogate modelling.  Proceedings of the Royal Society A, 463(2088):3251–3269, 2007.
  * [6] J. Fröhlich, C. P. Mellen, W. Rodi, L. Temmerman, and M. A. Leschziner.  Highly resolved large-eddy simulation of separated flow in a channel with streamwise periodic constrictions.  Journal of Fluid Mechanics, 526:19–66, 2005.
  * [7] W. Gao, W. Cheng, and R. Samtaney.  Large-eddy simulations of turbulent flow in a channel with streamwise periodic constrictions.  Journal of Fluid Mechanics, 900:A43, 2020.
  * [8] G. Geraci, M. S. Eldred, and G. Iaccarino.  A multifidelity multilevel Monte Carlo method for uncertainty propagation in aerospace applications.  In 19th AIAA Non-Deterministic Approaches Conference, 2017.
  * [9] M. B. Giles.  Multilevel Monte Carlo path simulation.  Operations Research, 56(3):607–617, 2008.
  * [10] L. Gilling, N. Sørensen, and L. Davidson.  Detached eddy simulations of an airfoil in turbulent inflow.  47th AIAA Aerospace Sciences Meeting, 24(AIAA 2009-270):1–13, 2009\.
  * [11] J. Goh, D. Bingham, J. P. Holloway, M. J. Grosskopf, C. C. Kuranz, and E. Rutter.  Prediction and computer model calibration using outputs from multifidelity simulators.  Technometrics, 55(4):501–512, 2013.
  * [12] R. B. Gramacy.  Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences.  Chapman Hall/CRC, Boca Raton, Florida, 2020.
  * [13] L. L. Gratiet and J. Garnier.  Recursive co-kriging model for design of computer experiments with multiple levels of fidelity.  International Journal for Uncertainty Quantification, 4(5):365–386, 2014.
  * [14] Z.-H. Han and S. Görtz.  Hierarchical kriging model for variable-fidelity surrogate modeling.  AIAA Journal, 50(9):1885–1896, 2012.
  * [15] D. Higdon, M. Kennedy, J. C. Cavendish, J. A. Cafeo, and R. D. Ryne.  Combining field data and computer simulations for calibration and prediction.  SIAM J. Sci. Comput., 26(2):448–466, 2004.
  * [16] M. D. Hoffman and A. Gelman.  The No-U-Turn Sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo.  Journal of Machine Learning Research, 15(47):1593–1623, 2014.
  * [17] K. Iwamoto, Y. Suzuki, and N. Kasagi.  Reynolds number effect on wall turbulence: toward effective feedback control.  International Journal of Heat and Fluid Flow, 23(5):678–689, 2002\.
  * [18] L. Jofre, G. Geraci, H. R. Fairbanks, A. Doostan, and G. Iaccarino.  Exploiting multi-fidelity strategies to quantify uncertainty in irradiated particle-laden turbulent flow.  Proceedings of the 7th European Conference on Computational Fluid Dynamics, pages 1–12, 2018.
  * [19] M. C. Kennedy and A. O’Hagan.  Predicting the output from a complex computer code when fast approximations are available.  Biometrika, 87(1):1–13, 2000.
  * [20] M. C. Kennedy and A. O’Hagan.  Bayesian calibration of computer models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63(3):425–464, 2001.
  * [21] J. Larsson, S. Kawai, J. Bodart, and I. Bermejo-Moreno.  Large eddy simulation with modeled wall-stress: recent progress and future directions.  Bulletin of the JSME, 3(1), 2016.
  * [22] M. Lee and R. D. Moser.  Direct numerical simulation of turbulent channel flow up to $\mathit{Re}_{{\it\tau}}\approx 5200$.  J. Fluid Mech., 774:395–415, 2015.
  * [23] B. Matern.  Spatial Variation.  Lecture Notes in Statistics. Springer New York, 1986.
  * [24] Y. Morita, S. Rezaeiravesh, N. Tabatabaei, R. Vinuesa, K. Fukagata, and P. Schlatter.  Applying Bayesian optimization with Gaussian process regression to computational fluid dynamics problems.  Journal of Computational Physics, 449:110788, 2022.
  * [25] L. W.-T. Ng and M. Eldred.  Multifidelity uncertainty quantification using non-intrusive polynomial chaos and stochastic collocation.  53rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, (AIAA 2012-1852):1–17, 2012.
  * [26] P. S. Palar, T. Tsuchiya, and G. T. Parks.  Multi-fidelity non-intrusive polynomial chaos based on regression.  Computer Methods in Applied Mechanics and Engineering, 305:579 – 606, 2016.
  * [27] R. Pasupathy, B. W. Schmeiser, M. R. Taaffe, and J. Wang.  Control-variate estimation using estimated control means.  IIE Transactions, 44(5):381–385, 2012.
  * [28] B. Peherstorfer, K. Willcox, and M. Gunzburger.  Survey of multifidelity methods in uncertainty propagation, inference, and optimization.  SIAM Review, 60(3):550–591, 2018.
  * [29] P. Perdikaris, D. Venturi, J. O. Royset, and G. E. Karniadakis.  Multi-fidelity modelling via recursive co-kriging and Gaussian-Markov random fields.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471(2179):20150018, 2015.
  * [30] C. E. Rasmussen and C. K. I. Williams.  Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning).  The MIT Press, 2005.
  * [31] W. C. Reynolds and G. W. M. F. Tiederman.  Stability of turbulent channel flow, with application to Malkus’s theory.  J. Fluid Mech., 27(part 2):253–272, 1967.
  * [32] S. Rezaeiravesh and M. Liefvendahl.  Effect of grid resolution on large eddy simulation of wall-bounded turbulence.  Physics of Fluids, 30(5):055106, 2018.
  * [33] S. Rezaeiravesh, R. Vinuesa, and P. Schlatter.  UQit: A python package for uncertainty quantification (uq) in computational fluid dynamics (cfd).  Journal of Open Source Software, 6(60):2871, 2021.
  * [34] S. Rezaeiravesh, R. Vinuesa, and P. Schlatter.  An uncertainty-quantification framework for assessing accuracy, sensitivity, and robustness in computational fluid dynamics.  Journal of Computational Science, 62:101688, 2022.
  * [35] P. Sagaut, S. Deck, and M. Terracol.  Multiscale and Multiresolution Approaches in Turbulence.  Imperial College Press, 2nd edition, 2013.
  * [36] J. Salvatier, T. V. Wiecki, and C. Fonnesbeck.  Probabilistic programming in Python using PyMC3.  PeerJ Computer Science, 2:e55, 2016.
  * [37] R. C. Smith.  Uncertainty Quantification: Theory, Implementation, and Applications.  Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2013.
  * [38] I. Sobol.  Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates.  Mathematics and Computers in Simulation, 55(1):271 – 280, 2001\.  The Second IMACS Seminar on Monte Carlo Methods.
  * [39] L. J. Voet, R. Ahlfeld, A. Gaymann, S. Laizet, and F. Montomoli.  A hybrid approach combining DNS and RANS simulations to quantify uncertainties in turbulence modelling.  Applied Mathematical Modelling, 89:885 – 906, 2021.
  * [40] H. Weller, G. Tabor, H. Jasak, and C. Fureby.  A tensorial approach to computational continuum mechanics using object-oriented techniques.  Computers in Physics, 12(6):620–631, 1998.
  * [41] D. C. Wilcox.  Turbulence Modeling for CFD.  DCW Industries, Inc. La Canada, California, 2006.
  * [42] H. Xiao, J.-L. Wi, S. Laizet, and L. Duan.  Flows over periodic hills of parameterized geometries: a dataset for data-driven turbulence modeling from direct simulations.  Computers and Fluids, 200:104431, 2020.
  * [43] D. Xiu and G. E. Karniadakis.  The Wiener–Askey polynomial chaos for stochastic differential equations.  SIAM Journal on Scientific Computing, 24(2):619–644, 2002.
  * [44] Y. Yamamoto and Y. Tsuji.  Numerical evidence of logarithmic regions in channel flow at $\mathrm{R}{\mathrm{e}}_{\tau}=8000$.  Phys. Rev. Fluids, 3:012602, 2018.