CEA, DES, IRESNE, DER A Marrel

B Iooss

V Chabridon

The ICSCREAM methodology: Identification of penalizing configurations in computer experiments using screening and metamodel -Applications in thermal-hydraulics

In the framework of risk assessment in nuclear accident analysis, best-estimate computer codes, associated to a probabilistic modeling of the uncertain input variables, are used to estimate safety margins. A first step in such uncertainty quantification studies is often to identify the critical configurations (or penalizing, in the sense of a prescribed safety margin) of several input parameters (called "scenario inputs"), under the uncertainty on the other input parameters. However, the large CPU-time cost of most of the computer codes used in nuclear engineering, as the ones related to thermal-hydraulic accident scenario simulations, involve to develop highly efficient strategies. This work focuses on machine learning algorithms by the way of the metamodel-based approach (i.e., a mathematical model which is fitted on a small-size sample of simulations). To achieve it with a very large number of inputs, a specific and original methodology, called ICSCREAM (Identification of penalizing Configurations using SCREening And Metamodel), is proposed. The screening of influential inputs is based on an advanced global sensitivity analysis tool (HSIC importance measures). A Gaussian process metamodel is then sequentially built and used to estimate, within a Bayesian framework, the conditional probabilities of exceeding a high-level threshold, according to the scenario inputs. The efficiency of this methodology is illustrated on two high-dimensional (around a hundred inputs) thermal-hydraulic industrial cases simulating an accident of primary coolant loss in a pressurized water reactor. For both use cases, the study focuses on the peak cladding temperature (PCT) and critical configurations are defined by exceeding the 90%-quantile of PCT. In both cases, the ICSCREAM methodology allows to estimate, by using only around one thousand of code simulations, the impact of the scenario inputs and their critical areas of values.

Introduction

In the framework of risk assessment in nuclear accident analysis, best-estimate computer codes are increasingly used to understand, model and predict physical phenomena and, ultimately, estimate safety margins. These codes (also known as "numerical simulators") usually take a large number of input parameters driving the phenomenon of interest or related to its physical and numerical modeling. The available information about some of these parameters is often limited or uncertain. The uncertainties come mainly from the lack of knowledge about the underlying physics and about the characterization of the input parameters of the model (e.g., due to the lack of experimental data) [START_REF] Baccou | SAPIUM: A generic framework for a practical and transparent quantification of thermal-hydraulic code model input uncertainty[END_REF]. Other additional sources of uncertainty can be considered such as the choice of a particular accident scenario. All these input parameters, and consequently the simulator output, are uncertain. In this context, it is essential to take the uncertainties tainting the results of computer simulations into account [START_REF] May | Probabilistic Methods for Evaluating Operational Transient Margins and Uncertainties[END_REF]. This major step for safety studies is referred to as "uncertainty quantification of numerical models" [START_REF] Smith | Uncertainty quantification[END_REF] in the statistical community and is called "Best-Estimate Plus Uncertainty" (BEPU) in nuclear safety analysis [START_REF] Wilson | Historical insights in the development of Best estimate Plus Uncertainty safety analysis[END_REF][START_REF] Sermer | Statistical Foundation for Decision Making in Nuclear Safety-Related Problems Using Best Estimate Plus Uncertainty Analyses[END_REF].

In this work, we focus on the identification of the penalizing configurations (corresponding to critical values of the output) of specific scenario inputs (i.e., defining the accident scenario considered), under the uncertainty of the other inputs. It is important to notice that the goal is not to perform a "worst-case" scenario analysis as often encountered in common nuclear safety practices. The idea here is to look for the inputs' domain leading to a high probability level that the physical variable of interest reaches some large values. Such an approach differs from finding the inputs' values which lead to the maximum of the variable of interest (worst-case analysis). In the uncertainty quantification literature, several families of methods can be found for this particular analysis: on the one hand, it is called "inversion" or "identification of an excursion set" [START_REF] Picheny | Adaptive Designs of Experiments for Accurate Approximation of a Target Region[END_REF][START_REF] Chevalier | Fast kriging-based stepwise uncertainty reduction with application to the identification of an excursion set[END_REF] when no uncertainty is considered; on the other hand, it is called "robust inversion" [START_REF] Richet | Inversion Algorithm for Civil Flood Defense Optimization: Application to Two-Dimensional Numerical Model of the Garonne River in France[END_REF] if the uncertainties about other inputs are considered. In the present paper, the proposed methodology rather belongs to the second type of methods. As mentioned in [START_REF] Richet | Inversion Algorithm for Civil Flood Defense Optimization: Application to Two-Dimensional Numerical Model of the Garonne River in France[END_REF], the real industrial application tackled in the present paper goes beyond standard canonical problems (e.g., purely optimization or inversion problems) and thus implies to develop a dedicated methodology. Our study is motivated and guided by the "Intermediate Break Loss Of Coolant Accident" (IB-LOCA) safety analysis, based on the numerical simulation of an accident of primary coolant loss in a pressurized water reactor [START_REF] Sanchez-Saez | Uncertainty analysis of large break loss of coolant accident in a pressurized water reactor using non-parametric methods[END_REF]. In the present paper, a realistic reactor-scale modeling of an IB-LOCA is considered [START_REF] Charignon | CathSBI, A new methodology for the revised French LOCA rules[END_REF][START_REF] Larget | How to bring conservatism to a BEPU analysis[END_REF]. This model, relying on the CATHARE2 code [START_REF] Geffraye | CATHARE2 V2.5_2: A single version for various applications[END_REF], enables to compute thermal-hydraulic transients with respect to a very high number (compared to previous simplified studies [START_REF] Iooss | Advanced methodology for uncertainty propagation in computer experiments with large number of inputs[END_REF]) of input uncertain parameters.

Standard approaches, related to our objectives and used in the nuclear industry (in particular on the IB-LOCA case), are for example:

• the inverse uncertainty quantification methods [START_REF] Mendizábal | Bayesian perspective in BEPU licensing analysis[END_REF][START_REF] Baccou | Development of good practice guidance for quantification of thermal-hydraulic code model input uncertainty[END_REF][START_REF] Damblin | Bayesian inference and non-linear extensions of the CIRCE method for quantifying the uncertainty of closure relationships integrated into thermal-hydraulic system codes[END_REF] to obtain the probabilistic distributions of some physical model inputs when experimental data are available,

• the Wilks method [START_REF] Wallis | Uncertainties and probabilities in nuclear reactor regulation[END_REF][START_REF] Petruzzi | Approaches, relevant topics, and internal method for uncertainty evaluation in predictions of thermal-hydraulic system codes[END_REF] to infer high quantiles of the model output variables of interest,

• the PLI method [START_REF] Larget | How to bring conservatism to a BEPU analysis[END_REF][START_REF] Iooss | BEPU robustness analysis via perturbed-law based sensitivity indices[END_REF] to detect the most influential inputs to be penalized amongst a large number of variables before applying the Wilks method.

• and finally the RIPS (Reduction of the Interval of variation of the Parameters of the Scenario) method, recently published in [START_REF] Lecoy | RIPS, a Statistical Method for Characterizing the Limiting Scenario in a BEPU Approach[END_REF], aims to characterize the limiting scenario in a BEPU approach.

In particular, RIPS allows to analyze the high-order (or low-order) quantiles of the output cumulative distribution function and determine, for each scenario input, the critical zone within its variation interval. A first issue of this method is that it relies on a quite subjective visual analysis, but its most important drawback is due to the intrinsic complexity in the tuning of the method.

Our goal is to provide a more robust and automatic methodology which allows to reduce the computational cost (in terms of number of code runs) compared to the RIPS approach. In uncertainty quantification studies, to solve the cost issue, a widely accepted approach consists in approximating the CPU-time expensive computer models by CPU-time inexpensive mathematical functions called "metamodels" (or "surrogate" models). These metamodels can be based, for instance, on polynomials, neural networks or Gaussian processes [START_REF] Fang | Design and modeling for computer experiments[END_REF]. The metamodel is built from a set of computer code simulations, and must be as representative as possible of the "true" code in the variable domain of the uncertain parameters while having good prediction capabilities. Nowadays, metamodels are extensively used in several engineering fields to solve industrial issues as it provides a multi-purpose tool [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF]: once fitted, the metamodel can be used, in conjunction with the costly computer code, to perform sensitivity analysis, as well as uncertainty propagation, optimization, or calibration studies see, e.g., [START_REF] Khuwaileh | Surrogate based model calibration for pressurized water reactor physics calculations[END_REF][START_REF] Wu | Kriging-based inverse uncertainty quantification of nuclear fuel performance code BISON fission gas release model using time series measurement data[END_REF][START_REF] Whyte | Surrogate model optimization of a "micro core" PWR fuel assembly arrangement using Deep Learning models[END_REF][START_REF] Christian | Dynamic PRA-Based Estimation of PWR Coping Time Using a Surrogate Model for Accident Tolerant Fuel[END_REF][START_REF] Puppo | Failure identification in a nuclear passive safety system by Monte Carlo simulation with adaptive Kriging[END_REF]).

An important issue is that the building process of the metamodel remains complex in the case of high-dimensional (e.g., typically several tens of inputs) numerical experiments. In order to build a metamodel in a efficient manner in such cases, [START_REF] Iooss | Advanced methodology for uncertainty propagation in computer experiments with large number of inputs[END_REF] proposed a methodology which combines several pre-existing advanced statistical tools: first, an initial space-filling design of experiments; second, a screening step (which aims at detecting influential inputs and non-influential ones) in order to reduce the dimension, and third, a specific building strategy (based on a sequential inclusion of variables) of a joint Gaussian process metamodel. Then, the resulting joint Gaussian process metamodel is used to accurately estimate high-order output quantiles. The efficiency of the methodology has been illustrated on a simplified IB-LOCA use case with 27 inputs and a total budget of 500 code simulations.

The objectives of the present study are different from those identified in [START_REF] Iooss | Advanced methodology for uncertainty propagation in computer experiments with large number of inputs[END_REF]. In the present paper, the goal is to perform an inversion and thus requires the new proposed statistical methodology, called ICSCREAM (pronounced "ice-cream"), for "Identification of penalizing Configurations using SCREening And Metamodel". Indeed, finding the penalizing configurations in input with respect to critical output values (e.g., typically, a high-order quantile) corresponds to identify specific areas of some specific inputs, while ignoring the uncertainty tainting the other ones. Moreover, since the considered use case is a reactor-scale model of an IB-LOCA, the number of input parameters is considerably larger than for standard mock-up cases. Parallel to that, the available simulation budget is fixed (around a thousands simulations). Therefore, a methodological challenge arises in the metamodel building process.

The paper is organized as follows. In Section 2, the general workflow of the IC-SCREAM methodology is first detailed. Two different IB-LOCA use cases and their corresponding datasets are then presented in Section 3. Thereafter, Sections 4 to 6 are dedicated to each step of the method: their theoretical foundations and their practical implementation are first detailed before presenting the results obtained on the two datasets. Section 7 provides some conclusions and prospects of this work. The Appendix (Section 7) provides an additional illustration of the methodology through an analytical test case.

ICSCREAM: a four-step-based methodology

For the sake of clarity, a few notations are introduced at this stage. Throughout the rest of this paper, the numerical model (i.e., computer code or simulator) is represented by the following input-output relationship:

M : X -→ Y X -→ Y = M (X) (1)
where the uncertain output variable Y and the d input parameters X = (X 1 , . . . , X d ) belong to some measurable spaces respectively denoted by Y and X ⊂ R d . As part of the probabilistic approach, the inputs are considered as random variables with probability distributions denoted by P X on X [START_REF] Helton | Uncertainty and sensitivity analysis in the presence of stochastic and subjective uncertainty[END_REF][START_REF] Oberkampf | Mathematical representation of uncertainty[END_REF]. Moreover, among the d uncertain inputs, the d pen scenario inputs (that need to be penalized) are denoted by X pen ⊂ X.

Notice that no general assumption about independence between the inputs is required here. Only the X pen are assumed to be independent from {X \ X pen }, i.e., from the remaining inputs. Dependency among the X pen has to be considered conjointly when computing the conditional probabilities for the identification of penalizing configurations (see Section 6). Then, it is supposed that only input-output observations (or realizations) of M are available. It is therefore assumed that we have a n-size sample of inputs and associated outputs denoted by (X s ,Y s ) where X s = {x (1) , . . . , x (n) } with

x (i) = x (i) 1 , . . . , x (i) d
denotes the matrix of n-size sample locations (also called the "experimental design"), and Y s = {y (1) , . . . , y (n) } the corresponding outputs observations with y (i) = M (x (i) ).

To tackle the problem of large input dimension and complex computer codes, the ICSCREAM methodology combines several statistical techniques which have been chosen for their relevance regarding the multiple constraints imposed by the problem of interest (i.e., both the dimensionality and the limited size of the learning sample). Figure 1 provides a general workflow of the main steps:

• Step 1: Design of experiments. Knowing the variation domain of the input variables X, a design of n numerical experiments is firstly performed to obtain the learning sample (X s ,Y s ). To apply statistical theory and methods, this sample must be chosen randomly, from Monte Carlo sampling techniques. The simplest choice is a crude Monte Carlo sample. However, "stratified" sampling techniques (such as Latin Hypercube sampling [START_REF] Loh | On Latin hypercube sampling[END_REF]) can also be used to ensure a more regular sampling of marginal distributions and a better convergence of statistical estimators. Another particularly relevant design strategy would be to use metamodel-based adaptive designs [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF]. This solution would allow to select, in a sequential way, new simulation points (in the inputs' domain) in order to improve the global accuracy of the metamodel or directly the estimation of penalizing configurations. However, this solution cannot be used in our case: we have to rely on a unique finite batch of random CATHARE simulations (performed at the beginning of the process). The objective is to provide an operational tool for the engineer without any extra interfacing with the simulation code and which is able to deal with a constrained simulation budget.

• Step 2: Preliminary screening and ranking with global and target HSIC measures. From the learning sample, a screening analysis is performed via statistical independence tests based on Hilbert-Schmidt Independence Criterion (HSIC) measures [START_REF] Lozzo | New improvements in the use of dependence measures for sensitivity analysis and screening[END_REF]. For this, HSIC are considered in global and target versions, "target" referring here to the area where the output exceeds a given critical value (i.e., a threshold). From the results of HSIC and target HSIC-based tests, the primary influential inputs (PII) and denoted X PII are identified. These PII are the most influential inputs regarding the output variability and can be ranked by decreasing influence. Another set of inputs of lesser (i.e., secondary) influence is gathered in X SII . Finally, all the remaining inputs are considered as global stochastic (i.e., unknown) inputs, denoted by X ε .

Unlike in a previous work [START_REF] Iooss | Advanced methodology for uncertainty propagation in computer experiments with large number of inputs[END_REF], HSIC are also used here in a goal-oriented version [START_REF] Marrel | Statistical developments for target and conditional sensitivity analysis: application on safety studies for nuclear reactor[END_REF] to tackle the objective of conditional probability estimation. Moreover, statistical independence tests based on HSIC are used to ensure a more rigorous interpretation in terms of ranking and screening, while taking into account the estimation error related to the limited sample size. Both global and target HSIC must then be combined into aggregated tests to provide a unique screening and ranking.

• Step 3: Building and validation of Gaussian process (Gp) metamodel. From the learning sample, a Gp metamodel is built to fit the simulator output Y , taking advantage of the results provided by Step 2. More precisely, one considers the ordered explanatory inputs X exp of the Gp metamodel. Basically, they consist in the X PII to which we add the inputs that have to be penalized (if not already selected in X PII ): X exp = {X PII ∪ X pen }. Note that the X SII are also included as explanatory variables but in a coarser way. Finally, the residual effect of the neglected inputs, merged into X ε , is captured using an additional noise effect. Concerning the estimation of all the Gp hyperparameters, similarly to [START_REF] Iooss | Advanced methodology for uncertainty propagation in computer experiments with large number of inputs[END_REF], a sequential process (which relies on the ranking deduced from the preliminary screening step) is used. More precisely, the information provided by HSIC-based tests is adapted in order to address the issue of the large input dimension. The inputs of lesser influence X SII , previously neglected in [START_REF] Iooss | Advanced methodology for uncertainty propagation in computer experiments with large number of inputs[END_REF], are now injected at the end of the sequential process to reduce the risk of mis-selection (which is common with such sparse data), while limiting the training efforts of the Gp metamodel. Finally, the accuracy and prediction capabilities of the Gp metamodel are controlled, either on a test sample (if available) or by cross-validation based on the learning sample.

• Step 4: Use of the metamodel for the identification of penalizing configurations. Once built, the Gp can be evaluated in prediction in a very intensive way (e.g., from 10 5 to 10 6 evaluations). Usually, a Gp metamodel can be used for a quantitative sensitivity analysis [START_REF] Veiga | Basics and trends in sensitivity analysis[END_REF] based on the output variance decomposition (e.g., estimation of Sobol' indices), as well as uncertainty propagation (e.g., estimation of high-/low-order quantiles or probabilities), leading in all cases to a large gain of computation time (see [START_REF] Iooss | Advanced methodology for uncertainty propagation in computer experiments with large number of inputs[END_REF] for an example in nuclear engineering).

Here, the Gp metamodel is used to estimate, within a Bayesian framework, the conditional probabilities of exceeding the critical quantile value, as a function of all the inputs that need to be penalized and ignoring the other uncertain ones. Thus, if the penalized inputs are independent, one will only have to compute a conditional probability per input and, in the case of dependent penalized inputs, this probability will turn into a probability hypersurface.

Again, the final purpose of the present paper differs from the one in [START_REF] Iooss | Advanced methodology for uncertainty propagation in computer experiments with large number of inputs[END_REF]: some modifications in the methodology were required, such as the use of target HSIC measures. Moreover, as it will be described in the next section, the datasets treated in the present work go far beyond by considering much more complex reactor-scale datasets with around hundred inputs (in contrast to the simplified mock-up use case with less than thirty inputs).

As stated before, the motivations for such a methodology are, before anything else, driven by an industrial application illustrated by two use cases described in the following section. However, it is important to notice that this methodology follows a logical modular sequence which makes it rather generic (in the sense that it can be applied to a large class of problems). To emphasize this, an illustration through an analytical test case is provided in Appendix (Section 7).

Description of the thermal-hydraulic use cases and datasets

In support of regulatory work and nuclear power plant design and operation, safety analysis considers the so-called "loss of coolant accident" which takes into account a double-ended guillotine break with a specific size piping rupture. More precisely, one considers here an intermediate break loss of coolant accident (IB-LOCA) [START_REF] Sanchez-Saez | Uncertainty analysis of large break loss of coolant accident in a pressurized water reactor using non-parametric methods[END_REF] in a three-loops pressurized water reactor from the French fleet1 . The output variable of interest is the maximal peak cladding temperature (PCT) of the hot rod during the accident transient. The reactor coolant system minimum mass inventory and the PCT are obtained shortly after the beginning of the accumulators' injection [START_REF] Charignon | CathSBI, A new methodology for the revised French LOCA rules[END_REF]. The IB-LOCA transients are simulated using the thermal-hydraulic system code (two phase flow six equations) CATHARE2 [START_REF] Geffraye | CATHARE2 V2.5_2: A single version for various applications[END_REF], jointly developed by CEA, EDF, Framatome and IRSN. Around a hundred inputs (vector X) are considered as uncertain and can be split into three different types [START_REF] Larget | How to bring conservatism to a BEPU analysis[END_REF]:

• Type 1: the boundary and initial conditions whose probability distributions are supposed to be known (as uniform or normal distributions);

• Type 2: the model parameters (e.g., models related to two-phase flow hydraulics, models associated to heat transfer and models describing the clad behavior). Their distributions can be obtained from experimental data, expert knowledge or recovered by assimilation of experimental databases [START_REF] Baccou | SAPIUM: A generic framework for a practical and transparent quantification of thermal-hydraulic code model input uncertainty[END_REF]. This leads to uniform, log-uniform, normal or log-normal distributions;

• Type 3: the scenario parameters (see their description below) which cover some variability between minimal m and maximal M bounds (they follow a uniform distribution that will be noted U[m, M] in the following).

All the inputs of types 1 and 2 are independent. Their detailed list, as well as their probability distribution and range of variations, are not given due to industrial confidentiality. The scenario inputs (type 3) are independent with those of the two other groups but some dependency might exist among them. Usually, the inputs corresponding to type 3 have to be taken at their worst-case values (corresponding to the maximal value which can be reached by the PCT) for safety demonstration [START_REF] Charignon | CathSBI, A new methodology for the revised French LOCA rules[END_REF][START_REF] Lecoy | RIPS, a Statistical Method for Characterizing the Limiting Scenario in a BEPU Approach[END_REF][START_REF] Larget | How to bring conservatism to a BEPU analysis[END_REF]. In practice, these worst-case values are unknown and only a domain of variation of each input is given. Thus, the idea is to find the penalizing values for these scenario inputs.

The following use cases are considered:

• Use-case #1 which corresponds to dataset IB-LOCA 1 with d = 96 uncertain inputs, including d pen = 2 inputs which have to be penalized: the size of the break, denoted by X 92 in the dataset (X 92 ∼ U[3, 4.2] inches), and the stopping time of the primary pumps, X 94 . X 94 follows a uniform distribution with a range of variation depending on the size of the break (globally included in [500, 1200] s). It is a small IB-LOCA case without "loss of offsite power" (LOOP) including a late main coolant pump coast down;

• Use-case #2 which corresponds to dataset IB-LOCA 2 with d = 97 uncertain inputs, including d pen = 10 independent scenario inputs: hot rod peaking factor elevation ( Note that the numbering of inputs is not the same for both use cases. Monte Carlo samples of n 1 = 889 and n 2 = 1000 CATHARE2 simulations are available for IB-LOCA 1 and IB-LOCA 2 use cases respectively. For both, the inputs are drawn according to their prior probability distributions. The histograms of the obtained values for the output of interest, namely the PCT, are given in Figure 2 (temperature is in • C). A kernel density estimation plot [START_REF] Parzen | On Estimation of a Probability Density Function and Mode[END_REF] is also added to provide a graphical illustration of the probability density function. Note that the number n of simulations is a compromise between the CPU time required for each simulation and the number of inputs. For uncertainty propagation and metamodel-building purpose, some rules of thumb propose to choose n at least as large as ten times the input dimension d [START_REF] Loeppky | Choosing the sample size of a computer experiment: A practical guide[END_REF][START_REF] Marrel | An efficient methodology for modeling complex computer codes with Gaussian processes[END_REF]. For the sake of clarity, one reminds that the β -order quantile (or β -quantile) is defined by q β = inf{y ∈ R | β ≤ F Y (y)}. From the learning sample, the empirical 90%-quantile of PCT is estimated to q (1) 0.9 = 673.18 • C (resp. q

X 1 ∼ U[2.4, 3.2] m), hot rod burn-up (X 14 ∼ U[515, 59000] MWj/t), discharge
(2) 0.9 = 671.85 • C) for IB-LOCA 1 (resp. IB-LOCA 2 ). As aforementioned, the ICSCREAM methodology aims to identify the values of the inputs X pen which yield to a high probability of exceeding this quantile. Figure 3 shows the convergence plots of the estimated empirical quantile, from bootstrap method (y-axes are scaled to the range of variation of PCT). In addition, the coefficients of variation of q 0.9 are estimated to be less than 1% (0.76% and 0.38% for IB-LOCA 1 and IB-LOCA 2 , respectively). The convergence of other statistics of PCT distribution is also verified with a coefficient of variation lower than 1% for the mean and around 2% for the standard deviation.

To illustrate the complexity of fitting the output according to such a large number of input, some scatter plots of the PCT (for IB-LOCA 2 dataset) with respect to some scenario inputs are displayed in Figure 4. In this case, only one scenario input, namely X 14 (hot rod burn-up), seems to have a clearly detectable influence on the PCT.

Theoretical and methodological details

An initial screening is performed directly from the learning sample in order to identify the primary and secondary influential inputs (resp. X PII and X SII in Figure 1). The objective is twofold: first, screening the inputs to reduce the dimension before the metamodeling step and second, ranking them by decreasing order of influence. To achieve it, we use the HSIC importance measures introduced by [START_REF] Gretton | Measuring statistical dependence with Hilbert-Schmidt norms[END_REF] and theoretically built upon cross-covariance operators in reproducing kernel Hilbert spaces (RKHS). In a nutshell, the underlying idea of HSIC measures is to capture the influence of an input X i on the output Y through the measure of their statistical dependence: the more "correlated" (in a specific sense that will be explained further) the two variables are, the more influential X i is. As highlighted by [START_REF] Veiga | Global sensitivity analysis with dependence measures[END_REF] or [START_REF] Lozzo | New improvements in the use of dependence measures for sensitivity analysis and screening[END_REF], HSIC measures show many advantages for global sensitivity analysis purposes [START_REF] Veiga | Basics and trends in sensitivity analysis[END_REF]. Some of them are briefly recalled hereafter. First, HSIC amounts to considering covariance between "feature" functions (or transformations) applied to each input-output couple of variables (i.e., an input X i and the output Y ). The set of candidate functions (possibly nonlinear) which can be applied is defined by a functional space and an associated kernel. Such a space can be of infinite dimension and thus allows to capture a very broad spectrum of forms of dependency between X i and Y . The HSIC, which is defined as the Hilbert-Schmidt (HS) norm of the cross-covariance operator (denoted by COV(•)), somehow "summarizes" the set of covariances between feature functions:

HSIC(X i ,Y ) = ||COV X i Y || HS = ∑ l,m |COV(u l (X i ), v m (Y ))| 2 ,
where (u l ) l≥0 and (v m ) m≥0 are orthonormal bases of the RKHS associated to X i and Y respectively. Furthermore, the so-called "kernel trick" allows to get rid of an explicit expression of the features (see [START_REF] Gretton | Measuring statistical dependence with Hilbert-Schmidt norms[END_REF] for more details). Thus, HSIC can be directly expressed as a linear combination of expected values of kernels and estimated in a very simple way at rather low cost (typically, a few hundred simulations). By denoting k i and k the kernels of the two RKHS respectively associated to X i and Y , it can be shown (see [START_REF] Gretton | Measuring statistical dependence with Hilbert-Schmidt norms[END_REF]) that the HSIC measure can be easily estimated by:

HSIC(X i ,Y ) = 1 n 2 Tr(L i HLH), (2)
where Tr(•) is the trace operator, L i and L are two Gram matrices defined by

L i = (k i (X (l) i , X (m) i )) 1≤l,m≤n and L = (k(Y (l) ,Y (m) )) 1≤l,m≤n with (X (m)
i ,Y (m) ) 1≤m≤n a n-sample of (X i ,Y ), and H = (δ lm -1/n) 1≤l,m≤n with δ lm the Kronecker operator.

Another key aspect relates to the choice of kernels. Characteristic kernels (such as, among others, the Gaussian one), allow to fully characterize the independence between the input X i and the output Y with HSIC. This leads to the following fundamental property: the nullity of HSIC is equivalent to independence between X i and Y . Therefore, one can build dedicated statistical independence tests for screening purposes [START_REF] Lozzo | New improvements in the use of dependence measures for sensitivity analysis and screening[END_REF]. The Gaussian kernel, widely used for real variables, is defined by 2 and parametrized by the bandwidth parameter λ , often set at λ = 1/σ 2 z with σ 2 z the empirical variance of the sample of Z. Thus, for a given input X i , statistical HSIC-based tests aim at testing the null hypothesis "H (i) 0 : X i and Y are independent", against its alternative "H (i) 1 : X i and Y are dependent". The significance level (i.e., the probability of rejecting the null hypothesis H 0 when it is true) of this test is hereinafter noted α test and usually set at 5% or 10%. Depending on the size n of the considered sample, several versions of HSIC-based tests are available: asymptotic versions (i.e., for large n) based on an approximation with a Gamma law [START_REF] Gretton | A kernel statistical test of independence[END_REF], spectral extensions and permutation-based versions [START_REF] Lozzo | New improvements in the use of dependence measures for sensitivity analysis and screening[END_REF] or adaptive strategies [START_REF] El Amri | Optimized HSIC-based tests for sensitivity analysis: application to thermal-hydraulic simulation of accidental scenario on nuclear reactor[END_REF] for non-asymptotic cases (i.e., small n). The p-value is the probability of obtaining HSIC values either equal to, or larger than, the observed HSIC (estimated from the learning sample) assuming that H 0 true (i.e., that X i and Y are independent). Such a p-value is used to decide whether to reject (or not) the null hypothesis. Beyond the screening task, the p-value of independence tests can be quantitatively interpreted for ranking the PII, since it can be viewed as a "margin" from independence. The lower the p-value, the stronger H (i) 0 is rejected and the higher the influence of X i . Furthermore, remember that the final objective of ICSCREAM is to identify the penalizing configurations and, more precisely, to accurately identify the critical input areas where the PCT exceeds q 0.9 (i.e., Y > q 0.9 ). Consequently, we consider an additional "target" sensitivity analysis based on so-called "target HSIC" (T-HSIC) indices [START_REF] Marrel | Statistical developments for target and conditional sensitivity analysis: application on safety studies for nuclear reactor[END_REF]. Applied here, T-HSIC and associated independence tests aim at measuring the influence of an input X i over the occurrence of the event {Y > q 0.9 }. T-HSIC is built with specific kernels and a weight function can be used in order to cope with the possible loss of information around the threshold (see [START_REF] Marrel | Statistical developments for target and conditional sensitivity analysis: application on safety studies for nuclear reactor[END_REF] for further details).

λ (z, z ) = exp -λ 2 (z -z )
Finally, both HSIC and T-HSIC tests' results are aggregated to extract a unique screening and ranking: heuristic choices [START_REF] Chabridon | Tools for global and target sensitivity analyses in the context of high-dimensional thermal-hydraulic numerical experiments[END_REF] or more robust alternatives (e.g., Bonferroni's correction, p-values combination as inspired from [START_REF] Heard | Choosing between methods of combining p-values[END_REF]) can be used. More details on Step 2 of the ICSCREAM methodology are given in [START_REF] Chabridon | Tools for global and target sensitivity analyses in the context of high-dimensional thermal-hydraulic numerical experiments[END_REF].

Application on IB-LOCA uses-cases

From the learning sample of the two IB-LOCA use cases, HSIC-independence tests are applied with a permutation-based approach for the estimation of p-values. The number of permutations is optimized using the algorithms proposed by [START_REF] El Amri | Optimized HSIC-based tests for sensitivity analysis: application to thermal-hydraulic simulation of accidental scenario on nuclear reactor[END_REF]. The results are given in Figures 5 and6. They illustrate the screening process based on p-values: variables whose p-values are below the level α test = 0.05 correspond to influential variables (red dots). In both use cases, a little bit more than twenty inputs are identified as influential (i.e., p-values < α test ) during the screening step following the aggregation of global and target HSIC-based tests. About five additional inputs of lesser influence are added in X SII (p-values in ]α test ; 2α test ]).

Concerning the ranking now, for IB-LOCA 1 , the two inputs to be penalized are the most influential: the stopping time of primary pumps (X 94 ) being the most influential one, followed by the break size (X 92 ). Then five other inputs are strongly influential, namely the upper plenum and core interfacial friction, core interfacial friction, the hot spot for the hot rod and two inputs relative to accumulators. A group of thirteen other variables of lower influence is also selected by global HSIC-tests. Similar results are obtained with T-HSIC-based tests, except that two additional inputs, namely the diphasic degradation law of pumps and the residual power, are selected as very influential.

For IB-LOCA 2 , a hard core of around sixteen variables are identified by all the tests with zero p-values. The most influential are the burn-up of the hot rod (X 14 ), the heat transfer coefficient of the wall-steam exchange of the hot rod (X 46 ) and the bubbly-slug in the core during blowdown phase, related to the interfacial friction (X 80 ). To these sixteen variables, five variables are added (p-values between 0 and 0.05). Regarding the final objective of penalizing the ten scenario parameters, only two of them, namely X 14 and, with a lesser influence, the liquid accumulator enthalpy (X 36 ), are part of the inputs selected by global and target sensitivity analysis. No dependency between the PCT and one of the eight other scenario parameters is detected. Finally, seven additional inputs are added in X SII .

Note that some convergence studies (with respect to the sample size) not presented here show that the total number of inputs in {X PII ∪ X SII } is relatively stable: inputs in X SII passing through the group of X PII as n (and so, the detection capability of the tests) increases. This suggests the gradual and ongoing convergence of the method, as well as its robustness. It also highlights the interest of considering this group of X SII when screening results are not converged yet. 5

Step 3: Gaussian process metamodeling

Parametric choices for the Gaussian process and estimation of the hyperparameters

The second step of the ICSCREAM methodology consists in building a metamodel, whose aim is to fit the simulator output Y (here, the PCT), based on the learning sample (X s ,Y s ). To do so, we use a similar approach than the one used in [START_REF] Iooss | Advanced methodology for uncertainty propagation in computer experiments with large number of inputs[END_REF], based on an homoscedastic (non-interpolating) Gp metamodel. The reader can refer to [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF] for a detailed review about Gp metamodel. As introduced in Section 2, the inputs variables are divided into three groups at the end of Step 2:

• the ordered explanatory inputs X exp : this group is made up of the X PII to which the inputs to be penalized X pen are added (if not selected in X PII yet). In this last case, inputs to be penalized are added at the end, in an arbitrary order, or guided by expert judgment on their supposed importance. We obtain X exp = {X PII ∪ X pen };

• a second group of inputs X SII of much less influence;

• the remaining and neglected inputs merged into X ε with X ε = {X \ {X exp ∪ X SII }}.

The output is thus redefined by Y = M ({X exp ∪ X SII }, X ε ) and the metamodeling process is focused on fitting the random variable Y |{X exp ∪X SII }2 . In other words, only the inputs in {X exp ∪ X SII } are considered as the explanatory inputs of the Gp. Basically, the Gp is built to approximate the expected value E(Y |{X exp ∪ X SII }). The residual effect of the other inputs (merged into X ε ) is captured using an additional "nugget" effect. Borrowed from geostatistics, a nugget effect assumes an additive white noise effect and relaxes the interpolation property of the Gp metamodel. It can be assumed to be either constant (homoscedastic) or dependent on the variables (heteroscedastic). Note that contrary to [START_REF] Iooss | Advanced methodology for uncertainty propagation in computer experiments with large number of inputs[END_REF], a simple Gp metamodel with an homoscedastic nugget effect is estimated, since fitting a joint Gp with heteroscedastic nugget is neither relevant nor realistic given the high dimensionality of the problem and the rather small sample size. Concerning the parametric choices for the Gp, a constant trend or a one-degree polynomial with a selection algorithm (e.g., the Lasso algorithm [START_REF] Efron | Least angle regression[END_REF]) is considered. Moreover, a stationary anisotropic Matérn 5/2 tensorized covariance function is used for the inputs in X PII while X SII are joint in a stationary isotropic Matérn 5/2 covariance function (still tensorized with the one of X PII ).

Note that, in order to have a robust and sufficiently automated methodology, some parameters should be fixed. For instance, we recommend in practice to use a Matérn covariance function which is quiet popular in the machine learning community as it covers a large spectrum of applications. Moreover, the resulting differentiability properties of the Gp metamodel (one or two times mean-square differentiable for Matérn 3/2 and 5/2, respectively) is a reasonable compromise given the high-dimensional input space, the sparsity of data and the potentially low regularity of the code output. Anyway, various covariance functions can be tested in order to select the best one from the analysis of validation metrics (as done here).

All the Gp hyperparameters (i.e., covariance parameters) are estimated by maximum likelihood estimation on the learning sample. More precisely, a sequential process which relies on the ranking deduced from the screening step is used (see [START_REF] Welch | Screening, predicting, and computer experiments[END_REF][START_REF] Marrel | An efficient methodology for modeling complex computer codes with Gaussian processes[END_REF][START_REF] Iooss | Advanced methodology for uncertainty propagation in computer experiments with large number of inputs[END_REF]). To summarize, at each iteration j, one more input of X exp is added to the covariance function, following the order of inputs in X exp . The ( j -1) hyperparameters' values estimated at the ( j -1)-th iteration are used as starting points for the optimization algorithm for the ( j -1) first hyperparameters. Once all the inputs of X exp have been included, the last iteration consists in adding X SII through the isotropic covariance function (see [START_REF] Iooss | Advanced methodology for uncertainty propagation in computer experiments with large number of inputs[END_REF] for the detailed sequential building process).

Once the Gp hyperparameters are estimated, the Gp is conditioned by the observations of the learning sample to obtain the so-called "Gp metamodel", i.e., the resulting conditional random process, which is still a Gp, denoted by Y Gp ({X exp ∪ X SII }). For each unobserved prediction point, it is therefore fully characterized by its mean and variance (see [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF] for an explicit expression of Gp mean and Gp covariance). The conditional mean, denoted by Y Gp ({X exp ∪ X SII }), is used as a predictor. The conditional variance, denoted by MSE[ Y Gp ({X exp ∪ X SII })], is also the mean squared error (MSE) of the predictor. This prediction variance is used to build a confidence interval around the prediction. Covariance between predictions is also available via the conditional covariance.

Validation of the Gp metamodel

The accuracy and prediction capabilities of the Gp metamodel are assessed by "K-fold cross-validation" [START_REF] Hastie | The elements of statistical learning: Data Mining, Inference, and Prediction[END_REF], given the limited budget of simulations. Note that the criteria presented in the following can obviously be calculated on a test basis if available. However, for the considered applications, no test sample is available. For K-fold crossvalidation, rule-of-thumb methods suggest to fix large values of K (typically 5, 10 or 20) since it is usually preferable to exploit a larger number of simulations for training purposes, at the expense of a loss of generalization of the estimated errors.

First, to quantify the accuracy of predictor, one can use the predictivity coefficient

Q 2 : Q 2 = 1 - ∑ n i=1 y (i) -Y (i) Gp,-i 2 ∑ n i=1 y (i) -1 n ∑ n i=1 y (i) 2 (3)
where y (i) and

Y (i)
Gp,-i are respectively the i-th observation of the learning sample and the corresponding prediction of the Gp metamodel built without y (i) . Q 2 corresponds to the coefficient of determination in prediction, computed by cross-validation on the learning sample. The closer to one the Q 2 , the better the accuracy of the metamodel. A plot of predicted values against observed values ( Y (i) vs. y (i) ) or a quantile-quantile plot can also be drawn. In the purpose of identifying the input area yielding to critical configurations, we also compute the rate of good prediction of {Y > q 0.9 }.

Second, the quality of the Gp prediction intervals should also be evaluated (intervals built from the prediction variance). For this, the histogram of the predicted standardized residuals can be plotted (a standard Gaussian distribution should be observed). The Predictive Variance Adequacy (PVA) is also computed: this criterion defined by [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyperparameters of Gaussian processes with model misspecification[END_REF], and more recently used in [START_REF] Demay | Model selection for Gaussian Process regression: an application with highlights on the model variance validation[END_REF], assesses whether the prediction errors are of the same order of the prediction variances or not:

PVA = log    1 n n ∑ i=1    ∑ n i=1 y (i) -Y (i) Gp,-i 2 MSE( Y (i) Gp,-i )       . ( 4
)
The smaller the PVA, the more reliable the prediction intervals. For example, a PVA around 0.2 (respectively 0.7) corresponds to prediction variances around 20% (respectively two times) too large or too small with respect to the prediction errors. Finally, we consider a graphical tool introduced by [START_REF] Marrel | Global sensitivity analysis of stochastic computer models with joint metamodels[END_REF]. This tool enables to compare the Gp prediction intervals of level α with respect to the proportions of observations that actually lie within these intervals. These proportions (i.e., the "observed" confidence intervals) can be visualized against the α-theoretical prediction intervals, for different values of α in [0, 1]. By definition, the more the points are located around the identity line, the better the adequacy is. This plot will be called "α-α plot" in the following. For the interested reader, a detailed interpretation of these numerous validation criteria of the Gp metamodel (as well as additional ones, not mentioned in the present work) can be found in [START_REF] Demay | Model selection for Gaussian Process regression: an application with highlights on the model variance validation[END_REF].

Application on IB-LOCA uses-cases

Some of the aforementioned diagnostic metrics have been computed with K = 10 folds for the cross-validation procedure and are plotted in Figures 7 and8 for the two IB-LOCA use cases. In both cases, the predictivity coefficient Q 2 is around 0.8 which can be considered as satisfactory and reasonable regarding both the final objective of the study (i.e., the uncertainty propagation) and the multiple constraints which have to be handled (limited size of the learning sample, high dimension in input, strong nonlinearity of the model). Note that, in the case of poor predictivity performance, the use of Gp-based iterative enrichment strategies can be relevant to add supplementary code simulations. Again, such a strategy is not possible in the context of the present study. Thus, for the two use cases, less than 20% of the output variability remains not explained by the Gp metamodel. More precisely, this includes both the inaccuracy of the Gp and the loss of information attributable to the inputs in X ε which have not been selected in the metamodel building process. Another possible reason for the part of unexplained output variability can be the behavior of the code itself. For instance, such a locally irregular behavior can be due to physical bifurcations or threshold effects. In our case, a convergence study, for which all the variables in the model were added sequentially (i.e., even the X ε ), has been performed. Such a study, not presented here for the sake of brevity, showed that no improvement of the predictivity has been detected. It is therefore reasonable to assume that the code behavior also plays a role in this remaining output variability. Finally, high rates of good classification for critical configurations are also obtained: 94% and 93% for IB-LOCA 1 and IB-LOCA 2 , respectively.

Note that building the Gp directly in dimension d = 97, without variable selection and sequential inclusion processes, leads to a poor estimation of the Gp (typically, the optimization procedure involved in the estimation of the Gp hyperparameters often fails) and yields to a Gp with a poor predictivity. This illustrates the practical interest of the proposed methodology. For comparison purposes, other metamodels such as robust linear regression (e.g., elastic net regularization of [START_REF] Zou | Regularization and Variable Selection via the Elastic Net[END_REF]), support vector machines [START_REF] Schölkopf | Learning with kernels: support vector machines, regularization, optimization, and beyond[END_REF] or regression trees [START_REF] Breiman | Classification and Regression Trees[END_REF] were tested. All of them yield poor predictivity with Q 2 lower than 0.5. Only a generalized additive model based on smoothing splines, namely the "Adaptive COmponent Selection and Shrinkage Operator" (ACOSSO [START_REF] Storlie | Surface estimation, variable selection, and the nonparametric oracle property[END_REF]), with a Q 2 around 0.7, showed competitive results (for IB-LOCA 2 only) with respect to the proposed Gp-based approach.

Concerning the quality of predictive variance and confidence intervals, low PVA values (0.15 and 0.02 for IB-LOCA 1 and IB-LOCA 2 , respectively) and α-α plots show that the Gp metamodel yields accurate confidence intervals in prediction, although being sometimes too conservative for central values of α. The Gp metamodel provides a reliable prediction error: the conditional distribution of Gp is accurate. In conclusion, the Gp metamodel can be confidently used in a Bayesian framework for uncertainty propagation in the next step. 6

Step 4: Application for the identification of penalizing configurations

Definition and computation of conditional probabilities

The final goal of the ICSCREAM methodology is to identify the penalizing configurations (i.e., corresponding to a maximal PCT) of the inputs of interest X pen , under the uncertainty of the other inputs. To do so, one can define, for any set of possible values of X pen , the probability of exceeding the critical value q 0.9 . This conditional probability can be estimated for a given value x pen ∈ X pen with the Gp-Bayesian approach by:

P(x pen ) = P[Y Gp (X exp , X SII ) > q 0.9 |X pen = x pen ] = 1 -E[1 Y Gp (X exp ,X SII )≤ q 0.9 |X pen = x pen ] = 1 -E[1 Y Gp ( X exp ,X pen )≤ q 0.9 |X pen = x pen ] = 1 -E E 1 Y Gp ( X exp ,X pen )≤ q 0.9 | X exp | X pen = x pen = 1 - X exp Φ   q 0.9 -Y Gp ( x exp , x pen ) MSE[ Y Gp ( x exp , x pen )]   dP X exp ( x exp ) (5)
where X exp = {X exp ∪ X SII } \ X pen denotes the Gp inputs deprived of X pen , X exp their domain of variation, dP X exp their probability density function, and Φ(•) the cumulative distribution function of the standard Gaussian distribution. Note that the independence between X exp and X pen is required and used to obtain the fourth line from the third one in Eq. ( 5). In the so-called "Gp-Bayesian" approach, the predictive distribution of the Gp is entirely propagated (including the Gp error prediction); it is strictly equivalent to a calculation by conditional simulations. This approach therefore differs from a "plug-in" approach where only the Gp predictor (i.e., the Gp mean) is considered. A comparison between the Gp-Bayesian and plug-in approaches is given for example in [START_REF] Iooss | Advanced methodology for uncertainty propagation in computer experiments with large number of inputs[END_REF] (for a quantile estimation problem).

In practice, this conditional probability can be computed for the different values of each input in X pen , as well as for any subgroup of l inputs (l = 2, . . . , dim(X pen )). Thus, it results in a one (or l)-dimensional function.

Results for IB-LOCA use cases

The conditional probabilities given by Eq. ( 5) are, in practice, computed by intensive crude Monte Carlo Gp-based simulations with respect to the couple (X 92 , X 94 ) in IB-LOCA 1 and for each input of X pen in IB-LOCA 2 . Here, we use 2 × 10 4 Monte Carlo simulations to estimate each probability. Note that, for a target probability to be estimated close to 0.1, such a budget corresponds to a coefficient of variation of 2% approximately. Consequently, it is more than enough (especially considering the accuracy of the metamodel itself). The satisfactory accuracy of the estimated probabilities can also be verified by a convergence graph as the number of Monte Carlo Gp simulations increases. The estimated probabilities are given in Figures 9 and10, respectively.

For the IB-LOCA 1 use case, the analysis of conditional probabilities reveals the strong interaction of the two scenario inputs (X 92 and X 94 ) in the occurrence of critical configurations. The worst cases (i.e., where the probability reaches maximal values) correspond to medium values of both inputs. Conditionally to a given break size (X 92 ), the probability function according to the stopping time of primary pumps (X 94 ) has a bell shape that reaches its maximum value for a stopping time which decreases linearly as the break size increases. Note that the worst configuration is obtained for a break size equal to 3.57 inches and a stopping time of the primary pumps of about 907.8 seconds, which leads to a probability of exceeding the quantile q 0.9 = 673.18 • C estimated at 0.55. These results can be enlightened by an analysis of the physical phenomenon. In a nutshell, the correlation between these two scenario inputs drives the degradation of the water inventory. The smaller the break size, the longer the pump will have to run for the same inventory degradation. As for the distinct left and right limits on the domain, they can also be explained. On the one hand, if X 92 < 3.3 inches, meaning that the break size is rather small, the water inventory does not degrade too much (whatever the primary pump does). This leads to a slow LOCA which can be contained by the protection systems which have enough time to intervene (hence, the net border). On the other hand, when the break size increases too much, the break tends to be prevailing and reduces the impact of the stop time of the primary pumps (hence, the fading area).

For the IB-LOCA 2 use case, amongst the ten scenario inputs, only two are significantly influential on the probability of exceeding 671.85 • C. The enthalpy of the liquid accumulator (X 36 ) has a moderate linear effect: the probability is increased by 0.03 for the lower bound of X 36 (penalizing value) and decreased by the same amount for the upper bound. This effect is consistent with the physical knowledge related to the impact of this input on the condensation effect during the safety injection. More interestingly, the strong, complex and non-monotonic effect of the burn-up is clearly illustrated: the probability of PCT is maximal for the lower bound ( P(x 14,min ) = 0.175) then decreases before increasing to reach a local maximum value around x 14 29000 Figure 9: IB-LOCA 1 use case: conditional probability P(X 92 , X 94 ) of exceeding q 0.9 , estimated from the Gp, with both surface and contour plot representations. MWd/t, almost as high as the global maximum value (the variation range of this input is [515; 59000]). It then decreases prior to increasing until another smaller peak, and then falls to a zero probability of exceeding q 0.9 , for the highest values of X 14 . The burn-up, which is an integral measure of how much energy is extracted from a primary nuclear fuel source, describes the fuel depletion as it supplies energy. When the fuel is too depleted (values of burn-up ≥ 50000 MWd/t), the residual power significantly decreases as the burn-up increases, yielding to low PCT values. On the contrary, for lower depletion rate (X 14 < 50000 MWd/t), the relationship between the burn-up and PCT is much more complex since this input controls a large number of other physical quantities related to the life of the hot rod. Therefore, the implication of the burn-up in several physical sub-models (as the pressure inside the cladding, the initial temperature of the fuel and the neutronic data), constituent of the global thermal-hydraulic model, explains the form of P(X 14 ): the first small bumps are due to fluctuations in fuel temperature over the life of the rod, while the large drop at the end is due to the decay of the neutronic data at the end of the rod lifetime.

From a statistical methodology point of view, it is really interesting to clearly highlight some behavior, which are known by physicists, but impossible to detect on the initial dataset by simple visualization (cf. scatter plots of Figure 4). The results obtained with the proposed ICSCREAM methodology are therefore consistent with physical knowledge and even go beyond this. Indeed, they enable to quantify the impact on the probability but also to precisely identify and justify the penalizing areas in input.

Conclusion and prospects

In the framework of risk assessment in nuclear accident analysis, it is essential to quantitatively assess the uncertainties tainting the results of best-estimate computer codes. Beyond the usual uncertainty propagation step, this paper has focused on identifying the penalizing configurations of some specific inputs (the scenario ones), under the uncertainty of the other inputs. This methodology, called ICSCREAM, was motivated by a study at the reactor-scale of an IB-LOCA scenario in a pressurized water reactor, with Figure 10: IB-LOCA 2 use case: conditional probabilities of exceeding q 0.9 for the ten inputs to be penalized X pen , estimated from the Gp on their respective ranges of variation [X min ; X max ]. The mean probability P[Y (X exp , X SII ) > q 0.9 ] is also estimated to 0.118 and plotted in black dotted line.

the thermal-hydraulic CATHARE2 code. In the considered use cases, around a hundred scalar inputs are uncertain and the penalizing values of two or ten of them must be identified. The output variable of interest is the peak cladding temperature (PCT) during the accident transient. The critical configurations are defined by PCT exceeding its quantile of probability 0.9. The fidelity and complexity of the numerical modeling, the limited budget of simulation and the very large number of uncertain inputs are real challenges that led to the development of a sophisticated and robust methodology based on advanced statistical tools.

Applied from a single Monte Carlo sample of CATHARE2 simulations, the IC-SCREAM methodology judiciously combines a step of sensitivity analysis to identify and rank the main influential inputs and to reduce the dimension, with a specific building procedure (based on a sequential inclusion of variables) of a Gaussian process (Gp) metamodel. The sensitivity analysis step relies on new statistical independence tests which aggregate information of global and target Hilbert-Schmidt independence criterion (HSIC) measures. From the sensitivity analysis results, the inputs are taken into account in a more or less fine way in the Gp metamodel, according to their supposed influence. The non-selected inputs are integrated in the variance of the Gp prediction. The robustness of the metamodel is therefore enhanced and its building is made possible in such a high-dimensional problem. Once built, the accuracy and prediction capabilities of the Gp metamodel are assessed: several graphical and quantitative tools are proposed for this purpose. Finally, the Gp metamodel can be used intensively for uncertainty propagation, estimation of conditional probabilities and for the inversion problem. The conditional probabilities of PCT exceeding the critical value are here estimated using the Gp within a Bayesian framework to take into account the uncertainty of the Gp prediction.

The efficiency of ICSCREAM methodology has been demonstrated on two IB-LOCA use cases of high complexity and from a learning sample of around a thousand simulations. For both cases, 80% of the PCT variance is explained by the Gp, with a high ability for identifying PCT critical areas. Although the metamodel could still be improved, it provides a reliable prediction error and accurate prediction intervals. This justifies its use in a Bayesian framework for the estimation of conditional probabilities of exceeding the 90%-quantile, according to each input (or group of inputs) that have to be penalized. This estimation, not directly feasible with the numerical thermalhydraulic model because of its computational cost, becomes now tractable via intensive simulations of the Gp metamodel. The analysis of the conditional probabilities has revealed the strong and non-monotonic individual influence of some scenario inputs for one use case, and the high interaction of the two scenario inputs for the other one.

ICSCREAM has been developed in a particular industrial context where the user cannot add new simulations to the current design of experiments. In the case where additional simulations could be available, a research track would be to develop dedicated adaptive sampling strategies. Following this idea, efficient Gp-based sequential enrichment schemes could be adapted to add code simulations in an efficient manner (see, e.g., the works of [START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF][START_REF] Chevalier | Fast kriging-based stepwise uncertainty reduction with application to the identification of an excursion set[END_REF][START_REF] Moustapha | Quantile-based optimization under uncertainties using adaptive Kriging surrogate models[END_REF][START_REF] Wang | Gaussian process based optimization algorithms with input uncertainty[END_REF][START_REF] Sabater | A Bayesian approach for quantile optimization problems with high-dimensional uncertainty sources[END_REF] which deal with various purposes). However, their direct use in such large dimensional problems still remains a challenge. Another important limitation of the ICSCREAM methodology is its use on large-size datasets, typically several thousands of simulations. The estimation and building of the Gp metamodel then become difficult to achieve and compute. To make the Gp applicable in this context, methods based on nested kriging [START_REF] Rullière | Nested Kriging predictions for datasets with large number of observations[END_REF] or variational inference [START_REF] Hensman | Gaussian Processes for Big Data[END_REF] could be studied. Finally, handling the "chaotic" code behavior (e.g., physical bifurcations or threshold effects) could be achieved by non-stationary Gp metamodels such as the treed Gp which allows a smart domain partitioning [START_REF] Gramacy | Bayesian Treed Gaussian Process Models With an Application to Computer Modeling[END_REF].

Appendix -Illustration on an analytical test case

Initially, the ICSCREAM methodology was directly developed to be applied to very high-dimensional real datasets. To further illustrate the robustness and efficiency of the method, we propose here an additional analytical test case, in lower dimension, but applying the same rules of thumb for the sample size. This case is designed to be "representative" of some features of the underlying physics at stake in an IB-LOCA model, in the sense that it exhibits interaction and non-monotonic conditional probabilities. This case is proposed to bring additional information to the previous applications by illustrating how a predictivity around 80% allows nevertheless to correctly estimate the conditional probabilities. Moreover, it can be seen as a tool for future research prospects, for example on adaptive sampling in large dimension.

The considered analytical model M f , inspired from the Friedman function [START_REF] Friedman | Multivariate Adaptive Regression Splines[END_REF] and defined in dimension d f = 20, is given by: M f (X) = a 1 sin 6π X 5/2 1 (X 2 -0.5) + a 2 (X 3 -0.5) 2 + a 3 X 4 + a 4 X 5 + r X6,...,X 15 [START_REF] Picheny | Adaptive Designs of Experiments for Accurate Approximation of a Target Region[END_REF] where r X 6 ,...,X 15 = The model depends only on the fifteen first inputs. The first term represents a strong and non monotonic interaction between the two first inputs. The second term is a quadratic function of X 3 while the other ones are linear. The inputs that have to be penalized are (X 1 , X 2 , X 3 , X 4 , X 5 ). The parameters for tuning the influence of the different inputs are chosen as follows: a = (5, 20, 8, 5, 1.5) . Under this parametrization, X 2 explains alone around 10% of the output variance, X 1 has no individual effect but its interaction effect with X 2 is strong (around 30% of the output variance). X 3 , X 4 and X 5 only have individual effects (no interaction) and explain, respectively, around 11%, 28% and 10% of the output variance. The effects of the ten remaining inputs (X 6 to X 15 ) represent around 11.5% of the output variance. Finally, we compute reference values (by intensive Monte Carlo simulations) of the 90%-quantile of the model output (approximately equal to q 0.9 = 14) and the one-dimensional mean effects as well as the one-dimensional conditional probabilities of exceeding q 0.9 , according to the inputs to be penalized (see Fig. 11).

Figure 11: M f analytical case: reference one-dimensional mean effects and conditional probabilities of exceeding q 0.9 for the five inputs to be penalized X pen .

To apply the ICSCREAM methodology, we consider a Monte Carlo sample of n = 200 evaluations of the M f model and the 90%-quantile is estimated from this learning sample ( q 0.9 = 13.8 is obtained here). Step 2 is performed with a level α test = 0.05 which yields to the selection of inputs X 1 to X 5 by HSIC measures and X 1 to X 6 by target HSIC ones. From an aggregation procedure based on Bonferroni's correction, X 1 to X 5 are selected in X PII and X 6 in X SII . From this selection, a Gp metamodel is built with the same parametric choices as in Section 5.3. Predictivity coefficients of Q 2 = 0.78 and Q 2 = 0.77 are obtained on a test basis and by cross-validation on the learning sample, respectively. Regarding the predictive variance, PVA is around 0.15. Diagnostic metrics of the Gp performance validation are illustrated in Figure 12. We observe that, despite the amount of unexplained variance (i.e., slightly more than 20%), the Gp prediction intervals are satisfactory. We are therefore in a configuration quite similar to those obtained on the two IB-LOCA use cases.

From the Gp metamodel, the conditional probabilities are estimated and provided in Figure 13. We observe a good approximation of the one-dimensional conditional probabilities (by comparing to the reference curves given in Fig. 11 right): first, the different behaviors are well captured; second, the penalizing areas of each input are clearly identified. Finally, a convergence study has also been performed: it shows the predictivity and Gp performance improvements as n increases (Q 2 around 0.85 and 0.94, for n = 400 and 600 respectively).

In the future, such a test case can be extended by adding, for instance, some locally irregular behavior or an interaction between an input to be penalized X A and another input X B (not to be penalized) whose global influence would be rather small on the output but significant on the conditional probability to X A . In this last case, the risk would be that X B would not be selected by HISC-based tests (e.g., if the sample size is too small), which would lead to a poor estimation of the conditional probability of X A . The challenge will be to be able to detect this type of error, perhaps by relying on convergence and stability plots of HSIC, or by building more powerful HSIC-based tests. Finally, the modeling and propagation of the "observed" error of the metamodel could be investigated and compared to the predictive error of the metamodel. Figure 13: M f analytical case: conditional probabilities of exceeding q 0.9 for the inputs to be penalized, estimated from the Gp.

Figure 1 :

 1 Figure 1: General workflow of the ICSCREAM methodology.

  line accumulator of the intact loop (X 34 ∼ U[800, 1900] m -4 ), discharge line accumulator of the broken loop (X 35 ∼ U[800, 1900] m -4 ), liquid accumulator enthalpy (X 36 ∼ U[33544, 213105] J/kg), temperature of safety injection (X 90 ∼ U[7, 50] • C), weight factor on flow rate of safety injection 1 (X 91 ∼ U[-1, 1]), tube plugging of steam generator of the intact loop (X 92 ∼ U[0, 0.09]), weight factor on flow rate of the emergency feedwater system of steam generator (X 93 ∼ U[-1, 1]), temperature of the emergency feedwater system of steam generator (X 94 ∼ U[7, 55] • C). It is a 15 IB-LOCA case with LOOP.

Figure 2 :

 2 Figure 2: Histogram of the PCT from the learning sample for IB-LOCA 1 (left) and IB-LOCA 2 (right) datasets. The estimated 90%-quantiles are indicated by a black dotted line.

Figure 3 :

 3 Figure 3: Convergence plots of the empirical estimator of 90%-quantile, according to the learning sample size n, for IB-LOCA 1 (left) and IB-LOCA 2 (right).

Figure 4 :

 4 Figure 4: IB-LOCA 2 : scatter plots of the learning sample with local polynomial regression of the PCT according to several scenario inputs.

  (a) Global-SA-oriented screening.(b) Target-SA-oriented screening.

Figure 5 :

 5 Figure 5: IB-LOCA 1 use case: p-values of HSIC and T-HSIC-based independence tests computed from the learning sample. The level α test = 5% is represented in black dotted line.

  (a) Global-SA-oriented screening. (b) Target-SA-oriented screening.

Figure 6 :

 6 Figure 6: IB-LOCA 2 use case: p-values of HSIC and T-HSIC-based independence tests computed from the learning sample. The level α test = 5% is represented in black dotted line.

  (a) Predicted vs. observed values.α-α plot.

Figure 7 :

 7 Figure 7: IB-LOCA 1 use case: diagnostics of Gp performance.

Figure 8 :

 8 Figure 8: IB-LOCA 2 use case: diagnostics of Gp performance.

a 5 √

 5 (∑ i=6...15 i 2 ) ∑ i=6...15√12i(X i -0.5) and X = (X 1 , . . . , X 20 ) are independent and uniform random variables on [0, 1].

  (a) Predicted vs. observed values. (b) α-α plot.

Figure 12 :

 12 Figure 12: M f analytical case: diagnostics of Gp performance.

The input data considered for this case, and thus the obtained PCT, do not correspond to actual industrial values.

Y |{X exp ∪ X SII } (i.e., Y knowing {X exp ∪ X SII }) is still a random variable as its value depends on the uncontrollable random vector X ε .

Acknowledgments

The writing of this paper was partly funded by the ANR project SAMOURAI (ANR-20-CE46-0013). We warmly thank our two reviewers whose remarks have greatly helped to improve the paper. We are also grateful to Vincent Larget, Mathieu Segond, Laurent Lefebvre and Jean-Christophe Lecoy for numerous discussions about this work. CATHARE2 code is developed under the collaborative NEPTUNE project, supported by CEA, EDF, Framatome and IRSN.