Bayesian inference and non-linear extensions of the CIRCE method for quantifying the uncertainty of closure relationships integrated into thermal-hydraulic system codes

TL;DR

This paper extends the CIRCE method for uncertainty quantification in thermal-hydraulic codes by incorporating Bayesian inference and non-linear modeling using Gaussian process emulators, enabling more accurate and efficient analysis.

Contribution

It introduces a Bayesian framework and Gaussian process emulators to improve the CIRCE method's handling of non-linearities and statistical uncertainties.

Findings

Bayesian approach provides posterior distributions of parameters.

GP emulators significantly reduce computational time.

Method applied to condensation closure relationships with successful results.

Abstract

Uncertainty Quantification of closure relationships integrated into thermal-hydraulic system codes is a critical prerequisite in applying the Best-Estimate Plus Uncertainty (BEPU) methodology for nuclear safety and licensing processes.The purpose of the CIRCE method is to estimate the (log)-Gaussian probability distribution of a multiplicative factor applied to a reference closure relationship in order to assess its uncertainty. Even though this method has been implemented with success in numerous physical scenarios, it can still suffer from substantial limitations such as the linearity assumption and the difficulty of properly taking into account the inherent statistical uncertainty. In the paper, we will extend the CIRCE method in two aspects. On the one hand, we adopt the Bayesian setting putting prior probability distributions on the parameters of the (log)-Gaussian distribution. The posterior distribution of the parameters is then computed with respect to an experimental database by means of Markov Chain Monte Carlo (MCMC) algorithms. On the other hand, we tackle the more general setting where the simulations do not move linearly against the multiplicative factor(s). MCMC algorithms then become time-prohibitive when the thermal-hydraulic simulations exceed a few minutes. This handicap is overcome by using Gaussian process (GP) emulators which can yield both reliable and fast predictions of the simulations. The GP-based MCMC algorithms will be applied to quantify the uncertainty of two condensation closure relationships at a safety injection with respect to a database of experimental tests. The thermal-hydraulic simulations will be run with the CATHARE 2 computer code.

Tables

Table 1. Table 1: Posterior mean estimates m ^ 1 subscript ^ 𝑚 1 \hat{m}_{1} and m ^ 2 subscript ^ 𝑚 2 \hat{m}_{2} , the corresponding confidence intervals and also the median estimates exp ⁡ m ^ 1 subscript ^ 𝑚 1 \exp{\hat{m}_{1}} and exp ⁡ m ^ 2 subscript ^ 𝑚 2 \exp{\hat{m}_{2}} of Λ 1 subscript Λ 1 \Lambda_{1} and Λ 2 subscript Λ 2 \Lambda_{2}

	${\hat{m}}_1$	$I{C}_{95\%}\left(𝔼\left[m_1|{𝐳}^f\right]\right)$	${\hat{m}}_2$	$I{C}_{95\%}(𝔼\left[m_2\right|{𝐳}^f]$)	$\exp {\hat{m}}_1$	$\exp {\hat{m}}_2$
Setting 1	$0.7602$	$[0.7518,0.7686]$	$-1.5071$	$[-1.5353,-1.4789]$	$2.1387$	$0.2216$
Setting 2	$0.6963$	$[0.6863,0.7064]$	$-1.3814$	$[-1.4145,-1.3484]$	$2.0063$	$0.2512$
Setting 3	$0.6959$	$[0.6876,0.7043]$	$-1.3706$	$[-1.4003,-1.3408]$	$2.0055$	$0.2540$

Table 2. Table 2: Posterior mean estimates of σ 1 2 superscript subscript 𝜎 1 2 \sigma_{1}^{2} and σ 2 2 superscript subscript 𝜎 2 2 \sigma_{2}^{2} and their corresponding confidence interval

	${\hat{\sigma }}_1^2$	$I{C}_{95\%}\left(𝔼\left[{\sigma }_1^2|{𝐳}^f\right]\right)$	${\hat{\sigma }}_2^2$	$I{C}_{95\%}(𝔼\left[{\sigma }_2^2\right|{𝐳}^f]$)
Setting 1	$0.1592$	$[0.1564,0.1621]$	$0.8883$	$[0.8564,0.9203]$
Setting 2	$0.1528$	$[0.1463,0.1593]$	$1.4785$	$[1.4305,1.5265]$
Setting 3	$0.1474$	$[0.1420,0.1528]$	$1.5365$	$[1.4877,1.5852]$

Table 3. Table 3: 95 % percent 95 95\% -fluctuation intervals of Λ 1 subscript Λ 1 \Lambda_{1} and Λ 2 subscript Λ 2 \Lambda_{2}

	$I{F}_{95\%}\left({\Lambda }_1\right)$	$I{F}_{95\%}\left({\Lambda }_2\right)$
Setting 1	$[0.82,4.80]$	$[0.02,1.59]$
Setting 2	$[0.73,4.44]$	$[0.01,2.93]$
Setting 3	$[0.76,4.39]$	$[0.01,3.03]$

Equations

$$M_{ref}(\mathbf{x})$$

$$M_{\lambda }(\mathbf{x})=\lambda \times M_{ref}(\mathbf{x}).$$

$${\mathbf{X}}^f=[{\mathbf{x}}_1^f,\cdots ,{\mathbf{x}}_n^f{]}^T\in {\mathcal{M}}_{n,q}(\mathbb{R})$$

$${\mathbf{z}}^f=(z_1^f,\cdots ,z_n^f{)}^T\in {\mathbb{R}}^n.$$

$$z_i^f=Y(M_{{\lambda }_i}({\mathbf{x}}_i^f))+{ϵ}_i$$

$${ϵ}_i\sim \mathcal{N}(0,{\sigma }_{{ϵ}_i}^2).$$

$$z_i^f=Y(M_{{\lambda }_i}({\mathbf{x}}_i^f))+b({\mathbf{x}}_i^f)+{ϵ}_i$$

$$M_{{\lambda }_i}({\mathbf{x}}_i^f):=(M_{{\lambda }_{i,1}}^1({\mathbf{x}}_i^f),\cdots ,M_{{\lambda }_{i,p}}^p({\mathbf{x}}_i^f))$$

$$z_i^f=Y_{{\lambda }_i}({\mathbf{x}}_i^f)+{ϵ}_i.$$

$${\Lambda }_j=\mathcal{N}(m_j,{\sigma }_j^2)\text{or}{\Lambda }_j=\mathcal{L}\mathcal{N}(m_j,{\sigma }_j^2),1\le j\le p.$$

$$\text{Cov}({\Lambda }_j,{\Lambda }_k)=0.$$

$$m=(m_1,\cdots ,m_p{)}^T\in {\mathbb{R}}^p$$

$${\sigma }^2=({\sigma }_1^2,\cdots ,{\sigma }_p^2{)}^T\in {\mathbb{R}}^{{+}^p}.$$

$$z_i^{\prime ,f}=h^T({\mathbf{x}}_i^f){\lambda }_i^{\prime }+{ϵ}_i$$

$$z_i^{\prime ,f}=h^T({\mathbf{x}}_i^f){\alpha }_i^{\prime }+{ϵ}_i$$

$$H=[h({\mathbf{x}}_1^f),\cdots ,h({\mathbf{x}}_n^f){]}^T\in {\mathcal{M}}_{n,p}(\mathbb{R}).$$

$$I{F}_{95\%}({\Lambda }_j)=[{\hat{m}}_j-1.96{\hat{\sigma }}_j,{\hat{m}}_j+1.96{\hat{\sigma }}_j]$$

$$I{F}_{95\%}({\Lambda }_j)=[\exp ({\hat{m}}_j-1.96{\hat{\sigma }}_j),\exp ({\hat{m}}_j+1.96{\hat{\sigma }}_j)].$$

$$\pi (b,{\sigma }^2|{\mathbf{X}}^f,{\mathbf{z}}^{\prime ,f})=\frac{\mathcal{L}({\mathbf{z}}^{\prime ,f}|{\mathbf{X}}^f,b,{\sigma }^2)\pi (b,{\sigma }^2)}{\mathcal{L}({\mathbf{z}}^{\prime ,f})}$$

$$\mathcal{L}({\mathbf{z}}^{\prime ,f}|{\mathbf{X}}^f,b,{\sigma }^2)={\int }_{{\mathbb{R}}^{n\times p}}{\mathcal{L}}^C({\mathbf{z}}^{\prime ,f},{\mathbf{A}}^{\prime }|{\mathbf{X}}^f,b,{\sigma }^2)d{\mathbf{A}}^{\prime }$$

$$\pi (b,{\sigma }^2,{\mathbf{A}}^{\prime }|{\mathbf{X}}^f,{\mathbf{z}}^{\prime ,f})\propto \mathcal{L}({\mathbf{z}}^{\prime ,f}|{\mathbf{A}}^{\prime },{\mathbf{X}}^f,b,{\sigma }^2)\pi ({\mathbf{A}}^{\prime },b,{\sigma }^2).$$

$$\pi ({\mathbf{A}}^{\prime },b,{\sigma }^2)=\pi ({\mathbf{A}}^{\prime }|b,{\sigma }^2)\pi (b,{\sigma }^2)$$

$${\mathbf{A}}^{\prime }(k+1)\sim {\mathbf{A}}^{\prime }|b(k),{\sigma }^2(k),{\mathbf{z}}^{\prime ,f},{\mathbf{X}}^f$$

$$b(k+1),{\sigma }^2(k+1)\sim b,{\sigma }^2|{\mathbf{A}}^{\prime }(k+1),{\mathbf{z}}^{\prime ,f},{\mathbf{X}}^f$$

$$\pi (b,{\sigma }^2)=\prod _{j=1}^p\pi (b_j,{\sigma }_j^2)=\prod _{j=1}^p\pi (b_j|{\sigma }_j^2)\pi ({\sigma }_j^2)$$

$$\pi(b_{j}|\sigma_{j}^{2})=\frac{\sqrt{a_{j}}}{\sqrt{2\pi}\sigma_{j}}\exp{\Big{\{}-\frac{a_{j}(b_{j}-\mu_{j})^{2}}{2\sigma_{j}^{2}}\Big{\}}}$$

$$\pi(\sigma_{j}^{2})=\frac{\gamma_{j}^{\psi_{j}}}{\Gamma(\psi_{j})}\sigma_{j}^{2(-\psi_{j}-1)}\exp{\Big{\{}-\frac{\gamma_{j}}{\sigma_{j}^{2}}\Big{\}}}.$$

$$(\hat{b},{\hat{\sigma }}^2)=\mathbb{E}[b,{\sigma }^2|{\mathbf{X}}^f,{\mathbf{z}}^{\prime ,f}].$$

$$(\hat{b},{\hat{\sigma }}^2)=\underset{b,{\sigma }^2}{\mathrm{argmax}}\pi (b,{\sigma }^2|{\mathbf{X}}^f,{\mathbf{z}}^{\prime ,f}).$$

$$\pi (\Lambda |{\mathbf{X}}^f,{\mathbf{z}}^{\prime ,f})=\int \pi (\Lambda |m,{\sigma }^2)\pi (m,{\sigma }^2|{\mathbf{X}}^f,{\mathbf{z}}^{\prime ,f})\text{d}m\text{d}{\sigma }^2.$$

Full text

The purpose of the CIRCE method is to estimate the (log)-Gaussian probability distribution of a multiplicative factor applied to a reference closure relationship in order to assess its uncertainty. Even though this method has been implemented with success in numerous physical scenarios, it can still suffer from substantial limitations such as the linearity assumption and the difficulty of properly taking into account the inherent statistical uncertainty. In the paper, we will extend the CIRCE method in two aspects. On the one hand, we adopt the Bayesian setting putting prior probability distributions on the parameters of the (log)-Gaussian distribution. The posterior distribution of the parameters is then computed with respect to an experimental database by means of Markov Chain Monte Carlo (MCMC) algorithms. On the other hand, we tackle the more general setting where the simulations do not move linearly against the multiplicative factor(s). MCMC algorithms then become time-prohibitive when the thermal-hydraulic simulations exceed a few minutes. This handicap is overcome by using Gaussian process (GP) emulators which can yield both reliable and fast predictions of the simulations.

The GP-based MCMC algorithms will be applied to quantify the uncertainty of two condensation closure relationships at a safety injection with respect to a database of experimental tests. The thermal-hydraulic simulations will be run with the CATHARE 2 computer code.

1 Introduction

The use of Best-Estimate system codes for nuclear applications in line with the so-called BEPU (Best-Estimate Plus Uncertainty) methodology has been encouraged from several years for reactor transient simulations and safety analysis. This includes code development as well as V&V (Verification and Validation) and UQ (Uncertainty Quantification) efforts. The paper written by D’Auria et al. (2012) gives an overview of the best practices for safety demonstration in the nuclear community and how the BEPU approach could meet them. As its name implies, the second stage of a BEPU methodology aims to quantify uncertainties related to the numerical simulation of a thermal-hydraulic transient under study. Such uncertainties can have multiple origins but it is commonly admitted that we can split them into two categories: those of numerical nature (out of the scope of the paper), and those related to the modeling of physical phenomena on which the code relies. The Verification step aims to assess the former, essentially by fixing code bugs and computing discretization errors induced by mesh size. Then, the Validation step aims to quantify the remaining gap between the simulations and the real system by means of some well chosen field experiments used as references (Bayarri et al., 2007). In the last two decades, the research dedicated to V&V and UQ activities has led to the development of exhaustive guidelines in the field of computational fluid dynamics (AIAA, 1998; ASME, 2009) which have been adapted for reactor safety licensing under the name of improved BEPU (Unal et al., 2011). The improved BEPU framework provides general requirements which help improving the seminal methods implemented by thermal-hydraulic experts such as CSAU (Boyack et al., 1990) and GRS approach (Glaeser, 2008). Both of these methods, developed respectively by the US Nuclear Regulatory Commission (NRC) and the German nuclear safety institution, propagate through the code the input uncertainties to the outputs of interest, such as the Peak Cladding Temperature (PCT). Their implementation for accidental transients at the full reactor scale has been the subject of the BEMUSE222Best Estimate Methods – Uncertainty and Sensitivity Evaluation program (De Crécy et al., 2008) where LB-LOCA333Large Break - Loss Of Coolant Accident analyses were performed to check whether both scaling and uncertainty propagation could be properly implemented. The probability distributions of uncertain input parameters were assumed to be known in BEMUSE because their determination was not the point of the project. This issue has actually been studied in two further international projects named PREMIUM444Post-BEMUSE REflood Models Input Uncertainty Methods (Skorek et al., 2019), following by the recent SAPIUM project555Systematic Approach for Input Uncertainty Quantification (Baccou et al., 2019). The former aimed to compare in a reflood scenario the performance of several methodologies for quantifying the uncertainty of closure relationships, including, in particular, that of CEA called CIRCE666In French the acronym “CIRCE” means Calcul des Incertitudes Relatives aux Corrélations Élémentaires and the Fast Fourier Transform-based method (Freixa et al., 2016). Unfortunately, the lack of consensus as well as the various practices between the participants raised the need of establishing guidelines and providing general recommendations on how this issue should be addressed. This was the goal of the SAPIUM project.

One of the main source of uncertainty is the one tainting the closure relationships integrated into the balance equations of thermal-hydraulic system codes. Such codes rely on a two phase flow six equations model, consisting of three balance equations for the liquid phase and three for the steam phase (mass, momentum and energy for each phase). The solution of such equations depends on closure relationships constructed from appropriate physical experiments called SETs (Separate Effect Tests). The SETs, often performed at a smaller scale than that of actual nuclear reactors, are instrumented in order to focus on some experimental Quantities of Interests (QoIs) as independently as possible of potential interactions with other phenomena. Then, based on a careful analysis of such tests, the closure relationships under study are established by thermal-hydraulic experts. At the end of this process, though, a model-form uncertainty can remain due to the lack of perfect knowledge of their mathematical expressions. This type of uncertainty can be quantified by the CIRCE method which applies a random multiplicative factor to each uncertain closure relationship (De Crécy and Bazin, 2001). Such factors are in fact modeled either by a Gaussian or log-Gaussian probability distribution whose parameters are estimated from the differences between the experimental QoIs and the corresponding simulations. This estimation stage requires a linear connection between the code outputs and the involved multiplicative factor(s).

In the CIRCE method, the parameters of the (log)-Gaussian distribution are estimated as the Maximum Likelihood Estimator (MLE), which is the configuration of the parameters at which the probability of observing the data is the highest. On the other side, the Bayesian paradigm is to represent uncertain parameters by random variables. First, a prior probability distribution is put on the uncertain parameters, then this distribution is converted into a posterior distribution through the Bayes formula (see for example Ghosh et al. (2007) for an introduction to Bayesian analysis). Various point estimates, such as the maximum a posteriori and the posterior mean can be computed as well as credible intervals. Besides, Bayesian inference is particularly suited for industrial contexts where the prior distribution can benefit from expert judgment, such as accounting for constraints between parameters, upper and lower bounds of variation or even guesses on the best fitting values. In the paper, we will present the Bayesian approach to estimate the parameters of the (log)-Gaussian distributions of the multiplicative factors. Then, the linear assumption of the CIRCE method will be dropped to tackle the case where the simulations move non linearly with respect to the multiplicative factors. However, although thermal-hydraulic simulations are moderately time consuming, intensive sampling required by Bayesian computation makes necessary the use of emulators. We have chosen Gaussian process (GP) emulators which can yield both reliable and fast predictions of the simulations.

Section $2$ recalls the foundations of the CIRCE method, and the statistical equation linked to it. Section $3$ deals with Bayesian estimation in the linear setting via Markov Chain Monte Carlo (MCMC) algorithms. Section $4$ goes ahead with the non linear setting and how it modifies MCMC sampling, now based on GP emulators. In Section $5$, the MCMC algorithms are applied to the COSI (Condensation at the Safety Injection) test facilities to quantify the uncertainty of two condensation models at a safety injection. Section $6$ gives conclusions and provides some ideas for future developments.

2 The CIRCE method

2.1 Motivation

The CIRCE method is devoted to Uncertainty Quantification (UQ) of physical models integrated into thermal-hydraulic system codes. Such models play the role of closure relationships of the balance equations on which thermal-hydraulic system codes are based. In a BEPU context, the implementation of the CIRCE method is motivated by the requirement of assessing the uncertainty of a few targeted quantities for safety, such as the Peak Cladding Temperature (PCT). This can be done by propagating through the code the uncertainty of the physical models to the PCT (if possible, both scaling and systematic biases should be taken into account as well). Prior to such a forward UQ stage, physical models are established by experts by means of one or several experimental database(s). The reference formulation of a physical model is written in the paper as

[TABLE]

with $\mathbf{x}$ being the vector of the involved thermal-hydraulic conditions such as pressure, temperature, geometrical properties, dimensionless numbers, void fraction, and so on. As such reference models are tainted by some degree of empiricism in their mathematical expressions, the UQ stage is compulsory. To this end, the CIRCE method applies a multiplicative factor $\lambda$ to each reference model of the thermal-hydraulic code. This allows to measure a relative model uncertainty, which is suited to deal with experimental databases having different orders of magnitude. A thermal-hydraulic simulation can then be run with respect to a perturbed model $M_{\lambda}(\mathbf{x})$ such that

[TABLE]

The more $\lambda$ moves away from $\lambda=1$, the more the engineer should question the accuracy of the reference model. Plausible bounds can sometimes be specified around $\lambda=1$ using expert knowledge. However, the way of representing the uncertainty between the lower and the upper bound often remains unclear. To do so, the CIRCE method adopts a probabilistic framework where the multiplicative factor is modeled by a (log)-Gaussian probability distribution $\Lambda$. When realizations $\lambda\thicksim\Lambda$ are not experimentally observed, the parameters of $\Lambda$ cannot be readily estimated by empirical mean and variance estimators. However, if the impact of $\Lambda$ on some thermal-hydraulic QoIs predicted by the code is strong enough, a statistical modeling of the mismatch between the experimental QoIs and the corresponding simulations at the reference model can be set to get back to the parameters of $\Lambda$. Let us consider, for example, an heat transfer model impacting a wall temperature in a pipe. The differences between the wall temperature measurements and the corresponding simulations could be used to estimate the (log)-Gaussian factor applied to the heat transfer model.

2.2 The statistical modeling

As explained previously, a database of experimental QoIs is needed to implement the CIRCE method. Throughout the paper, the superscript $f$ (for field) will refer to any experimental quantity. The set-ups of the database are gathered in the matrix $\mathbf{X}^{f}:$

[TABLE]

where each site $\mathbf{x}_{i}^{f}\in{\mathbb{R}}^{q}$ is made up with the physical conditions of the $i$-th experiment ($1\leq i\leq n$). The corresponding QoI at $\mathbf{x}_{i}^{f}$ is denoted by $z_{i}^{f}\in{\mathbb{R}}$. The vector gathering all these QoIs is denoted by

[TABLE]

Let $Y(.)$ be a thermal-hydraulic system code. The CIRCE method relies on a statistical equation which relates $\mathbf{z}^{f}$ to the corresponding simulations run at $\mathbf{X}^{f}$. The differences between the two are assumed to be realizations of $n$ independent zero-mean Gaussian random variables. Hence, for $1\leq i\leq n$,

[TABLE]

with $\lambda_{i}$ being missing (unobserved) and

[TABLE]

The random variable $\epsilon_{i}$ of Equation (6) is the experimental error (or observation error) whose variance $\sigma_{\epsilon_{i}}^{2}$ encompasses the measurement error and a residual variability between the computer code and the physical reality. Other sources of uncertainty could also contribute to the mismatch. In Kennedy and O’Hagan (2001) a function $b(\mathbf{x}_{i}^{f})$, called ”code discrepancy”, is explicitly taken into account between the physical system and the computer code, which leads to the equation

[TABLE]

Such a modeling has recently be applied to the nuclear field, but considering $\lambda$ as a vector of calibration parameters (Wu et al., 2018a, b) (see D for details on the difference between CIRCE and calibration methods). In the paper, we omit this term because the CIRCE method presumes that the gap between $\mathbf{z}^{f}$ and the simulations is mostly due to the uncertainty tainting the physical models. Moreover, the inclusion of such a code discrepancy can cause non identifiability problems from a statistical point of view (Brynjarsdottir and O’Hagan, 2014).

Equation (5) includes the case where the code outputs depend on several uncertain models, $\lambda_{i}$ being a $p$-dimensional sample of the random vector $\Lambda:=(\Lambda_{1},\cdots,\Lambda_{p})^{T}$. In such a case, we have

[TABLE]

where $\lambda_{i,j}$ is the unobserved realization applied to the $j$-th reference model $M^{j}_{ref}(.)$ ($1\leq j\leq p$).

Equation (5) does not mean a composition of the two functions $Y$ and $M_{\lambda_{i}}(.)$. In fact, thermal-hydraulic simulations result from an iterative process by which the balance equations are solved in both time and space. The reference model values are modified across iterations whereas the realization $\lambda\thicksim\Lambda$ is constant and should be considered as an extra input value. Hence, Equation (5) can be rewritten under a more readable way:

[TABLE]

The CIRCE method then relies on several important assumptions. First, all the experimental variances $\sigma^{2}_{\epsilon_{i}}$ are assumed known ($1\leq i\leq n$). Second, each component $\Lambda_{j}$ of $\Lambda$ follows either a Gaussian or a log-Gaussian probability distribution$:$

[TABLE]

The next assumption is not compulsory to run CIRCE, but highly recommended$:$ the $p$ variables $\Lambda_{j}$ are independent from one another. For $1\leq j\neq k\leq p$, we set

[TABLE]

The number of parameters to be estimated is therefore limited to $2p$ instead of $2p+p(p-1)/2$ if cross-covariance terms were estimated as well. In the sequel, we will use concise notations for the $p$ means and variances of $\Lambda$ to be estimated$:$

[TABLE]

and

[TABLE]

2.3 The estimation stage

The estimation of $(m,\sigma^{2})$ requires a linear approximation of the code outputs at a nominal value $\lambda^{\star}$ using a first order Taylor approximation. It is quite frequent that thermal-hydraulic simulations move linearly with respect to the values of the multiplicative factors in a large vicinity of $\lambda^{\star}$. For each input site $\mathbf{x}_{i}^{f}\in\mathbf{X}^{f}$, Equation (5) is then substituted by

[TABLE]

where

$\lambda_{i}^{\prime}:=\lambda_{i}-\lambda^{\star}\thicksim\mathcal{N}(m-\lambda^{\star},\Sigma:=\textrm{diag}(\sigma^{2}))$ is the missing model sample shifted of $\lambda^{\star}$,

• 2.

$h(\mathbf{x}_{i}^{f})\in{\mathbb{R}}^{p}$ is the vector of partial derivatives of the code at $\mathbf{x}_{i}^{f}$,

• 3.

$z_{i}^{\prime,f}$ is the difference between the experimental QoI at $\mathbf{x}_{i}^{f}$ and the corresponding simulation at $\lambda^{\star}$.

The location $\lambda^{\star}$ is chosen by default equal to $\mathbf{1}_{p}$, meaning that the linear approximation is constructed at the reference model.

In the thermal-hydraulic field, such linear approximations are often more accurate with respect to $\alpha=\log{\lambda}$. Equation (5) is then replaced around $\alpha^{\star}=\log{(\lambda^{\star})}$ by

[TABLE]

where

$\alpha_{i}^{\prime}=\alpha_{i}-\alpha^{\star}\thicksim\mathcal{N}(m-\alpha^{\star},\Sigma:=\textrm{diag}(\sigma^{2}))$ is the missing model log-sample shifted of $\alpha^{\star}$,

• 2.

$h(\mathbf{x}_{i}^{f})\in{\mathbb{R}}^{p}$ is the vector of partial derivatives of the code at $\mathbf{x}_{i}^{f}$ (with respect to $\alpha=\alpha_{\star}$).

The location $\alpha^{\star}$ is chosen by default equal to $\mathbf{0}_{p}=\log{\mathbf{1}_{p}}$ meaning that the log-linear approximations are still constructed at the reference model.

Both Models (14) and (15) are statistically identifiable if and only if the rank of the corresponding matrix of the partial derivatives matrix is equal to $p$. This matrix is denoted by

[TABLE]

A small condition number is also expected to ensure reliability of estimates (Belsley et al., 1980). The Maximum Likelihood Estimate (MLE) of $(m,\sigma^{2})$, denoted by $(\hat{m},\hat{\sigma}^{2})$, is then computed by means of a variant of the Expectation-Maximization (EM) algorithm called ECME (A. P. Dempster and Rubin, 1977; Celeux et al., 2010).

The last stage of CIRCE, but not the least, is to choose between the Gaussian distribution $\mathcal{N}(\hat{m},\hat{\sigma}^{2})$ and the log-Gaussian distribution $\mathcal{LN}(\hat{m},\hat{\sigma}^{2})$ where $(\hat{m},\hat{\sigma}^{2})$ are the estimates of Models (14) and (15) respectively. The choice goes in favor of the distribution making the best linear approximation of the code outputs in the vicinity of $\lambda^{\star}$ or $\alpha^{\star}$ respectively. If the Gaussian distribution is chosen for $\Lambda_{j}$, then the $95\%$-fluctuation interval of the multiplicative factor is

[TABLE]

Otherwise, if the log-Gaussian distribution is chosen instead, then the $95\%$-fluctuation interval of the multiplicative factor is

[TABLE]

Either of these intervals is valid if the linearity assumption is well verified across it. Eventually, the probability distribution $\Lambda$ should be propagated through the code to check whether the resulting simulations envelop the physical experiments $\mathbf{z}^{f}$. This is called Envelop calculations in the CIRCE guidelines. Also important is to compare the predicted residuals with the Gaussian density to ensure they broadly match each other (such validation steps of the CIRCE method are out of the scope of the paper).

Remark** 1**

The statistical uncertainty is assessed with confidence intervals centered on $(\hat{m},\hat{\sigma}^{2})$ by using either bootstrapping, or the Fisher information matrix. However, this uncertainty is not taken into account in the computation of the fluctuation intervals.

In the next section, we lay out the Bayesian version of the CIRCE method where the uncertain parameters $m$ and $\sigma^{2}$ are now modeled by random variables. A prior distribution is first put on $(m,\sigma^{2})$, then the posterior distribution of $(m,\sigma^{2})$ is computed using the Bayes formula.

3 The Bayesian inference of the CIRCE method

3.1 General principle

Without loss of generality, we present the Bayesian equations related to Model (15) where $\Lambda$ is specified as a log-Gaussian probability distribution. As the code is linearized at $\alpha^{\star}=\log{\lambda^{\star}}$, we deal with Bayesian estimation of $b:=m-\alpha^{\star}$ and $\sigma^{2}$.

Let $\mathbf{A}^{\prime}=[\alpha_{1}^{\prime},\cdots,\alpha_{n}^{\prime}]^{T}\in\mathcal{M}_{n,p}({\mathbb{R}})$ and $\mathbf{z}^{\prime,f}=(z_{1}^{\prime,f},\cdots,z_{n}^{\prime,f})^{T}$ be respectively the matrix of the $\alpha^{\star}$-shifted missing values of the multiplicative factors and the vector of the $\alpha^{\star}$-shifted experimental QoIs. Bayesian estimation of $b$ and $\sigma^{2}$ begins by specifying a joint prior probability density $\pi(b,\sigma^{2})$. Then, applying the Bayes formula gives the posterior density

[TABLE]

where the numerator is the product between the likelihood related to Equation (15) and the prior distribution $\pi(b,\sigma^{2})$, while the denominator is the marginal likelihood of $\mathbf{z}^{\prime,f}$ playing the role of normalization constant. The likelihood in Equation (19) can be written as

[TABLE]

where $\mathcal{L}^{C}$ is the complete likelihood including the missing data $\mathbf{A}^{\prime}$. By taking into account these data as unknown quantities as well, the Bayes formula gives

[TABLE]

The prior distribution in Equation (21) can be expanded on

[TABLE]

where $\pi(\mathbf{A}^{\prime}|b,\sigma^{2})$ is the density of the product between $n$ Gaussian distributions having same mean $b$ and same variance $\sigma^{2}$. Equation (21) shows a hierarchical Bayesian modeling where no distinction is made between the missing data $\mathbf{A}^{\prime}$ and the parameters $(b,\sigma^{2})$, both being unknown and modeled by random variables. As the computation of the denominator in Equation (19) is challenging, the posterior distribution (21) cannot be readily computed. Instead, an efficient way to generate samples from this distribution is to run Markov Chain Monte Carlo (MCMC) algorithms relying on the numerator only.

3.2 Computation of the full conditional posterior distributions

In such a linear setting, MCMC methods can be used to generate samples from the posterior distribution (21). MCMC methods are able to generate correlated samples of intractable probability distributions. The algorithm we present below is referred to as the blocked Gibbs sampler, which was originally called substitution sampling (Gelfand and Smith, 1990). Its sampling scheme follows the prior structure in Equation (22), but now conditional on experimental data $(\mathbf{z}^{\prime,f},\mathbf{X}^{f})$.

Blocked Gibbs sampler

Set starting values: $(b(0),\sigma^{2}(0))$ 2. 2.

Choose a number $N_{mcmc}$ of iterations, then for $1\leq k\leq N_{mcmc}$, generate in a loop:

[TABLE]

[TABLE]

If the conjugate Gaussian-inverse gamma prior distribution is put on $(b,\sigma^{2})$, then the distributions (23) and (24) are analytically tractable (see A for the proof). Such a prior is defined by

[TABLE]

where for $1\leq j\leq p$

[TABLE]

and

[TABLE]

The Greek letter $\Gamma(.)$ refers to the Euler gamma function. The hyperparameters $\psi_{j},\gamma_{j},\mu_{j},a_{j}$ are set depending on the amount of prior information on $(b,\sigma^{2})$. If a data-dominated posterior distribution is expected, typically in absence of prior information on $(b,\sigma^{2})$, the values $\mu_{j}=0$, $a_{j}=\epsilon$, $\psi_{j}=\epsilon$ and $\gamma_{j}=\epsilon$ with $\epsilon$ being close to [math] should be used (Spiegelhalter et al., 2004). However, a non negligible impact of $\epsilon$ on the shape of the posterior distribution may remain when $\mathbf{z}^{\prime,f}$ is of moderate size (Gelman, 2006).

The blocked Gibbs algorithm generates samples $\{\mathbf{A}^{\prime}(k),(b(k),\sigma^{2}(k))\}_{k}$ by a Markov chain process whose stationary distribution matches the posterior distribution (21). However, the practical convergence of such MCMC algorithms is often not obvious and need be checked (see C for a look on this issue). Posterior samples $\{(b(k),\sigma^{2}(k))\}_{k}$ can be used for computing point estimates of the posterior distribution of $(b,\sigma^{2})$ as well as posterior credible intervals. The posterior mean is

[TABLE]

Another much used point estimate is the Maximum A Posteriori (MAP)$:$

[TABLE]

Eventually, posterior samples $\{m(k)\}_{k}$ can be obtained from the posterior samples $\{b(k)\}_{k}$ by moving them from $\alpha^{\star}$.

Remark** 2**

Other MCMC algorithms than the blocked Gibbs sampler can be implemented to compute the posterior distribution of $(b,\sigma^{2})$ (see A for statistical details).

3.3 Bayesian fluctuation intervals

A nice property of the Bayesian framework is that the statistical uncertainty of $(m,\sigma^{2})$ can be taken into account in computing the fluctuation intervals of the multiplicative factors, which makes them more conservative than those computed by the standard version of CIRCE. Such intervals are derived from the posterior predictive distribution of $\Lambda=(\Lambda_{1},\cdots,\Lambda_{p})$ whose density is computed as

[TABLE]

where

[TABLE]

with $\pi(\Lambda_{j}|m_{j},\sigma_{j}^{2})$ being the log-Gaussian density of $\Lambda_{j}$ ($1\leq j\leq p$).

The predictive density $\pi(\Lambda|\mathbf{X}^{f},\mathbf{z}^{\prime,f})$ has no closed-form but Monte Carlo simulation can be used to generate samples of it. Let $N$ be the number of samples we want to generate. Then, in a loop, for $1\leq k\leq N$ :

pick up a posterior sample $m(k),\sigma^{2}(k)\thicksim m,\sigma^{2}|\mathbf{X}^{f},\mathbf{z}^{\prime,f}$, 2. 2.

generate a realization $\lambda(k)\thicksim\Lambda|m(k),\sigma^{2}(k)$.

Such a simulation scheme yields samples $\lambda_{j}(1),\lambda_{j}(2),\cdots,\lambda_{j}(N)$ to be used for kernel density estimation of (30) as well as computing empirical quantiles of $\Lambda_{j}$ ($1\leq j\leq p$). Let $q_{0.025}(\Lambda_{j})$ and $q_{0.975}(\Lambda_{j})$ be respectively the $2.5\%$ and $97.5\%$ empirical quantiles of $\Lambda_{j}$. The resulting $95\%$-fluctuation interval of $\Lambda_{j}$ is equal to

[TABLE]

3.4 Main limitation of the linear approach

On the one hand, the implementation of CIRCE (and its Bayesian version) can be subject to a strong user effect in constructing the matrix $H$ of the partial derivatives in Equation (16). As thermal-hydraulic system codes tend to switch from a flow regime (or thermal regime) to another, the values of the finite differences can be strongly impacted by the selected increment. On the other hand, when the size of $\mathbf{z}^{f}$ is moderate ($30\leq n\leq 100$), the resulting fluctuation intervals are likely to be very spread out. Hence, the linear approximations may fail to match the code outputs for realizations $\alpha$ (resp. $\lambda$) not close enough to $\alpha^{\star}$ (resp. $\lambda^{\star}$). Both of these limitations motivate us to tackle the general setting where the code outputs move non linearly against $\alpha$.

4 The non linear setting

Let $\mathbf{A}=[\alpha_{1},\cdots,\alpha_{n}]^{T}\in\mathcal{M}_{n,p}({\mathbb{R}})$ be the matrix of the missing log-samples of the multiplicative factors. As there is no longer linearization of the code outputs, we then revert to the estimation of $m$. The Bayes formula applied to Equation (9) with $\lambda=\exp{\alpha}$ gives

[TABLE]

where $\pi(m,\sigma^{2})$ can still be chosen as a Gaussian-inverse gamma prior density. Then, a blocked Gibbs sampler can be run in a similar way to the linear case.

Blocked Gibbs sampler (non-linear version)

Set starting values: $(m(0),\sigma^{2}(0))$ 2. 2.

Choose a number $N_{mcmc}$ of iterations, then for $1\leq k\leq N_{mcmc}$ generate:

[TABLE]

[TABLE]

The conditional distribution (35) is unchanged compared to the linear case (except we estimate $m$ instead of $b$) whereas the probability density of (34) is now written as

[TABLE]

which depends on the code outputs $Y_{\exp{\alpha_{i}}}(\mathbf{x}_{i}^{f})$ for $1\leq i\leq n$. This distribution has no longer closed-form and then MCMC methods such as a Metropolis-Hastings (MH) algorithm can be used for sampling from it (Metropolis et al., 1953) (see B.2).

The more time consuming the simulations are, the more time consuming the blocked Gibbs sampler will be. Because several ten of thousands of iterations are commonly required for this sampler to converge to the posterior distribution, then running it in a reasonable time, say, up to one day, is no longer feasible when the simulations require a couple of minutes or more. This handicap can be circumvented by replacing the code outputs with emulator predictions in Equation (36). We will use the Gaussian process (GP) emulator which can deliver fast predictions of time demanding simulations as well as uncertainty bounds around them (Sacks et al., 1989; Santner et al., 2003).

4.1 The Gaussian Process (GP) emulator

Let the computer code output $Y_{\lambda}(\mathbf{x}_{i}^{f})$ be a function of $\lambda\in\mathcal{L}$ with $\mathcal{L}\subset{\mathbb{R}}^{p}$. When simulations are not instantaneous, the code can be considered as an unknown quantity. Hence, from a Bayesian point of view, we can put assume the functional prior

[TABLE]

as the sum of a deterministic and a stochastic part being a second order stationary Gaussian Process (GP). The positions of $\lambda$ and $\mathbf{x}_{i}^{f}$ are reversed in Equation (37) to point out that the GP is a random function of $\lambda\in\mathcal{L}$ at a fixed site $\mathbf{x}_{i}^{f}$. Let $\mathbf{D}_{M}$ be a design of $M$ locations of the multiplicative factors$:$

[TABLE]

The corresponding set of simulations run at $\mathbf{x}_{i}^{f}$ is denoted by $Y_{\mathbf{x}_{i}^{f}}(\mathbf{D}_{M})$. The prediction of the code output at a new location $\lambda\in\mathcal{L}$ is given by the random process conditional on $Y_{\mathbf{x}_{i}^{f}}(\mathbf{D}_{M})$, which is still Gaussian$:$

[TABLE]

where the two terms in the brackets are the predictive mean and variance respectively. Their mathematical expressions are given in B.1. Equation (39) can be referred to as the posterior GP (Currin et al., 1991). The mean function in Equation (37) is commonly assumed as either constant

[TABLE]

or a linear function as the product between a vector of regression functions and regression parameters $\beta_{i}:=\beta_{\mathbf{x}_{i}^{f}}\in{\mathbb{R}}^{m}$ of regression parameters$:$

[TABLE]

The covariance function of a second order stationary GP is decomposed into a scale parameter $\sigma_{i}^{2}:=\sigma^{2}(\mathbf{x}_{i}^{f})\in{\mathbb{R}}^{+}$ multiplied by a correlation function indexed by a vector of correlation lengths l_{i}:=l(\mathbf{x}_{i}^{f})\in{\mathbb{R}}^{p}$$:

[TABLE]

The most commonly used correlation functions are those of Matérn type due to their flexibility. Let us define $\delta\lambda=|\lambda-\lambda^{\prime}|$ with $\lambda,\lambda^{\prime}\in\mathcal{L}\subset{\mathbb{R}}$ and $\Gamma(.)$ and $K_{\nu}(.)$ being the Gamma and modified Bessel function of the second kind respectively. The one dimensional Matérn function is written as

[TABLE]

with $\nu$ being a hyperparameter setting the regularity of the GP’s trajectories. This function quantifies to which extend the code outputs $Y_{\mathbf{x}_{i}^{f}}(\lambda)$ and to $Y_{\mathbf{x}_{i}^{f}}(\lambda^{\prime})$ are related to each other according to the deviation $\delta\lambda$. The higher $\nu$, the smoother the GP trajectories. When $p\geq 2$, a multidimensional Matérn function can rely on the deviation $\delta\lambda$ below$:$

[TABLE]

with $\lambda=(\lambda_{1},\cdots,\lambda_{p})^{T}$ and $l_{i}=(l_{i,1},\cdots,l_{i,p})^{T}$ being respectively the vectors of the $p$ deviations in each dimension and $p$ correlation lengths. Then, an isotropic Matérn function is constructed by replacing $\delta\lambda$ by (44) into Equation (43), which gives

[TABLE]

In real problems, $\nu$ should be properly selected according to the smoothness of the code outputs. The hyperparameters $\beta_{i}$, $\sigma_{i}^{2}$ and $l_{i}$ will be estimated from $Y_{\mathbf{x}_{i}^{f}}(\mathbf{D}_{M})$ as part of the theorem given in the next section.

4.2 The GP-based MH within Gibbs sampler

The following theorem rewrites the probability density (36) by using GP emulation.

Theorem 1

Let us assume that one GP emulator is constructed per site $\mathbf{x}_{i}^{f}\in\mathbf{X}^{f}$ for $1\leq i\leq n$. The whole set of simulations used for fitting all the GP emulators is denoted by

[TABLE]

Based on every GP at $\mathbf{x}_{i}^{f}\in\mathbf{X}^{f}$ given by Equation (39), the so-called GP-based approximation of the probability density (36) is written as

[TABLE]

with $\hat{\theta}_{i}:=(\hat{\beta}_{i},\hat{\sigma}_{i}^{2},\hat{l}_{i})$ being a point estimate of the hyperparameters of the predictive mean and variance.

Proof 1

A Bayesian modular approach is implemented in the spirit of Bayarri et al. (2007) (see B.3 for details).

This theorem provides an approximation of the density (36) in which the thermal-hydraulic simulations are substituted for the predictive means and variances of the posterior GPs. The interest of using the probability density (47) in comparison to (36) is that the computer code outputs are replaced by fast-to-evaluate Gaussian predictions. The density (47) should be interpreted as the natural GP-based emulation of (36) because, on the one hand, the code outputs are simply replaced by the predictive means of the GPs and, on the other, the predictive variances are added up to the experimental variances so that the potential departure between the predictive means and the actual code outputs is taken into account. The mismatch between $\mathbf{z}^{f}$ and the predictive means of the GPs then consists of both the experimental and the emulator uncertainties.

The MH step is now much sped up in sampling (47) instead of (36). The resulting GP-based MH within Gibbs algorithm consists of sampling the two conditional distributions (47) and (35) one after the other. The more the GP emulators can yield accurate and reliable predictions of the thermal-hydraulic simulations, the closer (47) is expected to (36) and, therefore, the closer the GP-based posterior distribution is expected to the exact posterior distribution (33).

Remark** 3**

If the predictive variances of the GP emulators are larger than or comparable to the experimental variances $\sigma_{\epsilon_{i}}^{2}$, then the GP-based posterior distribution can hardly match the exact posterior distribution.

In the next section, we present the COndensation at the Safety Injections (COSI) experimental database from which the uncertainty of two condensation models at a safety injection will be quantified.

5 Application to COSI tests

5.1 Description

The COSI database is established from a SET facility which reproduces the cold leg and the emergency core cooling injection of a nuclear power plant with a power and volumic scale of 1/100. This facility allows to characterize two heat transfer phenomena to which the two multiplicative factors $\Lambda_{1}$ and $\Lambda_{2}$ are applied$:$

$\Lambda_{1}$ is applied to the condensation created on the free surface far from the injection (Janicot and Bestion, 1993),

• 2.

$\Lambda_{2}$ is applied to the condensation due to the turbulent mixing in the vicinity of the injection jet (Gaillard et al., 2015).

The COSI facility is composed of$:$

• $\bullet$

a cold leg,

• $\bullet$

an emergency core cooling injection,

• $\bullet$

a downcomer,

• $\bullet$

a boiler to regulate the pressure in the test facility.

We have considered steady-state tests. Various thermal-hydraulic input conditions $\{\mathbf{x}^{f}_{i}\in{\mathbb{R}}^{4}\}_{i}$ have been tested, such that$:$

• $\bullet$

water height in the cold leg is between 0 and 60%,

• $\bullet$

injection temperature is between 20 and 80 degrees Celsius,

• $\bullet$

pressure in the test facility is between 20 and 70 bars,

• $\bullet$

injection flow rate is between 0.1 and 0.6 kg/s.

The instrumentation has been composed by thermocouples in the cold leg as well as a boiler flow rate measuring device to follow the temperature evolution in the liquid layer and the condensation flow rate in the facility respectively.

The vector $\mathbf{z}^{f}$ is composed by $n=50$ condensation flow rate values with respect to various conditions $\mathbf{x}_{i}^{f}$ ($1\leq i\leq n$). The standard deviation of $z_{i}^{f}$ is equal to $\sigma_{\epsilon_{i}}=10^{-4}$ for all experiments. The corresponding simulations run with the CATHARE 2 code need a couple of minutes.

5.2 Preliminary estimations using the standard CIRCE method

Both Gaussian and log-Gaussian parametrizations were compared to each other (see Equation 14 and 15 respectively). The latter gave the best approximation of the code outputs and then $\Lambda_{1}$ and $\Lambda_{2}$ are log-Gaussian distributions. Applying the CIRCE method yielded the following biases$:$

[TABLE]

As these values are far away from [math], the linear approximations made at $\alpha^{\star}=(0,0)$ were poor in the vicinity of $(\hat{b}_{1},\hat{b}_{2})$. A way of improving the robustness of the CIRCE method is to make successive linearizations until small biases are obtained. This is referred to as the iterative CIRCE method (Nouy and de Crécy, 2017). For the COSI tests, we carried out four iterations, the last being at $\alpha^{\star}=(\log{1.73},\log{0.48})=(0.55,-0.73)$. The biases computed by CIRCE at this $\alpha^{\star}$ were equal to

[TABLE]

Equation (18) then gave the resulting $95\%$-fluctuation intervals of $\Lambda_{1}$ and \Lambda_{2}$$:

[TABLE]

5.3 Bayesian estimation using the linear approximations

By using a diffuse Gaussian-inverse gamma prior distribution with $\epsilon=10^{-2}$, the blocked Gibbs sampler presented in Section 3 was run through an executable C++ macro. We generated $3\times 10^{6}$ MCMC samples $\{(b(k),\sigma^{2}(k))\}$ in some minutes, then stored into a TTree object of the TDataServer class of the Uranie platform (Blanchard et al., 2019).

The Geweke testing was applied to determine the number of iterations, called burn-in period before the chain really converges to the posterior distribution (see C). As the four $p$ values for $b_{1}$, $b_{2}$, $\sigma^{2}_{1}$ and $\sigma^{2}_{2}$ were all above $0.05$, no evidence of non convergence was detected in the first part at this level of significance. As a result, we removed the first $10\%$ samples and kept the remaining $90\%$, i.e. $2.7\times 10^{6}$ samples. Then, we got back to the posterior distribution of $(m_{1},m_{2})$ by moving of $\alpha^{\star}=(0.55,-0.73)$ the posterior samples of $(b_{1},b_{2})$. The corresponding posterior means of $m_{1}$ and $m_{2}$ were equal to

[TABLE]

implying

[TABLE]

Neither $\exp{(\hat{m}_{1})}$ nor $\exp{(\hat{m}_{2})}$ is close to $1$. The log-factor $\log{\Lambda_{1}}$ is thus positively biased whereas $\log{\Lambda}_{2}$ is negatively. The posterior means of $\sigma^{2}_{1}$ and $\sigma^{2}_{2}$ were equal to

[TABLE]

The effective sample size (e.s.s.) values were computed from the posterior samples of $m_{1}$, $m_{2}$, $\sigma_{1}^{2}$, $\sigma_{2}^{2}$ (see C). The e.s.s. values allowed us to compute the $95\%$ confidence intervals of these parameters$:$

[TABLE]

and

[TABLE]

Even though the autocorrelation between consecutive samples was moderate, only every $100$-th sample was picked up to present the histograms of $(\exp{(m_{1})},\exp{(m_{2})})$ and $(\sigma_{1}^{2},\sigma_{2}^{2})$, on Figures 1 and 2 respectively. The strong negative dependence between $m_{1}$ and $m_{2}$ means that a larger median for $\Lambda_{1}$ can be counterbalanced by a smaller one for $\Lambda_{2}$, which is in line with simulation of the COSI tests.

To finish, the predictive distribution of $\Lambda=(\Lambda_{1},\Lambda_{2})$ was sampled by the Monte Carlo algorithm explained in Section 3.3 in order to deduce the $95\%$-fluctuation intervals of $\Lambda_{1}$ and \Lambda_{2}$$:

[TABLE]

These intervals are more spread out than those given by Equation (50). This is because the uncertainty of $(m,\sigma^{2})$ is fully taken into account in the Bayesian computation. However, we put into question their validity because the accordance between the linear approximations and the actual code outputs was quite poor in the tails of the intervals (56) (see Section 3.4). This motivated to redo such a statistical analysis via the non linear setting to give more credence to the fluctuation intervals.

5.4 Bayesian estimation in the non linear setting

As seen in Section 4, the GP-based MH within Gibbs sampler requires fitting one GP emulator per input site $\mathbf{x}_{i}^{f}\in\mathbf{X}^{f}$ by means of the learning simulations $Y_{\mathbf{x}_{i}^{f}}(\mathbf{D}_{M})$ ($1\leq i\leq 50$).

The method we applied to construct $\mathbf{D}_{M}$ as well as a testing design $\mathbf{D}^{test}_{M}$ consists of the three steps below$:$

construction of a design $\mathbf{D}_{2M}$ as a Space Filling Design (SFD) with $2M$ locations $\lambda\in\mathcal{L}$. The method maximinSA_LHS from the Dice Design R library was used for doing so (Dupuy et al., 2015); 2. 2.

construction of $\mathbf{D}_{M}$ by using the wspDesign method from the same library. This method can take $M$ locations out of $\mathbf{D}_{2M}$ so that the $M$ remaining ones form a SFD of $\mathcal{L}$; 3. 3.

construction of a testing design $\mathbf{D}_{M}^{test}=\mathbf{D}\setminus\mathbf{D}_{M}$ to check the reliability of GP predictions.

[TABLE]

The choice of $\mathcal{L}$ depends on the values taken by the two multiplicative factors. Based on physical expertise on the two condensation models we chose $\mathcal{L}:=(0,5)\times(0,5)$. Figure 3 shows both the learning and testing designs with $M=500$ locations for each.

Every GP emulator at $\mathbf{x}_{i}^{f}$ ($1\leq i\leq 50$) was constructed with a constant mean (40) and the Matérn covariance function (45). Several values of $\nu$ were tested comprising $\nu=5/2$, $\nu=3/2$ and $\nu=1/2$. The Modeler module of the URANIE platform was used to estimate the vector $\theta_{i}$ of hyperparameters by maximizing the likelihood with respect to the learning simulations (see B.1). Then, both the predictive mean and variance in Equation (39) can be evaluated at any $\lambda\in\mathcal{L}$. The goodness-of-fit of every GP was assessed through the predictive coefficient $Q^{2}$ which measures the part of variance of the true code outputs explained by the predictive means, and also through the RMSE (Root Mean Squared Error) as a measure of bias in the predictions$:$

[TABLE]

with $\bar{Y}^{test}_{\mathbf{x}_{i}^{f}}$ being the mean of the $M$ simulations run over the testing design $\mathbf{D}^{test}_{M}$. Same formula can be applied to $\mathbf{D}_{M}$ using leave-one-out (loo) predictions$:$

[TABLE]

where

$\mu^{M-s}_{\mathbf{x}_{i}^{f},\hat{\theta}_{i}}(.)$ is the predictive mean computed with respect to the $M$ learning simulations except the $s$-th,

• 2.

$\bar{Y}_{\mathbf{x}_{i}^{f}}$ is the mean of the $M$ simulations run over $\mathbf{D}_{M}$.

The closer the $Q^{2}$ values are to $1$, the better the GP predictive means are. Almost all the values of $Q_{loo}^{2}(\mathbf{x}_{i}^{f})$ and $Q_{test}^{2}(\mathbf{x}_{i}^{f})$ for $1\leq i\leq 50$ are above $0.99$ regardless of the value of $\nu$. Then, Figure 4 presents the values of $RMSE_{loo}(\mathbf{x}_{i}^{f})$ and $RMSE_{test}(\mathbf{x}_{i}^{f})$. We can see that the values corresponding to $\nu=1/2$ are twice as large as those corresponding to $\nu=3/2$ and $\nu=5/2$. The latter $\nu$ are thus the most promising for point predictions.

In a second stage, the standardized prediction errors (s.p.e.) were computed to evaluate the reliability of the covariance function. The value of $\nu$ should be judged as acceptable either if the empirical distribution of the s.p.e. matches that of the standard Gaussian distribution $\Phi$ or if most of the s.p.e. are closer to [math] than expected. The latter is when the predictive variances have been overestimated. For each GP emulator at $\mathbf{x}_{i}^{f}\in\mathbf{X}^{f}$, we calculated a coverage measure giving the percentage of the s.p.e. covered by the $95\%$ bilateral interval of $\Phi$. Figure 5 displays box plots of these percentages over $\mathbf{X}^{f}$. Overall, the best agreement with $95\%$ is achieved with $\nu=3/2$ whereas $\nu=5/2$ and $\nu=1/2$ show respectively under-coverage and over-coverage. We can conclude that the predictive variances of the GP emulators fitted with $\nu=5/2$ are underestimated (under conservative predictions) and those fitted with $\nu=1/2$ are overestimated (conservative predictions).

Several GP-based MH within Gibbs samplers (see Section 4.2) were run to assess the possible impact of $\nu$ on the shape of the resulting GP-based posterior distribution. We considered three settings$:$

Setting $1$: The GP-based MH within Gibbs algorithm relies on GP emulators constructed with the value $\nu=5/2$ for every $\mathbf{x}_{i}^{f}\in\mathbf{X}^{f}$; 2. 2.

Setting $2$: The GP-based MH within Gibbs algorithm relies on GP emulators constructed with the value $\nu=3/2$ for every $\mathbf{x}_{i}^{f}\in\mathbf{X}^{f}$; 3. 3.

Setting $3$: similar to Setting 2, but the value $\nu=1/2$ is used for sites $\mathbf{x}_{i}^{f}\in\mathbf{X}^{f}$ where the coverage measure with $\nu=3/2$ is below $95\%$. This is a conservative option to avoid any underestimated predictive variances.

The Gaussian-inverse gamma prior (25) was used with $\epsilon=10^{-2}$ like in the linear case. For each setting above, $3\times 10^{6}$ MCMC samples were generated by the GP-based MH within Gibbs algorithm. As the autocorrelation was stronger than in the linear case, illustrated by a poor mixing of the samples, four chains were run in parallel from different starting locations. In Setting $1$, only the second chain passed the Geweke testing at level $5\%$. In Setting $2$, the third and fourth chains passed it. In Setting $3$, all the chains passed it except the fourth one. The Gelman and Rubin’s statistic was also computed to ensure that the four chains actually converged to the posterior distribution that they have in common, as it should in theory (see C).

Among the chains that passed the Geweke testing, we selected in each setting the one having the highest e.s.s. after cutting off the burn-in period. Tables 1 and 2 present respectively the posterior means of $m_{1}$, $m_{2}$ and $\sigma_{1}^{2}$, $\sigma_{2}^{2}$ along with the corresponding confidence intervals at $95\%$. We have expected that both $\sigma_{1}^{2}$ and $\sigma_{2}^{2}$ were higher in Settings $2$ and $3$ than in Setting $1$ because the predictive variances of the GP emulators are larger in those two settings. Actually, Table 2 shows that only the posterior mean of $\sigma^{2}_{2}$ is increased. Table 3 presents the resulting fluctuation intervals computed in each setting. The fluctuation interval of $\Lambda_{2}$ is then much more spread out in both Settings $2$ and $3$ than in Setting $1$ whereas that of $\Lambda_{1}$ does not change much (and even a little bit narrower). A preliminary sensitivity analysis of the COSI database showed that $\Lambda_{2}$ is applied to the condensation model that most impacts the total condensation flow rate. This factor is thus more affected by changes into the GP-based posterior distribution (47).

Since the GP emulators in Setting $1$ yield underestimated predictive variances, the fluctuation intervals of Setting $2$ and, even more, Setting $3$ should be preferred over that of Setting $1$. This is even more critical within the BEPU framework where conservative fluctuation intervals are wanted.

The predictive variances of the GP emulators in Settings $2$ and $3$ are on average about equal to $10^{-7}$ while the variance of $z_{i}^{f}$ is equal to $\sigma_{\epsilon_{i}}^{2}=10^{-8}$ ($1\leq i\leq 50$). Therefore, a significant gap may still exist between the fluctuation intervals reported in Table 3 and those that would have been computed from the exact posterior distribution if the CATHARE 2 simulations had been instantaneous (or at least much faster) (see Remark 3).

6 Conclusions

6.1 Summary

The paper has laid out a Bayesian methodology to quantify the uncertainty of the physical models which close the balance equations of thermal-hydraulic system codes. This work should be received as the Bayesian implementation of the CIRCE method in that the joint posterior probability distribution of $(m,\sigma^{2})$ is computed instead of a point estimate $(\hat{m},\hat{\sigma}^{2})$. The interest in doing so is that the statistical uncertainty of $(\hat{m},\hat{\sigma}^{2})$ can be taken into account in computing the fluctuation intervals of the multiplicative factors applied to the physical models. The fluctuation intervals are then expected to be more spread out than those derived from the standard CIRCE method, which is particularly relevant for BEPU studies where some degree of conservatism is required in assessing uncertainties. In addition, we have tackled the extended Bayesian framework where the linear assumption of the code outputs with respect to the factors has been dropped.

An efficient MCMC algorithm, known as the blocked Gibbs sampler, has been implemented both in the linear and non linear settings to generate samples of the posterior distribution of $(m,\sigma^{2})$. As the cost of this algorithm strongly increases with the time duration of the thermal-hydraulic simulations, a Gaussian process(GP)-based blocked Gibbs sampler has been implemented in the non linear setting. Unlike other surrogate models such as neural networks, polynomial chaos expansion, the GP emulator is particularly suited for computer experiments because it can interpolate the learning simulations as well as providing Gaussian predictions of the code outputs at any other realization of the multiplicative factors. The GP-based blocked Gibbs sampler relies on a conditional posterior distribution of the missing (log)-values of the multiplicative factors in which the predictive means and variances of the GP emulators replace the expensive simulations. This sampler has been applied to the COSI tests to quantify the uncertainty of two condensation models at a safety injection.

A great attention has been paid on the choice of the covariance function so that the actual smoothness of the code outputs is captured as best as possible. Several Matérn covariance functions have been tested on the COSI simulations. The smaller the smoothness parameter $\nu$, the less smooth the Gaussian process trajectories. The values $\nu=1/2$ and $\nu=3/2$ have yielded much more robust predictions of the condensation flow rate than the value $\nu=5/2$ for which the predictive variances have been underestimated. As a consequence, the fluctuation intervals of the multiplicative factors computed with $\nu=5/2$ (Setting $1$) have been different than those computed with a more suited $\nu=3/2$ (Setting $2$) or even a more conservative option (Setting $3$).

The more the GP predictions match the true values of the code outputs, the more the GP-based posterior distribution is expected to be close to the exact posterior distribution (33) in some sense, for example, according to the Kullback-Leibler divergence. We could surely prove such a convergence property in the same way than done in Damblin et al. (2018). The proof would rely on the theory of kernel interpolation (Schaback, 1995). From a practical point of view, however, when the number $M$ of learning simulations is moderate, it remains challenging to guaranty that the actual gap between both of those distributions is negligible or even not that large.

6.2 Future works

According to the previous discussion, the improvement of the accuracy of the GP emulators should be aimed for. This can simply be done by increasing the number of learning simulations so that the predictive variances are significantly down in comparison with the variance of the experimental error. However, the computational cost of the GP-based Gibbs algorithm would increase (see B.4). We could suggest two ways for facing this limitation$:$ evaluating the GP emulators in parallel at every iteration of the Gibbs sampler and speeding up GP evaluations for large datasets (Gramacy and Apley, 2015; Rullière et al., 2017).

Another tough question concerns the choice of the prior distribution for the variances $\sigma^{2}$ of the multiplicative factors. If few informative inverse-gamma priors are put on such parameters, as we did in the paper, then the significant prior mass near [math] does not vanish in the posterior distribution. This can lead to a Bayesian procedure that is not as much as data-dominated as expected, especially in small data situations. Some possibly more suited prior distributions are promoted in (Gelman, 2006) and should be tested.

Finally, because the posterior samples generated by the GP-based blocked Gibbs algorithm are highly correlated, other proposal densities for the Metropolis-Hastings algorithm could be tested or even advanced samplers such a Delayed Rejection Adaptive Metropolis (DRAM) could be implemented (Haario et al., 2006). In this way, we should improve the mixing speed of the Markov chain and, therefore, the precision of point estimates as well as resulting fluctuation intervals of the multiplicative factors.

7 Acknowledgments

This work has been funded by the tripartite project devoted to Uncertainty Quantification, consisting of French Alternative Energies and Atomic Energy Commission (CEA), Electricity of France (EDF) and Framatome (FRA). The authors thank the URANIE team for guidance in using the software and careful reading of the paper, as well as Lucia Sargentini and Alberto Ghione for fruitful discussions on the COSI tests. Finally, we would like to thank the two anonymous reviewers who help improving the content of the paper significantly.