{'fraction_non_alphanumeric': 0.029686880697582244, 'fraction_numerical': 0.006004756242568371, 'mean_word_length': 4.939383239171375, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Development and Realization of Validation Benchmarks

In the field of modeling, the word validation refers to simple comparisons between model outputs and experimental data. Usually, this comparison constitutes plotting the model results against data on the same axes to provide a visual assessment of agreement or lack thereof. However, there are a number of concerns with such naive comparisons. First, these comparisons tend to provide qualitative rather than quantitative assessments and are clearly insufficient for making decisions regarding model validity. Second, they often disregard or only partly account for existing uncertainties in the experimental observations or the model input parameters. Third, such comparisons can not reveal whether the model is appropriate for the intended purposes, as they mainly focus on the agreement in the observable quantities. These pitfalls give rise to the need for an uncertainty-aware framework that includes a validation metric. This metric shall provide a measure for comparison of the system response quantities of an experiment with the ones from a computational model, while accounting for uncertainties in both. To address this need, we have developed a statistical framework that incorporates a probabilistic modeling technique using a fully Bayesian approach. A Bayesian perspective yields an optimal bias-variance trade-off against the experimental data and provide an integrative metric for model validation that incorporates parameter and conceptual uncertainty. Additionally, to accelerate the analysis for computationally demanding flow and transport models in porous media, the framework is equipped with a model reduction technique, namely Bayesian Sparse Polynomial Chaos Expansion. We demonstrate the capabilities of the aforementioned Bayesian validation framework by applying it to an application for validation as well as uncertainty quantification of fluid flow in fractured porous media.


Abstract
In the field of modeling, the word validation refers to simple comparisons between model outputs and experimental data. Usually, this comparison constitutes plotting the model results against data on the same axes to provide a visual assessment of agreement or lack thereof. While comparisons between model and data are at the heart of any validation procedure, there are a number of concerns with such naive comparisons. First, these comparisons tend to provide qualitative rather than quantitative assessments and are clearly insufficient as a basis for making decisions regarding model validity. Second, naive comparisons often disregard or only partly account for existing uncertainties in the experimental observations or the model input parameters. Third, such comparisons can not reveal whether the model is appropriate for the intended purposes, as they mainly focus on the agreement in the observable quantities.
These pitfalls give rise to the need for an uncertainty-aware framework that includes a validation metric. This metric shall provide a measure for comparison of the system response quantities of an experiment with the ones from a computational model, while accounting for uncertainties in both, in a rigorous way. To address this need, we have developed a statistical framework that incorporates a probabilistic modeling technique using a fully Bayesian approach. A Bayesian perspective on a validation task yields an optimal bias-variance trade-off against the experimental data and provide an integrative metric for model validation that incorporates parameter and conceptual uncertainty. Additionally, to accelerate the validation process for computationally demanding flow and transport models in porous media, the framework is equipped with a model reduction technique, namely Bayesian Sparse Polynomial Chaos Expansion. We demonstrate the capabilities of the aforementioned Bayesian validation framework by applying it to an application for validation as well as uncertainty quantification of fluid flow in fractured porous media.

Introduction
Over the last century, the computational modeling in the field of porous media has witnessed tremendous improvement. After decades of development, the state-of-the-art simulators are now capable of solving coupled partial differential equations governing the complex subsurface multiphase flow system within a practically large spatial and temporal domain. Given the importance of computational modeling, assessment of the reliability of models in light of the purpose of a given simulation is of paramount importance. For this reliability assessment, there exists various complementary measures, such as unit testing of individual components of the computational model, comparisons with analytical solutions or other computational models.
While assessing the simulation results with aforementioned measures could hint us towards the correctness of the model, measuring the computational model performance against experimental data can enhance the confidence in the models. However, comparisons with experimental data can be very cumbersome, since the experiment in question has probably not been designed to meet the objective of the validation of a computational model. Moreover, disagreement with experimental data can have other reasons than the deficiency of the simulation code, such as an insufficient description of the experimental set-up or considerable uncertainties associated with the measurements. On the other hand, pure code intercomparison studies cannot ensure that a successfully participating model indeed maps the reality.
The goal of this PhD project is to formulate and conduct benchmarks, which assist in the uncertainty-aware validation of several computational models developed within the Collaborative Research Center 1313. With the intended validation benchmarks developed in this project, we plan to overcome the forgoing shortcomings.

Validation of Physical Models in the Presence of Uncertainty
In the field of modeling, it is fairly common to use the word validation to refer to simple comparisons between model outputs and experimental data. Usually this comparison constitutes plotting the model results against data on the same axes to provide visual assessment of agreement or lack thereof. According to Moser et al. [24], while comparisons between model and data are at the heart of any validation procedure, there are a number of concerns with such naive comparisons. First, these comparisons tend to provide qualitative rather than quantitative assessments. Such qualitative assessments are often essential and informative. However, they are clearly insufficient as a basis for making decisions regarding model validity. Second, naive comparisons often disregard or only partly account for existing uncertainties in the experimental observations or the model input parameters. Without accounting for these uncertainties, it is not possible to appropriately determine whether the model and data agree. Third, such comparisons can not reveal whether the model is appropriate for the intended purposes, as they mainly focus on the agreement in the observable quantities.
These pitfalls of straightforward but naive comparisons gives rise to the need for a framework that includes a validation metric. This metric must measure the system response quantities of an experiment with the ones from a computational model, while accounting for uncertainties in both in a rigorous way. This fact is widely recognized, particularly in the statistics community, and there are a number of possible approaches, but here we employ a statistical framework. We have developed such framework by means of a Bayesian validation framework that incorporates parameter and conceptual uncertainty. Incorporation of a fully Bayesian approach yields an optimal bias-variance trade-off against the experimental data and provide an integrative quantity for model validation. Additionally, in order to guarantee the feasibility of the Bayesian validation framework for computationally expensive models, we accelerate the computations for expensive models via model reduction techniques.
The remaining of this report is structured as follows: Chapter 2 introduces components of the Bayesian validation framework. This is followed by Chapter 3, which elucidates the contribution of a surrogate modeling to offset the computational cost. Lastly, the application of the Bayesian validation model to flow simulation models for fractured porous media is presented.

Validation from Bayesian Perspective
Bayesian epistemology offers a robust framework for characterizing scientific inference since its simple concept lies in the fact that rational belief comes in degrees that can be measured in terms of probabilities. Moreover, Bayesian epistemology has resulted in the useful elucidation of notions such as confirmation. Thus, it is proven to form a viable method for data-driven validation of computer simulations and can provide a solid basis for a sound evaluation of computer simulations.
Bayesian epistemology dates back to ideas by Thomas Bayes and was reinforced by many scientists in the twentieth century, such as Bruno de Finetti, Frank P. Ramsey, and Leonard Savage, among others. Bayesianist believes that trust comes in degrees that are measured in terms of probabilities, or probability densities when dealing with characteristics with a continuous range of values. They generally tend to proceed with the following three steps. First, they formulate plausible hypotheses related to a simulation. Second, they consider rational degrees of belief in these hypotheses. Thirdly, they apply the bayesian principles. Following these steps yields posterior probabilities that reveal the degree of trust, one shall rationally invest in the hypotheses.
What are the relevant hypotheses regarding the validation of a simulation? These could be building up or losing trust in simulations. In terms of validation, the hypothesis is whether the model is able to satisfactorily represent the real system of interest. Moreover, for assessment of physical phenomena in question, several representations, i.e. models, might exist with different approaches and assumptions to analyze the occurring processes.

Bayes' theorem
The Bayesian approach to validation extensively exploits Bayes' theorem. This theorem is a combination of traditional probabilities and statistics. Let us assume A and B are two events. The conditional probability of event A given that we know that the event B has occurred, P (A|B), can be cast as: where P (B) denotes the probability of event B, which poses positive value, and P (A ∩ B) signifies the probability that both A and B occurred. Given that P (A ∩ B) is equal to P (B ∩ A), Eq. (2.1) can be recast as follows: Bayes' theorem can simply be obtained by dividing the right-hand side of Eq. (2.2) by P (B), which gives: Bayes' theorem, as expressed in Eq. (2.3), connects the conditional probability P (B|A) to other conditional probability P (A|B). P (A) is the prior probability and P (A|B) stands for posterior probability, which are two basic concepts specific to Bayesian methods. While the posterior is the conditional probability based on the event B, the prior is unconditioned probability which is used to integrate the expert knowledge or previous experience.

Bayesian inference
Bayesian statistics makes use of Bayes theorem, described earlier, to fit a statistical model to the problem at hand. This is achieved by updating the prior knowledge on hyper-parameters Θ, defined by the conditional probability Θ ∼ P (Θ), with possibly few observation data points. These hyper-parameters are treated as random variables, whose subjective definition must disclose the available information before any measurement of the quantity of interest, Y, is performed.
Bayesian theorem in the context of the statistical inference can be recast as the following: where P (θ|Y) denotes the posterior distribution of the hyper-parameters, p(Y|θ) the Assuming that a measured data set of Y = (Y 1 , ..., Y N ) T with the independent realization is available, the probability of observing the data can be defined by the likelihood function as the following: where Σ denotes the covariance matrix, which includes the measurement error as well as all other error sources. P (Y), in Eq. (2.4), is a normalization factor, which ensures that the posterior probability sums up to one. It is also known as marginal likelihood or evidence and can take the following form: The posterior distribution in Eq. (2.4) provides a summary of the inferred information regarding the hyper-parameters after updating the prior knowledge with the observed data. However, the practical computation of the posterior distribution P (θ|Y) is nothing but trivial. There exists analytical expressions only for special choices of the prior distribution P (θ). For more involved cases, one approach is to use sampling methods, such as Monte Carlo, or Markov Chain Monte Carlo sampling methods to approximate the posterior distribution.
One advantage of using the Bayesian inference in the task of model validation is that the uncertainty on hyper-paramters θ can be incorporated into prior and posterior evaluations of quantities of interest y by so-called predictive distributions. The prior predictive distribution can be obtained by the following expression: which is average of the conditional distribution of the prior distribution over the prior distribution P (θ). Similarly, the posterior predictive distribution P pred (y|Y) can be computed by averaging the conditional distribution P (y|θ) over the posterior parameter distribution P (θ|Y), defined in Eq.

Model comparison in Bayesian framework
To compare alternative computational models with possibly distinct conceptual models, various strategies have been suggested in the literature. The benefits of this comparison are twofold. First, this evaluates their strengths and weaknesses. Second, their predictive ability is assessed. Hoeting et al. [14] proposed Bayesian model averaging (BMA) as a formal statistical approach, which allows comparing alternative conceptual models, testing their adequacy, combining their predictions into a more robust output estimate, and quantifying the contribution of conceptual uncertainty to the overall prediction uncertainty.
The BMA method is grounded on Bayes' theorem, which, as mentioned earlier, combines a prior belief about the efficacy of each model with its performance in replicating a common measurement data set. BMA can be regarded as a Bayesian hypothesis testing framework, combining the idea of classical hypothesis testing with the ability to examine multiple alternative models against each other in a probabilistic manner. It returns model weights that represent posterior probabilities for each model to be the most appropriate one from the set of proposed competing models. Additionally, the computed weights can provide a ranking and a quantitative comparison of the competing models. Bayes' theorem closely follows the principle of parsimony or Occam's razor [1], in that the posterior model weights offer a compromise between model complexity and goodness of fit, also known as the bias-variance trade-off [9].

Bayesian Model Averaging Framework
Let us consider that N m plausible, competing models M k are available. The posterior predictive distribution of a quantity of interest θ in Eq. (2.8) given the vector of observed data Y can be expressed as: where P (M k |Y) being discrete posterior model weights.The weights can be interpreted as the Bayesian probability of the individual models to be the best representation of the system from the pool of competing models. The model weights (posterior probabilities of models) are given by Bayes' theorem, which can be recast for a set of M k competing models as: where P (M k ) denotes the prior probability, also known as the subjective model credibility that model M k could be the the most plausible model in the set of models before any comparison with observed data have been made. Hoeting et al. [14] proposed that a "reasonable, neutral choice" could be equally likely priors, i.e. P (M k ) = 1/N m , in case of paucity of prior knowledge regarding the merit of the different models under study. The denominator in Eq.(2.10) is the normalizing constant of the posterior distribution of the models and can simply be obtained by determination of the individual weights. Since all model weights are normalized by the same constant, this normalizing factor could even be neglected. Thus, the ranking of the individual models against each other can be represented by the proportionality: p(Y|M k ) expresses the Bayesian model evidence (BME) term, as introduced in Section 2.1.2. BME is also referred to as marginal likelihood or prior predictive, since it quantifies the likelihood of the observed data based on the prior distribution of the parameters. The BME term can be estimated by integration over the full parameter space Θ k , which is known as Bayesian integral by Kass and Raftery [17], and is expressed as: with θ k being the parameter vector of model M k with the dimension of N p,k . Θ k denotes the parameter space of model M k , and P (θ k |M k ) is the corresponding prior distribution. The likelihood or probability of the parameter set θ k of model M k to have generated the measurement data set is represented by p(Y|M k , θ k ) in Eq. (2.12). For more details on the properties of BME and a comparison of available techniques to evaluate this term, the reader is referred to Schöniger et al. [30]. Here,we perform a brute-force MC integration over each model's parameter and input space to obtain the BME values.

Bayesian hypothesis testing
The Bayesian hypothesis testing framework was first introduced by Jeffreys [16]. The key component in this framework is the Bayes factor, whose extensive description in the context of practical applications is given by the review paper of Kass and Raftery [17]. They define the Bayes factor as the ratio of BME for two alternative models.
The Bayes factor is regarded as a measure for significance in Bayesian hypothesis testing. It quantifies the evidence (literally, as in Bayesian model evidence) of hypothesis M k against the null-hypothesis M l . Jeffreys provided a rule of thumb in his book, Theory of probability [16] for the interpretation of Bayes factor values on a log 10-scale. The grades of evidence is summarized in Table 2.1.  [16] log 10 (BF ) Interpretation 1 -3 anecdotal evidence 3 -10 substantial evidence 10 -100 strong evidence > 100 decisive evidence Following this suggestion, a Bayes factor which lies between 1 and 3 indicates evidence in favor of M k that is "not worth more than a bare mention", a factor of up to 10 represents "substantial" evidence, and a factor between 10 and 100 can be regarded "strong" evidence. Finally, a Bayes factor greater than 100 admits "decisive" evidence, i.e., it can be used as a threshold to reject models based on poor performance in comparison to the best performing model in the set.

Theoretical Upper Limit for Model Performance
Bayes Factor in the context of Bayesian hypothesis testing provides a performance comparison of pairwise competing models. However, in a validation benchmark task, we are also interested in comparing their performance to the best achievable performance. Schöneger et al. [29] argued that this theoretical upper limit for model performance exists when the measurement data set has noise. They propose that this limit can be established via determining a distribution of BME for a so-called theoretically optimal model (TOM), which is also dubbed as a sure-thing hypothesis by Jaynes [15]. They define the observed data set as TOM, as it gives an exact fit with zero bias while having a minimum number of parameters, i.e. exactly zero, which is equivalent to zero variance. Stated differently, the TOM indicates the expected best possible performance in presence of measurement error. This theoretically optimal distribution of BME in presence of measurement noise can be defined as the distribution of likelihoods of the observed data set given the perturbed data sets. Assuming that measurement 2.2 Model comparison in Bayesian framework 9 errors follow a Gaussian distribution and are independent and identically distributed, the TOM performance (shown as log-BME) has a distribution of the weighted sum of normal squared residuals. Consequently, this distribution can be defined by the chi-square distribution [13] as: with k being the number of degrees of freedom, which is equal to the size of the observed data set N s .

Uncertainty Propagation in Bayesian Analysis
A classic Bayesian analysis, discussed in the previous chapter, requires propagation of the parametric uncertainty through the given computational model. This task is also known as the uncertainty propagation (UP). Typically, a significant number of model evaluations are needed in order to yield convergent statistics. In practice, however, the computational complexity of the underlying computational model, as well as the total available computational budget severely restrict the number of evaluations that one can actually carry out. In such situations, the estimates produced by the Bayesian analysis lack sufficient trust, as the limited number of model evaluations can yield additional uncertainty.
The most common approach when dealing with expensive models is to replace them with easy-to-evaluate surrogates. Simply put, one evaluates the model on a set of design points and then strives to establish an accurate relationship between the response surface and the design points. Then, the original computational model can be substituted by its surrogate in the Bayesian analysis. The polynomial chaos expansion (PCE) is one of the most rigorous approach to UP, thanks to its strong mathematical basis and ability to provide functional representations of stochastic quantities. However, the accuracy of the prediction of these surrogate models, trained with only a handful of simulations is debatable. This argument is rooted in the fact that the surrogates do not attempt to quantify the epistemic uncertainty associated with their predictions.
The goal of this chapter is to highlight how a surrogate model using a PCE can be constructed for computationally intensive models with as few simulations as the computational budget allows. It is also shown how the Bayesian formalism can be materialized by employing the concept of PCE to account for the uncertainty in surrogate's predictions. Moreover, we introduce a set of sequential adaptive sampling strategies, in which one attempts to augment the initial design in an iterative manner. Doing so, interesting regions in the parameter space are properly explored, avoiding the waste of computational resources as opposed to the so-called one-shot designs. These regions are more likely to provide valuable information on the behavior of the original model responses.

Polynomial Chaos Expansion
In a probabilistic framework, uncertainties in input parameters are modeled via random variables. These input uncertainties can be investigated using a Polynomial chaos expansion (PCE). This method provides the means to develop an approximation to the map between inputs and the quantities of interest (QoI). This mapping is both computationally tractable and sufficiently accurate. The main idea of PCE is to expand a QoI with a finite variance in a suitably built basis of multivariate polynomials that are orthogonal to the joint probability density functions of the inputs. It is worth noting that the random variables are assumed to be statistically independent or may be linearly correlated. The linear correlation can be handled by adequate linear transformation [26].
A PCE is a linear regression that includes linear combinations of a fixed set of nonlinear functions with respect to the input variables, known as basis functions (Section 3.2.1). The PCE of the random variable Y can be cast as the following: where Ψ α (X) represents multivariate polynomials orthogonal with respect to f X and α denotes a multi-index that represents the components of the multivariate polynomials Ψ α . The c α ∈ R are the corresponding coefficients (coordinates). For practical reasons, the sum in Eq. (3.1) needs to be truncated to a finite sum, by introducing the truncated polynomial chaos expansion: where A ⊂ N M denotes the set of selected multi-indices of the multivariate polynomials. A standard truncation scheme can be defined as all polynomials in the M input variables of total degree less or equal to p: For more other truncation schemes, the reader is referred to [23].

Polynomial basis functions
The multivariate polynomials Ψ α (X) are tensor product of the univariate polynomials: The univariate orthonormal polynomials ψ (i) α i (x i ) must satisfy the following: where i represents the input variable with respect to which they are orthogonal as well as the corresponding polynomial family, j and k the corresponding polynomial degree, f X i (x i ) is the ith-input marginal distribution and δ jk is the Kronecker delta. The classical families of univariate orthonormal polynomials are given for reference in [36]. For detailed description of each of this classical families, the reader is referred to [40].
The calculation of polynomial basis via the classical families is grounded in the fact that an exact knowledge about the probability density functions is available. However, the information about the distribution of the data is distinctly restricted in engineering applications, most importantly when environmental influences or natural phenomena are of interest or when prediction is involved. For instance, the material properties of subsurface reservoirs are not readily available to shed light on their distribution. Oladyshkin & Nowak [26] demonstrate that statistical moments are the only source of information that is propagated in all polynomial expansion-based stochastic approaches. The author leverage this fact to propose an arbitrary polynomial chaos expansion (aPCE), that can operate with probability measures that may be implicitly and incompletely defined via their statistical moments. Using aPCE, one can build the orthonormal polynomials even in the absence of the exact probability density function f X (x).

Calculation of the coefficients using Bayesian sparse learning
Recently, sparse learning algorithms in the context of linear modeling has received attention. These algorithms set many weights to zero in the estimator predictor function. Sparsity is an attractive concept, which offers elegant complexity control, over-fitting control, feature extraction, the potential for characterization of meaningful input variables along with the practical computational speed and compactness. As stated earlier, PCEs also belong to linear regression models and employing the concept of sparsity There are many mathematical approaches when dealing with a regression problem that lead to a sparse solution. These approaches have led to the emerge of numerous sparse solvers in the compressed sensing (e.g. [2]), as well as in the sparse PCE. The proposed solvers in the context of PCE can be categorized into convex optimization solvers, greedy methods, iteratively re-weighted methods, and Bayesian sparse learning, which is also known as compressive sensing. For further details on different solvers in each category, the reader is referred to the comprehensive survey of Lüthen et al. [21]. Here, we employ sparsity within a Bayesian framework by a Bayesian sparse learning method to be able to provide a probabilistic prediction, i.e. a prediction with the associated uncertainty. This prediction uncertainty can be used as the expected error occurring when replacing the original computational model by a possibly less accurate surrogate.
In Bayesian sparse learning, one imposes a sparsity-inducing prior on the coefficients (weights) of the predictors (Ψ α i in the expansion (3.2)), whose parameters are considered to be random variables with a hyperprior. Then, the posterior of the weights are inferred e.g. using a fast marginal likelihood maximization algorithm [39]. This learning process leads to extremely sparse inferred predictors, since they yield relatively few non-zero c i parameters. That means a significant number of the predictors give posterior distributions centered at zero.
Let the target variable be Y, which is given by a deterministic function y(X, c) with an additive Gaussian noise which reads: where is a zero mean Gaussian random variable with precision (inverse variance) β. Hence, the equation can be cast as: where this is a function of the parameter c and β.
We introduce a Gaussian prior distribution over the parameter vector c by giving each of the weight parameters c i a separate hyper-parameter α i . Thus, the prior of the weights reads as: where α = {α i , ..., α M } denotes the precision of the prior over its associated weight parameter c. The form of prior is ultimately responsible for the sparsity properties of the model (for more details, see [38]). The posterior distribution, conditioned on the model responses, is given by combining the likelihood in Eq. (3.8) and the prior Eq. in (3.9) according to Bayes' rule. This posterior, given α, can take the form: which is Gaussian N (c|µ, Σ) with: where Ψ is the design matrix of the size N × M with elements Ψ ni = ψ i (x n ), and A = diag(α i ). The values of α and β can be determined via type-II maximum likelihood [3], also known evidence approximation in the machine learning literature [11,22].

Prediction with Bayesian sparse PCE
Having found values α * and β * for the hyperparameters that maximize the marginal likelihood, one can evaluate the predictive distribution over Y for a new input x by: (3.12)

On the accuracy of the surrogate model 15
The predictive mean is given by Eq.(3.11) with c set to the posterior mean µ, and the variance of the predictive distribution is given by: where Σ is calculated by Eq.(3.11) in which α and β set to their optimized values α * and β * .

On the accuracy of the surrogate model
Once the surrogate model is constructed, its accuracy and predictive capability need to be assessed. In case, a N v independent set of inputs and outputs, a.k.a a validation set (x (1) , M(x (1) )), ..., (x (j) , M(x (j) )) , is available in addition to a training set for training the surrogate model, the validation error can be computed as: V al ) denotes the mean of model responses for the validation set. Since computation of the aforementioned error requires a large number of model evaluations, it is only computationally tractable for simple analytical functions. To avoid additional model evaluations for assessing the accuracy of the surrogate model, an error based on the already evaluated ED is more desirable.
One common approach is the leave-one-out error (LOO), proposed by [8,35] explicitly introduced for PCE. This error, denoted by LOO , is composed of rebuilding N surrogate models in sequential (M P C\i ), using the original experimental design excluding i-th set (X \ x (i) ). Then the prediction error at the excluded set (x (i) ) is computed. For more details, see [5]. Blatman [4] shows that calculating N independent surrogates is not needed, when using the linear superimposition of orthogonal terms, which is the case for PCE. Alternatively, the error can be calculated analytically from a single surrogate based on all sets in the ED using the following equation: where h i is the i-th diagonal entry of the experimental matrix A(A T A) −1 A T . 16

Experimental design
The computational cost of construction of a surrogate model and its accuracy crucially depends on the number of required evaluations of the computationally expensive forward model on the so-called experimental design (ED). ED is a set of training samples from the joint distribution of the input parameters. Properly designed ED has proved vital for simultaneous reduction of the effect of noise and bias errors which can raise the confidence in the task of Bayesian analysis. This has motivated researchers to examine assorted strategies for constructing the sample set {Ψ (i) } N i=1 beyond the standard MC sampling. In this context, the influence of different experimental designs on predictions have been adequately addressed in the literature [34,10,33,28,7,12].
In general, the sampling approaches can be categorized into two groups: classical sampling and sequential sampling. The common practice in classical (one-shot) sampling is to choose the experimental design P grounded only in the information that is available prior to any model evaluation, e.g. existence of noise, relevance of the input parameters, measurement precision. Then, the computational model is evaluated on the selected samples in the ED, and the surrogate model is finally created. This approach is also known as one-shot approach, as all the sample points in the ED are specified at once and no later evaluations of additional samples are made. This is a quite challenging task, since the determination of an optimal sample size is hindered by lack of prior knowledge about the model behavior.
To tackle this problem, flexible sequential sampling strategies have been proposed, which sequentially determine the samples in the design using the information from previous iterations. The sequential sampling approaches can be grouped into two categories: space-filling sequential sampling and sequential adaptive sampling. Spacefilling approaches make sure that the generated samples cover over the entire domain evenly. These sampling approaches are usually developed from some one-shot sampling criteria by selecting the points in a sequential manner [7]. However, the adaptive sequential sampling, also known as adaptive sampling and active learning in machine learning [31], makes more informed choices of samples via the surrogate model itself or data that it learns from, and hence, achieves better performance with fewer points than the space-filling sampling, resulting in saving the simulation cost of expensive computational models [20].
Here, we adopt a sequential adaptive sampling experimental design (SAED). Firstly, an initial batch of samples is selected via a one-shot experimental design. This design can be produced by a common DOE approach, such as Latin Hyper-cube sampling or random sampling. Next, the model is evaluated provided the previously selected samples. Then, the surrogate model is trained to construct a relationship between the ED and the quantities of interest. After construction of the surrogate model, its accuracy is estimated using an error metric, e.g. Leave-one-out error in Eq. (3.15). Given that the initial ED is chosen to be small, the estimated error metric most probably suggests that the enrichment of ED is needed until either the accuracy requirement is met or the maximum allowed total number of runs is reached. Algorithm 1 summarizes a typical sequential experimental design method. The selection of these additional samples is made by an active learning sampling strategy. Some of these strategies will be explained later in detail. Finally, a new surrogate model is built using all the data gathered thus far, and the model accuracy is estimated again. If either the prescribed accuracy level of the surrogate model or total number of samples is still not reached, the entire sample selection process is repeated. Through sequential selection of samples, more information is available to improve sampling, compared to classical design of experiments [6]. The ultimate goal of this algorithm is to reduce the overall number of samples, as evaluating the samples (running the simulation) is the dominant cost in the entire surrogate modelling process. Essentially, non-optimal designs require more resources to make inferences on the features of interest with the same (or less) level of reward that an optimal design would. Hence, when the computational bottleneck is the evaluation of the QoI for any given realization, the additional computational cost of constructing an optimal design is justifiable.
In what follows, the learning strategies (design criteria) for sequential enrichment of the ED are introduced. These criteria can be implemented for emulators being trained for model calibration, i.e. the measurement data is available and can assist in the task of Bayesian model validation. In Section 3.4.1, a Bayesian optimal design is introduced for sequential learning method. The objective is to identify a new design, among some prospective designs, which maximizes the expected utility. Moreover, these strategies are investigated with an analytical example with 10 input parameters.

Bayesian active learning
The Bayesian framework provides a principled approach for incorporating prior information and/or uncertainties concerning the statistical model via a utility function which encapsulates the experimental goals. Given the current model with design d, the optimal next design point d * is the one which maximizes the expected utility function U (d) over the design space D with respect to the future data y and model parameters θ: One needs to consider every possible observation that could be obtained from an experiment with each design, and then evaluate the relative likelihoods and statistical values of these observations. Ultimately, the design that maximizes the expected utility is selected as the optimal design at each stage of experimentation. In the following sections, we will first show how three choices of utility functions will be introduced. These utilities lead to valid measures of information gain, model evidence and information entropy.

Bayesian Inference for Bayesian sparse PCE
Bayesian sparse PCE provides predictions as a mean value µ y (θ, x, y, z, t) and standard deviation σ y (θ, x, y, z, t), as discussed in Section 3.2.3. Therefore, initial knowledge on model response y(θ, x, y, z, t) in each point of space (x, y, z) and time t for the given exploration parameter set d + from the parameter space D is encoded in the Gaussian prior probability distribution N (µ y (d + , x, y, z, t), σ y (d + , x, y, z, t)). Thus, the prior probability distribution of model response y(θ, x, y, z, t) for the given parameter set d + is forming response space Y that is a multivariate Gaussian, denoted as N d + (µ y , σ y ). According to the Bayesian framework (Section 2.1.2), we can obtain a posterior probability distribution p d + (y|Y) of the model response for the given parameter set d + , incorporating the observed data Y: where the term p d + (Y|y) is the likelihood function that quantifies how well the surrogate model predictions y(d + , x, y, z, t) drawn from the multivariate Gaussian N d + (µ y , σ y ) match the observed data Y and the term p d + (Y) denotes the Bayesian model evidence value for the given parameter set d + .
Assuming independent and Gaussian distributed measurement errors, the likelihood function p d + (Y|y) can be written as: (3.18) where y(d + , x, y, z, t) comes from N ωe (µ y , σ y ) and N * denotes the number of measurement points.

Model evidence based utility
As discussed in Section 2.2, Bayesian model evidence (BME) can be regarded as a metric to rank competing models. Here, we leverage this property of BME to identify the next suitable training point in the parameter space during the sequential design. In each iteration, we compute the BME value for prospective design point d + using the following expression: BME BODE in Equation (3.19) can be approximated by: where the term on the right hand side denotes the expected value E N d + (µy,σy) of the likelihood p d + (Y|y) over the prior N d + (µ y , σ y ) provided by the surrogate's prediction. Consequently, the next training point for the surrogate model, i.e. d * BODE ∈ D, can be identified by maximizing the model evidence BME BODE , one can find the next training point :

Information gain utility
A utility function based on mutual information is known as one of the most widely used Bayesian design criteria, which is based on relative entropy. This utility function includes Kullback-Leibler divergence [19] and seeks to maximize the expected information gain in moving from the multivariate Gaussian prior N d + (µ y , σ y ) to the posterior p d + (y|Y) during the learning procedure.
Formally, the relative entropy DKL BODE [p d + (y|Y), N d + (µ y , σ y )] can be defined for each 3.4 Experimental design 20 candidate sampling point d + from the parameter space D as following: Following [27], one can avoid multidimensional integration in Eq. (3.22) by: Therefore, the optimization problem to select the next training point take the following form: It is worth mentioning that Eq. (3.24) depends not only on BME BODE values from Eq. (3.23), but also on the cross entropy represented by term E p d + (y|Y) (ln [p d + (Y|y)]). This term reflects how informative the likelihood is (see details in [27]). Moreover, we obtain the last term via a rejection sampling technique using the evaluations from the already trained surrogate model.

Information entropy-based utility
Another utility for selection of the next training point in a BODE has its root in information entropy [32] and is often used in machine learning. Here, we aim at reducing the expected information loss during the sequential design. The information entropy H BODE [p d + (y|Y)] to asses information loss for each parameter set d + as the candidate for next training point can be computed by: According to [27], information entropy in Eq. (3.25) can be written as following: We obtain all terms in Eq. (3.26) using prior-bases or posterior-bases sampling on surrogate model's prediction, avoiding any multidimensional integration using methods such as rejecting sampling. Therefore, our optimization problem takes the following form:
To assess the prediction accuracy of y(θ, t), in Eq. (3.28) comparing to the synthetic observed data y * , we use the likelihood function in (2.5), assuming independent and Gaussian distributed error of σ = 2. Due to the random nature of the discussed sampling techniques, the results presented in the following are obtained by running 60 independent replications. As convergence criteria, we monitor the difference of Bayesian model evidence and Kullback-Leibler divergence defined as below with their reference values: (3.29) Figure 3.1 reveals that with only 150 runs for a highly non-linear problem with the parameter space of 10 dimensions, the BME and DKL values converge to the reference value, obtained by running a million of model evaluations. It can also be noticed that the information gain utility (DKL) and information entropy based utility (infoEntropy) have reached an acceptable convergence after only 100 runs.  Flow in porous media is often characterized by very strong heterogeneities,in particular fractures, whose influence is important for the understanding of the overall systems' behavior in many natural and technical applications. There are many competing model concepts to represent the flow in fractured porous media. These models vary drastically in their level of geometric details, simplifications and abstraction. In this chapter, we apply the Bayesian Validation Framework, introduced earlier, to set up a benchmark for two models of flow simulation in fractured porous media.

Experimental setup
Two different cases were considered in the experiments. In the first case, the fractures are connected from one side of the sample to the other (

Forward models
Two model variants for flow simulation in fractured porous media have been investigated in this study. The first model (B01) employs a phase-field fraction representation, while the second (B03) makes use of a sharp fracture model. The model B01 simulates the solid fluid interaction within the porous medium under the consideration of Biot's Theory of consolidation. Moreover, the flow of the fluid in the porous medium is described by Darcy's law where the permeability within the cracks is increased by an additional permeability tensor. It models Poiseuille-type flow within the cracks and is derived by the lubrication theory. The solution is obtained by the finite element method in FEAP [37].
The Model B03 follows a discrete fracture network approach. All fractures are geometrically resolved conforming to the mesh. A mixed-dimensional model is used for Darcy flow both in the bulk porous medium and the fractures. The permeability of the fractures is determined by the Poiseuille approximation using the hydraulic diameter. The solution is obtained by a Finite-Volume method implemented in DuMu x [18].  Uncertain parameters and their ranges for the model validation analysis are summarized in Table 4.1.

Error quantification
Bayesian statistics updates the prior belief on the model by comparing the model responses with the measured data. It can include not only the errors and uncertainties in the observed data, but it can also take into consideration other sources of errors, such as numerical errors and surrogate model's error. These errors can be considered in the covariance matrix Σ in Eq. (2.5), assuming that they are follow a normal distribution.
In the following, we will shed light on three sources of error used in the validation benchmark of the flow simulation models for fractured porous media.

Measurement errors
We use the standard deviation of the measured value for each measurement sensor as the experimental error. These values are extracted for two cases and for different flow rates for the calibration and validation, separately. The pressure measurement's distribution for two different cases for the calibration step (i.e. with the flow rate 0.2 [ml/min]) are shown in Figure 4.3.

Numerical errors
Here, we only investigate the discretization error that originates from a certain choice of the meshing size. Following [25], we take a heuristic approach to quantify this error is to make a comparison between different mesh spacing and to make use of a generalized Richardson extrapolation to estimate the error. The Richardson extrapolation takes the following form: f k denotes the exact solution to the discrete equation on a mesh with spacing h k (known),f stands for the exact solution to the original PDE (unknown). g p is the error term coefficient andp indicates the observed order of accuracy. Here, we seek the error with an order of one. Thus, the unknowns, i.e.f and g p , can be easily determined via a least square method.

Surrogate model's errors
As discussed earlier, we replace the computationally intensive model in the Bayesian framework to offset the computational cost. By doing so, a new source of error is introduced, namely surrogate modeling's error. In this study, we set the limit for computational budget to 200 simulation runs for each model. The prediction error of the surrogate model, introduced in Section 3.2.3, is taken into accounts in the calculation of the likelihoods in Eq. (2.5).

Bayesian model validation
For the validation within the Bayesian approach, the following steps are taken. First, a surrogate model is trained based on the simulation results, obtained by the original computational models, here the pressure readings at 9 sensors shown in Figure  4.1. Then, the Bayesian updating is performed in that the prior knowledge on the distribution of uncertain parameters is revised by comparing the model outputs with the measurement, using Eq. (2.4). This can be done either by a brute-force Monte Carlo sampling (a.k.a rejection sampling) or a more sophisticated methods, for instance Markov-Chain Monte Carlo. Afterwards, the resulting posterior distribution is used to compare the model with a new observed data set. In what follows, some preliminary results will be presented.

Preliminary results
As discussed in Section 2.2, Bayesian model evidence (BME) provides a reasonable metric for assessment of the model performance. BME is sometimes also referred to as marginal likelihood because it quantifies the likelihood of the model to have produced the observed data averaged over the complete prior parameter space. To also consider the uncertainty associated with this value and its influence on our analysis, the measured data has been perturbed with a Gaussian noise.
Moreover, we use the distribution of BME for TOM, following [30], computed by a chisquare distribution, as upper limit for model performance, as defined in Section 2.2.3. Figure 4.4 shows the distribution of the BME for the models B01 and B03, as well as The BME distributions for the competing models with the perturbed data set the upper limit performance.These values are obtained by varying the measurement data set to account for the uncertainty associated with BME. Model B01 shows slightly better performance, as its BME distribution is closer to that of TOM.
One approach to further compare models is hypothesis testing in a Bayesian setting. This can be done by a so-called Bayes factor, that has been presented in Section 2.2.2. The Bayes factor is a measure for significance in Bayesian hypothesis testing. It quantifies the evidence of hypothesis that M l is the data generating model against a null-hypothesis. The null-hypothesis can be defined as one model to be the best model among the set of models under investigation. Figure 4.5 illustrates the distribution of Bayes factor of two models in log-scale. The significance levels, suggested by Jeffreys [16] are also marked here with vertical solid lines. One can interpret that there is no  The Bayesian factor for two scenarios: (a) Model B03 performs better as the data-generating model (b) theoretically optimal model has the best performance substantial evidence to favor B01 against B03, as shown in Figure 4.5a. We can also compare the performance of each model with the upper limit for performance, provided by the theoretically optimal model, in short "TOM". In Figure 4.5b, the distribution of the Bayes factor includes negative values, that means the evidence is in favor of the null hypothesis instead of the actual hypothesis. Thus, it can be concluded that there is substantial evidence in favor of the TOM against Model B01.