Edgar Jaber
email: edgar.jaber@edf.fr

Vincent Chabridon

Emmanuel Remy

Michael Baudin

Didier Lucor

Mathilde Mougeot

Bertrand Iooss

Sensitivity Analyses of a Multi-Physics Long-Term Clogging Model For Steam Generators

Keywords: Global Sensitivity Analysis, Polynomial Chaos Expansions, Sobol' indices, Hilbert-Schmidt Independence Criterion, Clogging, Steam Generators

Long-term operation of nuclear steam generators can result in the occurrence of clogging, a deposition phenomenon that may increase the risk of mechanical and vibration loadings on tube bundles and internal structures as well as potentially affecting their response to hypothetical accidental transients. To manage and prevent this issue, a robust maintenance program that requires a fine understanding of the underlying physics is essential. This study focuses on the utilization of a clogging simulation code developed by EDF R&D. This numerical tool employs specific physical models to simulate the kinetics of clogging and generates time dependent clogging rate profiles for particular steam generators. However, certain parameters in this code are subject to uncertainties. To address these uncertainties, Monte Carlo simulations are conducted to assess the distribution of the clogging rate. Subsequently, polynomial chaos expansions are used in order to build a metamodel while time-dependent Sobol' indices are computed to understand the impact of the random input parameters throughout the whole operating time. Comparisons are made with a previous published study and additional Hilbert-Schmidt independence criterion sensitivity indices are computed. Key input-output dependencies are exhibited in the different chemical conditionings and new behavior patterns in high-pH regimes are uncovered by the sensitivity analysis. These findings contribute to a better understanding of the clogging phenomenon while opening future lines of modeling research and helping in robustifying maintenance planning.

Introduction

Steam generators (SGs) are heat exchangers within pressurized water nuclear reactors (PWRs). PWRs are composed of two distinct water circuits responsible for the heat exchange. In this system, water from the primary circuit initially passes through the pressure vessel where it is heated due to nuclear reactions in the core. After this, it proceeds to the SG where it transfers its heat to the water in the secondary loops (mainly, three or four loops, hence three or four SGs). The resulting steam exits the SGs through their upper openings and then flows through turbines to generate electrical power. In this process, the primary fluid circulates within a bundle of Ushaped tubes, which are stabilized by tube-support plates (TSPs) (a general SG scheme and TSP is shown in Fig. 1). Adjacent to these tubes, the secondary fluid flows through the holes of the TSP and undergoes vaporization due to the heat from the primary fluid in the tubes. After prolonged periods of operation, a clogging phenomenon may appear in the flow-holes of the TSPs of some SGs. Clogging refers to a reduction in the flow area within the TSPs due to the accumulation of iron oxide deposits. These deposits primarily result from the corrosion of the secondary circuit, after which the oxides are carried by the secondary flow. This phenomenon has diverse impacts on the SGs, in particular it causes a localized redistribution of flow among the TSP flow-holes, elevating the risk of vibration, SG tube ruptures, and hampering the SGs' response to operational changes. To address this issue, the clogging rate (denoted as τ c ) is evaluated during PWR outages using non-destructive video examinations (an example is shown in Fig. 1). The clogging rate is an average quantity that gives a measure of how much the holes of a TSP are clogged. Subsequently, preventive maintenance consisting in chemical cleanings procedures are carried out to reduce the clogging rate and the secondary fluid pH can be supervised by chemical conditioning.

To achieve better control of this phenomenon and better plan the maintenance, EDF R&D has developed a physico-chemical clogging model and numerical simulation chain known as "THYC-Puffer-DEPOTHYC" [START_REF] Feng | Tube support plates blockage of PWR steam generators: thermalhydraulics and chemical modeling[END_REF][START_REF] Prusek | Modélisation et simulation numérique du colmatage à l'échelle du sous-canal dans les générateurs de vapeur[END_REF], that will be called "THYC-Puffer-DEPO" (TPD) in the rest of this paper for the sake of clarity (and to avoid any confusion between the THYC computer code and the DEPOTHYC one). This model is calibrated with respect to the real in-situ data by determining the optimal value of a scalar calibration parameter named a v that will be detailed later (see equation ( 2)). This computational chain is a spatio-temporal simulator that takes a large number of inputs in order to fully model an entire SG and all the U-shaped tubes as shown in Fig. 1. However, among all these inputs, most of them are related to geometric quantities describing the SG and are thus purely deterministic, while only a few of them are tainted by uncertainties, either due to physical variability, or because of lack of knowledge about these physical quantities.

In accordance with the generic probabilistic framework for the treatment of uncertainties (De Rocquigny, Devictor, & Tarantola, 2008), input probability distributions have been elicitepd by expert judgment in order to quantify and propagate these uncertainties. Thus, a first study has been proposed by [START_REF] Lefebvre | Improving the Predictivity of a Steam Generator Clogging Numerical Model by Global Sensitivity Analysis and Bayesian Calibration Techniques[END_REF]. In this work, several techniques have been applied, from uncertainty propagation through Monte Carlo sampling, surrogate modeling based on an artificial neural network, global sensitivity analysis using Sobol' indices [START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF][START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF] and Bayesian calibration of the a v parameter to enhance and update the sensitivity analysis' results. This first work enabled to open the path to some sound and relevant thoughts about the computational model and the input physical variables involved. However, this previous study [START_REF] Lefebvre | Improving the Predictivity of a Steam Generator Clogging Numerical Model by Global Sensitivity Analysis and Bayesian Calibration Techniques[END_REF] was conducted on an earlier version of the computational chain (named "COLMATHYC") that did not account for the long temporal evolution of clogging.

The primary objective of the present paper is to adopt a rather similar approach as the one proposed in [START_REF] Lefebvre | Improving the Predictivity of a Steam Generator Clogging Numerical Model by Global Sensitivity Analysis and Bayesian Calibration Techniques[END_REF] (without going up to the Bayesian calibration step in the present study) while tackling several major differences with this previous work: firstly, a new version of the computational chain providing a long time-dependent evolution of the clogging rate; secondly, to try to incorporate the chemical conditioning (which is a key information arising from the operating strategy and that may drastically affect the clogging dynamics); and thirdly, to develop a fully "given-data" methodology in order to estimate complementary sensitivity measures that might be useful for a better understanding of the clogging phenomenon, its better monitoring and a more adapted predictive maintenance. All in all, these uncertainty studies, coupled with new future data, could help enhancing the whole computational chain [START_REF] Cheng | Machine learning with data assimilation and uncertainty quantification for dynamical systems: a review[END_REF]. The term "given-data" here means that the whole methodology relies on a single input-output dataset of Monte Carlo simulations. Such a framework is often encountered in industry as such complex simulation codes are expensive-toevaluate and thus require to be deployed on high performance computing facilities (clusters) in order to make studies be tractable. Therefore, the main difficulty is to collect the maximum amount of relevant information (e.g., in order to to build an efficient surrogate model or to perform sensitivity analysis) using this finite sample and dedicated estimators. The aim of global sensitivity analysis (GSA) is to study how the output variability can be aportioned to the variability of the inputs [START_REF] Da Veiga | Kernel-based ANOVA decomposition and Shapley effects -Application to global sensitivity analysis[END_REF][START_REF] Saltelli | Global Sensitivity Analysis[END_REF]. These two core topics, i.e surrogate modeling and GSA, play a central role in the UQ framework and have been widely studied in the literature. For instance, one can mention the works of [START_REF] Heuveline | A hybrid generalized polynomial chaos method for stochastic dynamical systems[END_REF] for the use of PCE applied to stochastic dynamical systems and [START_REF] Schöbi | Polynomial-Chaos-based Kriging[END_REF] for an efficient combination between PCE and Kriging in order to take advantage of both types of surrogate models. Regarding GSA, much effort has been put on providing robust estimators of Sobol' indices [START_REF] Janon | Uncertainties assessment in global sensitivity indices estimation from metamodels[END_REF][START_REF] Lamboni | Global sensitivity analysis: an efficient numerical method for approximating the total sensitivity index[END_REF][START_REF] Puy | A comprehensive comparison of total-order estimators for global sensitivity analysis[END_REF]. The present work tries to exploit some rather standard tools based on polynomial chaos expansions (PCE) based surrogate model in order to estimate time-dependent Sobol' indices, but also more recent kernel-based indices based on the so-called Hilbert-Schmidt independence criterion (HSIC) indices [START_REF] Da Veiga | Global Sensitivity Analysis with Dependence Measures[END_REF][START_REF] De Lozzo | Sensitivity analysis with dependence and variance-based measures for spatio-temporal numerical simulators[END_REF][START_REF] Gretton | Measuring statistical dependence with Hilbert-Schmidt norms[END_REF][START_REF] Marrel | Statistical developments for target and conditional sensitivity analysis: application on safety studies for nuclear reactor[END_REF] whose formulation is purely dedicated to the given-data framework.

The present paper is organized as follows. In Section 2, the physical and numerical modeling of the clogging phenomenon is briefly presented as well as the probabilistic modeling of the input uncertainties. Numerical results obtained from Monte Carlo simulations are presented and analyzed as a starting point of the proposed study. Then, Section 3 presents the PCE surrogate model and its ability, in a time-dependent context, to easily provide Sobol' indices while putting an emphasis on the robust and accurate design of such a metamodel. Section 4 proposes to apply various types of HSIC-based indices (global, target and conditional) in order to capture some complementary information about time-dependent sensitivities. All the numerical results are discussed and analyzed regarding the current available knowledge in the literature. Finally, Section 5 draws some conclusions about the proposed tools and gives some perspectives for future work.

2 Physics and numerics of the clogging phenomenon in steam generators

Overview of the physical modeling

The clogging model proposed and described in [START_REF] Feng | Tube support plates blockage of PWR steam generators: thermalhydraulics and chemical modeling[END_REF][START_REF] Prusek | Modélisation et simulation numérique du colmatage à l'échelle du sous-canal dans les générateurs de vapeur[END_REF] is decomposed into three main steps. These steps are naturally associated with the different time scales in the clogging phenomenon. First, the two-phase flow of the secondary fluid inside the SG is described by a system of conservation laws for mass, momentum and energy. A stationary state of this system is then computed. Among these stationary thermal hydraulic quantities, there is ρ m , C ℓ and U ℓ which are respectively, the density of the liquid-gas mixture, the liquid quality and the velocity of the liquid phase. The liquid phase velocity field serves as a flux to a system of two transport equations governing the mass fractions of iron oxide Γ p , Γ s in its different chemical states (solid particles and soluble particles, respectively). To these transport equations, a solubility map of iron-oxide Γ max s inside the SG is provided as well as initial conditions Γ p (0), Γ s (0). Finally the third phenomenon is a local balance equation for the time variation of the mass of iron oxide m c with fluxes Φ p , Φ s linked respectively to the vena contracta and the flashing mechanisms [START_REF] Prusek | Modélisation et simulation numérique du colmatage à l'échelle du sous-canal dans les générateurs de vapeur[END_REF].

The clogging model boils down to the following system of mixed partial and ordinary differential equations (with the two initial conditions on the last line):

         ∂ t Γ p + U ℓ • ∇ Γ p = f p (Φ p , Φ s , Γ s , Γ max s ) ∂ t Γ s + U ℓ • ∇ Γ s = f s (Φ p , Φ s , Γ s , Γ max s ) dmc dt = Φ p + Φ s Γ s (0, .) = Γ s (0), Γ p (0, .) = Γ p (0) , (1)
where, for u = {s, p}, one has Γ u = ρ m C ℓ Γ u and where f u are essentially additive functions of their variables. As for the particle flux related to the vena contracta mechanism which is of interest here, it is given by [START_REF] Prusek | Modélisation et simulation numérique du colmatage à l'échelle du sous-canal dans les générateurs de vapeur[END_REF]:

Φ p = a v k v (ρ p -ρ ℓ ) U 2 z d 2 p µ ℓ Γ p , (2)
where a v is the calibration non-dimensional parameter of the vena contracta mechanism, k v is a geometric parameter of the flow-hole, ρ p , ρ ℓ are the densities of particles and of the liquid phase respectively, U z := U m • e z is the vertical component of the homogeneous two-phase mixture velocity, µ ℓ is the dynamic viscosity and d p is the diameter of the particles. Once the system above is solved, the mass m c is multiplied by the volumic mass of iron oxide to obtain the volume V c . Finally the average clogging rate τ c is obtained through the following empirical correlation:

τ c = α (1 -exp{-βV c }) . (3)

Overview of the numerical model

For the numerical model description, computer codes are denoted by functions g * :

x * → g * (x * ) where * stands for the name of the code and x * is its vector of deterministic input variables. In case the input vector is random, capital bold letters X * will be used. The output will sometimes be written Y * = g * (X * ).

The stationary state of the conservation laws for the two-phase mixture is obtained using the THYC software developed by EDF R&D [START_REF] David | Three Dimensional Thermal-Hydraulic Simulation In Steam Generators With THYC Exchangers Code -Application[END_REF], denoted by g THYC in the following. It is based on a numerical finite-volume scheme using a porous-medium approach to provide results at the component scale. Once the thermal hydraulic state Y THYC is computed, the g DEPO code is used to solve the transport-balance equations. The output variable of interest is a clogging rate increment between two instants t < t ′ :

∆ t ′ t τ c := τ c (t ′ ) -τ c (t) . (4)
Moreover, to better control the clogging rate, the pH of the secondary circuit may be modified by addition of basic organic solvents. This action is named chemical conditioning and the solubility maps of iron oxide have to be computed accordingly.

This is done with the code g Puffer that provides the iron oxide solubility maps Γ max s for a given SG geometry, thermal hydraulic state and chemical conditioning.

For long periods of simulation (e.g., simulating more than 10 years of the SG operating time), the stationary state for the thermal hydraulic quantities may be not physically valid anymore. Indeed the clogging deposits build additional headlosses that need to be taken into account. This is why a coupling has to be made between the three preceding codes which gives rise to the global architecture of the THYC-Puffer-DEPO simulation chain presented in Fig. 2. The code is initialized at iteration j = 0 by computing the stationary thermal hydraulic state y 0 THYC . The code g THYC is supposed to have no uncertainties on its input variables. The vector x 0 THYC contains SG geometry parameters and different physical correlations that are considered well known. After this, the code g Puffer provides a solubility map y Puffer for a given chemical conditioning and SG geometry in x 0 Puffer . Again, g Puffer is supposed to have deterministic input variables. Once the solubility map y Puffer is built, the input vector for the deposit module g DEPO is the concatenated outputs of the preceding codes as well as the uncertain input variables X DEPO that will be described further. The latter contains more than just uncertain variables but for the sake of clarity, it is assumed that it contains only the inputs related to the uncertainty analysis performed. The deposit module g DEPO will calculate clogging increments ∆ j+k j+k-1 τ c between times t j+k-1 and t j+k after which several conditions are checked. If the total amount of clogging between t j and t j+k is greater than 30%, then the head-losses have to be updated and the stationary thermal hydraulic states recalculated. The inputs X j THYC are updated and some of its components become random given the uncertainty on the output of g DEPO . If there is a modification of the chemical conditioning at time t χ , then there is a recalculation of the solubility maps. If t j+k corresponds to a maintenance cleaning date of the SG, then the clogging rate is lowered according to the type of cleaning performed and g THYC is relaunched. Finally, if t j+k corresponds to the ending time of the simulation, t f , the code provides the variation curve of the clogging rate τ c . For such a computational chain, a typical simulation horizon is 50 years (i.e., the typical operating duration of a SG without considering any other issue leading to a potential replacement). Finally, a unitary call of THYC-Puffer-DEPO takes about 5 hours on an high-performance computing infrastructure.

Uncertainty quantification and propagation

Among the numerous input variables of the code, the ones retained by experts for the uncertainty quantification treatment consists of d = 7 independent input variables:

X DEPO = (α, β, ϵ e , ϵ c , d p , Γ p (0), a v ) .
(5)

As defined earlier, this vector does not contain all g DEPO input variables but only gathers the uncertain ones. In comparison with the study proposed in [START_REF] Lefebvre | Improving the Predictivity of a Steam Generator Clogging Numerical Model by Global Sensitivity Analysis and Bayesian Calibration Techniques[END_REF], the solubility map Γ max s is not considered here. The reason for this is that the previous version of the code -COLMATHYC -did not use g Puffer for computing the solubility map and the clogging modeling has evolved since. The probabilistic modeling of these input variables is provided in Table 1. The input probability distributions as well as supports have been set by expert judgment and are in line with the previous study [START_REF] Lefebvre | Improving the Predictivity of a Steam Generator Clogging Numerical Model by Global Sensitivity Analysis and Bayesian Calibration Techniques[END_REF]. The probability laws elicitepd are either Gaussian (X i ∼ N (µ, σ 2 ) with mean µ and variance σ 2 ) or triangular (X i ∼ T (a, b, c) with mean b on the interval [a, c]). More information concerning the physics of the variables can be found in [START_REF] Prusek | Modélisation et simulation numérique du colmatage à l'échelle du sous-canal dans les générateurs de vapeur[END_REF]. Table 1 Probabilistic modeling of uncertain input variables.

In order to perform the uncertainty propagation through the THYC-Puffer-DEPO chain, a coupling was built with the OpenTURNS (Open source initiative for the Treatment of Uncertainties, Risks ' Crude Monte Carlo simulations were executed by drawing n = 10 3 samples according to the previously mentioned joint probability distribution of X DEPO defined in Table 1. In order to leverage the computational cost, these simulations were launched in parallel on a cluster, in order to simulate a t f = 50 years period (approximately, one week on high-performance computing infrastructure for the whole computational process). The resulting trajectories are shown in Fig. 3. The remaining sections of this paper focus on the results obtained for the hot leg (HL) at the uppermost section of the SG in z max (see the SG scheme in Fig. 1). This choice stems from the fact that the top portion of the HL is known to be the most clogged within the SG and is the primary area inspected by video examinations and therefore used for maintenance planning.

From the simulations in Fig. 3, several remarks can be drawn. Firstly, in the case of this particular SG, two distinct cleaning procedures have been executed all along the operating duration time of the SG: one is a curative treatment while the other one is a preventive cleaning. Typically, operators are aware that preventive cleaning is less effective in reducing clogging rates compared to curative cleaning. Secondly, one can see that the secondary fluid has undergone three distinct chemical conditionings at two different pH levels. These conditionings are denoted by χ i , i ∈ {1, 2}, representing the two main chemical species, while the pH is represented as a binary variable with low or high as outcomes. A notable observation is that these various chemical conditionings yield differences in the clogging kinetics. This disparity becomes obvious in the data presented in Fig. 4, which illustrates the distributions of a clogging delta at different times in specific chemical conditionings. The findings reveal a certain deceleration of the simulated clogging kinetics when employing a more basic chemical conditioning. This fact is supported by operational observations. In addition to the time-dependent analysis, some scatter plots illustrating the relationships between inputs and output in the rank-space are provided in Fig. 5, at three distinct time points: t 0 < t 1 < t 2 , which correspond to the three different chemical conditionings.

τ c (t 0 ) τ c (t 1 ) α τ c (t 2 ) β e c d p Γ p (0) a v
Fig. 5 Scatter plots in the rank-space for the three different chemical conditionings.

The dependency structures observed in these scatter plots closely resemble those presented in [START_REF] Lefebvre | Improving the Predictivity of a Steam Generator Clogging Numerical Model by Global Sensitivity Analysis and Bayesian Calibration Techniques[END_REF]. Specifically, a nonlinear correlation is observed between the diameter of iron-oxide particles d p and the clogging rate τ c at all the time points. In the physical space, this nonlinear correlation reveals that the particles responsible for high-clogging rates are the ones of intermediate diameter at all chemical conditionings and clogging levels. Additionally, a nonlinear correlation is better seen within at time point t 2 (corresponding to the χ 2 with high-pH regime) between the clogging rate porosity ϵ c and τ c . It is important to note that this particular correlation is a novel finding compared to the previous study in [START_REF] Lefebvre | Improving the Predictivity of a Steam Generator Clogging Numerical Model by Global Sensitivity Analysis and Bayesian Calibration Techniques[END_REF] where the porosity of fouling and clogging were considered to be the same (i.e., ϵ e = ϵ c = ϵ). Thus, the sensitivity analysis results found it to be non influential. It is worth pointing out that in this prior study, the pH did not play any important role since the solubility maps were constant on all the SG mesh and distributed according to the uncertainty law chosen. In the more advanced version of the clogging computational chain used in our work, the Puffer code provides a mesh-varying solubility map of iron oxide according to the different chemical conditionings.

Finally, linear correlations are observed between the clogging rate and the initial condition of particle mass-fraction Γ p (0) as well as for the calibration parameter a v . Consistent with the previous study, no significant correlation is identified with the empirical correlation parameters α and β, as well as the fouling porosity rate ϵ e . These few preliminary visual observations on the input-output dependency-structures will be more accurately and quantitatively studied with the results of the given-data sensitivity analysis in the following sections.

3 Variance-based sensitivity analysis using polynomial chaos expansions

Time-dependent surrogate modeling via polynomial chaos expansions

Given the time-costly nature of the computational chain and the main objective of performing global sensitivity analysis, a polynomial chaos expansions(PCE) metamodel is chosen as a suitable candidate regarding these purposes. Note that, if any other well-known family of metamodels would have been suitable too (e.g., Gaussian processes, artificial neural networks, support vector machines, ...), the use of a PCE metamodel here is motivated by the low-dimensional problem, the relatively low-size training sample (n = 10 3 ) and the robust functional guarantees of these types of approximations (Le Maître & Knio, 2010). On top of that, the ability of this metamodel to provide Sobol' indices as a by-product is a non-negligible advantage when one wishes to perform given-data global sensitivity analysis. Similar methodology can be found in [START_REF] El Garroussi | Tackling random fields non-linearities with unsupervised clustering of polynomial chaos expansion in latent space: application to global sensitivity analysis of river flooding[END_REF][START_REF] Méndez Rojano | Uncertainty quantification of a thrombosis model considering the clotting assay PFA-100®[END_REF].

The clogging rate output of the TPD code is now written as a generic function g:

g : (t, z, X) ∈ [0, t f ] × [0, z max ] × X → g(t, z, X) ∈ [0, 100] , (6)
where X := X DEP ∈ X = X 1 × . . . × X d ⊂ R d ∼ P X is the vector of inputs in (5), X is the cartesian product of the support of the distributions in Table 1 and P X the product measure of these laws (since the input variables are assumed to be independent). For a choice of inputs X 0 , the output is a function (t, z) → g(t, z, X 0 ). As already mentioned, all the simulation results will be considered at the specific altitude z max which corresponds to the top of the highest TSP in the HL of the SG. Therefore the output is a random vector of dimension equal to the time discretization N :

(g(t 1 , z max , X 0 ), . . . , g(t N , z max , X 0 )) ∈ R N , (7)
where in this case t N = t f . It is assumed the vector model lies in L 2 P X (X ; R N ). A Hilbert basis {φ α } α∈N d made of tensorized orthonormal polynomials is then used [START_REF] Xiu | The Wiener-Askey Polynomial Chaos for Stochastic Differential Equations[END_REF]. For Gaussian distributions, Hermite polynomials are used and for triangular distributions, since the orthonormal family is not analytically explicit, an adaptive Stieltjes algorithm (Gautschi, 2004, Ch .2, Sec. 2.2.3.1) is used for building it (to be more specific, an adaptive version of the initial algorithm based on a Gauss-Kronrod quadrature rule is implemented in OpenTURNS). At each time step, the PCE of highest degree p ∈ N is constructed as:

g(X) = α∈J g α φ α (X) , (8)
where

g α = g 1 α , . . . , g N α ∈ R N for all α ∈ J ⊆ {α ∈ N d , d i=1 |α i | ≤ p}.
The index set J = J 1 ∪ . . . ∪ J N is determined by performing, for each k ∈ {1, . . . , N }, the following two-step method [START_REF] Blatman | Anisotropic parcimonious polynomial chaos expansions based on the sparsity-of-effects principle[END_REF], [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF]):

• A hyperbolic enumeration rule of quasi-norm 0 < q ≤ 1 is chosen. This method allows to enhance the selection of marginal indices in the sparse basis for a choice of small quasi-norm; • The coefficients G k := {g k α } α∈J k are calculated by solving the following leastsquare problem using the least-angle regression method [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF] :

G k = arg min g k ∈R |J k | E   g k (X) - α∈J k g k α φ α (X) 2   . (9)
Once the coefficients have been computed at each time-step, the index set is built by interweaving the basis terms.

In the metamodel-building procedure, the quality of the metamodel is measured through the use of the predictivity coefficient Q 2 (Fekhari, Iooss, Muré, Pronzato, & Rendas, 2023). Since the problem is time-dependent, this coefficient is estimated at each time step, using a test dataset D test of size m < n at every time discretization (t 1 , . . . , t N ) over [0, t f ]. Typically, the output trajectories presented in Fig. 3 are discretized in N = 75 timesteps. Thus, one has:

∀k ∈ {1, . . . , N }, Q 2 (t k ) = 1 - m j=1 g(t k , X test j ) -g(t k , X test j ) 2 Var g(t k , X test j ) , Q 2 = 1 N N k=1 Q 2 (t k ) , (10
) with Q 2 the time-averaged predictivity coefficient. Consequently, the hyperparameters q and p are optimized while maximizing the evolution of Q 2 .

To select the best hyperparameters, the variability of Q 2 with respect to the permutation of the training/test samples is studied. From the available data batch, 75% is used for training and 25% for testing while 5 permutation attempts are made. The q-norm is studied by fixing a maximal degree of 2. This q-norm being the number of basis elements retained, it is obvious that, for a full basis, the metamodel could overfit, hence the fall in Q 2 observed in Fig. 6. From these results, and in order to minimize the computation time, a q-norm of 0.5 is chosen as a good compromise (Q 2 ≈ 0.89), since increasing it does not really improve Q 2 .

Then, the degree p is optimized and a satisfactory level of Q 2 is obtained for a degree 4 at a q-norm of 0.5 as shown in Fig. 7 (Q 2 ≈ 0.97). Finally, by considering the evolution of the sequence Q 2 (t k ) for t k ∈ {t 1 , . . . , t N }, one can clearly see a variation of the coefficient through time, with a decrease in the regions near the chemical cleaning discontinuities. However, globally, the PCE metamodel obtained can be confidently used for performing sensitivity analysis.

Q 2 (%)
Fig. 7 Boxplots for degree validation with q = 0.5.

Time-dependent Sobol' indices based on polynomial chaos expansion

Reminders on the theory

Assuming that, for every k ∈ {1, . . . , N }, the output component g k ∈ L 2 P X (X , R), one has access to the unique "functional analysis of variance" (FANOVA) decomposition [START_REF] Hoeffding | A class of statistics with asymptotically normal distribution[END_REF][START_REF] Sobol | Sensitivity analysis for non-linear mathematical models[END_REF][START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF] such that:

g k (X) = γ⊆{1,...,d} g k γ (X γ1 , . . . , X γ d ) = γ⊆{1,...,d} g k γ (X γ ) . ( 11
)
For every γ = (γ 1 , . . . , γ d ), with |γ| ∈ {1, . . . , 2 d -1}, the Sobol' index of order |γ| at time t k is defined by the ratio of the variance of g k γ over the total output variance:

S γ (t k ) = Var g k γ (X γ ) Var (g k (X)) . ( 12
)
While several sampling-based techniques are available in the literature to estimate these indices (see, e.g., (Da Veiga et al., 2021, Ch. 4) for a review), the PCE framework enables their direct computation as a simple by-product as proposed in [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF]. By setting:

J γ = α ∈ N d , |α| ≤ p : α ℓ > 0, ∀ℓ ∈ γ & α ℓ = 0, else , (13)
the terms in ( 13) can be regrouped and identified as the unique factors in the FANOVA decomposition such that:

g k γ = α∈Jγ g k α φ α (X) . ( 14
)
Sobol' indices are then directly assessed with simple calculations. Using the orthonormality of the basis terms, the following expression is obtained:

S γ (t k ) = α∈Jγ (g k α ) 2 |α|≤p (g k α ) 2 . ( 15
)
In what follows, special attention is given to the first-order index S i and the totalorder index S T i for i ∈ {1, . . . , d}, as well as the global interaction index S * , defined by: 16)

S i (t k ) = Var g k i (X i ) Var (g k (X)) , S T i (t k ) = γ|i⊂γ S γ (t k ) and S * (t k ) = 1 - d i=1 S i (t k ). (
The first-order and total-order indices correspond respectively to the influence of a specific input X i on the output variance and the total effect of input X i taking into account its interactions with the other inputs. Moreover, it will be of interest to plot the non-normalized first order indices t k → Var g k i (X i ) .

Analysis of numerical results on the THYC-Puffer-DEPO chain

In Fig. 8, the temporal evolution of the variance of every input is presented. Likewise, in 9 temporal evolutions of various Sobol' indices computed from the PCE metamodel are provided, regarding the three chemical conditionings which are recalled in the background of the plots. Before delving into the details, several initial observations merit consideration. Notably, alongside the instances of chemical cleaning, there are discernible spikes in the first and total order indices. It is imperative to emphasize that these spikes are a numerical artifact stemming from the construction of the metamodel and its susceptibility to discontinuities within the data. Another remark is that the proposed study corroborates the prior results obtained in [START_REF] Lefebvre | Improving the Predictivity of a Steam Generator Clogging Numerical Model by Global Sensitivity Analysis and Bayesian Calibration Techniques[END_REF]. The predominant portion of the output variance is primarily influenced by a subpart of the inputs, mainly d p , a v and Γ p (0). It is worth noting that, in contrast to the previous study, the solubility maps Γ max interactions terms have a negligible impact all along time evolution. This result is consistent with the total-order indices (see Fig. 9, on the right) where the indices seem to be almost equal to the first-order indices. As a consequence, one can say that the TPD code behaves here, conditionally to the input probabilistic model given in Table 1, almost as an additive function of the input variables. Moreover, this refined version of the code affords the opportunity to discern the impact of chemical conditioning on the hierarchy of influential variables. Specifically, in the χ 2 with high-pH regime, there is a notable transfer of variance contribution from d p and Γ p (0) to ϵ c , while the influence of a v remains relatively constant across the three chemical conditionings. The porosity rate of clogging deposits, denoted by ϵ c , quantifies the proportion of free space volume to the total volume occupied by the iron oxide deposit. It serves as a measure of how effectively the iron oxide particles are aggregated within the system. It is conventionally assumed to remain constant throughout the SG and is considered independent of the operational conditions. Our discovery presents an intriguing avenue for various physical interpretations. Iron oxides inherently possess charged particles that establish crystalline bonds, primarily Van der Waals interactions, based on their electronic structure. Altering the hydrogen ion charge of the solution can have a profound impact on this electronic structure. This phenomenon is notably exemplified in colloidal aggregation [START_REF] Liu | Aqueous aggregation behavior of citric acid coated magnetite nanoparticles: Effects of pH, cations, anions, and humic acid[END_REF]. In light of this, it can be inferred that under high-pH conditions, the iron oxide deposits tend to exhibit a greater degree of aggregation. Consequently, this translates into an increased significance of uncertainty in clogging porosity on the system's output. Furthermore, the stiffening or rigidification of clogging deposits in a high-pH regime could potentially account for the observed deceleration in clogging kinetics. These findings shed light on the intricate interplay between chemical conditions, particle aggregation, and clogging behavior within the system.

4 Kernel-based given-data sensitivity analysis 4.1 The Hilbert-Schmidt independence criterion for global sensitivity analysis

Reminders on the theory

Another method for performing sensitivity analysis is through the use of kernel methods like the Hilbert-Schmidt independence criterion (HSIC). Introduced in the machine learning community by [START_REF] Gretton | Measuring statistical dependence with Hilbert-Schmidt norms[END_REF], its use as a powerful importance measure in sensitivity analysis is due to [START_REF] Da Veiga | Global Sensitivity Analysis with Dependence Measures[END_REF]. It allows to uncover deep dependency structures between each input variable and the output quantity of interest. The underlying idea is to compare the joint probability distribution P Xi,Y with the product of distributions P Xi P Y . To do so, it makes use of generalized operators in Reproducing Kernel Hilbert Spaces (RKHSs) [START_REF] Berlinet | Reproducing Kernel Hilbert Space in Probability and Statistics[END_REF]. In (7), the output space is denoted by Y = Y 1 ×. . .×Y N = g(t 1 , X )×. . .×g(t N , X ) ⊂ R N where N is the time discretization. Let i ∈ {1, . . . , d} and k ∈ {1, . . . , N } be an input parameter index and a time-output index. Let F i , G k be two RKHSs of respective kernels κ i , κ k that are assumed to be characteristic [START_REF] Sriperumbudur | Universality, Characteristic Kernels and RKHS Embedding of Measures[END_REF]). On the space

F i × G k consider the kernel ν ik = κ i ⊗ κ k , defined for all (X, X ′ ) ∈ X 2 i , (Y, Y ′ ) ∈ Y 2 k by: ν ik ((X, Y ), (X ′ , Y ′ )) = κ i (X, X ′ )κ k (Y, Y ′ ) . (17)
To alleviate the notations for Y ik = g(t k , X i ), one will consider the pair (X i , Y ik ) which will be written as (X i , Y k ). Finally, the generalized covariance operator is defined as:

C XiY k = E Xi,Y k [ν ik ((X i , Y k ), .))] -E Xi E Y k [ν ik ((X i , Y k ), .))] , (18)
and the HSIC measure between the input X i and the output Y k is defined as the Hilbert-Schmidt norm of the generalized covariance operator:

HSIC(X i , Y k ) = ∥C XiY k ∥ HS . (19)
One theoretical result which strongly motivated the use of HSIC indices for global sensitivity analysis is that [START_REF] Da Veiga | Global Sensitivity Analysis with Dependence Measures[END_REF][START_REF] Gretton | Measuring statistical dependence with Hilbert-Schmidt norms[END_REF]:

HSIC(X i , Y k ) = 0 ⇐⇒ X i ⊥ Y k . (20)
Therefore, assessing independence between an input and an output can be achieved through the use of HSIC indices. From a statistical viewpoint, one needs to estimate this quantity from a dataset. Various estimators have been derived, either based on U-statistics or V-statistics (Da Veiga et al., 2021, Ch. 6). For instance, the following trace-based estimator is often used in practice:

HSIC(X i , Y k ) = 1 n 2 Tr (L i HL k H) , (21)
where L i , L k are the Gram-matrices on the samples (x

(p) i , y (p)
k ) 1≤p≤n , thus given by

(L i ) pq = κ i (x (p) i , y (q) k ) and (L k ) pq = κ k (x (p) i , y (q)
k ) and H is a shift matrix defined by (H) pq = δ p,q -1 n pq for all 1 ≤ p, q ≤ n. For a given input-output pair (X i , Y k ), a null-hypothesis can be derived:

H 0 : HSIC(X i , Y k ) = 0. ( 22
)
Let S T = n × HSIC(X, Y ) be the test statistic and S T,obs be the observed statistic.

The corresponding p-value associated to H 0 is then defined by:

p val = P( S T ≥ S T,obs |H 0 ). (23)
A large p-value is a sign of independence and in our applications a standard limit value p val = 0.05 is chosen above which independence achieves statistical significance. It is important to emphasize that it is not a measure of dependency but its non-nullity is a potential indicator of dependence and of influence on the output. Compared to the Sobol' indices, there is no underlying FANOVA decomposition for the standard formulation of HSIC indices as used in the present paper (see [START_REF] Da Veiga | Kernel-based ANOVA decomposition and Shapley effects -Application to global sensitivity analysis[END_REF] for a deeper discussion). Usually a normalized index is used [START_REF] Marrel | Statistical developments for target and conditional sensitivity analysis: application on safety studies for nuclear reactor[END_REF] but the latter does not sum to one. It is defined by:

R 2 HSIC (X i , Y k ) = HSIC(X i , Y k ) (HSIC(X i , X i )HSIC(Y k , Y k )) . ( 24
)
Deriving such FANOVA methods for HSIC indices is an ongoing topic of research [START_REF] Da Veiga | Kernel-based ANOVA decomposition and Shapley effects -Application to global sensitivity analysis[END_REF]. The HSIC index only yields a global analysis and is not sensitive to the interactions between the inputs (as with Sobol' indices). The use of HSIC here can be seen as complementary to the previous analysis in Section 3.

Analysis of time-dependent global HSIC indices

Fig. 10 shows the corresponding R 2 HSIC indices and the associated p-values of the statistical tests, carried out on the TPD Monte Carlo sample used previously. The sensitivities to the discontinuities are not present anymore. This is due to the given-data estimation without relying on any approximation of the function. The main conclusions from the Sobol' analysis here persist, namely that the influential variables are a v , d p and Γ p (0) in all chemical conditionings and ϵ c in the high-pH-χ 2 conditioning. Moreover the ranking of the influential variables for every fixed time is also the same as in the previous analysis. Here the use of the entire dataset allows to compute the p-value with an asymptotic estimator. All the influential variables are well under the threshold value. It can be concluded that with this complementary computation, the same results hold and therefore the robustness and relevance of the sensitivity analysis is achieved.

Target HSIC indices

Reminders on the theory

With the help of HSIC indices, one can also easily access the target indices. In other words, this is about uncovering which indices contribute the most to the outputs belonging to a given critical region C [START_REF] Marrel | Statistical developments for target and conditional sensitivity analysis: application on safety studies for nuclear reactor[END_REF]. It makes use of a filter function f C : Y → R + that selects the simulations falling in C by transforming the output distribution. Usually a hard-thresholding with an indicator function 1 C is possible but this can forget outputs that are "close" to the region of interest and also leads to a notable decrease of the available data for estimation. This can be overcome using a modified covariance kernel as explained in [START_REF] Marrel | Statistical developments for target and conditional sensitivity analysis: application on safety studies for nuclear reactor[END_REF]. The filter function used is defined by:

Y j = f C (Y j ) = exp - inf y∈C ∥Y j -y∥ s σ Yj ∈ [0, 1] , (25)
where s = 1/5 and σ Yj is the estimated standard deviation of the output sample at t j . The hyperparameter s may seem arbitrary but is the one retained for this type of analysis [START_REF] Marrel | Statistical developments for target and conditional sensitivity analysis: application on safety studies for nuclear reactor[END_REF]. It can be optimized accordingly but this is out of the scope of this work. The V-statistic estimator (20) can still be used with the modified variable to yield the target HSIC indices:

T-HSIC(X i , Y j ) = HSIC(X i , Y j ), (26)
and the corresponding p-value associated to the hypothesis test H 0 can still be used.

Analysis of time-dependent target HSIC indices

A threshold of 70% is chosen as a target region for the clogging simulation dataset. Indeed this highly-clogged regime is critical for the operation of the SG. The associated variation of T-HSIC indices amount is displayed in Fig. 11. For a certain amount of time, the T-HSIC index is negligible because there is no relevant data on the output set. Indeed all the outputs are weighted to zero because they are still far from the target region. Once the target set size increases, one can see the influential variables appear but the statistical test is not yet significant. Here a permutation-based estimator is used for computing the p-value, given the scarcity of the data. Indeed, this is linked to the number of critical simulations that appear progressively. From this local analysis of the highly-clogged simulations, it becomes manifest that the clogging porosity ϵ c becomes the most influential variable. It should be stressed that this was apparent in the previous analyses in Figs. 9 and 10 in the long-time high-pH χ 2 conditioning. Indeed it constitutes a new finding that is here quantitatively established with the help of the local target indices.

Conditional HSIC indices

Reminders on the theory

Conditional sensitivity analysis aims at identifying the influential input variables knowing that the output falls within a certain critical region C. Reformulating this allows to uncover the variables driving the event {Y ∈ C}. This can be achieved with the help of the filter function f C , by a standard procedure [START_REF] Marrel | Statistical developments for target and conditional sensitivity analysis: application on safety studies for nuclear reactor[END_REF]. The conditioning formula gives:

Y j |{Y j ∈ C} ∼ P f Yj , s.t ∀ measurable A ⊂ Y, P f Yj (A) = A f C (y) E [f C (Y j )]
dP Yj (y). ( 27)

By conditioning, the underlying input space is also modified; thus, for all i ∈ {1, . . . , d}, X i |{Y j ∈ C} ∼ P f Xi defined similarly. The joint distribution involved in the general covariance operator ( 18) is changed accordingly (X i , Y j ) ∼ P f Xi,Yj . Moreover a Vstatistic estimator similar to (21) can be derived.

Analysis of time-dependent conditional HSIC indices

The identical threshold value as in the target analysis is here maintained. The same progression in target set sizes persists due to the utilization of the identical weight function (25). Notably, a significant distinction arises, as there is a noticeable absence of a null C-HSIC index throughout the entire simulation duration. This phenomenon can be attributed to the influence of conditioning on the input variables; however, the resultant effect lacks interpretability when considered in the context of sensitivity analysis. This phenomenon is evident in the uniform distribution of HSIC "shares" at the initial stages of the upper plots in Fig. 12. Once this conditioning effect stabilizes, the ranking of C-HSIC at the 70% threshold remains consistent, and the previous findings regarding the clogging porosity rate are reaffirmed. Following the implementation of preventive cleaning and under high-pH conditions in the χ 2 chemical conditioning, this variable emerges as the most influential, while the previously dominant variable, d p , becomes negligible, as further corroborated by the p-value plot.

Conclusion and perspectives

In this study, the complex clogging rate phenomenon observed in the steam generators of pressurized water reactors is studied with the help of a panel of surrogate modeling and sensitivity analysis tools. The complex multiphysics THYC-Puffer-DEPOTHYC clogging code, developed by EDF R&D, is used and a complete uncertainty quantification and propagation is performed using Monte Carlo simulations. As a result, a n-size sample (n = 10 3 ) of output trajectories is obtained for a specific steam generator while simulating the time evolution of the clogging rate over 50 operating years. The analysis of the results enable to recover qualitative results that are in accordance with a previous published study [START_REF] Lefebvre | Improving the Predictivity of a Steam Generator Clogging Numerical Model by Global Sensitivity Analysis and Bayesian Calibration Techniques[END_REF] which tried to tackle the same phenomenon but with a different computational chain. However, the present study enabled to highlight new results in terms of inputs-output influence. These qualitative results are corroborated by quantitative results obtained via time-dependent Sobol' indices obtained as a by-product of a high-fidelity metamodel based on polynomial chaos expansions. Through these indices, one can understand that the interactions between inputs are relatively weak, leading to a quasi-additive structure of the model. The provided time-dependent Sobol' indices enable to identify and rank the most influential inputs. In addition, the use of another class of indices based on the Hilbert-Schmidt Independence Criterion (HSIC) enable to offer a refined (in terms of the nature of influence, i.e going beyond the output variance) and complementary insight about the input-output relationships. Such indices are versatile since they can offer three types of analyses (global, target and conditional). All in all, the results are still in coherence with the previous ones obtained with the variance-based importance measures, but they shed new light on the ranking in the different output critical domains, with respect to the evolution in time. As a perspective, some improvements could be achieved through the use of other class of metamodels (see, e.g [START_REF] Schöbi | Polynomial-Chaos-based Kriging[END_REF]) and other sensitivity indices could be investigated (see, e.g [START_REF] Lüthen | Global sensitivity analysis using derivative-based sparse Poincaré Chaos Expansions[END_REF]).

Fig. 1

 1 Fig. 1 SG scheme, TSP and example of video examination during an PWR outage (© EDF).

Fig. 2

 2 Fig. 2 THYC-Puffer-DEPO code architecture.

  First empirical correlation parameter in (3) N (101.6, 4.0) β Second empirical correlation parameter in (3) N (0.0233, 0.0005) ϵe Porosities of the fouling deposits T (0.2, 0.3, 0.5) ϵc Porosities of the clogging deposits T (0.01, 0.05, 0.3) dp Iron oxide particle diameter (m) T (0.5, 5.0, 10.0) × 10 -6 Γp(0) Initial data for solid particles mass fraction transport equation T (1.0, 4.5, 8.0) × 10 -9 av Calibration parameter of the vena contracta physical mechanism T (0.1, 7.8, 12) × 10 -4

Fig. 4

 4 Fig. 4 Densities of clogging kinetics for various chemical regimes.

Fig. 6

 6 Fig.6Boxplots for the q-norm validation with p = 2.

sFig. 8

 8 Fig. 8 Evolution of the variance contribution of each uncertain input.

Fig. 9

 9 Fig. 9 Evolution of first-order, interactions (left) and total-order (right) time-dependent Sobol' indices.

  Fig. 10 Time variation of the estimated HSIC indices and p-values associated to the statistical test of independence.

Fig. 11

 11 Fig. 11 Time variation of the estimated, normalized T-HSIC indices, the associated target set size and the p-values associated to the statistical test of independence.

Fig. 12

 12 Fig. 12 Time variation of the estimated, normalized C-HSIC indices and p-values associated to the statistical test of independence.

  N Statistics, see https://openturns.github.io/www/)

			1000 clogging trajectories -HL at z max .
	χ 1 , low pH
	χ 1 , high pH
	χ 2 , high pH
	Mean trajectory
	t 0 , t 1 , t 2	τ c
	Curative cleaning
	Preventive cleaning
		2000	4000	6000	8000	10000	12000	14000	16000
					t (d)
	Fig. 3 Trajectories obtained from Monte Carlo simulations
	(with curative and preventive chemical cleanings).

software, which is an open source C++/Python library developed by a consortium of industrial partners and public research institutions (Airbus Group, EDF, IMACS, ONERA and Phimeca)

[START_REF] Baudin | Open TURNS: An industrial software for uncertainty quantification in simulation[END_REF]

.

Acknowledgments

The authors would like to thank several people: Dr. Qingqing Feng, Dr. Stéphane Pujet and Morgane Garo Sail (EDF R&D / MFEE Department) for fruitful discussions about the clogging computational model, Jérôme Delplace (EDF Nuclear Division) for discussions about the industrial aspects as well as Dr. Régis Lebrun (Airbus) for fruitful discussions about OpenTURNS and the "Adaptive Stieltjes Algorithm" used in this work. This work is part of a PhD program funded by the French National Association for Technological Research (ANRT) under Grant n o 2022/1412.